Nonexistence of motility induced phase separation transition in one dimension

We introduce and study a model of hardcore particles obeying run-and-tumble dynamics on a one-dimensional lattice, where particles run in either +ve or -ve $x$-direction with an effective speed $v$ and tumble (change their direction of motion) with a constant rate $\omega.$ We show that the coarse-grained dynamics of the system can be mapped to a beads-in-urn model called misanthrope process where particles are identified as urns and vacancies as beads that hop to a neighbouring urn situated in the direction opposite to the current. The hop rate, same as the magnitude of the current, depends on the total number of beads present in the departure and the arrival urn; we calculate it analytically and show that it does not satisfy the criteria required for a phase separation transition. Tumbling is generally detrimental to the stability of jamming; thus, our results for this restricted tumbling model strongly suggest that motility induced phase separation transition can not occur in one dimension.

An important class of nonequilibrium systems is that of active matter systems (AMS) [1] where the individual constituents are self-propelled; instances of such systems include bird flocks [2], bacterial colonies [3], photophoretic colloidal suspensions [4] and actin filaments [5] etc. They exhibit a number or interesting features like large number-fluctuations [1], clustering and pattern formation [4]. A major area of interest in the study of AMS has been the so-called motility-induced phase separation (MIPS) [6][7][8][9][10][11][12][13] which refers to spatially separated high and low density regimes. Such aggregation or clustering of particles has been observed experimentally in many active matter systems [4]. Relevance of the aggregation process has also been proposed as a mechanism of formation bacterial biofilms [3], which are sources of infection.
Occurrence of MIPS relies on an argument that effective velocity of active particles decrease in crowded or high density regions formed either by explicit dependence of local density or merely by exclusion. Naturally such a slowing down of movement further increases the density of particles and gives rise to a feedback loop allowing the stable high density (liquid-like) regions to form and coexist with a low density (gas-like) phase elsewhere. MIPS has been widely investigated in simulations and apparent phase separation has been observed. Theoretical investigations of this phenomenon have thus far concentrated on continuum models [8][9][10] where motility parameters, such as particle flux or velocity are characterized as functions of the coarse-grained local density [6,7]. Lattice models of active particles have been studied in one and two dimensions numerically [14][15][16] with run and tumble particles (RTPs). RTPs move at a fixed speed along the direction of their orientation (a run) until they tumble and change their orientation. In one dimension (1D), the two orientations (say, ±) are usually referred to as the internal degrees of the particle (spin), which flips with a certain rate. Analytical studies of these lattice models are limited. Thompson et. al. [13] have introduced a model of self propelled particles with RTP dynamics; in 1D. These models exhibit inhomogeneous density profiles when particle velocities depend on their position. Recently Slowman et. al. [17,18] have obtained an exact solution for two RTPs and found jamming induced attraction between the particles of the opposite spins, which indicates that, for many particle systems, a phase separated state might originate from these attractive interactions. Later, Dandekar et. al. [19] have obtained a mean-field solution of RTPs in 1D which turned out to be a good approximation when tumbling rate is large.
An element of surprise in the formation of a phase separated state without any explicit attractive interaction has generated much excitement to the study of MIPS and raised questions about the stability of such states in 1D in absence of any explicit interaction or spatial potential. Recent works have added to the doubt by showing that MIPS phase transition in 2D belongs to the Ising universality class [20] which does not have a counterpart in one dimension. In this article we argue and show explicitly using 1D lattice models of RTPs that indeed MIPS transition can not occur in 1D; the inhomogeneous states observed in numerical simulations and in hydrodynamic models are only long lived transient states.
First we introduce a generic model of hardcore RTPs in 1D with a restricted tumbling dynamics and show that its coarse-grained dynamics can be mapped to a beadsin-urn model, namely a misanthrope process [21] where beads hop to their neighbouring urn, situated in the opposite direction of the particle current, with a rate same as the magnitude of current. The functional form of hop rate is determined from the exact steady state results of the model with only two RTPs. To determine if MIPS transition is possible, we use the following criterion. If a system of hardcore particles phase separates as its density ρ crosses a threshold ρ * then the maximum density at which it remains homogeneous is ρ * . Since systems with homogeneous densities are well described in the grand canonical ensemble (GCE) by a unique chemical potential µ (or fugacity z = e µ ), we argue that phase separation transition is possible in a system when its density in GCE attains a maximum value ρ * = Max[ρ(z)] which is less than unity (the density of a fully occupied lattice). Nonexistence of MIPS transition in restricted tumbling model would imply that MIPS can not occur in any other RTP model in 1D where tumbling occurs more frequently.
The restricted tumbling model: We introduce a generic model of RTPs on an one dimensional periodic lattice with sites labeled by i = 1, 2, . . . L. The sites are either empty (represented by τ i = 0) or occupied by at most one RTP τ i = ± having orientation (spin) ±. Particles follow a run dynamics, where RTPs move forward or backward with rates p ± and q ± respectively. Along with this, they can tumble and change their spin with rate ω as follows, Tumbling is restricted here in the sense that only those particles which are assisted from right by other particles can tumble their direction. This restriction helps us getting an approximate steady state of the system without tampering the main aim: the proposition that a stable MIPS state can not be sustained in 1D. Since frequent tumbling of particles helps the system to clear jamming, a proof of nonexistence of MIPS in our model necessarily guarantees its nonexistence in any other model that has more liberal tumbling dynamics. Hereafter we refer to the model following dynamics (1) and (2) as restricted tumbling model (RTM). Although RTM is defined for generic rates (p ± , q ± ) we study the case p ± = q ∓ where the run dynamics exhibit a symmetry transformation, namely simultaneous interchange of parity (left right) and spin (+ −), that keeps the dynamics invariant. This symmetry was present for both run-and tumble-dynamics in 1D lattice models studied earlier [17,19]. When p ± = q ∓ , it is also ensured that in the limit when lattice spacing vanishes [22], a single particle dynamics of RTM reduces to that of a RTP moving in continuum space with same speed v = p − − q − = q + − p + along +ve and −ve x-directions. Note that, under parity transformation (left right) the tumbling dynamics of our model is modified as tumbling now occurs for only those particles which are assisted by other particles from left. But, for p ± = q ∓ , a left-assisted tumbling dynamics leads to the same steady state as the right-assisted tumbling. This can be verified easily from the exact mapping of these models to the corresponding beads-in-urn models (see later discussions).
A special case of RTM with p + = α = q − , p − = 0 = q + and unrestricted tumbling dynamics + ω ω − was studied earlier by Slowman et. al. [17] and an exact steady state solution was obtained for a system of two RTPs. It turned out that these two particles experience an effective attractive interaction in the steady state when their spins are opposite; it is envisaged that this attraction might be the source of MIPS states observed in corresponding hydrodynamic models. In comparison, in Eq. (2) we have dropped one of the transition +0 ω ω −0; as a consequence, particles do not tumble if they are not assisted by a right neighbour.
Mapping to beads-in-urn model: Any microscopic configurations {τ i } of RTM can be viewed as urns containing beads -each particle is an urn that contains beads which are uninterrupted sequence of 0s (vacancies) to the right of the particle (as described in Fig. 1(a)). The spin ± of the particle is termed as the internal degree of the urn. Thus we have a beads-in-urn model of N urns indexed by k = 1, 2, . . . N, each carrying an internal degree σ k = ± and m k = 0, 1, 2 . . . beads. The dynamics (1) and (2) now translate to hopping of a bead from urn k to k + 1 (k − 1) with rate q σ k+1 (p σ k ), and flipping of internal degrees σ k → −σ k with rate ωδ m k ,0 . The total number of beads N k=1 m k = L − N ≡ M is conserved by the dynamics. Like particle density ρ = N L , the bead density η = M N = 1−ρ ρ is also conserved. Note that in this beads-in-urn model the internal degrees of the urns can flip only when they are empty; this restriction forces k-th urn either to transfer a bead (when m k > 0)or to change the internal degrees (when m k = 0) and help us getting an exact steady state. It is easy to see that a left-assisted tumbling dynamics with same rate ω will also map to the same beads-in-urn dynamics when particles are identified as urns containing number of beads same as the consecutive vacancies to their left and the hope rates are p ± = q ∓ .
The mapping of RTM to beads-in-urn model is exact but its steady-state could not be obtained analytically. We proceed to develop a coarse-grained picture. In the steady state of the urn model, the local bead current J (summed over ± degrees) effectively transports the beads from one urn to its neighbour situated along the direction of total current. Since hop-rates (q σ k+1 , p σ k ) in the original beads-in-urn model were dependent on spins of neighbouring urns it is expected that the local bead current must depend on the number of beads present in neighbouring urns, i.e. J ≡ J(m k , m k+1 ). This current can be set as the effective hop-rate of a coarse-grained model where urns lose their internal degrees and a single bead hops from urn k to (k + 1) with rate u(m k , m k+1 ) = J(m k , m k+1 ); rightward hopping (k to (k + 1)) is considered assuming that the current is flowing in +ve x-direction. Thus, in this coarse-grained picture (see Fig 1(b)), all urns are equivalent (as they lose their internal degrees) and the hop-rate depends on the number of beads present in the departure and the arrival urn; such a process is called a misanthrope process (MAP) [21].
In fact, mapping of hardcore particle systems to urn model with an exact or effective coarse-grained dynamics, similar to the dynamics of a zero range process (ZRP) [23] are quite reliable and have helped researchers [24] earlier to establish non-existence of phase separation transition in certain lattice models [25] where rigorous numerical simulations have exhibited apparent phase separated states. It also helped in predicting true phase separation transition in many other models [24,26,27]. In contrast, mapping to that of misanthrope process, that we introduce here, provides a better coarse-grained picture as steady-state correlation between neighbouring urns are retained here. The bead-current J(m k , m k+1 ) flowing across the urns can be computed from numerical simulations (will be discussed later), but that does not help us to compute ρ(z) in grand canonical ensemble. To calculate ρ(z) we need functional form of J(m 1 , m 2 ) which can be calculated exactly using matrix product ansatz (MPA) [28] for a system of two urns containing M number of beads (i.e., L = M + 2), each one following the dynamics described in Fig. 1(a).
For urn models, a matrix product steady state (MPSS) can be obtained following Ref. [29]. We now consider RTM model, which is mapped exactly to the urn model described in Fig. 1(a). The steady state probability of a generic configuration {σ k m k }, where k th urn (spin σ k ) has m k beads, is given by a matrix product ansatz, where matrix X σ k (m k ) represents the k th urn having internal degree σ k and m k beads. The δ-function here ensures that the total number of beads M are conserved.
These matrices are constrained to follow a matrix algebra so that P ({σ k m k }) defined above must satisfy the steady state condition dP dt = 0 for the dynamics in Fig  1(a). We find (see Appendix) that for N = 2, matrices X σ (m) have a 2 × 2 representation (for any ω > 0), The steady state probabilities of two urns containing m 1 , m 2 beads are then, P σ1σ2 (m 1 , with m 2 = M − m 1 . Thus, the average local current carried by the beads when the two urns have (m 1 , m 2 ) particles is For RTPs, which need to satisfy the condition p ± = q ∓ , where v = p − − q − = q + − p + and γ = p+ p− (as in Eq. (4)). Note that J(m 1 , m 2 ) depends only on the sum of its arguments, i.e., J(m 1 , m 2 ) ≡ J(m 1 + m 2 ). We will now set J(m 1 + m 2 ) as the hop-rate of beads in the coarsegrained model, i. e., u(m k , m k+1 ) = J(m k +m k+1 ). This urn model is a misanthrope process where hop-rate is a function of total number of beads present in the departure and the arrival site. It turns out that the steady state of this specific misanthrope process has a factorized form, The grand partition function with a fugacity y that controls the total number of beads M N k=1 m k is In RTM, both N, M = N k=1 m k vary keeping the system size L fixed. To account for that we introduce another fugacity z, so that the new partition function is,  (1) and (2) (equivalently an urn model described in Fig. 1(a)). In each case, statistical averaging is done for more than 10 7 samples.
which gives rise to N = z ∂ ∂z ln Z(z, y) and M = y ∂ ∂y ln Z(z, y). We now set N + M ≡ L to obtain z in terms of y, z = L (1+L)F (y)+yF (y) . Then, The maximum value of the RTP density, obtained when y → 0, ρ * = 1 (fully occupied lattice). Thus the system remains homogeneous for any density 0 < ρ < 1 and it can not phase separate (following the criterion we discussed). The above argument is based on a coarse-grained picture where the hop rate u(m, n) ≡ u(m + n) is taken same as the average local current of beads. In the following we employ a method to calculate J(.) numerically from Monte Carlo simulations of the model and compare it with Eq. (6).
To simulate the dynamics we must set p ± = q ∓ required for the system to have a valid RTP dynamics, which gives γ = p+ p− in Eq. 4. Without loss of generality we can set p − = 1 = q + , by choosing a suitable time unit; then, p + = q − = γ and the speed of RTPs v = q + − p + = 1 − γ. We also consider γ ≤ 1 (γ > 1 case can be explored directly by using left/right and +/− symmetry). From Eq. (6), J(m) = v Qm (1 − γ m ), which has an asymptotic form (for large m), This implies that u(m) −1 is a linear function of m −1 with slope c(1 − γ) −1 and y-intercept (1 − γ) −1 , which we verify from the Monte Carlo simulations of the urn model (Fig 1(a)). For a given value of γ, ρ, ω first we allow the system to relax for a long time starting from a random initial configuration. The system may take a very long time to reach a true phase separated state when it exists, but the hoping dynamics in the coarsening regime given by u(m 1 , m 2 ) = J(m 1 , m 2 ) can predict, well in advance, if the system is approaching towards a inhomogeneous (MIPS) or a homogeneous state.
In the coarsening regime we consider a large time interval and calculate (F r (m 1 + m 2 ), F l (m 1 + m 2 )), the number of times beads move to (right, left) when the departure and arrival urns have exactly m 1 and m 2 beads respectively (internal degree of the urns are ignored). Also, we keep track of F (m 1 + m 2 ), the number of jump-events attempted during that interval. Clearly, u(m) = (F r (m) − F l (m))/F (m). In Fig. 2(a) Fig. 2(b), p(m) exhibits exponential distributions that match very well with the prediction when ω is large. As ω → 0 the exponential feature remains persistent but the value of y differs substantially from the theoretical value 1 − ρ. This is because ergodicity is broken at ω = 0; the system there falls into one of the fully jammed (or absorbing) configuration and remains there.
Essentially, the coarse-grained picture turns out to be a good description of the RTP model as p(m) decays exponentially for large m as predicted -rest of the details are less relevant because an exponential form of p(m) is enough to assure that the fugacity in GCE can always be tuned to secure any desired particle density 0 < ρ < 1. Such a system can not support any stable MIPS phase and settles to form a homogeneous density profile for all ω > 0, γ ≥ 0.
The above conclusion can also be obtained from using an approximate matrix product steady state (MPSS). Matrix representations (4), that provides exact MPSS exclusively for N = 2, are also excellent approximations for larger N (justified in the Appendix). With these matrices, for N > 2, the grand partition function Z(z, y) and density ρ(y) are given by Eqs. (17) and (18) Clearly, the maximum density that can be achieved in GCE by tuning y is ρ * = 1 (when y = 0) and thus, this RTP model can not undergo a phase separation transition at any ρ < 1. One can safely extend these results for restricted tumbling dynamics to other RTP models where tumbling occurs more frequently; this is because tumbling is generally detrimental to the stability of MIPS. Our conclusions are consistent with the recent results [20] that MIPS transition in 2D belongs to the Ising universality class that does not have an one dimensional analogue.
In summary, we show that phase separation of free hardcore-RTPs with constant run and tumble rates is not possible in 1D. One may however add some crucial features which are known to enhance or freshly produce phase separated states of passive particles, like invoking explicit attractive interaction [24] or making tumbling rates to decrease with L (so that it vanishes in the thermodynamic limit) [30] or explicitly forcing the run dynamics to depend on (and reduce substantially with increase of) local particle density [23] or adding impurities [31]. Then a phase separation transition may occur, but will it keep its charm and glory to be identified as the motility induced phase separation, particulary when the transition is anyway expected for similar system of passive particles (without motility)? Recently Kourbane-Houssene et. al. [32] have introduced a RTP model where the difference of run-rates (or effective velocity) are taken proportional to 1 L and the tumbling rate is proportional to 1 L 2 (downplayed by a factor 1/L compared to the run rates); using an exact coarse-grained hydrodynamic description they show that a homogeneous phase in 1D loses its stability in certain parameter regimes. Another way might be to use strongly biased tumbling rates where, say, + → − occurs much more frequently than − → +. In this case a phase separation transition occurs [30] when q ± = 0, where the dynamics of RTM reduces to that of a two species exclusion process [33]. Its extension to small q ± 0, is a RTP model (having a good continuum limit) and it is reasonable to assume that the phase separation features may also survive there. Yet another possibility is to introduce defects. Recent studies [34] have shown that a jammed phase does exist in RTM like models with defects. More investigations are required in all these directions to confirm if RTP models in 1D can phase separate.

APPENDIX
The dynamics (1) and (2) of RTM can be mapped exactly to an urn model described in Fig. 1(a) where beads hop from site k to site k + 1 (or site k − 1) with rates q σ k+1 (or p σ k ) respectively. The probability density of a generic configuration {σ k m k } evolves following the Master equation, = −(p σ k + q σ k+1 )P (. . . , σ k−1 m k−1 , σ k m k , σ k+1 m k+1 , . . . ) + q σ k P (. . . , σ k−1 m k−1 + 1, σ k m k − 1, σ k+1 m k+1 , . . . ) + p σ k+1 P (. . . , σ k−1 m k−1 , σ k m k − 1, σ k+1 m k+1 + 1, . . . ) − ωδ m k ,0 P (. . . , σ k−1 m k−1 , σ k m k , σ k+1 m k+1 , . . . ) + ωδ m k ,0 P (. . . , σ k−1 m k−1 , −σ k m k , σ k+1 m k+1 , . . . ) (12) where first three terms in the right hand side corresponds to the run dynamics and the rest describes tumbling at a generic site k. In the steady state d dt P ({σ k m k }) must vanish; this, along with the matrix product ansatz (3) We now introduce some suitable choice of auxiliary ma-tricesX σ k ,σ k+1 (m k , m k+1 ), yet to be determined along with X σ k (m k ), so that both k H R k and k H T k vanish separately; one such cancellation scheme for H R k is, We find that a choiceX σ,σ (m, n) = h σσ X σ (m)X σ (n) with some scalar parameter h σσ does satisfy the steady state condition with 2 × 2 matrices when γ = p++q− p−+q+ , h +− = 0 = h −+ and These matrices also satisfy the condition k Tr[H T k ] = 0 set by the tumbling dynamics because X σ (0)X σ (m) = X σ (m) for all σ, σ , m. The only troubling part is that h σσ s depend implicitly on m, n violating the assumption that they are constants. This implicit dependence of h ++ and h −− on m, n drops out when (i) q ± = 0 (all particles move in the same direction), (ii)γ = 1 (which sets the speed of RTPs v = 1 − γ = 0 when p ± = q ∓ ). In both cases we have an exact MPSS, but neither of these cases constitutes the scenario of MIPS. Yet another case is N = 2 where matrices given by Eq. (15) leads to an exact MPSS. This is because the cancellation scheme in Eq. (13) acts on product of three consecutive matrices which are not present when N = 2; thus, one can make h σσ independent of m, n by setting safely h σσ = 0 for all σ, σ . Steady state probabilities for N = 2 is given by Eq. (5). Now we proceed for larger N and get an approximate MPSS while dependence of h σσ on m, n are ignored and both h ++ and h −− are taken as q σ (1 − γ σ ) ∀m, n ≥ 0. We will see that the matrices (15) provide a MPSS which are an excellent approximation to the exact ones. The canonical partition function of the system is Note that F (y) N acts as the partition function of the system when N is fixed. From Z(z, y) = To verify if MPSS obtained here is indeed a good approximation let us calculate and compare from Monte Carlo simulations, the steady state values of η + , the average number of beads per + urn and ρ + , the fraction of urns having internal degree +, Since simulations are done at some specific L, N, we can use F (y) N as the partition function of the system; thus p + (m) = y m /F (y) and p − (m) = γ m y m /F (y) and, Using density-fugacity relation (18), both η + and ρ + can be obtained for different ρ. In Fig. 3 we plot η + and ρ + as a function of γ (dashed lines), for different ρ in the range (0.1, 0.9), along with those obtained from the Monte Carlo simulations of the model (solid lines). They match quite well for all γ < 1, indicating that, the approximate MPSS describes the RTP model very well.