Interaction vs inhomogeneity in a periodic TASEP

We study the non-equilibrium steady states in a totally asymmetric simple exclusion process with periodic boundary conditions, also incorporating (i) an extra (nearest-neighbour) repulsive interaction and (ii) hopping rates characterized by a smooth spatial inhomogeneity. We make use of a generalized mean-field approach (at the level of nearest-neighbour pair clusters), in combination with kinetic Monte Carlo simulations. It turns out that the so-called shock phase can exhibit a lot of qualitatively different subphases, including multiple-shock phases, and a minimal-current shock phase. We argue that the resulting, considerably rich phase diagram should be relatively insensitive to minor details of either interaction or spatial inhomogeneity. As a consequence, we also expect that our results help elucidate the nature of shock subphases detected in previous studies.


Introduction
The totally-asymmetric simple-exclusion process (TASEP) is an elementary stochastictransport model, with a variety of applications to out-of-equilibrium physics, ranging from biological microsystems to vehicular traffic [1,2], and intriguing theoretical connections (Kardar-Parisi-Zhang universality, random matrix theory) [3,4].For such reasons, and since it is simple enough that it can be treated analytically up to a certain extent, this model has inspired a lot of fundamental studies in non-equilibrium statistical mechanics.The model is usually defined on a linear chain, whose nodes can be occupied by at most one particle, and where each particle can hop to the adjacent node, provided the latter is empty: totally asymmetric means that hopping can occur in one direction only.Over the years, many variants of the model have been investigated, with the aim of incorporating extra features, which may be relevant in different physical contexts.For example, a first important issue is the distinction between periodic or open boundary conditions.In the latter case, if the chain is coupled to particle reservoirs that inject and extract particles at opposite ends, the system exhibits different non-equilibrium steady states and a number of phase transitions among them, controlled by injection and extraction rates (boundary-induced phase transitions) [5].This noticeable case is still exactly solvable, and the steady-state solution dates back to the 1990s [6,7,8,9].
Among other interesting issues, that of spatial inhomogeneity of hopping rates has been frequently considered in the literature.In case of periodic boundary conditions, it turns out that several different types of inhomogeneities (a unique slower rate [10,11,12], randomly distributed rates [13,14], smoothly varying rates [15,16,17]) give rise to the same kind of phenomena, namely, steady states featuring domain walls (or shocks), which separate regions at different densities.Similar shocks can also emerge if the system is allowed to adsorb and/or desorb particles in its bulk (Langmuir kinetics) with proper size scaling [18,19,20,21,22,23].Note that, in such cases, the rates are spatially uniform, but the density profile is not (even in the steady state), due to violation of the particle-number conservation constraint.A different type of extra feature, that can be incorporated in the model, is some kind of interaction beyond the basic exclusion constraint, in the simplest case a nearest-neighbour (NN) interaction [24,25,26,27,28,29,30,31,32,33,34].Such interactions can give rise to many different non-equilibrium steady states, and therefore very complex phase diagrams (in particular, up to 7 phases for strong enough repulsive interaction and open boundary conditions) [25,26,27,28,29,30,31].
Of course, it may also be of interest to investigate models incorporating more than one of the above features, as the latter may likely occur together in physical systems.For instance, [35] investigates a TASEP-like model with both hoppingrate inhomogeneity and Langmuir kinetics, which the authors argue to be minimal "ingredients" for modeling motor-protein motion along microtubules in living cells.Still motivated by biological transport phenomena, as well as by methodological issues, [36] investigates the role of NN interactions in a TASEP with Langmuir kinetics, though Three regimes observed in [38], as a function of the (mean) particle density: (a) no interaction or weak interaction; (b) mildly strong interaction; (c) very strong interaction.Phase tags are explained in the text.
without inhomogeneities.Furthermore, [37,38] investigate the interplay between NN interactions and rate inhomogeneities, in the presence of open [37] or periodic [38] boundary conditions.In particular, the recent work by Pal and Gupta [38] (periodic boundary conditions) points out that a NN repulsive interaction, in combination with a very simple rate inhomogeneity (represented by a regular function with a single minimum and a single maximum) gives rise to a seemingly non-trivial phenomenology.The presence of three different regimes is revealed, as the repulsion strength increases (see figure 1).In case of weak repulsion (figure 1a), the behavior of the system is qualitatively equivalent to that of a system without NN interaction [16], with a single phase at intermediate densities, characterized by the presence of a shock in the density profile and maximal current (S MC phase), and by two phases, at high and low density, characterized by a spatially-modulated density profile but without shock, and lessthan-maximal current (HD and LD phases).This regime corresponds to a (local) current-density relation featuring a unique "central maximum", at density 1/2.In case of stronger repulsion, where it is known that the current-density relation develops a "central minimum" at density 1/2 and two "side maxima", the other two regimes emerge.In one of them, at "mildly strong" interaction (figure 1b), the S MC phase is characterized by two shocks, one of which (the "descending" one) localized at the minimum-rate position.Furthermore, the separation regions, between this phase and the HD and LD phases, are characterized by a single shock, somewhat similar to the "ordinary" shock phase (i.e. that of a model without NN interaction), but with apparently non-maximal current.The nature of these phases (called "S phases") remains not completely clarified, since in this regime there appear large discrepancies between theory and simulations, and the analysis relies mainly on the latter, which are unfortunately affected by significant finite-size effects.In the other regime, at "very strong" interaction (figure 1c), the presence of yet another S phase is observed at densities close to 1/2, with a minimum of the current between the two maximal-current S MC regions.
In this paper we aim at elucidating the phenomenology that emerges from the work described above, which we find quite interesting, especially since, as previously discussed, it originates from the interplay of just two simple extra "ingredients" of the model, compared to ordinary TASEP (rate inhomogeneity and repulsive NN interaction).To do this, we are going to study a model in principle identical to that of Pal and Gupta, but with a different constraint on the rates, which we will call the KLS condition (since it was originally introduced by Katz, Lebowitz and Spohn [24,25,26]).The underlying idea is explained below.On the one hand, we expect that the overall physical behaviour remains qualitatively unaffected by the change in the rate constraint, specifically because it can be seen that the modified model undergoes the same kind of "transition" (from unimodal to bimodal) in the current-density relation.On the other hand, it is known that the steady state of the model with KLS condition, in the simpler case of homogeneous rates, coincides with the equilibrium state of a onedimensional lattice gas with NN interaction (a detailed proof is given by Dierl, Einax and Maass [29]).As a consequence, a generalized mean-field theory including NN pair correlations, which we denote as pair approximation (PA), is exact in this case.Thus, for our system, characterized by the assumption of smoothly-varying rates, we expect that the PA remains very accurate.Indeed, all the numerical evidences collected in this work suggest that it exactly reproduces the smooth sections of the density profiles, so that we can obtain a very detailed analysis of the various phases, practically without finite-size effects.The paper is organized as follows.In sections 2 and 3 we introduce the model and the PA theory, respectively.Section 4 describes, in terms of density profiles, the different regimes occurring in the non-equilibrium steady states of the model.The results are made more systematic in section 5, which presents and explains the analytical phase diagrams, and the way they have been worked out.Section 6 is devoted to numerical simulations, in order to check the accuracy of the PA.We recap our findings and draw some conclusions in section 7. The technical details are reported in four appendices.

