Dynamical Transitions in a One-Dimensional Katz–Lebowitz–Spohn Model

Dynamical transitions, already found in the high- and low-density phases of the Totally Asymmetric Simple Exclusion Process and a couple of its generalizations, are singularities in the rate of relaxation towards the Non-Equilibrium Stationary State (NESS), which do not correspond to any transition in the NESS itself. We investigate dynamical transitions in the one-dimensional Katz–Lebowitz–Spohn model, a further generalization of the Totally Asymmetric Simple Exclusion Process where the hopping rate depends on the occupation state of the 2 nodes adjacent to the nodes affected by the hop. Following previous work, we choose Glauber rates and bulk-adapted boundary conditions. In particular, we consider a value of the repulsion which parameterizes the Glauber rates such that the fundamental diagram of the model exhibits 2 maxima and a minimum, and the NESS phase diagram is especially rich. We provide evidence, based on pair approximation, domain wall theory and exact finite size results, that dynamical transitions also occur in the one-dimensional Katz–Lebowitz–Spohn model, and discuss 2 new phenomena which are peculiar to this model.


Introduction
The Totally Asymmetric Simple Exclusion Process (TASEP, see References [1,2] for recent reviews), originally inspired by biological traffic problems [3,4], has now gained a prominent role in non-equilibrium statistical physics, becoming a paradigmatic minimal model, with a status analogous to that of the Ising model in equilibrium statistical physics. It can be simply generalized by including additional processes and/or interactions, obtaining models with a small number of parameters that exhibit rich stationary state phase diagrams, as generalized Ising models do in the equilibrium case.
The TASEP is defined on a one-dimensional lattice, whose nodes can be empty or occupied by one particle. Particles can hop only in one direction and the physics depends strongly on the boundary conditions. The most interesting situation is obtained with open boundary conditions [5], when the lattice is connected to particle reservoirs at both ends. The non-equilibrium stationary state (NESS) of the model is then determined by the particle densities of the boundary reservoirs.
Exact results are available for the NESS phase diagram [6][7][8][9], which can be rationalized in terms of the so-called theory of boundary-induced phase transitions [10], based on certain maximum and minimum useful in our investigation of the AS model [19]. Finally, the extrapolation of exact finite size results was first used in Reference [24], with very accurate results which could have led to an early discovery of the dynamical transition. Again, we found this technique useful in our study of the AS model [19].

