Arrested States in Persistent Active Matter: Gelation without Attraction

We explore phase separation and kinetic arrest in a model active colloidal system consisting of self-propelled, hard-core particles with nonconvex shapes. The passive limit of the model, namely cross-shaped particles on a square lattice, exhibits a first-order transition from a fluid phase to a solid phase with increasing density. Quenches into the two-phase coexistence region exhibit an aging regime. The nonconvex shape of the particles eases jamming in the passive system and leads to strong inhibition of rotations of the active particles. Using numerical simulations and analytical modeling, we quantify the nonequilibrium phase behavior as a function of density and activity. If we view activity as the analog of attraction strength, the phase diagram exhibits strong similarities to that of attractive colloids, exhibiting both aging, glassy states and gel-like arrested states. The two types of dynamically arrested states, glasses and gels, are distinguished by the appearance of density heterogenities in the latter. In the infinitely persistent limit, we show that a coarse-grained model based on the asymmetric exclusion process quantitatively predicts the density profiles of the gel states. The predictions remain qualitatively valid for finite rotation rates. Using these results, we classify the activity-driven phases and identify the boundaries separating them.


I. INTRODUCTION
Active matter, constituted of particles that convert ambient energy to directed motion, has emerged as an important class of non-equilibrium systems with examples ranging from bacterial suspensions to synthetic colloids. Being driven out of equilibrium at microscopic scales, the collective dynamics of these systems are far richer [1][2][3][4][5][6][7] than thermal systems, which are bound by fluctuation-dissipation relations.
A particular collective behavior that has been widely studied is motility-induced phase separation (MIPS) [8][9][10][11][12][13]. MIPS, a kinetic phenomenon, is striking in its similarity to equilibrium phase separation such as in passive colloids with attractive interactions [14]. The nonequilibrium phases and transitions between them while having exact analogs in equilibrium systems [8][9][10][11][12][13] exhibit anomalous fluctuations that can be traced back to their non-equilibrium nature. The universality of MIPS has led to the proposition that activity mimics attraction [10]. Under certain conditions, the non-Brownian random walks representing active-particle dynamics can be mapped onto systems with detailed balance [15].
In addition to phase separation, passive colloids exhibit dynamically arrested phases in the form of glasses and gels. These two types of disordered, amorphous solids have distinct structural and dynamical signatures. The glass transition occurs in both repulsive and attractive colloids at packing fractions close to random close packing and is structurally homogeneous on large length scales. Gelation in attractive colloids leads to strongly heterogeneous states with fluid-like regions coexisting with an arrested, percolated, dense phase [14]. Active analogs of the glass transition [16][17][18][19][20][21][22][23][24] and jamming [25,26] have been explored extensively in active Brownian particles (ABPs) interacting via repulsive po-tentials. A review of the emergent behavior of active particles in crowded environments appears in [27]. Recent work on an extreme limit of ABPs with long persistence time of their self-propulsion direction has revealed fluctuations in the dense limit that are qualitatively different from those at short persistence times [28]. In this limit of long, yet finite persistence times, clustering and heterogeneous dynamics analogous to passive gels has been observed [29] lending further credence to the idea that activity can act as an effective attractive interaction. Similarly, soft disks with infinitely persistent active motion along quenched random directions also display a transition to an absorbing jammed phase above a critical density [30]. Other studies have focussed on the variation of persistence times in such systems using the static fluid structure as well as non-equilibrium velocity correlations [19][20][21][22]. The activity-induced change in the effective attraction in these systems depends on the microscopic details of the particles as well as their persistence times, and therefore activity may enhance or suppress glassy dynamics.
In this paper, we explore the two well-known paradigms of dynamical arrest in passive colloids, gelation and glass formation, in a lattice model of ABPs with purely repulsive interactions but with a non-convex shape that can interlock and hinder rotations. Using numerical simulations and a coarse-grained model based on a mapping to an asymmetric simple exclusion process (ASEP) [33], we classify the activity-induced phases and construct a non-equilibrium phase diagram. A novel feature of this phase diagram is the appearance of a phase with coexisting voids and solids separated by an interface wetted by an active fluid. The coarse-grained model quantitatively predicts the width of this interface, which combined with conservation laws leads to predictions of the non-equilibrium phase boundaries separating the ar- The color-bar represents the stationary time for each particle on log scale, log 10 (τi(tmax)) (see Section III). Unoccupied lattice sites are colored white. Particles that have not taken a single step through the duration of the simulation are colored darkest. The phases we find by varying the density ρ and activity ∆v are (a) steady state passive fluid, (b) passive fluid-solid/aging glass states, (c) finiteactivity states resembling the aging regime of the passive system, (d)-(f) active phase with void-solid coexistence (see Section V). There is a progression from majority non-arrested states (d), to majority arrested states (e)-(f) with increasing density. Panel (g) displays a color map of the fraction of arrested states (see Section IV). Locations of the snapshots shown in (a)-(f) are marked on the phase diagram with pink circles. ρrcp in the figure denotes the random-close-packing density of hard crosses: the maximal density that can be reached via the RSAD process [31]. The morphology exhibited by the snapshots at high densities (e)-(f) are reminiscent of high-density gels in attractive colloids, while (c)-(d) resemble the gel-bubble states at low densities [32]. rested states. The main result, summarized in Fig. 1, is that activity triggers arrest into a percolated phase of immobile particles akin to a gel in attractive colloids. At high densities, the transition is from an aging glass, whereas at low densities it is from an active fluid. To our knowledge, this model provides the first realization of an activity-induced transition from a repulsive glass to a gel.
As described in detail below, we study an active lattice gas model [11,34] of hard-core, cross-shaped particles on a square lattice, which is also referred to as the N3 model since each cross prevents the occupation of the first, second and third neighbors of its central square, as shown in Fig. 2. In the passive limit, this is the simplest lattice-gas that exhibits a finite-density first-order transition from a fluid phase to a sublattice-ordered phase with tenfold symmetry [35]. The sublattice-ordered states can be further grouped into right-handed and left-handed chiral order. In continuum, experiments have demonstrated the emergence of long-range chiral order in crystals of crossshaped particles [36]. As the density is quenched into the two-phase coexistence region [31,37], one observes a crossover from a simple fluid to a slowly coarsening or ag-ing regime in which concentrated immobile clusters with local crystalline order emerge. There is evidence for the existence of a glass transition [38] in this passive system in the form of diverging timescales and the appearance of dynamical heterogeneities [39,40]. This paper is organized as follows. Section II describes our model and the simulation methodology. In Section III, we quantify the spatially heterogeneous dynamics that is visible in the snapshots shown in Fig. 1. Next, in Section IV, we classify the states into two categories, arrested and non-arrested, based on measurements of the mean-squared-displacements (MSDs) of the particles. The MSD measurements also distinguish between aging, glassy states, and gel-like arrested states. This classification is used to construct the phase diagram shown in Fig. 1 (g). In Section V, we present a coarse-grained model of the dynamics that leads to a prediction of the density profile. We compare these results to the density profiles obtained in our numerical simulations in Section VI, and construct a non-equilibrium phase diagram that delineates states based on the density profiles of the arrested states. This classification connects the dynamical signatures of arrest shown in Section IV to phase sepa-  2. (a) Nearest neighbor labels and hopping rates for our active hard-cross model on a square lattice. The particles can perform thermal moves in any of the four directions with a rate v0 (black arrows), and active moves with a rate v0 + ∆v along the direction of their orientation (red arrow). (b) The cross at the central site (green) is prevented from rotating by the presence of neighboring crosses occupying the fourth and fifth nearest-neighbor sites.
ration. Lastly, in Section VII, we consider the effects of finite persistence times. The appendices provide further details of the passive, glassy dynamics, and discuss finite size effects.