The model
The model we consider is defined on a linear chain of L nodes, labelled as i = 0, . . ., L − 1, with periodic boundary conditions.Each node can be empty or occupied by a single particle.We introduce occupation-number variables n t i , taking value 0 or 1 if node i at time t is respectively empty or occupied.Each particle can hop in only one direction (conventionally from node i to node i + 1), provided the destination node is empty.The rate of hopping from node i (to i + 1) may depend on both the position i and on the i − 1 and i + 2 node configurations, so that it will be denoted as w i (n i−1 , n i+2 ).This type of dependence can be regarded as a NN interaction.Node indices are always understood modulo L, according to periodic boundary conditions.
To begin with, we shall consider the most general case where, for each node i, we have 4 different hopping rates (associated with all possible occupancy states of the 2 nodes i − 1 and i + 2), with a fully arbitrary dependence on the position i.The PA theory will be first developed in such a general case.Subsequently, we shall take some more restrictive hypotheses, along the lines of [38] and previous literature [27,28,29,30,36,37].In particular, we shall assume the following expressions for the rates, whose meaning is discussed below: The first restrictive assumption contained in (1) is that the spatial dependence of the rates is fully incorporated in a unique function λ(x), being independent of the occupancy states of the forward and backward nodes.As a consequence, the ratios between rates associated to different forward/backward occupancies are specified only by the position-independent coefficients p, q, r, s.In the following we will often call λ(x) the rate modulation function.The second important assumption is the continuity of the rate modulation function, meaning in practice that the spatial dependence of the rates is "smooth", as it occurs only through the scaled position variable x i/L (in order to respect boundary conditions, we can assume that λ(x) is periodic with period 1).Furthermore, as anticipated in the introduction, we shall assume that the rate modulation function is characterized by a unique minimum and a unique maximum (over its period).In particular we shall consider a sinusoidal function, taking values in a positive interval [λ min , λ max ], that is With this particular choice (irrelevant to the qualitative behaviour of the model) the rate modulation function turns out to be completely specified by the only 2 parameters λ min and λ max .Throughout our analysis, we shall usually consider one of them (for instance λ min ) fixed as a reference value, varying the ratio λ max /λ min , which we shall call rate modulation ratio.We shall also introduce some restrictive hypotheses on the p, q, r, s parameters.As in Pal and Gupta's work [38], we shall first assume p = s, which turns out to induce a particle-hole symmetry.‡ In fact, one of the four parameters can be set equal to 1 without loss of generality, so the imposed condition will actually be leaving q and r as free parameters.Actually, we find it more convenient to define and thence to consider v and q as independent parameters, with r fixed by (4).We shall call v the interaction parameter, as it discriminates the attractive case (q > r, or v > 1) ‡ Throughout the paper, we shall often use the term "particle-hole symmetry" to denote the mirror symmetry of the current-density relation, although this is just one aspect of the symmetry.
from the repulsive one (q < r, or v < 1).It can be seen that the situation studied in [38] corresponds to adding the simple condition q r = p s , which, together with the symmetry hypothesis (3) and definition (4), leads to This assumption is based on a thermodynamic consistency argument, and was previously proposed by Dierl, Maass and Einax [27], so we shall denote it as the DME condition.
Conversely, as anticipated in the introduction, we will assume the KLS condition, namely which, still in combination with (3) and ( 4), entails It is worth noting that, even though all the physical results that we report in the article have been obtained with the KLS condition, most of the related analytical theory is developed without it, in particular retaining only the symmetry assumption (3) and keeping both v and q as free parameters.