Model and Methods
The 1D KLS model is defined on a linear chain of N nodes, labelled by i = 1, . . . , N. Each node can be empty or occupied by a single particle. We introduce occupation number variables n t i taking value 0 (respectively 1) if node i is empty (resp. occupied) at time t. A particle can hop in one direction only (say rightward, or from node i to node i + 1), provided the destination node is empty. The hopping rate from node i to node i + 1 depends on the configuration of the nodes i − 1 and i + 2 and will be denoted by Γ i,i+1 (n i−1 , n i+2 ). Following References [12,13] we shall use the Glauber rates It can be shown [13] that, with this choice of the hopping rates the NESS probability distribution in the bulk (e.g., in a model with periodic boundary conditions, with a given particle density which will be independent of time) is the equilibrium Boltzmann distribution of a 1D lattice gas model with repulsive (for V > 0) nearest-neighbour interaction V and thermal energy k B T = β −1 = 1, that is with a dimensionless Hamiltonian where µ is a chemical potential, conjugate to the particle density, and the time index in the occupation number variables has been dropped, since we are dealing with a steady state. A slightly more general result actually holds [13]: the NESS is described by the above Hamiltonian if the hopping rates satisfy Γ(0, 1) + Γ(1, 0) = Γ(0, 0) + Γ(1, 1) .
These conditions are clearly satisfied by the Glauber rates Equation (1), but other choices are possible, for instance the AS model, where Γ(0, n i+2 ) = Γ(1, n i+2 ). It is important to stress here that the NESS distribution does not satisfy detailed balance with respect to the dynamics of the 1D KLS model, so the NESS is truly out of equilibrium. The relationship between the bulk NESS and the 1D lattice gas has important consequences, since the equilibrium distribution of a 1D lattice gas (or Ising) model factors in a simple way into a product of local marginals, and the PA becomes exact [25]. In order to be more specific it is useful to introduce some notation for marginal distributions and expectation values. In the following we will denote the marginal distribution for a cluster of nodes starting at i, at time t, by We will also denote the local densities and correlations by respectively. In the NESS we will drop the time index (e.g., the local density will be ρ i ) and in a bulk NESS we will also drop the position index (e.g., the bulk density will be ρ). The local densities and correlations can be used to parameterize the 1-node and 2-node marginals as follows The factorization property of the 1D lattice gas, and hence of a bulk NESS of our model, can now be expressed by writing that k-node marginals (k ≥ 3) factor according to that is into a product of 2-node and 1-node marginals (see e.g., Reference [25] for a derivation of this property for an equilibrium 1D model with nearest-neighbour interactions, which includes the Hamiltonian Equation (2)). The PA [18,19], as well as the MCAK [12,13] and related techniques [26][27][28][29][30], is based on the assumption that the factorization Equation (14) holds at any time, and this explains why these techniques can give exact results for a bulk NESS.
In terms of entropy, we can introduce 1-node and 2-node cluster entropies (in units such that k B = 1) Of course, cluster entropies can be defined analogously for clusters of any length. In particular, for a k-node cluster, Equation (14) implies that where we have also introduced the 2-node entropy cumulants [25,31] A simple consequence of Equation (17) is that the 3-node cumulants [25,31] S l,l+1,l+2 = S l,l+1,l+2 − S l,l+1 − S l+1,l+2 + S l+1 , together with all the higher order ones, vanish in a bulk NESS of our model. Considering the PA (and other techniques sharing the same factoring assumption), we can say that it approximates the kinetics by neglecting 3-node (and higher order) entropy cumulants.
Using the factorization Equation (14) we can write kinetic equations for 1-node and 2-node marginals, that is for ρ t i and φ t i , at the PA level. Introducing the probability current and marginalizing the full master equation of the model, one obtainṡ The above equations, derived in the Appendix A, are exact but not closed in that they depend on 4-node marginals. The PA is indeed a possible closure scheme, using the factorization Equation (14) as an approximation: Before considering the kinetics according to the above equations, it is necessary to discuss bulk solutions for the NESS in some detail, and then specify the boundary conditions. In order to find a bulk NESS solution, that is a solution being invariant with respect to both time and position, we In order to reduce the above equation to its simplest possible form we introduce the bulk correlator We can then write the 1-node and 2-node marginals as Equation (25) can then be rewritten as 1 which establishes the relationship between ρ and η characterizing the bulk steady state, and in particular it allows us to determine η as a function of ρ as We can now compute the bulk current from Equation (20). Using Equations (25)-(31) we obtain where we have also defined With the Glauber rates Equation (1) we obtain a = 1/2 and b = 0 for the present model. Equation (34) can also be used for the AS model, in which case one obtains a = 0 and b = 1. For pure TASEP instead one has a = b = 0 and η = 1.
The bulk current-density relation (fundamental diagram) Equation (34) is illustrated in Figure 1 as a function of both the repulsion V and the density ρ. It must be stressed that this result has already been obtained in References [12,13] (using the MCAK), where the current was also shown to exhibit a single maximum at density 1/2 for V < V * = 2 ln 3 and 2 maxima, symmetric with respect to a minimum at density 1/2, for V > V * . In the following we consider a system with open boundary conditions, coupled to 2 particle reservoirs, characterized by particle densities ρ L (the left reservoir) and ρ R (the right one). Specifying the coupling between the system and the reservoirs, that is the rates of particle injection from the left reservoir and extraction to the right reservoir, is a delicate issue, and several choices are possible. We will use the so-called bulk-adapted boundary conditions [11][12][13]26], which make the theory of boundary-induced phase transitions applicable. These conditions can be defined in such a way that they would yield bulk NESS solutions in semi-infinite systems. In order to fix ideas, consider the left boundary: the injection rate at node 1 will depend on the occupation of node 2. We would like to determine this injection rate in such a way that a semi-infinite system with nodes i = 1, 2, . . ., coupled to a left reservoir of density ρ L , exhibits a bulk NESS of the type discussed above, with density ρ L . This can be achieved if the left reservoir is equivalent to another semi-infinite system (with nodes at i = 0, −1, . . .) whose NESS is a bulk one with density ρ L . Within this framework, the time evolution of our system, at the PA level, is still described by Equations (21) (for i = 1, . . . , N) and (22) (for i = 1, . . . , N − 1), provided we define the 4-node marginals P t i [k10n] pertaining sites "outside the system" (namely, for i = −1, 0, N − 2, N − 1) in a suitable way. In particular, we have to generalize Equation (23) (which holds for i = 1, . . . , N − 3) to the following one where we have definedP where in turn P L [·] and P R [·] denote the bulk probabilities (Equations (27)-(31), together with Equation (33)) evaluated for ρ = ρ L and ρ = ρ R , respectively. From the above equations one can derive those reported in Reference [13] for the boundary rates, if Glauber rates are chosen according to Equation (1). With the above choice for the coupling between system and reservoirs, the 1D KLS model is now fully specified, and we can discuss how our methods, namely the PA, the mDWT and the extrapolation of finite size results, can be applied to characterize the NESS and the relaxation towards it.
The PA is based on Equations (21) (for i = 1, . . . , N), (22) (for i = 1, . . . , N − 1) and (36) (for i = −1, . . . , N − 1), which can be rewritten as a whole by introducing the (2N The NESS x = (ρ 1 , φ 1 , . . . , ρ N−1 , φ N−1 , ρ N ) will be given by the condition f (x) = 0, and relaxation near the NESS will be described by the relaxation matrix M, with elements In particular, the smallest eigenvalue λ 1 of M is the slowest relaxation rate, that is the inverse of the longest relaxation time, which, for our investigation of the dynamical transition, is the quantity we are mainly interested in. It is important to stress that the PA gives exact results for the bulk NESS solutions, but not for the relaxation rate λ 1 . It is therefore important to obtain estimates of λ 1 from independent techniques. In this work we will use the mDWT and the extrapolation of exact finite size results, as we did in the case of the AS model [19].
The mDWT was proposed in Reference [16], on the basis of a comparison between the standard DWT result for the relaxation rate of pure TASEP and the exact one. In the DWT, the relaxation rate is given by where D L,R = J(ρ L,R )/|ρ R − ρ L |. It turns out that the DWT result is exact in the slow phase of pure TASEP, and the dynamical transition corresponds to a maximum of the DWT rate. In the mDWT, which is exact by construction for pure TASEP, one takes the DWT result in the slow phase and the maximum rate in the fast phase. Based on our study of the AS model [19], we do not expect the mDWT to remain exact for models with interactions, but just to provide an independent qualitative confirmation of the PA results. Finally, the third independent approach we will use to estimate the relaxation rate of our model will be the extrapolation of finite size results, along the lines of References [24,32], where very accurate results were obtained for pure TASEP. The smallest nonzero eigenvalue of the transition matrix of the full master equation is numerically computed for a set of small N values, and the N → ∞ limit is then extrapolated using the Bulirsch-Stoer (BST) algorithm [33,34]. The parameter ω, characterizing the leading term in the expected size dependence, has been set at 2 (as in our study of the AS model [19]), based on the exactly known finite-size behaviour of the relaxation rate for pure TASEP [14,16], after verifying numerically that this value gives near-optimal results according to the criterion proposed in Reference [34].

Results
In the present section we compute the relaxation rate of the 1D KLS model with bulk-adapted boundary conditions, for V = 2V * (that is, in the case of the dashed line in Figure 1), as a function of the reservoir densities. We choose this particular value because in this case the NESS phase diagram, reported for convenience in Figure 2, has already been determined in Reference [13] by means of extremal current principles which turn out to be exact in the case of bulk-adapted boundary conditions. In Reference [13], it has been shown indeed that, as far as the NESS is concerned, this case exhibits all relevant new features of the model. This NESS phase diagram exhibits 7 different phases and a symmetry with respect to the diagonal ρ L + ρ R = 1 due to the particle-hole symmetry of the model. In the portion of the phase diagram corresponding to each phase we report an identifier of the phase (I to VII) and the value of the bulk density. This is indeed a quantity which characterizes the NESS because, in the thermodynamical limit N → ∞, the local density tends to the bulk density over the whole lattice, except in boundary regions whose size remains finite as N → ∞. Phases I and VI are LD phases, whose bulk density ρ L is determined by the left reservoir. Symmetrically, phases III and V are HD phases, whose bulk density ρ R is determined by the right reservoir. In addition, there are the 2 maximal current phases II and VII, whose bulk densities ρ M1 and ρ M2 correspond to the maxima of the fundamental diagram along the dashed line in Figure 1, and the minimal current phase IV, whose bulk density 1/2 corresponds to a minimum in the fundamental diagram. Among the various NESS transitions separating these phases, it is important to observe here that LD and HD phases are always separated by discontinuous (in the bulk density) transitions, denoted by solid lines in Figure 2. Exactly at these transitions, the 2 phases coexist. In the following we will compute the relaxation rate λ 1 , using the PA, the mDWT and extrapolation of finite size exact results, in the HD phases. The corresponding results for the LD phases can be obtained by symmetry, and λ 1 vanishes at the discontinuous transitions between LD and HD phases. In the maximal and minimal current phases instead the relaxation is not exponential, so the relaxation rate is not defined. More precisely, the relaxation is power-law in these phases, and the relaxation rate vanishes upon approaching these phases from LD or HD phases.
We will be especially interested in singularities in the relaxation rate, representing dynamical transitions, which separate regions where λ 1 depends only on the bulk density (that is, on just one reservoir density) from regions where it depends on both reservoir densities. In order to exemplify the various possible transitions, we shall consider three cases in detail, corresponding to the thin solid lines in Figure 2, that is ρ R = 0.75 (HD phase V), 0.98 (HD phase V) and 0.45 (HD phase III).
Before turning to a detailed analysis, let us anticipate our main results. For ρ R = 0.75 a situation similar to pure TASEP is found, that is a slow phase close to the coexistence line, separated from a fast phase by a dynamical transition. For ρ R = 0.98, due to the proximity of the minimal current phase, new features appear, suggesting the possibility of a new dynamical phase and additional dynamical transitions. Finally, for ρ R = 0.45, we clearly find 2 slow phases close to the coexistence lines. According to the PA, these slow phases are separated by 2 dynamical transitions from a central phase which is a TASEP-like fast phase, although this picture is not fully confirmed by BST results, maybe due to small system sizes.
As previously mentioned, we start our analysis with the simplest case, that is ρ R = 0.75, in HD phase V. The relaxation rate λ 1 as a function of ρ L is reported in Figure 3. The PA results (black lines) clearly suggest that, in the thermodynamical limit N → ∞, 2 phases can be distinguished. For small ρ L (slow phase), close to the (discontinuous) transition with the LD phase VI, λ 1 depends on ρ L , starting from 0 at the transition and then increasing. Notice that in this phase the size dependence is very weak, in particular the asymptotic value is approached with an exponential decay in N. For large ρ L (fast phase), λ 1 becomes independent of ρ L (while the size dependence can be appreciated, the approach to the asymptotic value is ∼ N −2 ) and a singularity, corresponding to a dynamical transition, can be clearly identified. This is the same behaviour observed in pure TASEP [16][17][18]. We can approximately locate the dynamical transition by extrapolating, using the BST algorithm with ω = 2, the relaxation rate in the fast phase for ρ L = 1 and by looking for the point where the N = 400 results reaches this value. The corresponding estimate of the transition point, computed for ρ R ∈ (ρ M2 , 1), will be reported in the dynamical phase diagram with a blue line. In Figure 3 we also report the results of the DWT (dashed red line) and the mDWT (solid red line). As in pure TASEP, the DWT result for λ 1 starts correctly from 0 at the transition with the LD phase VI, then increases with ρ L until it reaches a maximum. In the pure TASEP case, the behaviour after the maximum was shown [16] to be incorrect, and the mDWT was constructed by replacing the right portion of the curve with a constant, equal to the maximum. Notice that the mDWT estimate is much smaller than the PA one (as already observed in the case of the AS model [19]), but the location of the dynamical transition (solid red line in Figure 6) is rather close to the PA one.
We also report the results of the BST extrapolation (grey dots) of finite size exact results. The relaxation rate is computed, in a numerically exact way, for N = 4 ÷ 24. The BST algorithm with ω = 2 is then applied to the subsequences N = 4 ÷ 24 to N = 17 ÷ 24. The results seem consistent with the qualitative behaviour found in the other techniques and suggest that the PA (respectively the mDWT) overestimate (resp. underestimate) the relaxation rate.
We now move on to the case ρ R = 0.98 (Figure 4), which is still in the HD phase V but, due to the proximity of the minimal current phase IV, exhibits a richer behaviour. Both the PA and DWT estimates behave in the same way as above for sufficiently small and sufficiently large ρ L . There is however an intermediate region, close to the minimal current phase IV, where new features appear. The PA estimate becomes nearly constant, whereas the DWT one exhibits another maximum, followed by a minimum. One might be tempted to think that 2 more singularities, corresponding to additional dynamical transitions, appear, and the apparently unphysical non-monotonical behaviour of the DWT could be fixed in the mDWT by replacing the oscillating part of the curve with a constant equal to the minimum (solid red line in Figure 4). The mDWT estimates of these additional dynamical transitions are reported in Figure 6 with dashed lines, and the fact that they start from 3-phase points, like the transition discussed above (and like those found in pure TASEP and in the AS model), is encouraging. On the other hand, the PA rate shows a very weak size dependence here, as it usually happens in the slow phase where the rate depends on ρ L , so the evidence is contrasting. The BST extrapolation of finite size results is again qualitatively consistent with the PA and mDWT, which in this case both overestimate the relaxation rate for most ρ L values. In the inset we can see that the mDWT seems very accurate (based on the comparison with the BST results) for ρ L < 0.05, and then deviates from the BST results. This may suggest that a dynamical transition occurs near ρ L = 0.05, giving some support to the scenario with additional dynamical transitions, though this issue is far from being settled. Finally, we switch to HD phase III and we consider the case ρ R = 0.45 ( Figure 5). Here the PA rates seem to suggest the existence of 2 dynamical transitions, separating a central fast phase from 2 slow phases, close to the LD phases I and VI. The 2 maxima exhibited by the DWT estimate are also consistent with this scenario. The approximate location of these dynamical transitions (blue lines in Figure 6) can be found in the PA following the criterion outlined above. It is instead more tricky to find a heuristic modification of the DWT, because the location and height of the 2 maxima depend in a non-trivial way on ρ R . In particular, for ρ R < 0.43 the highest maximum becomes the right one. A criterion which leads to results qualitatively consistent with the PA would be to choose the constant value of the rate in the fast phase equal to the lowest maximum, but this would lead to an apparently unphysical singularity in the dynamical transition lines close to ρ R = 0.43, where the 2 maxima exchange their role. For this reason we do not report the transition lines estimated according to this criterion in Figure 6. The BST results seem to confirm the DWT ones close to the LD phases I and VI, but not in the central region, supporting the scenario with 2 slow phases separated by a central fast phase. Let us note indeed that in the central region the BST results are affected by a very high noise level (in some cases giving extrapolated results out of the figure range), which suggests they are not reliable at all and that longer sequences (by now beyond our current computational resources) would be needed to obtain convergence.
The full dynamical phase diagram, including the dynamical transitions in the LD phases, obtained by symmetry, is reported in Figure 6.

Discussion
Dynamical transitions are singularities in the relaxation (towards a NESS) rate, which do not correspond to any transition in the NESS itself. They have been already found in the TASEP [16] and, with approximate methods, in both the AS model [19] and the TASEP with Langmuir kinetics [20]. Dynamical transitions occur in those NESS phases where the bulk density equals the density of a boundary reservoir (HD and LD phases) and separate a fast phase, where the relaxation rate λ 1 depends only on the bulk density (suggesting that the slowest relaxation mode is a bulk one), from slow phases, where λ 1 depends on the density of both boundary reservoirs (suggesting that the slowest relaxation mode is boundary-dependent). The terms fast and slow were introduced in Reference [19] on the basis of the observation that, keeping fixed the bulk density (that is the density ρ R of the right reservoir in the HD phases) and varying the density of the other reservoir (ρ L in the HD phases), the relaxation rate takes its maximum value in the fast phase. In the NESS phase diagram the slow phases are typically located close to HD-LD coexistence (discontinuous transition) lines, where λ 1 vanishes. In a slow phase then the relaxation rate grows from 0 at the coexistence line to its maximum value, which depends on the bulk density, at the dynamical transition line.
We have investigated dynamical transitions in the 1D KLS model, which generalizes the TASEP by assuming that the hopping rate depends on the occupation state of the 2 nodes adjacent to the nodes affected by the hop. The 1D KLS model is more general than the AS model, where the dependence is only on the adjacent node in the direction of motion, and this dependence reflects in a richer fundamental diagram, which can exhibit 2 maxima and a minimum, and consequently in a richer NESS phase diagram. Following Reference [13] we have chosen Glauber rates and bulk-adapted boundary conditions. We have used the same techniques as in our previous work on the AS model [19], namely the PA, the DWT and its heuristic modification (mDWT), and the BST extrapolation of exact finite size results. The use of the PA is well justified, since with Glauber rates and bulk-adapted boundary conditions many PA results for the NESS are exact.
Our results, summarized in the dynamical phase diagram in Figure 6, seem to confirm the general picture outlined above, with slow phases close to coexistence lines, separated from fast phases by dynamical transitions. In addition, the 1D KLS model seems to exhibit certain qualitatively new phenomena. In the HD phase V we find evidence of additional phase transitions, possibly related to the appearance of a new dynamical phase close to the minimal current phase IV (by symmetry, this also applies to the LD phase I). This evidence is however only partial, and needs to be verified by additional investigations. Moreover, in the HD phase III, the nature of the central phase, a fast phase according to the PA, is not clearly confirmed by the BST results, and further investigations are also needed in this case (again by symmetry, this also applies to the LD phase VI). In both cases, further work may involve increasing the largest lattice sizes that we used in the exact finite size approach and/or moving in the parameter space of the model. In particular, varying the repulsion V in the range V > V * can be expected to shrink, or enlarge, the portions of the phase diagram occupied by the various phases, possibly making the analysis of the dynamical transition numerically simpler in certain cases. The results for V ≤ V * can instead be expected to be qualitatively similar to pure TASEP (V = 0), hence less interesting from the point of view of the present study.
Finally, it would also be interesting to check how the overall picture is affected by a change in the boundary conditions. The bulk-adapted boundary conditions that we have used here, following Reference [13], are not the only possibility. For instance, in Reference [13], a different choice was considered, named equilibrated-bath couplings, which leads to a significantly different NESS phase diagram. It is reasonable to expect that the dynamical transitions would also be strongly affected by such a change in the boundary conditions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
This appendix is devoted to a formal derivation of the time-evolution Equations (21) and (22) for the local densities and nearest-neighbor correlations, respectively. Let us start from the generic master equatioṅ where of course P t [n] ≡ P t [n 1 . . . n N ] denotes the probability of the lattice configuration n ∈ {0, 1} N at time t and W[m → n] denotes the transition rate from m to n. The matrix of transition rates can be expanded as a sum over the elementary processes occurring in the system, namely injection, extraction and hopping. In formulae we can write where we assume in particular that = 1, . . . , N − 1 corresponds to hopping from node to + 1, whereas = 0 and = N correspond respectively to injection (at node 1) and extraction (from node N). One key assumption of the model under investigation is that the hopping rate Equation (1) depends on the configurations of nearby nodes (namely, − 1 and + 2 for the to + 1 hopping process). Given that such nodes actually exist only for = 2, . . . , N − 2, in these cases we can write where δ denotes a Kronecker delta and a bar over occupation numbers denotes the boolean complement. Note that, apart from the hopping rate Γ, the remainder of the above expression is meant to ensure that W [m → n] takes a nonzero value only when the following conditions are verified: Equations (A4) and (A5) together mean that, in the "departure" configuration m, node is occupied and node + 1 is empty, whereas the opposite happens in the "arrival" configuration n. Equation (A6) means that all the other nodes do not change their occupation state from m to n.

Appendix A.1. Local Density Equation
The local densities, Equation (6), can be written in terms of the full probability distribution function P t [n] as where it is understood that the sum runs over all possible lattice configurations n. Taking the time derivative of Equation (A7) and using the expression ofṖ t [n] obtained by Equation (A1) in combination with Equation (A2), by simple algebra we can obtaiṅ Let us now consider the inner sum in Equation (A8). Taking into account Equation (A3), we can observe that in the elementary summation term (n i − m i ) W [m → n] there always appears a factor (n i − m i ) δ(m i , n i ) ≡ 0 except when i = or i = + 1. As a consequence, the sum over in Equation (A8) turns out to contain just two terms, namely = i and = i − 1, which we analyze below.
which can be obtained using Equation (A3) to express the elementary transition rate W i [m → n] and observing that, since m i and n i are boolean variables, m i 2 ≡ m i and n i n i ≡ 0. In the other case, = i − 1, by analogous arguments we get In the end, plugging Equation (A9) and Equation (A10) into Equation (A8) (and remembering that the sum over contains only these two terms), we obtaiṅ By marginalization, this last equation is equivalent to Equation (21). Let us stress the fact that the above derivation is correct as long as the elementary transition-rate matrices W i and W i−1 take the form of Equation (A3), which is the case (according to the previous discussion) for i = 3, . . . , N − 2. The remaining cases, namely i = 1, 2, N − 1, N, require the definition of boundary rates (including injection and extraction rates), and will be discussed below in the last paragraph.

Appendix A.2. Local Correlation Equation
The local nearest-neighbor correlations, Equation (7), can be written in terms of P t [n] as where as usual the sum runs over all possible lattice configurations n. The calculation proceeds in close analogy with the density equation, leading tȯ Let us then consider the inner sum in Equation (A13). Still taking into account Equation (A3), in this case we can observe that in the elementary summation term (n i n i+1 − m i m i+1 ) W [m → n] there appears a factor (n i n i+1 − m i m i+1 ) δ(m i , n i ) δ(m i+1 , n i+1 ) ≡ 0, except when i = or i = ± 1. As a consequence, the sum over in Equation (A13) turns out to contain (at most) three terms, namely = i and = i ∓ 1, which we analyze below.
which as usual is obtained using Equation (A3) to express the elementary transition rate W i [m → n]. This term can be argued to vanish, by observing that m i+1 m i+1 ≡ 0 and n i n i ≡ 0.
In the end, plugging Equation (A15) and Equation (A16) into Equation (A13) (and recalling that the sum over contains only these two terms, as we have argued that even the term = i vanishes), we obtaiṅ By marginalization, this last equation is equivalent to Equation (22). Even in this case, let us stress the fact that the derivation is correct as long as W i , W i−1 and W i+1 can be written as in Equation (A3), which is the case for i = 3, . . . , N − 3. The remaining cases i = 1, 2, N − 2, N − 1 depend on the definition of boundary rates, which we discuss in the next paragraph.

Appendix A.3. Bulk-Adapted Boundary Conditions
As already mentioned above, the reason why the previous derivation leaves out some boundary equations is that the rates of certain processes have not been defined yet. Such processes are in particular: injection, extraction, and hopping from node 1 to 2 and from node N − 1 to N. As far as hopping is concerned, a special definition is required because in the two cases, respectively, the previous or next node does not exist, so that the rate Γ would remain undetermined. The definition of all these rates is important as it corresponds to fixing the boundary conditions or, in other words, coupling the system to the boundary reservoirs. In order to do so, one possibility is to introduce some extra "virtual" nodes, namely i = −1, 0 at the left boundary and i = N + 1, N + 2 at the right boundary, and to reinterpret injection and extraction as hopping processes, respectively from 0 to 1 and from N to N + 1, all with the same rate Γ (in this framework, nodes i = −1 and i = N + 2 are just needed to ensure that all hoppings possess both a previous and a next node). To complete this alternative description, one also needs to define an extended probability distribution function incorporating the virtual nodes, which we shall denote as P t [n −1 n 0 n 1 . . . n N n N+1 n N+2 ], to be used to compute the averages in Equations (A11) and (A17) even in the boundary cases. One possibility for defining the extended distribution, which we have considered in this paper and which is usually denoted as bulk-adapted [13], amounts to building up a continuation of the "ordinary" distribution P t [n] ≡ P t [n 1 . . . n N ] on both the left and right boundary, in a pairwise factorized form like that of Equation (14), which we have previously mentioned to be exact in a bulk steady state. In particular, define P t [n −1 n 0 n 1 . . . n N n N+1 n N+2 ] ≡ P L [n −1 n 0 ] P L [n 0 n 1 ] P L [n 0 ] P L [n 1 ] P t [n] P R [n N n N+1 ] P R [n N+1 n N+2 ] where, as in the main text, P L [·] and P R [·] denote the bulk probabilities (Equations (27)-(31), together with Equation (33)) evaluated respectively for ρ = ρ L and ρ = ρ R . The basic intuition is that, if the steady state is characterized by a bulk region on either the left or the right side of the system (which we recall to be the case in LD and HD phases, respectively), such a bulk region smoothly continues outside the system. From a formal point of view, this choice for the boundary conditions also has the advantage that time-evolution equations for all local densities ρ t i (for i = 1, . . . , N) and all local correlations φ t i (for i = 1, . . . , N − 1) take the same form of Equations (A11) and (A17), and therefore Equations (21) and (22).
Let us finally observe that, in the framework of the PA, the ordinary distribution P t [n] is itself expressed in the (approximate) pairwise-factorized form As a consequence, plugging Equation (A19) into Equation (A18), we can write P t [n −1 n 0 n 1 . . . n N n N+1 n N+2 ] = ∏ N+1 where the pair and node probabilitiesP t i [·] are precisely those defined in Equations (37) and (38). The 4-node marginals can then be written as in Equation (36).