II. MODEL AND SIMULATIONS
We study a model of hard-core particles on the square lattice with exclusion up to third nearest-neighbors [35,38,[41][42][43]. Each particle can be represented as a hard, cross-shaped object occupying five lattice sites, see Fig. 2. The highest possible density of ρ = 0.2 is achieved for a perfectly ordered arrangement of crosses. This equilibrium model is the simplest hard-core exclusion model on a lattice that exhibits a first-order transition from a fluid to solid phase, with coexistence of fluid at density ρ fluid ≈ 0.16 with a crystalline solid at ρ solid ≈ 0.19 [35,38,[41][42][43]. There are 10 distinct sublattice orderings possible for the crystalline packings. The competition between these phases leads to frustration at high densities, and indeed this model displays a glass transition at higher densities [31,39,40]. By adding activity we can, therefore, study active analogs of the glass transition and gelation in passive colloids [14,32,44].
In the active generalization of the model [45], each particle is assigned an active direction which can point along any of the four lattice directions 0, π 2 , π, 3π 2 . The particles perform active Brownian walks on the lattice with a rate v 0 + ∆v along the active direction, where v 0 is the "thermal" hopping rate along each of the four lattice directions. Each particle can change its active direction by ± π 2 with a rotation rate D R . The thermal diffusion coefficient, D T = a 2 v 0 , where a is the lattice spacing, is set to unity in our simulations. Since there is no energy scale in this model, the only role of temperature is to set the magnitude of the diffusion coefficient, which simply fixes the unit of time. The expected self-propulsion velocity for a single cross in the dilute limit is v p = a∆v [45], and hence the translational Peclet number is given by P e t = vpa D T = ∆v v0 . Since we fix v 0 = 1, we use ∆v to represent P e t .
The active dynamics we prescribe for hard crosses are identical to those implemented in simulations of MIPS for a simple-exclusion lattice gas model [11]: squares on a square lattice, which do not exhibit an equilibrium phase transition or glassy dynamics. Further, we consider rotation of the active direction and require a rigid rotation of the whole cross. Consequently, rotations are disallowed for crosses that have neighboring crosses which occupy either the fourth or fifth nearest-neighbor site, see Fig. 2 (b). We note that unlike simulations of continuum active dynamics, D T cannot be set to zero because of a kinetic trap which only exists for random walks on a lattice [11]. This is especially true for non-rotating active particles: without thermal moves to free them, non-rotating active particles become immediately trapped upon collision [46,47].
We use a continuous time, rejection-free, kinetic Monte Carlo algorithm to implement the active dynamics [48][49][50]. All allowed events in the system are assigned a rate, and the relative weight of each rate determines the probability for the event to occur. Time proceeds by randomly selecting an event, and then advancing the clock by an interval − log(r)/R, where r is a uniform random variable, and R is the sum of all non-zero rates for the allowed events at a give time t. The time increments after each event are exponentially distributed with mean ∆t = 1 R . This algorithm is especially efficient for simulations at large densities, where most moves are disallowed by the excluded volume constraint.
The initial states of the system are prepared using a Random Sequential Adsorption and Diffusion (RSAD) process [37] which can generate disordered packings up to a maximum density of ρ rcp = 0.1717... [38], corresponding to the random-close-packing density for hard crosses. We study a range of global densities between ρ = 0.10 and ρ = 0.17, and a range of activity values ∆v from 0 to 1. Note that even though ∆v is a dimensionless activity, it may take values larger than unity [45].The system domain is a two-dimensional square box of linear length L, periodic boundary conditions, and a fixed total number of particles N = ρL 2 . Unless otherwise stated, we present results for L = 450. The longest simulation time is set to t = 2 × 10 6 , which is much larger than the α-relaxation time at ρ = 0.16, τ α ≈ 100 (see Appendix A).
In this work we focus primarily on the limit D R → 0, where the persistence of the active motion becomes infinite. In this limit the self-propulsion directions of the particles are "quenched" random variables, and we assign these uniformly with an equal number of the particles having an active direction along each of the four lattice directions. Studying active matter in this limit provides us with the opportunity to probe the strongest depar-tures from equilibrium systems [28], and as we show in Sections V and VI, allows us to construct a coarse-grained model and obtain quantitative estimates of the density profiles observed in the arrested states. Because of the non-convex shape of the particles, we expect that even for finite rotation rates (D R > 0), neighboring crosses should inhibit one another from rotating [45]. Therefore including rotational locking is important for modeling cases in which non-convex self-propelled particles must reorient their body axis in order to change their direction of motion. In the presence of this rotational locking effect, small but finite D R leads to qualitatively similar results as the infinite persistence-time limit, including overall global arrest due to percolating gel-like structures, as shown in Section VII.

III. DYNAMICAL HETEROGENEITY AND ACTIVITY-INDUCED AGING
Our measure of dynamical heterogeneity is based on the definition of "stationary times" for each cross. The stationary time τ i (t) at the observation time, t, is defined to be the time that cross i has spent at its currently occupied site, r i , namely τ i (t) = t − t i , where t i is the time at which particle i arrived at r i . Distributions of τ i (t) provide a quantitative measure of the spatial heterogeneity of the dynamics, and is closely related to the distribution of persistence times used to analyze the glass transition in kinetically constrained models [51][52][53]. In the aging regime and at large activities, the distributions of τ i depend explicitly on time, therefore, spatial configurations of τ i (t max ) are used to construct the color bar in Fig. 1. The distributions of τ i (t max ) are shown in Fig. 3.
We first discuss the nature of dynamical heterogeneities in the passive system. The adsorption process of the RSAD protocol generates, at time t = 0, a uniform configuration with density ρ in which crystalline order is minimized. In the passive system, the diffusion of the crosses at t > 0 leads to equilibrium configurations for ρ ≤ 0.1625. As shown in Fig. 1 (a), these equilibrium states have very little dynamical heterogeneity. The probability distribution of the stationary times, P (log 10 τ i (t max )) has a single peak with an increasingly broad tail at large τ i as ρ → 0.1625 ( Fig.  3 (a)). For ρ > 0.1625, the equilibration process is interrupted by an aging/coarsening process in which clusters of increasingly immobile particles grow at an exceedingly slow rate, which prevents the system from reaching a time-translational invariant state. The onset of aging is indicated by the appearance of a bimodal distribution ( Fig. 3 (a)) and distinct clusters of particles with τ i (t max ) t max appearing within a background of particles with 10 1 < τ i (t max ) < 10 4 ( Fig. 1 (b)). Bimodal distributions of persistence times have been used to identify dynamical heterogeneities in several glass-forming kinetically constrained models [51,54]. We show in Appendix A that standard measures such as the self-intermediate scattering function also indicate the onset of aging at ρ ≈ 0.1650. At non-zero activity, P (log 10 τ i (t max )) can develop a bimodal structure at densities lower than 0.165, as seen in (Figs. 3 (b)-(d)). Our results, therefore, indicate that activity leads to a lowering of the onset-density for glassy dynamics, in agreement with the results of continuum active glass-forming models [18,55]. As in the passive, aging system, these bimodal distributions are associated with states that do not have time-translational invariance. However, in the active fluids we observe two distinct classes of such states. At densities ρ ≥ 0.16, the states ( Fig. 1 (c)) resemble the passive aging fluid ( Fig.  1 (b)) with growing clusters of immobile, solid-like regions suspended in a fluid. Increasing the activity in this density regime leads to states with percolated clusters of immobile particles as seen in Figs. 1 (f). At densities lower than this regime of activity-induced aging, the appearance of a bimodal distribution in P (log 10 τ i (t max )) with increasing activity is accompanied by a clear spatial separation of the particles into "voids" and dense regions accommodating the most immobile particles (Fig 1 (d)-(e)). In this regime, we observe large variations in the final structures from one simulation run to another.
The variance of the stationary times, Q(t) = (τ i (t) − τ i (t) ) 2 , provides a global measure of the time evolution of spatial heterogeneity in our dynamics. Fig. 4 shows the ensemble averaged time series, Q(t) , of Q(t) at different densities and activities. In the passive sys- tem, Q(t) rapidly reaches a small steady state value in the non-aging regime of densities ( Fig. 4 (a)). Both the equilibration time and the magnitude of the steadystate dynamical heterogeneity grow with density, until for ρ > 0.1625, the equilibration time surpasses the maximum simulation time. In this aging regime, Q(t) grows indefinitely. Activity drives strong growth of Q(t) at lower densities. However, as shown in Figs there is a delay time, τ active , before which the system dynamics match the passive system. To quantify this delay time, we define τ active to be the time at which the derivative of the stationary-time variance becomes greater than a small threshold, d Q(t) /dt > . For any activity the deriviative is positive for t > 0, but it may remain small for some time: we find = 10 2 gives a robust signature for the approximate time when the growth rate of Q(t) first becomes significantly different from zero, see Fig. 4. We have checked that the choice of does not significantly change our results other than providing a scale. This delay time increases with decreasing activity or density. At very weak activities, therefore, we observe the passive behavior of Q(t) since τ active increases beyond our maximum simulation time. For instance, at ρ = 0.1600, Fig 4 (c) shows that the growing dynamical heterogeneity only appears for activities ∆v ≥ 0.025. In attractive colloids, the long-time behavior of MSD's has been used to analyze the gel-glass transition [44]. In the next section, we quantify the dynamical arrest of the active states using this measure. The phase diagram shown in Fig. 1 (g) is constructed from these measurements. The states with percolated regions of immobile particles, seen in Figs. 1 (e)-(f), are dynamically arrested gel-like states [32]. The non-arrested states ( Fig. 1 (d)) resemble fluids with suspended "beads" of gels [32]. In the arrested states, the long time behavior of Q(t), which is the analog of the zero-wavevector, four-point susceptibility, χ 4 (q = 0, t) [56], resembles that observed at small wave vectors in a model of chemical gelation [57]. A striking feature of the arrested and non-arrested states at high activities is the appearance of "voids". We show below that the appearance of these voids is a manifestation of an extreme form of MIPS for these non-rotating particles that can be understood from ASEP dynamics. Physical gels arise from arrested phase separation [14], the gel-like states in our active lattice gas similarly seem to arise from arrested MIPS.

IV. COLLECTIVE ARREST, PRESENCE OF "ABSORBING STATES"
In this section, we explore the appearance of arrested states through measurements of the MSD of individual particles. The MSD of all particles (i = 1, ..., N ) in the system is defined as In order to determine if states are arrested, we study the log derivative of this MSD for each system, defined as which has been used for studying dynamics in a model of gelation [44]. This observable is an instantaneous exponent that indicates whether the MSD is diffusive (γ = 1), ballistic (γ = 2), or arrested (γ = 0). If there is percolation of a dynamically arrested phase, then γ(t) → 0 at long times. Since the percolated phase acts as a solid boundary for the fluid-like particles, this is expected: for a random walk confined within a box, it is well known that the MSD crosses over from diffusive growth R 2 ∼ t to a flat plateau R 2 ∼ t 0 once the walker reaches the walls. Fig. 5 shows MSD and γ(t) measurements for individual runs at large and intermediate activity values, ∆v = 0.20 and ∆v = 0.07. We can construct a meanfield model for the behavior of γ(t) in the non-arrested states by using known results about the dynamics of a single non-rotating active tracer moving in a background of passive particles with density ρ [45]. For a non-rotating active tracer moving on a lattice at ρ = 0, the MSD is given by ∆r 2 i (t) = D 0 (4+∆v)t+D 2 0 ∆v 2 t 2 [11]. This form may be generalized to higher densities by introducing a density-dependent diffusion coefficient D(ρ), giving Note that this requires a measurement of D(ρ) for the passive lattice gas. For a single active tracer, therefore, γ(t) displays a smooth crossover from 1 to 2 at a characteristic timescale, which we can estimate from γ(t * ) = 3/2, yielding t * = (4 + ∆v)/(∆v 2 D(ρ)). This crossover time increases as the density is increased or activity is decreased. This is the same trend as exhibited by τ active . Thus, any non-trivial collective behavior arising from the activity appears once the ballistic motion takes over. We define arrested runs to be those for which γ(t max ) < 0.5. For strong activities all runs fall into the arrested class, whereas at lower values of ∆v only a finite fraction of runs become arrested. When considering the ensemble of possible dynamical trajectories, we observe two different types of activity-driven behavior: arrested states for which γ(t) → 0 at long times, and active states with γ(t) fluctuating around a value of 2. It is clear from Fig. 5 that the differentiation between arrested and nonarrested runs emerges only at times longer than t * .
The behavior of γ(t) offers the clearest contrast between activity-induced arrest, as seen in our model, and the attraction-induced gelation seen in passive colloids. In passive colloids, diffusing particles become arrested either to indicate a glass or gel transition. In contrast, as seen from Fig. 5, it is the persistent motion, indicated by the ballistic behavior, that leads to arrest in the active hard crosses. The coarse-grained model for density inhomogeneity that we present in the next section is consistent with this picture.
In passive systems, low-density gels can exhibit subdiffusive behavior at intermediate times but γ(t) asymptotes to unity since there are always particles that have finite mobility and can diffuse [44]. In fact, the MSD exhibits only weak signatures of gelation in passive colloids [14,58,59]. The behavior we observe is akin to colloids trapped in a porous environment [60]. Since the arrested states (Figs. 1 (e)-(f)) are characterized by a percolating network of immobile particles with no gaps between them, all the active crosses are effectively trapped and thus γ(t) → 0 at long times. A similar suppression of diffusion was demonstrated in the 2D kinetically constrained model studied in [61], and also in the context of diffusion through biological anisotropic fibrous environments [62]. In the non-arrested states ( Fig. 1 (d)), there are pores through which the particles can escape and thus they exhibit sub-diffusive behavior at intermediate times, as seen in Fig. 5, but γ(t) → 2 at long times as the particles recover their persistent, ballistic motion.
The color map of arrested and non-arrested states shown in Fig. 1 (g) was constructed from a measurement of the probability of becoming arrested at different values of (ρ, ∆v). Between 10 to 35 runs were conducted at each set of parameter values, and the color bar indicates the fraction of those runs for which γ(t max ) < 0.5. Fig. 1 (g) shows that larger activity is required to generate arrested states at lower densities. For states in the range (ρ > 0.16, ∆v < 0.03), none of the states be- come arrested within our simulation time since the time t * needed for the active particles to exhibit ballistic motion is dramatically slowed down by both the small value of the diffusion coefficient D(ρ), as well as the very weak activity ∆v. Consequently, the states in this range of densities resemble those found in the passive system.

V. COARSE GRAINED DENSITY PROFILES AND HYDROSTATIC LENGTHSCALE
Since dynamic differentiation emerges between different realizations of the simulations over a range of activity and density values, a further question arises about how these different classes of states differ structurally at long times. For our system of non-rotating, infinitely persistent active crosses, the morphology of the arrested states (see Figs. 1 and 6) depict voids coexisting with a solidlike (ρ solid ≈ 0.19), highly immobile system spanning network, and a "wetting" layer of particles with intermediate mobility and density. The appearance of voids seem to be a feature of MIPS in the zero-rotation limit [28,63].
Starting from a random orientation of particles, as the dynamics of the system progresses, clusters of particles emerge with opposing active orientations, i.e. oriented towards each other, forming a solid-like high density region. At large times, mobile particles have an average active drift towards fluid-solid interfaces, as the particles oriented away from an interface have had sufficient time to travel to the other boundaries in the system. This ac-

Increasing activity
Increasing density FIG. 6. Contour maps of the coarse-grained density profile of arrested states at the longest simulation time tmax = 2 × 10 6 . We used a coarse-graining box size with area L 2 /900, so that each box encloses 15 2 lattice sites. The colorbar represents the local density which varies from 0 ≤ ρ ≤ 0.2. The arrows represent the gradient of the density field, the lowest density point for each figure has been shifted to the center of the frame (making use of the periodic boundary conditions). Rows show fixed activity ∆v = 0.20, 0.09, 0.05, and columns show densities increasing left to right ρ = 0.13, 0.14, 0.15, 0.16. Density fluctuations for a typical non-arrested active liquid state (ρ = 0.12, ∆v = 0.040) are shown in the lower left corner. As the density is increased at a fixed activity, the width of the solid phase increases, and the active liquid fills the void region until it disappears. We observe that the width of the interface region is controlled solely by the activity ∆v.
tive flux towards the interface causes an increasing density in the vicinity of the solid, giving rise to a diffusive current away from the solid. In order to model this process, we coarse grain the system to construct a spatially varying density field ρ(x, y). We consider a 1D section of the system perpendicular to an interface (oriented in the y-direction for convenience) between an active fluid and the solid, giving rise to a linear density profile ρ(x). The exclusion due to particles in adjacent rows as well as their lateral diffusion give rise to correlations, which we ignore for large enough coarse graining blocks. This preference for biased motion perpendicular to the interface can also be modeled using the well-known ASEP model which incorporates both the hard-core exclusion along with diffusion and biased motion. The steady state density profile in the arrested states can then be derived from a hydrodynamic treatment of the ASEP [33], as we show below.
The ASEP is a paradigmatic model where many exact statements can be made regarding the coarse grained dynamics in a non-equilibrium system. We can therefore, through this mapping, write equations for the coarse grained densities appearing at late times in our active crosses system. In steady state, the density ρ(x) is independent of time, and there is no net particle current between the different coarse grained blocks. There are two components to this current determined by the density profile ρ(x): (1) A diffusive (or thermal) current J T arising due to the spatial variations in density, which to lowest order is (2) An active current J A proportional to the density of particles. However, if neighboring blocks are at high densities, this current decreases due to exclusion. Once again to lowest order we have where α is an as yet undetermined proportionality constant. These are essentially the mean field currents in ASEP [64].
In steady state, the net currents are zero. Hence, combining Eqs. (4) and (5), we obtain: The only two homogeneous solutions to this equation are voids with ρ = 0 and the solid state with ρ = ρ solid . Eq. (6) We note that this "hydrostatic lengthscale" diverges in the limit of zero activity. This divergence as the activity is decreased is shown in Fig. 8 along with the theoretical prediction from Eq. (7) showing near perfect agreement with the ∆v −1 decay. Note that Eq. (7) does not involve the global density of the system, and the correlations between rows in our two-dimensional lattice can provide non-trivial corrections to the derived behavior for larger densities.

VI. NON-EQUILIBRIUM PHASE DIAGRAM
As seen in Fig. 6, two types of profiles are observed in the arrested states: coexistence of an "active liquid" with solid regions (bottom row), and a solid network, punctuated by voids that have a characteristic size, and an active liquid interface separating the two. In addition, there are non-arrested, active liquid states with density fluctuations of amplitude much smaller than the solid density. In the previous section, we showed that our theory correctly predicts the variation of the width of the interface. In this section, we extend our analysis to construct a non-equilibrium phase diagram of the arrested states, and provide a theory for the emergent length scales characterizing the voids and the solid regions.
We can derive phase boundaries between the three types of states by considering the conditions that must be satisfied at a given density and activity to create each of these configurations. For the arrested states at strong activities, the dense immobile solid (ρ solid ≈ 0.19) is bordered by an active liquid interface that can be fit by a linear profile (as an approximation to the sigmoidal solutions of Eq. (6)) with slope m = dρ dx and width ξ, such that ξm = (ρ solid − ρ void ) = 0.19. We use these linear density profiles (as shown in Fig. 7) for all states, including the active liquid, in the computation of the phase diagram. The total length of the system is fixed at L, and the total number of particles in the system is conserved as N = ρL.
The first condition needed to create a solid region along with an active liquid interface is that there must be enough mass available in the system to populate both these regions (as in Fig. 7 (a)). The total area under the trapezoid must conserve the total number of particles in system at a given global density ρ, yielding N = ρL = (ρ solid )(ξ + l solid ). (8) Since the interface width is determined by the activity (Eq. (7)), the active liquid region can only contain a fixed activity-dependent mass, and consequently the width of the solid region depends on both the density and the activity as The solid first appears at the point where l solid = 0, yielding the equation for the solid-active liquid phase boundary ∆v = D αL Below this activity value, the interface width is larger than L/2, therefore we expect the active liquid phase to contain all the mass in the system. So far we have only imposed conservation of mass on the system. To determine when voids first appear, we must also consider whether the shape of the trapezoidal solid-active liquid profile can fit within the total space of the system. Since none of the regions overlap, their lengths sum to the system size L. We therefore have (see Fig. 7 Since the thickness of the solid region is constrained by conservation of mass, Eq. (9), we can solve for how much space remains available for the void region, To determine when a void region may first appear, we find the point l void = 0. This yields the second phase boundary between states containing voids and those without voids ∆v = D αL Below this activity the active liquid is confined within a space smaller than its preferred width ξ. Phase boundaries based on the 1D linear profile analysis, presented above, are shown in Fig. 9 (a). In Fig. 9 (b) we show that a numerical classification of the three types of non-equilibrium phases agrees qualitatively with the 1D theory. It is straightforward to extend the treatment developed in this section to a derivation of the phase boundaries based on 2D density profiles. We find that the qualitative features of the resulting phase diagram do not change as compared to the 1D case studied here. The states in the phase identified as the solid-active liquid, ). We have used the observed value D α = 10, ρ solid = 0.19, and the simulated system size L = 450. (b) Phase diagram obtained from the numerically sampled phase space of the active lattice gas, including the passive and the active states. Each point is classified as either pure active liquid, solid + active liquid, or solid + active liquid along with void regions. Open symbols denote the transition from the passive liquid to the aging glassy regime. These classifications are based on ensemble-average density profiles. The topology of this phase diagram is the same as that observed in attractive colloids [14,32,65]. The black region denotes densities larger than ρ solid = 0.19. Note that the predicted phase diagram shown in (a) is valid only at finite activities and is expected to accurately capture the physics only at large enough activities where activity-induced correlations overwhelm the correlations intrinsic to the passive hard crosses at ρ ≈ ρrcp. This is because, as discussed in the text, the coarse-grained model includes only simple exclusion, not the extended N3 exclusion.
based on the density profiles, are dynamically the aging glassy states shown in Fig. 1. The arrested, gel-like states, shown in Fig. 1 are in the solid-active liquid-void region of Fig. 9.
In Appendix C we show that the variation of the phase For the smaller activity value, the non-rotating crosses show phase separation into solid and void-like regions, with an intervening active liquid region (c). In contrast, at DR = 0.005, there is a simpler fluid-solid coexistence indicated by a clearly bimodal distribution (d).
boundaries with system size is consistent with the predictions of the theory derived above. In the infinite system size limit, the ASEP analysis implies that voids can be accommodated at any density and activity, as the phase boundaries become lines with infinite slopes at ρ = 0 and ρ = ρ solid . This feature is a consequence of the simple exclusion process, which cannot describe the non-trivial correlations in the passive, hard-cross system (c. f. Fig.  9 (b)). We are currently exploring avenues for incorporating these effects, possibly through the construction of a large-deviation function that incorporates the physics in of void-solid coexistence of persistent hard crosses, encapsulated in the ASEP model, and the glass transition physics of passive hard crosses.

VII. EFFECT OF FINITE ROTATION RATE
In this section we discuss the effect of a finite but small rotation rate on the arrested phase separation observed for non-rotating hard crosses (D R = 0) derived above. We perform simulations in two regimes D R = 0.005, ∆v = 0.1, and D R = 0.005, ∆v = 0.5. It should be noted that D R is the attempted rotation rate. The non-convex shape of the hard-crosses leads to strong rotational locking that leads to a smaller effective rotation rate, which also decreases as the density increases.
We find that at large activities ∆v ≥ 0.5, a percolated, globally arrested solid network, similar to the ones observed at D R = 0, is able to form in the same regime of densities as shown in Fig. 1. However, there are structural differences between the arrested states at zero and non-zero D R . Specifically, the voids observed at D R = 0 are replaced by a low-density, gaseous fluid at D R = 0.005 as seen in Fig. 20. The fluid phase still exhibits an inhomogeneous density profile with the highest densities occurring at the interface with the solid. In addition, a finite rotation rate seems to create a rounding of the facets separating the solid from the fluid. At lower activities, the fluid density becomes homogeneous, as seen in Fig. 21. The morphology of these states resemble those observed in the active-aging regime (Fig. 1  (c)) rather than the percolating arrested solid observed at D R = 0, seen in Fig. 21 (a).
The ensemble-averaged density distributions shown in Fig. 10, demonstrate that at D R = 0.005, the "void" peak disappears. In addition, there is a clear signature of two-phase coexistence between a fluid and a solid, with little dependence of the peak positions on the global density [66]. With increasing activity, the peak marking the fluid density shifts significantly to lower densities whereas the solid peak stays roughly pinned ρ ≈ 0.19. Since the ensemble-averaged density distributions do not depend on the global density at a fixed activity (Fig. 21), a binodal line can be constructed from the positions of the peaks, as shown in Fig. 11. Interestingly, this phase diagram looks different from that of spherical ABPs [66], and remarkably similar to the one observed in activelydriven dumbbell shaped particles, which are also nonconvex [12].
We can analyze the stability of the D R = 0 phase diagram to small but non-zero D R . On the lattice, since the particles can only have discrete orientations, we can compute the mean free path of the particles in the persistent direction. Since the density of rotation events in time have a gap distribution D R exp(−D R t), and the average time between succesful rotation events scales 1/D R [45], the average length traveled without a rotation move, is given by ξ R ∼ (∆v)(1/D R ). In order for the rotational diffusion to affect the motion of the particles, there must be a sufficient time for the particles to move before reaching the dense regions and become a part of the wetting fluid. Therefore, when the average length traveled in the persistence direction is comparable to the size of the voids l void , the phases appearing will no longer correspond to the zero rotation limit. This provides a crossover value D * R = ∆v/l void , above which the rotational diffusion takes on significant effect. Since the size of the voids can be estimated from Eq. (12), we can estimate the crossover rotation rate below which voids are expected to appear in the system For the values of ∆v used in the simulations, and using the values D/α = 10 and L = 450 deduced from the data, we estimate the crossover value for density ρ = 0.1 to be D R = 0.005 for ∆v = 0.1, and D R = 0.009 for ∆v = 0.5. Since the crossover value of D R decreases with ∆v, we expect that as the activity is decreased, the qualitative features change from those observed in arrested states at zero rotation, and begin to resemble a MIPS between a low and high density fluid, as is observed in Fig. 21.
To summarize, we find that the most significant difference between hard crosses with D R = 0 and D R = 0.005 is the disappearance of voids, which also leads to changes in the phase diagram. In earlier studies, the appearance of voids has been associated with active particles with infinitely persistent self-propulsion directions [63].

VIII. DISCUSSION
In this paper we have studied a lattice gas of active particles with a non-convex shape that leads to strong inhibition of rotations. We showed that unlike the usual MIPS observed in active particles with large rotational diffusion, highly persistent particles are able to arrest the phase separation leading to states that are characterized structurally by a void-solid coexistence with a fluid of relatively mobile particles "wetting" the void-solid interface. The voids have a characteristic size that depends on density and activity in contrast to coarsening, as they would if the phase separation was not arrested. Dynamically, the glassy dynamics characterizing the aging regime of the passive system transitions to complete arrest in the phase-separated regime. In addition, activity enlarges the density range of the aging regime. We also show that adding moderate rates of rotational diffusion does not change this picture qualitatively.
By appealing to the dynamics of a single active tracer in a passive background of hard crosses [45], we showed that the dynamic differentiation between arrested and non-arrested states appears at times longer than that needed for the tracer dynamics to cross over from diffusive to ballistic. Thus, the gel-like arrest, in contrast to the glassy caging dynamics of the passive system, is driven by activity. The morphology of the arrested states, however, closely resembles gelation or arrested phase separation in passive colloids with attractive interactions [32,65]. Our analysis shows a clear crossover from the passive glassy dynamics to gel-like arrest at high densities [65].
We used the fact that our persistently active dynamics effectively causes one-dimensional motion against an interface, to map the late-time dynamics in the arrested states onto the Asymmetric Exclusion Process (ASEP). This microscopic mapping of the long time behavior of the system to a well-known lattice model allowed us to invoke well-established coarse graining procedures which we used to describe the non-trivial collective behaviour observed in our system. Building on this understanding, we used the predicted length scale to map out a non-equilibrium phase diagram that predicts non-trivial phases, a non-trivial topology, as well as non-trivial finite size scaling. While the position of the phase boundaries are not in quantitative agreement with our numerical results, the topology matches the numerical results over most of the phase space. The ASEP-based predictions fail at low activity and high densities because it takes into account only simple exclusion, and thus does not incorporate the shape-induced frustration and glassy dynamics in the passive system at high enough densities.
In this paper we make concrete predictions about the phase behavior and the morphology of arrested states when rotation of the active direction is strongly inhibited. This should occur naturally in collections of active entities with non-convex shapes [67]. Numerical simulations have explored MIPS [12] and glassy dynamics [56] in active dumbbells. In these studies, the mechanism of decorrelation of the active direction is, however, related to thermal diffusion [56] and not independently controlled. It would be interesting to explore dynamical arrest in simulations of rigid non-convex shapes where the rotation rate is affected by rotational locking. Experimental investigation of phase separation and arrest in collections of active colloidal particles with non-convex shapes offer possibilities of testing our theoretical predictions. We are not aware of any such experimental investigation, however, extending studies such as the Brownian dynamics of hard crosses [36] to include self propulsion seem feasible.

ACKNOWLEDGMENTS
We thank Mustansir Barma, Pinaki Chaudhuri, Chandan Dasgupta, Madan Rao, Abhishek Dhar, Eli Eisenberg and Andras Libal for helpful discussions. The work of CM and BC has been supported by the Brandeis MR-SEC (NSF-DMR 1420382). CM was hosted as graduate fellow at KITP, with support from the National Science Foundation under Grant No. NSF PHY-1748958 and the Heising-Simons Foundation. This work was also partially supported by the Israel Science Foundation Grant No. 968/16 and by a grant from the United States-Israel Binational Science Foundation.

Appendix A: Dynamics of Passive Hard Crosses
In this appendix, we provide a brief description of our analysis of slow dynamics and the "glass transition" in the passive hard cross system (∆v = 0) [31,39,40]. This analysis was performed to establish a baseline for the dynamics in the absence of activity. The initial states were created through a quench into the two-phase coexistence region using the Random Sequential Adsorption and Diffusion (RSAD) protocol [37]. The hard-crosses were then evolved according to the dynamics described in Section II in the main text, and two standard measures were used to probe the glassy dynamics [68,69]: (i) the self-intermediate scattering function (ISF) and the (ii) mean-squared-displacement (MSD) of individual crosses. The slowest relaxation occurs at the wave-vector with magnitude q at which the static structure factor has a peak. For N particles, the ISF is defined as where the wave vector q has magnitude q. The MSD is defined as: In a supercooled liquid, both the ISF and MSD should respect time-translational invariance, and should be independent of the waiting time, t w . From Figs. 12 and 13, this expectation is met for densities ρ < 0.1625. In this supercooled regime, we can fit the long-time decays of the ISF (Eq. A1), averaged over a range of t w from 200,000 to 800,000, to a stretched exponential form [70,71], F q (∆t) ∝ exp(−(∆t/τ α ) β ), as shown in Fig. 14  (a). The τ α extracted from this fit increases by two orders of magnitude over the density range 0.12 to 0.16 ( Fig. 14 (b)). Within our simulation time window of t max = 2 × 10 6 , passive hard crosses do not reach a time-translationally invariant state for densities ≥ 0.165: both the ISF and the MSD depend on t w , at these densities. We have established that in the aging regime the ISF exhibits a t/t w scaling, as shown in Fig. 15. Scalings of this form are characteristic of the aging regime in glassy systems [72][73][74].

Appendix B: Ensemble Averaged Density Distributions
The standard evidence for MIPS is based on measurements of ensemble-averaged density distributions. In this appendix, we present numerical results for these distributions over the full range of densities and activities. These ensemble averages include both arrested and nonarrested states. As in standard equilibrium phase separation, MIPS phase boundaries are drawn based on coexisting peaks in the density distributions, whose positions do not depend on the global density but which trade intensity as the global density is varied [66]. In our system, this expectation is met at low activities (Fig. 16) where we observe a continuation of the two-phase coexistence characteristic of the passive hard-crosses.
At higher activities, as seen in Fig. 17, the density distribution shows more structure. The ensemble averaged density distributions for the passive hard crosses exhibit a fluid-solid coexistence between ρ ≈ 0.16 and 0.19. In the presence of activity, the fluid peak at ρ ≈ 0.16 becomes broader and ultimately develops a peak at very low densities, characteristic of the voids seen in the arrested states. Before the void peak emerges clearly, in the solid-active liquid regime of the phase diagram, the positions of these low-density peaks are observed to depend both on activity and global density (c. f. Fig. 17 (c)-(d)). This observation is consistent with the fact that the phase separation is "arrested". As we have shown using the ASEP-based coarse-grained theory, the arrest of this phase separation leads to a wetting layer of width ξ within which the active liquid is confined. The density of particles in this active layer is a function of both activity and global density, as observed from the snapshots presented in the bottom row of Fig. 6.

Appendix C: Finite-Size Dependence of Phase Boundaries
The coarse-grained theory, based on ASEP that we have used to construct the non-equilibrium phase diagram [32] of arrested states at D R = 0, predicts that the phase boundaries depend explicitly on the system size L. In the infinite-system size limit, voids can form at any density and activity, according to this theory. We have performed limited studies of the size dependence of the phenomena we observe with the sole intent of checking whether the results of numerical simulations agree qualitatively with the predictions of the theory. In Fig 18 we compare ensemble averaged density distributions for system sizes L = 225, 450, and L = 600 at four different values of ∆v at ρ = 0.1500. The results show that the activity at which phase separation commences decreases with increasing L. For example, at ∆v = 0.030, the system is a homogeneous fluid at L = 225, 450, but there is fluid-solid coexistence for L = 600. Similarly, at ∆v = 0.05, the active-liquid solid phase is not observed at L = 225.
As mentioned above, the reason for the strong finite size effects is the ability of the system at D R = 0 to form voids at any activity and density. This feature is validated in the simulations, as seen in Fig 18, which shows regions that appear as low-density fluctuations in small systems transition to voids at larger L. These effects are in turn a consequence of the emergence of a single length scale, ξ, which depends only on ∆v, characterizing the density variation. When the system size is larger than ξ, voids open up to ensure mass conservation. Conversely, squeezing the system down to sizes smaller than ξ forces the interfaces to overlap. The ensemble averaged density distributions reflect this overlap through the development of a broad tail on the low density side of the solid peak in the the regime of solid-active liquid coexistence in finite systems. Since this phase disappears in the thermodynamic limit, it cannot be characterized by thermodynamic measures such as binodals derived from the peaks of the density distributions. As seen in Fig.  17, the shape of the distribution that interpolates between the void and the solid peaks depends on the global density and the system size (Fig. 18).   Color bars denote the stationary times of each cross, log 10 (τi). For this large activity value ∆v = 0.5, the percolating gel-like solid network can still form for active crosses with a finite rotation rate, suggesting that the rotational-locking mechanism of crosses greatly facilitates the global arrest of the system. For finite rotation rates, the rectangular empty void regions found for non-rotating crosses become filled with a finite density gas, and the boundaries of the low-density region become more rounded. The term "gas" here refers to very dilute regions and "liquid" is used to the describe higher densities (no significant physical difference is meant to be implied).