The analytical methods
Let us now present the PA, in the most general case introduced above.It is useful to define first some notation for marginal distributions and expectation values.Let us denote the marginal distribution at time t, for a cluster of consecutive nodes starting at i, by Moreover, let us define specific symbols for the local densities and NN correlations, respectively: As we are dealing with binary random variables, the latter quantities completely specify 1-node and 2-node marginals (for NN pairs), and can be used to parameterize them as follows: and Now, from the master equation one can derive time-evolution equations for ρ t i and φ t i , through a marginalization procedure (a detailed derivation is given in [31]).As far as the local densities are concerned, we obtain a typical continuity equation, namely where J t i represents the probability current from i to i + 1 at time t.Such current can be written as a sum of 4 contributions, one for each possible combination of backward and forward occupation states, in formulae where Regarding local correlations, we have the following equation In the end we can see that the time-derivatives of ρ t i and φ t i can be written in terms of 4-node marginals, so the resulting time-evolution equations, though exact, are not closed.A possible closure scheme is naturally provided by the PA [39], also known as Bethe approximation [40] or 2-cluster approximation [2].In the specific case, such approximation reads where 1-node (site) and 2-node (NN pair) marginals can be expressed in terms of ρ t i and φ t i through equations ( 12) and ( 13).An analogous technique has already been applied to the NN-interacting TASEP, for instance in [36] (for a system with Langmuir kinetics) and [37] (with site-dependent hopping rates).In these papers, the resulting equations have been used as an intermediate step to determine a continuum limit, leading to slightly different methods, respectively denoted as cluster mean-field and correlated cluster mean-field.Conversely, in this paper we adopt a simpler strategy, as done in [31], that is, we perform a direct time-integration of the discrete system (at finite L), obtaining the steady state as a long-time limit of the numerical procedure (more precisely, we define the limit by requiring that the magnitude of all time derivatives stays below a certain threshold).This procedure does not need a considerable computational power, so we can easily reach quite large sizes (all the results reported in this section refer to L = 10000), such that the density profiles obtained are (in most cases) practically indistinguishable from the continuum limit.Apart from this difference, our approach is equivalent to the cluster mean-field theory of [38], since both are based on the same pairwise factorization (18).We nonetheless prefer to retain the old-fashioned term PA, in order to avoid confusion with analogous approximation strategies, that take into account higher-order clusters [2,41].
In this work we shall also make use of the current-density relation, in a form derived from the PA theory.This relation, which by definition refers to the continuum limit, is precisely a function returning the value of the current given that of the local density.Such a function is completely specified by the local values of the hopping rates, which play the role of parameters.Collectively denoting by w an array of 4 possible rates (associated with the different occupancy states, as described above), we can write the current-density relation as Since the rates w depend on the position, equation ( 19) establishes a constraint to the spatial variations of the density ρ, such as to impose that the value of the current J remains fixed.Now, let us observe that, if all the hopping rates are multiplied by the same factor, say λ, the current gets multiplied by the same factor, in formulae In mathematical terms we may say that the function F is homogeneous (of degree 1) with respect to the parameters.This property follows from a simple physical argument, but it can also be verified explicitly (see Appendix A).In our case it is specially important, since, as introduced above, the position dependence of the rates is given by a common prefactor, namely the rate modulation function λ(x).As a consequence, in order to determine the possible density profiles, we need to invert a single current-density function (i.e. the one characterized by the position-independent parameters p, q, r, s), according to equation For the model considered here, the inversion can be done analytically (see Appendix D).
In the following, we shall usually call reduced current the (position-dependent) quantity defined as As known, the current-density function is not invertible in the narrow sense or, in other words, its inverse is a multi-valued function.Therefore, equation (21) does not determine a unique density profile, but a (small) set of possible profiles.The actual profile can be determined by imposing that the integral average of the function ρ(x) is equal to the actual mean particle density in the system, say ρ, which can thus be regarded as a control parameter.In formulae Since, according to equation ( 21), the profile ρ(x) depends not only on the parameters but also on the value of the current J , equation ( 23) determines the latter as a function of ρ.It can be seen that, in certain ranges of ρ values, there does not exist a single continuous profile, being solution of (21).In such a case, the actual stationary profile is made up of continuous sections, each one being solution of ( 21), separated by discontinuities (shocks).Let us observe that, even in case of a single shock, equation ( 23) alone is unable to determine both resulting unknowns (namely current and shock position) simultaneously.However, it can be seen that in this situation the system dynamics determines the value of the current according to an extremal principle (maximum or possibly minimum), analogous to the one stated in [25,26].Even with this extra condition, the shock position may not be completely determined (more so when there are more than one shocks), and in this case some shock-stability criterion must come into play.We will also observe that, when a particular relation among the parameters occurs, the shock positions remain anyway indeterminate, and that this physically corresponds to the onset of two delocalized (but synchronized) shocks, in analogy to what has been reported by Banerjee and Basu (in a model with rate modulation function characterized by two equivalent absolute minima) [16].
In the framework of the PA theory, the current-density relation for the model considered here has already been obtained, in explicit analytical form, in various previous works [30,31,38] (note that most of such papers actually deal with homogeneous rates, but this fact turns out to be irrelevant, due to the "smoothness" assumption).Here however we report (in Appendix A) a slightly more general derivation, which holds even in the absence of the symmetry assumption (3) (albeit we do not analyze this case further).As mentioned in the introduction, the physically relevant fact is that the current-vs-density diagram undergoes a transition from a unimodal shape, for attractive or weakly repulsive interactions, to a bimodal one, for stronger repulsive interactions.Taking v and q as free parameters, the transition line between the two regimes, shown in figure 2, is defined by the following equation Also the latter equation has already appeared in the literature §, but for the sake of completeness we report a derivation of it in Appendix B. In the same figure 2 we report lines representing the DME and KLS conditions, respectively expressed by equations (6a) and (8a).By comparing the latter equation with (24) and eliminating q we obtain the transition value for the interaction parameter, namely v = 1/3.As anticipated in the introduction, both aforementioned lines turn out to cross the transition line, and one can also observe that they are tangent to each other at point q = v = 1 (i.e. for vanishing NN interaction).From now on we will always assume that equations ( 8) hold, and that accordingly the steady-state dynamics of the system are specified only by the interaction parameter v, obviously besides the rate modulation ratio λ max /λ min .§ Equation ( 26) in [38] precisely corresponds to our equation ( 24), whereas equation ( 17) in [30] (seemingly quite different) actually includes some extra nonphysical solutions.In figure 3 we plot the (reduced) current-density relation J = F p,q,r,s (ρ), evaluated in the hypotheses described above, namely (3) and ( 8), for a given value of the interaction parameter v in the bimodality region.The explicit expression, parameterized by v only: can be obtained by plugging (A.10), (A.16) and (A.20) into (D.1), of course along with the KLS condition (8a).As previously mentioned, we can see that such relation is symmetric under the ρ → 1 − ρ transformation (mirror symmetry), as a consequence of (3), and is characterized by a convex region around a "central minimum" at density ρ • = 1/2 and two "side maxima" at densities ρ ± * = 1/2 ± ∆.The function values corresponding to the (local) minimum and the (absolute) maxima, which we respectively denote by J • and J * , are derived analytically in Appendix C. Plugging the KLS condition (8a) into (C.3) and (C.5), we respectively get We will see that the knowledge of these values allows one to determine, analytically, the transition between the two different regimes emerging in the strong interaction region.Figure 3 also illustrates the meaning of equation ( 21), which, given the reduced current J(x), allows one to determine the possible density profiles, by inverting the current-density function (details in Appendix D).It can be seen that, apart from degenerate cases, there can be either 2 or 4 possible density solutions, respectively for 0 < J(x) < J • or J • < J(x) < J * .Due to the mirror symmetry, they turn out to be pairwise complementary, that is, according to the notation in figure 3,

Density profiles
In this section we first describe the results that we have obtained at finite size, by the numerical PA method, introduced in the previous one.As mentioned above, at the size considered (L = 10000) the smooth sections of the density profiles turn out to be practically indistinguishable from the analytical solutions, obtained in the continuum limit.We consider two different combinations for the relevant model parameters, representative of the two different regimes that we have observed.In particular, keeping the rate modulation ratio fixed at λ max /λ min = 1.5, we take two different values of the interaction parameter, namely v = 0.15 and v = 0.10 (recall that smaller v means stronger repulsion).We respectively denote the two regimes as multi-shock (MS) and small-shock (SS), for reasons that will be immediately clear from the description.The MS regime is illustrated in figure 4. For low enough values of the mean density ρ, the system exhibits a smooth profile, with a density maximum at position x min = 1/2, corresponding to the minimum of the rate modulation function (2).We denote this phase as the LD (low-density) phase.Let us now imagine to increase ρ.At a certain value (roughly ρ ≈ 0.20), the steady-state profile begins to develop a shock, initially placed at x min , with vanishing amplitude, and progressively moving backward (i.e.towards smaller x values, opposite to the particle flux), with increasing amplitude.We denote the latter phase, characterized by 1 shock, as S1L (where L stands for low density).Such a displacement of the shock position for increasing ρ goes on until it reaches a peculiar position, say x • .After that, a further increase of ρ no longer affects this first shock, but it induces the onset of a second shock, initially placed at x • (with vanishing amplitude) and progressively moving forward (with increasing amplitude).The transition to this new phase, characterized by 2 shocks (and denoted as S2L), occurs roughly at ρ ≈ 0.27.Once again, the displacement of the second shock goes on until it reaches another peculiar position, in this case coinciding with x min .Upon further increasing ρ, even this shock remains "locked", while a third shock sets on, initially placed at x min (with vanishing amplitude) and subsequently moving backward (with increasing amplitude).
The transition to the latter 3-shock phase (denoted as S3L) occurs roughly at ρ ≈ 0.34.Furthermore, around ρ ≈ 0.41, the third shock joins the first one at x • , giving rise to a single shock (with double amplitude), which still moves backward for increasing ρ.Note that the shock placed at x min is still there, so that we have another 2-shock phase (denoted as S2), which extends over a range of average densities including the halffilling value ρ = 0.5.Needless to remark that all observed shocks are properties of the steady state, so what we call "displacements" do not happen over time but are actually changes in the steady state of the system, in response to changes in its mean density ρ.For ρ > 0.5 the behaviour of the system can be argued from the particle-hole symmetry.
In particular we can see that transforming the average density as ρ → 1 − ρ entails the following transformation for the density profiles As a consequence, it is possible to identify the high-density counterparts of the various shock phases S1L, S2L, S3L, which we respectively denote as S1H, S2H, S3H.Of course, also a smooth high-density phase exists, which we denote as HD.As far as the current is concerned, it can be observed that it keeps a constant maximum value in all the shock phases, and that the value computed at finite size is practically indistinguishable from the continuum-limit value where we recall that J * denotes the maximum of the reduced current-density relation (figure 3).Looking at figure 3, one can also argue that the x • position, introduced above, is the one where the intermediate-density solutions become degenerate, namely , and the reduced current takes value J(x • ) = J • (with a slight abuse of language, we will call this the critical point).As a consequence, from equations ( 22) and ( 30), with J = J max , we can write the equation for x • as where we recall that J • and J * are known analytically from ( 26) and ( 27) (see also Appendix C).Taking into account the specific form (2) of the rate modulation function λ(x), we obtain Let us note, in figure 4, that the shock we claim to be placed at x • , in the S2L and S3L phases, is in fact slightly displaced to the left, and the same goes for the S3L/S2 transition, where the two-shock merging appears to occur at x slightly lower than x • .Even though we do not display the results, we have a quite clear evidence that both such discrepancies are finite-size effects, since they tend to vanish upon increasing the number of nodes L. In fact, it is not surprising that finite size effects are especially relevant in the neighbourhood of the critical point, where the derivative of ρ(x) diverges.
Let us now observe that, if the rate modulation ratio is not large enough, specifically then equation (31) has no more solutions, meaning that no critical point can exist.Keeping the left-hand side fixed, the above condition can also occur if the right-hand side becomes exceedingly large, as a consequence of a stronger repulsive interaction.
In both cases, a very different transition scenario takes place (SS regime), which we illustrate in figure 5.For low values of the mean density ρ, the system still exhibits a smooth profile (LD phase), which undergoes a transition to the 1-shock phase (S1L) upon increasing ρ (the transition occurs roughly at ρ ≈ 0.19).As in the MS regime, the shock is initially placed at x min = 1/2, with vanishing amplitude, and, upon further increasing ρ, it progressively moves backward with increasing amplitude.Here begins the difference with respect to the MS regime.Due to the lack of a critical point, the shock is allowed to go through the entire system, so that in figure 5 we see it exiting our "observation window" at x = 0 and reentering at x = 1.At this point, corresponding to the maximum of the rate modulation function (2), the shock takes maximum amplitude, and subsequently proceeds backward, with decreasing amplitude.For ρ ≈ 0.40, it finally comes back to x min and disappears, giving rise to another phase featuring a smooth profile, which we denote as LD'.In the whole S1L phase the current takes the constant maximal value (30) (apart from the usual vanishing discrepancy, due to finite size), whereas in the LD' phase the current turns out to decrease upon increasing ρ.Note that the latter is a typical feature of a high-density phase, even though in this case it takes place in a low-density range, i.e. ρ < 0.5.This is clearly related to the presence of side maxima (specifically the one at lower density) in the current-density relation.Upon further increasing the mean density, roughly above ρ ≈ 0.46, we can observe the onset of another shock, which in this case originates at x = 0, that is, at the maximum of the rate modulation function.The shock moves forward on increasing ρ, it goes trough the whole system, with maximum amplitude at x min , and finally disappears at its initial position (x = 1 in our view).We denote this last 1-shock phase as S1', in order to distinguish it from the "normal" 1-shock phase (S1), occurring when the current-density relation is unimodal (that is, in case of weak [38] or lacking [16] NN interaction).A first peculiar feature of the S1' phase is the descending shock (i.e.going from higher to lower density), which does not appear in any other 1-shock phase (including S1L and S1H).Another one is that the current is constant (i.e.independent of ρ) but minimal, specifically taking value where we recall that J • denotes the minimum of the current-density relation (figure 3).
In fact, the onset of a minimal-current phase is related to the presence of the central minimum in the current-density relation, as recognized for instance in [25,26], and it has actually been observed in several similar models.For ρ > 0.5 the behaviour of the system is still constrained by the particle-hole symmetry, so that, as done in the MS regime, it is possible to identify the high-density counterparts of the phases encountered at ρ < 0.5, which we respectively denote as HD, S1H, HD' (the S1' phase extends roughly up to ρ ≈ 0.54).The symmetry also entails that in the HD' phase the current increases upon increasing ρ, which is actually a typical feature of a low-density phase.
We conclude this section by discussing the fact that, as already mentioned in section 3, in principle it is possible to completely determine the density profiles in the continuum limit, even without making use of the numerical PA method, but using only the analytical solutions of equation (21).Figures 4 and 5 clearly confirm that, in the shock phases, the smooth sections of the density profiles always match with the analytical profiles evaluated at maximal or minimal current.In all the observed phases, including the MS regime, there is always only one shock whose position, say x s , depends on the mean density of the system, while any other shock is "locked".To determine x s , we need to solve equation ( 23) (the integral must be done numerically), for a suitable profile ρ(x), being the union of known continuous-profile sections, with a discontinuity at x s .We refer to section 6 for more details about the specific equations for the various phases.Here we just discuss the fact that such equations generally have two possible solutions, and how the unphysical one can be detected.In our case, the matter is specially simple, since the smooth profiles are characterized by a mirror symmetry with respect to x min = 1/2 (following from the same property of the rate modulation function), which is broken by shock profiles.As a consequence, if ρ(x) is a solution of equation (23), the "mirrored" profile ρ(1 − x) is also a solution, and therefore a shock placed in x s in the former profile, is shifted to 1 − x s in the latter.Furthermore, if the shock is "ascending" in the former profile, it becomes "descending" in the latter (and vice versa), and consequently one is stable and the other unstable, so that ultimately only one of the two profiles can be the stationary one.Let us recall that the stability of a shock depends on the propagation velocity of small perturbations of the density profile, i.e. on the derivative of the current-density relation (the so-called collective velocity) [25,26,42].In particular, stability occurs if the propagations on the two sides of the shock are opposite and directed towards the shock itself (see in particular equation ( 9) in [26]).By comparing figures 4 and 5 with the current-density relation in figure 3, it can be verified that the profiles obtained by the numerical PA method always satisfy the stability criterion.The aforementioned "locked shocks", observed in 2-and 3-shock phases, are partly an exception, and can be considered marginally stable, as they involve at least one point of the current-density function with zero derivative.In any case, it can be seen that the numerical PA method is convenient also in this respect, since, by following the dynamical evolution of the system, it "automatically" determines the stable stationary state.

Phase diagrams
In this section we summarize some results which have already been outlined in the previous one, but which can be described more precisely and exhaustively in terms of phase diagrams.First, we have seen that the interesting behaviours (i.e. the two regimes tagged MS and SS) occur when the current-density relation is bimodal, for a sufficiently strong repulsive interaction (i.e. for small enough values of the interaction parameter: v < 1/3).For v > 1/3 the phenomenology is qualitatively equivalent to that of a system without NN interaction, where a single shock phase (S1) occurs [16].We call this the largeshock (LS) regime, since the (unique) shock, occurring between a low-density (LD) and a high-density (HD) phase, can assume greater amplitudes than those observed in the MS and SS regimes.The transition between the latter two regimes is also controlled by another parameter, namely the rate-modulation ratio λ max /λ min .The transition line is determined by (33) taken as an equality, where the right-hand side depends on v according to (26) and (27) (see also Appendix C), thence The resulting phase diagram is shown in figure 6, where for clarity we have also marked the specific points analyzed in the previous section.By varying the mean density ρ, in the previous section we pointed out two different sequences of phase transitions (characterizing the MS and SS regimes) and we approximately identified the ρ values at which the transitions take place.Taking into account the analytical profiles obtained from the continuum PA theory (section 3), it is possible to precisely determine these threshold values, and also their evolution, in response to arbitrary variations of the control parameters.Let us first consider the MS regime and in particular the density profiles occurring at the 4 transitions, observed in the previous section.We can see that, limited to the (x • , x min ) interval, such profiles correspond respectively to the 4 possible analytical solutions, which we have denoted as ρ − 1 (x), ρ − 2 (x), ρ + 2 (x), ρ + 1 (x) (in increasing order of density).Moreover, in the remaining intervals [0, x • ) and (x min , 1) (which, taking into account the periodicity of the system, are effectively equivalent to a single interval (x min − 1, x • )), all the transition profiles match with the lowest-density solution ρ − 1 (x).Thus, generically denoting with ρA/B the mean density at which the transition between phases A and B takes place, we can write where an appropriate subscript reminds us that the profiles are always calculated at maximum current.As far as the SS regime is concerned, it can be seen that the profiles occuring at the 3 observed transitions correspond entirely to single analytical solutions, which however can be at either maximum or minimum current.In particular, with the same notation used above, we have where we note that, still due to periodicity, (37a) is equivalent to (36a).All the above integrals can be solved numerically with negligible computational effort, just because the integrand functions are all known analytically.Keeping the rate-modulation ratio fixed and varying the interaction parameter, we obtain the phase diagram shown in figure 7, where again we have marked (thin dotted lines) the two cases studied in the previous section.Note that the diagram has been completed by symmetry, adding the complementary phases that can be observed in the high density range (ρ > 1/2).On the other hand, keeping the interaction parameter fixed and varying the rate-modulation ratio, we obtain the phase diagram shown in figure 8. Let us finally observe that, in both figures 7 and 8, almost all transitions (namely, all those denoted by solid lines) are continuous.In fact they are characterized by the appearance or disappearance of a In the first equation, splitting into two integrals is unnecessary, but it better highlights the analogy with subsequent ones.shock, whose amplitude does not exhibit any discontinuity, as a function of the control parameters.The only exception is the transition denoted by the dashed line, which occurs simultaneously with the crossover between the MS and SS regimes.We refer to the next section for a closer discussion of this "crossover regime".By now we only mention the fact that, in this peculiar situation, the critical point x • becomes degenerate with its complementary 1 − x • (we shall denote such phenomenon as coalescence), which entails that the shock occurring at that point in the S2L and S3L phases (and similarly in S2H and S3H) is no longer locked.As a result, the mean-density constraint equation ( 23) is no longer sufficient to determine the locations of all shocks.It is reasonable to expect that physically this corresponds to the onset of so-called delocalized and synchronized domain walls [16], i.e. two shocks whose positions are indeed stochastic processes (random walks), yet "rigidly" bound to each other because of equation ( 23).

Numerical simulations
The main purpose of this section is to provide quite robust evidence that, under the KLS condition, the PA theory in the continuum limit describes the behavior of the system in a very accurate (plausibly exact) manner, thus also confirming the accuracy (or even the exactness) of the phase diagrams shown in the previous section.To do this we are going to compare the density profiles obtained from kinetic Monte Carlo (KMC) simulations (with increasing sizes: L = 2000, 5000, 10000) with those obtained from the analytical theory (in the L → ∞ limit).We focus on the shock phases, since in the smooth phases an even better matching is obtained, at relatively small sizes, of the order of L = 500, 1000.The simulations are carried out using the well-established Gillespie algorithm.As usual, each simulation is first run for a "settling time" t set , to ensure that the system relaxes to steady state, after which averages are computed over an "averaging time" t ave .These characteristic times have been empirically adjusted, leading us to choose (at the largest system size considered) t set = 2 × 10 5 and t ave = 10 6 (the time unit is fixed by λ min = 1).
As far as the analytical profiles are concerned, we recall that just one shock position varies with the mean density, and that such position must be determined according to the scheme outlined at the end of section 4, i.e. taking into account both the density equation ( 23) and the stability criterion.As previously mentioned, the results of the numerical PA method help us to select one of two possible profile shapes, on which the only constraint equation ( 23) is left to be imposed.The equations can be conveniently written, making use of the mean-density transition values ( 36) and (37), and evaluating the corrections that occur at a given mean density ρ as a function of the shock position.In particular, for the different phases observed in the MS regime we have x min x min x min where we see that the generic unknown x A , denoting the shock position in a given phase A, appears solely as lower bound of the integration interval.Regarding the SS regime, in the S1L phase the equation is equivalent to (38a), whereas in the S1' phase we have The results, relating to the MS and SS regimes, are shown respectively in figures 9 and 10.First of all, one observes that, as previously mentioned, the smooth sections of the simulated profiles match perfectly with the analytical solutions.Also, in the shock phases, two considerably different situations are observed, depending on whether only one shock or more than one are there.In the former case, i.e. in phases S1L and S1' (figures 9 and 10), there is a clear convergence of the simulation results towards the analytical ones, as the system size increases.Conversely, in the latter case, i.e. in phases S2L, S3L and S2 (figure 9), although the expected trend is still quite evident, extremely relevant finite-size effects occur, such that, even at the maximum size considered, some shocks are still quite far from the theoretical prediction.As already noticed in section 4 about the numerical PA results, the critical point seems to play a key role in these effects.Indeed, especially in phases S2L and S3L, the largest discrepancy between theory and simulation is observed for the shock lying in the vicinity of this point, whereas the similar discrepancy for the nearby (rightward) shock can be interpreted just as a side effect of the former, mediated by the mean-density constraint (23).Needless to recall that, due to symmetry, all the above remarks could be analogously repeated for the respective high-density phases.Further plausible evidence of exactness for the PA theory is obtained by analyzing the so-called fundamental diagram [2], i.e. the current J as a function of the mean density ρ.In the framework of the theory, it is convenient to compute the latter as a function of the former, by integrating the appropriate analytical density profiles, evaluated at a given current J .In figure 11 we report the results, for the previously considered values of the parameters (representing the MS and SS regimes), and again we make a comparison with the corresponding simulations (due to the mirror symmetry, we only consider one half of the diagram, for ρ ∈ [0, 1/2]).We note that the increasing part of the theoretical curves corresponds to the LD phase, while the decreasing part corresponds to the LD' phase.The former occurs in both regimes, for J ∈ [0, J max ], whereas the latter occurs only in the SS regime, for J ∈ [J min , J max ].In the two cases, the mean density as a function of the current is respectively obtained as ρLD (J ) = Conversely, the plateau regions, where the current is independent of the mean density, correspond to the various shock phases (at maximal or minimal current).By the way we note that, evaluating (40) at J max or possibly J min , we recover the mean-density values at the transitions with the shock phases, previously calculated by (37), i.e.
In figure 11 these values are explicitly tagged.Regarding the simulations, we have studied just a few values of the mean density in a very accurate way.It can be seen that in the smooth phases the agreement is excellent already at the smallest size considered (L = 2000), such that the deviation with respect to the theoretical prediction is not even detectable at the scale of the figure.In the shock phases there is a small residual discrepancy, which however shows a clear tendency to reduce for increasing size.
In the remainder of this section, we try to give a somewhat more detailed description of what happens at the transition between the MS and SS regimes, which we previously called the crossover regime.From the viewpoint of the continuum limit (in the PA theory), we have already noticed that in this situation the two so-called critical points x • and 1 − x • , occurring in the MS regime, degenerate into a single "coalescence point" (in the graphs we see this point split at the left and right edges: x = 0, 1).As a result, the shock placed at the coalescence point in the S2L and S3L phases is no longer locked, which gives rise to an extra degree of freedom in the shock positions.It can therefore be expected that the S2L and S3L phases in the crossover regime are characterized by 2 shocks whose positions vary over time with a random motion (plausibly with a slower kinetic than that of particle hopping), but with a constraint on the relative position imposed by equation ( 23) (that is -physically speaking -by particle conservation).As already mentioned above, we speak in this case of delocalized and synchronized shocks [16].In order to visualize this phenomenon, we carry out KMC simulations in which, after waiting the usual settling time, we average over a certain number of (disjoint) successive time windows of duration t ave = 10 5 .Figure 12 shows some results, corresponding to the usual rate modulation ratio value λ max /λ min = 1.5 (so one can still refer to the phase diagram in figure 7).The expected phenomenon is clearly visible, that is, by superimposing density profiles computed on different time windows, shocks are observed in different positions.In this regard, it should be noted that the choice of time windows shorter than usual does not only have the role of limiting the required computational effort, but also that of better highlighting the shock dynamics itself (too short a window would result in a noisy profile, whereas a very long window would Crossover regime for λ max /λ min = 1.5 and different mean-density values: ρ = 0.25 (S1L phase), 0.35 (S1L-S2L crossover), 0.40 (S1L-S3L crossover), 0.45 (S1'-S3L crossover).Solid lines denote KMC simulation results for L = 10000 and v = 0.1405, 0.1410, 0.1415, 0.1420 (only for ρ = 0.40, 0.45), 0.1425 (only for ρ = 0.45), averaging over different time windows (see the text).Dotted lines denote all solutions of (21) at v ≈ 0.1392 (theoretical crossover value) and maximal current J max = λ min J * .
show a smooth profile, without shocks).Of course, the synchronization effect cannot be appreciated from the figure, but it can actually be verified, albeit quite roughly, by examining the different profiles separately.
Let us now describe in detail the four cases displayed in figure 12.At low meandensity values (specifically at ρ = 0.25), it can be seen that the density profile never reaches the coalescence point, so the system is completely unaffected by crossover.The profile is characterized by a single very stable shock, qualitatively equivalent to that of the S1L phase (see figures 9 and 10).Actually, in the phase diagram of figure 7, the point corresponding to this case is located inside the S1L phase region.At higher mean-density values (specifically at ρ = 0.35), we begin to observe 2 shocks (which for simplicity we will call "upper" and "lower", relating to density values) and the delocalization phenomenon, that is, shock positions vary, within a certain range, depending on the time window considered.At one end of the range, the lower shock is positioned at the coalescence point (x = 0, 1) and the upper one at some position x ∈ (0, 1/2).The resulting profile is therefore completely analogous to that of the S2L phase (see figure 9).However, as previously mentioned, in this crossover situation the lower shock is no longer locked, i.e. it can move backwards (in the figure we see it reentering from the opposite side, at x = 1), up to an extreme position x ∈ (1/2, 1), which corresponds to the upper shock getting back to the coalescence point (x = 0), with vanishing amplitude.This opposite extreme situation is therefore equivalent to the S1L phase.In the phase diagram in figure 7, the corresponding point is located along the S1L-S2L section of the crossover line.At even higher mean-density values (ρ = 0.40, 0.45) one can observe up to 3 shocks (say "upper", "intermediate", and "lower").For both cases, one of the extremes of the delocalization range corresponds to the lower shock positioned at the coalescence point, with an overall profile equivalent to that of the S3L phase (see figure 9).As already seen above, however, the lower shock can move backwards, and this induces (due to the usual conservation constraint) a forward displacement of the upper shock, while, at least initially, the intermediate-shock position remains locked at x = 1/2.If the displacement of the lower shock goes beyond a certain limit, the upper shock can come to vanish at the x = 1/2 position, after which the intermediate shock gets "unlocked" and can proceed backwards.In this situation the profile is equivalent to that of the S2L phase.Now, if the mean density is not too high (e.g.ρ = 0.40), the intermediate shock can even reach the coalescence point (x = 0), with vanishing amplitude, when the lower shock is still in a position x ∈ (1/2, 1), so that the extremal profile is still that of the S1L phase (in the diagram of figure 7 we are along the S1L-S3L section of the crossover line).If, on the other hand, the mean density exceeds a certain limit (e.g.ρ = 0.45), it can be seen that the intermediate shock can never reach the coalescence point, but at most a position x ∈ (0, 1/2), with nonzero amplitude.This corresponds to the fact that the lower shock proceeds backward until it vanishes at position x = 1/2, so that the resulting extremal profile is that of the S1' phase (see figure 10).Still with reference to figure 7, in this last case we are indeed along the S1'-S3L section of the crossover line.In general, we can say that the situations illustrated in figure 12 exhaust the types of delocalized shocks, that can be observed in this model.Another general fact, which is worth noting, is that the whole crossover phenomenology, described above, is observed in the simulations over a (narrow) interval of interaction parameter values, not at a single special value, as predicted by the theory in the continuum limit.Furthermore, the interval does not even cover the theoretical crossover value (v ≈ 0.1392), determined by equation (35).In fact, the simulations carried out at the theoretical v value still exhibit extremely stable shocks, qualitatively indistinguishable from those of the SS regime, shown in figure 10.We have reason to believe that this discrepancy too, like others observed previously, is a finite-size effect, since it can be verified that both the amplitude of the aforementioned interval and the distance from the theoretical value, which are anyway small, also show a clear decreasing tendency, upon increasing size.

Conclusions
In this paper we have studied a TASEP with periodic boundary conditions, hopping rates characterized by a smooth spatial modulation (with a unique global minimum), and a NN repulsive interaction.The model is in principle analogous to that considered in a previous study [38], yet with the interaction parameters characterized by a different constraint, denoted as the KLS condition.Such a constraint makes the model more amenable to analytical treatment, specifically by means of the so-called PA theory (a generalized mean-field theory, based on a NN pair cluster).Actually, we conjecture that, in combination with the KLS condition, the analytical theory is asymptotically exact, except in the vicinity of shocks.Although we cannot rigorously prove this statement, we provide considerable numerical evidences.
Our main contribution is to have elucidated the nature of so-called S phases, namely, certain non-equilibrium steady states, that the cited previous study [38] had detected as qualitatively distinguished by the more "universal" smooth and shock phases [16].In particular, we have pointed out that such phases exhibit different features, depending on whether one considers intermediate or strong repulsion regimes.In the intermediate regime, it turns out that the S phases can be regarded as subphases of the maximalcurrent phase, and can also be distinguished into several further subphases, with density profiles displaying up to 3 shocks (MS regime).On the other hand, in the strong regime it turns out that the S phases can be either maximal-or minimal-current phases, but always displaying just 1 shock (SS regime).Furthermore, we have pointed out the existence of a crossover regime, emerging in between the aforementioned regimes and featuring delocalized and syncronized shocks [16].The analytical PA theory has allowed us to work out the whole phase diagram of the model in great detail.As a consequence of the above conjecture, we claim that such phase diagram may be exact as well, under the KLS assumption, just qualitatively correct otherwise.In fact, we have argued that the whole, considerably rich phenomenology emerging in this model can be traced back to the onset of a bimodal current-density relation, so that the choice of the KLS condition is relevant only at a quantitative level.Thus, we claim that our results can describe in a qualitatively correct way the physics of a whole class of similar models, differing only by the absence of the KLS condition, and including in particular the one studied in [38].In the course of our work we have indeed collected several evidences of such a claim, that have not been reported.Here, in figure 13, we report just one significant example, relating to the analogous model with DME condition (i.e. the one considered in [38]), in the intermediate-repulsion (MS) regime.We can see in particular that the PA theory quantitatively fails, especially in predicting the location of the critical points.However, the sequence of density profiles, upon incresing the mean particle density, turns out to be qualitatively equivalent to that obtained with the KLS condition in the same regime (figure 9).
Let us finally note that the present work naturally opens the way to the investigation of a number of related models, for which the analytical theory developed here could be applied without relevant changes.Among such extensions, which are of course beyond the scope of this paper, we first have in mind the very same model with open boundary conditions ( [17] considers a closely related model, though without NN interaction).We also include the possibility of relaxing some of the restrictive assumptions about hopping rates, in particular (1) and (3).Relaxing (1) means to consider that the interaction itself may feature a spatial modulation, whereas relaxing (3) amounts to remove the mirror symmetry of the current-density relation.The former issue is mainly methodological (specifically to test whether the PA theory is still reliable in such a case), the latter is meant to make the model slightly more realistic as a model for vehicular traffic [2,32].It would also be interesting to carry out a deeper investigation about the random motion of shocks in the crossover regime, possibly developing an ad-hoc theory based on a Fokker-Planck equation, along the lines of [16] and [43].Furthermore, one could investigate to what extent the observed phenomena can be affected by any concurrent processes, such as Langmuir kinetics [21,35] and/or local stochastic resetting [44,45,46,47].We are considering these types of questions as possible topics for future work.The latter equality follows from (A.10), which states that a given (admissible) standardcurrent value has two possible corresponding densities, with unit sum, precisely Now, in terms of the reciprocal variable 1/η, the correlator-density relation (A.13) reads clearly analogous to (A.10).Consequently, also this equation has in principle two solutions for 1/η, with unit sum, but we easily argue that only the larger one satisfies (A.17), with equality occurring precisely for v = 0.For more clarity, in figure A1 we display the admissible region in the plane I vs 1/η, along with some instances of the correlator-density relation (A.19), obtained for different values of v.We note in particular that, excluding the borderline case I = 0, η = 1, the regions corresponding to v < 1 or v > 1 are respectively characterized by η > 1 or η < 1. Remembering the correlator definition (A.9), this physically means that we have, respectively, a larger or smaller probability to find an empty node, if we know that one of its neighbors is occupied.These are clearly the fingerprints of a repulsive or attractive interaction, respectively.Some results from KMC simulations, reported as well in figure A1, confirm that the PA theory in the continuum limit provides extremely accurate (possibly exact) results in regions characterized by smooth profiles.On the other hand, a definite deviation corresponds to a shock region, where as usual the results are highly sizedependent.By the way, we also note that in the shock region the correlator η can become smaller than 1, even in case of repulsive interaction, and that the effect is more pronounced for larger system size.

Appendix B. Transition to bimodality
Making use of the continuum PA theory, in this appendix we determine the transition where, as described in section 3, the current-density relation becomes bimodal, with the central minimum at density ρ • = 1/2 and the two side maxima at densities ρ * = 1/2 ± ∆.We understand that this calculation refers to a "particle-hole-symmetric" case, where (3) holds, and the appropriate expression for the (reduced) current is (A.where we see that, apart from parameters, the current is expressed as a function of the only correlator variable η.We also recall that, as stated by (A.20), η does not depend on ρ explicitly, but only through the standard current I.As a consequence, the derivative can be evaluated as We can now make the following observations.In order for the extra stationary points to be placed in the physically meaningful interval (0 < ρ * < 1), also being distinct from the "normal" stationary point (ρ * = 1/2), the required condition in terms of standard current is, according to (A.10), 0 < I * < 1/4.The latter condition corresponds, through equation (A.20), to different conditions for the correlator, depending on whether v < 1 or v > 1.Altogether, we can write Furthermore, analyzing equation (B.9), one can argue that (B.10) maps to respective conditions for r, namely Let us observe that the right-hand side of the above inequalities turns out to be either larger than 1 (for v < 1) or negative (for v > 1).As a consequence, in both cases (B.11) inherently satisfies the realness condition for the correlator expression (B.9), which is of course 1/r < 1.Now, plugging (A.16a) into (B.11),we can rephrase the conditions in terms of the basic parameters v and q, obtaining 1 − q where we immediately note that the (1 − v) term can be simplified.Since this term is respectively positive for v < 1 or negative for v > 1, in the latter case (only) the inequality sign gets reversed, so that in the end we have a unique inequality, which can be written as Following the above derivation, one can also argue that (B.13) with an equality sign corresponds to the borderline case I * = 1/4, ρ * = 1/2, which represents, in the currentdensity relation, a degeneration of the side maxima into the "normal" central maximum.Such a behavior makes reason to regard this phenomenon as a sort of second-order phase transition.Let us finally remark that in principle we have correctly taken into account both possibilities v < 1 or v > 1, but in the end we deduce that, for physically meaningful (i.e.positive) q values, (B.13) can be verified only for v < 1, that is for repulsive interactions.By means of (B.11), we have previously argued that this case also implies r > 1.This last remark will come in handy in the next section.
immediately obtained from (C.3) and (C.5).Note that the right-hand side of equation (C.7) can easily be verified to take value 1 at the bimodality transition (24).where we have defined a symbol for the first-order coefficient, explicitly denoting the dependence on J, that is ψ(J) a + J/r − (1 − v 2 )(a − J/r) 2 .(D.4) Note that in (D.3) the dependence on the model parameters is fully incorporated into the ψ(J) function.In particular, the fact that the zeroth-order coefficient can be written as J ψ(0) can easily be argued from (D.4), taking into account the parameter identity (B.3).For completeness, we can also determine ψ(J) in terms of only the elementary parameters v and q, by plugging into (D.4) the appropriate expressions for r and a, namely equations (A.16).After some algebra we obtain Figure 1.Three regimes observed in[38], as a function of the (mean) particle density: (a) no interaction or weak interaction; (b) mildly strong interaction; (c) very strong interaction.Phase tags are explained in the text.

Figure 2 .Figure 3 .
Figure 2. Parameter plane (q vs v): the shaded region is the one characterized by a bimodal current-density relation.The thin dashed and solid lines respectively denote the DME and KLS conditions (see the text).

Figure 4 .
Figure 4. MS regime: λ max /λ min = 1.5, v = 0.15.Different line types denote density profiles, computed by the PA method at finite but large size (L = 10000), for different values of the mean density ρ (see the legends).Thin dotted lines denote the analytical solutions (21) at maximal current J max = λ min J * .

Figure 6 .
Figure 6.Reciprocal rate-modulation ratio (λ min /λ max ) vs interaction parameter (v) phase diagram.Tags denote different regimes, described in the text.The transition line between the MS and SS regimes (dashed line) is defined by equation(35).The shaded region is the one characterized by a bimodal current-density relation.Cross symbols mark the particular cases analyzed in the previous section (λ min /λ max = 2/3, v = 0.10, 0.15).

Figure 7 .
Figure 7.Phase diagram at constant rate-modulation ratio (λ max /λ min = 1.5): mean density ρ vs interaction parameter v. Solid lines denote continuous transitions.A dashed line denotes the crossover regime between MS and SS (see the text).Phase tags are also explained in the text.Thin dotted lines mark the particular cases analyzed in the previous section (v = 0.10, 0.15).
Figure A1.The admissible region in the plane I vs 1/η is delimited by a thin dash-dotted line.Thin solid lines represent the correlator-density (A.19), for different values of v (the analytical continuations outside the admissible region are denoted by dotted lines).Thicker lines denote KMC simulation results for v = 0.4, λ max /λ min = 1.5, ρ = 0.5 and L = 2000 (dashed line), 5000 (dash-dotted line), 10000 (solid line).

( i )
As discussed in Appendix A, the η vs I relation (A.20) turns out to be strictly monotonic (excluding the trivial case v = 1), so we always have dη dI = 0 .(B.6) (ii) From the standard-current definition (A.10), we have dI dρ = 1 − 2ρ , (B.7) which entails that ρ • = 1/2 is always a stationary point for the current-density relation.(iii) From equation (B.4) we immediately get dJ dη = r/η 2 − (r − 1) 1 − v 2 , (B.8) from which we argue that an extra stationary point may appear for η = η * .19)allows us to obtain the corresponding standard-current value I * , which in turn corresponds to 2 possible densities ρ * , determined by (A.18).
Appendix D. Inverting the current-density relationIn this last appendix we show how to derive the standard current I as a function of the reduced current J and of the model parameters, once again in the framework of the continuum PA theory.The resulting equation, in combination with (22) and (A.18), allows us to determine density profiles analytically, at given values of the physical current J .The calculation still refers to a particle-hole-symmetric case, for which the appropriate expression of the current is (B.1), equivalent toJ = r a − η (a − I) .(D.1)Assuming to fix J, the above equation represents a linear relation between the standard current I and the reciprocal correlator 1/η, that is1 η = a − I a − J/r .(D.2)Another equation relating the above two quantities is provided by (A.19), the so-called correlator-density relation.Plugging the former into the latter, we immediately obtain a quadratic equation for I, whose coefficients depend on the reduced current J, besides model parameters.This equation can be written as I 2 − I ψ(J) + J ψ(0) = 0 , (D.3)