PAPER • OPEN ACCESS

# Vectorial active matter on the lattice: polar condensates and nematic filaments

To cite this article: Josué Manik Nava-Sedeño et al 2023 *New J. Phys.* **25** 123046

View the [article online](#) for updates and enhancements.

You may also like

- Computational approach to the structure spectral analysis of GERT-network models of mobile object monitoring systems
I V Kovalev, D I Kovalev, N A Testoyedov et al.

- Design and implementation of flight test parameter calibration data management system
Xue Bai, Ning Tang and Wei Li

- Mobile refrigeration system for precool and warm up of superconducting magnets
S K Gandla and R C Longsworth
<figure><img src="image_1.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

12 July 2023

REVISED

19 November 2023

ACCEPTED FOR PUBLICATION

12 December 2023

PUBLISHED

27 December 2023

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_2.png" /></figure>

PAPER

# Vectorial active matter on the lattice: polar condensates and nematic filaments

Josué Manik Nava-Sedeño¹,²,* ID, Haralampos Hatzikirou³ ID, Anja Voß-Böhme¹,⁴ ID, Lutz Brusch¹, Andreas Deutsch¹,* and Fernando Peruani²,³,* ID

¹ Technische Universität Dresden, Center for Information Services and High Performance Computing, Nöthnitzer Straße 46, 01062 Dresden, Germany

² Department of Mathematics, Faculty of Sciences, Universidad Nacional Autónoma de México, Circuito Exterior, Ciudad Universitaria, 04510 Mexico City, Mexico

³ Mathematics Department, Khalifa University, PO Box 127788, Abu Dhabi, United Arab Emirates

⁴ Fakultät Informatik/Mathematik, Hochschule für Technik und Wirtschaft, Dresden, Germany

⁵ Laboratoire de Physique Théorique et Modélisation, UMR 8089, CY Cergy Paris Université, 95302 Cergy-Pontoise, France

* Authors to whom any correspondence should be addressed.

E-mail: manikns@ciencias.unam.mx, andreas.deutsch@tu-dresden.de and fernando.peruani@cyu.fr

**Keywords:** velocity alignment, swarming, nematic order, cellular automaton

Supplementary material for this article is available online

## Abstract

We introduce a novel lattice-gas cellular automaton (LGCA) for compressible vectorial active matter with polar and nematic velocity alignment. Interactions are, by construction, zero-range. For polar alignment, we show the system undergoes a phase transition that promotes aggregation with strong resemblance to the classic zero-range process. We find that above a critical point, the states of a macroscopic fraction of the particles in the system coalesce into the same state, sharing the same position and momentum (polar condensate). For nematic alignment, the system also exhibits condensation, but there exist fundamental differences: a macroscopic fraction of the particles in the system collapses into a filament, where particles possess only two possible momenta. Furthermore, we derive hydrodynamic equations for the active LGCA model to understand the phase transitions and condensation that undergoes the system. We also show that generically the discrete lattice symmetries—e.g. of a square or hexagonal lattice—affect drastically the emergent large-scale properties of on-lattice active systems. The study puts in evidence that aligning active matter on the lattice displays new behavior, including phase transitions to states that share similarities to condensation models.

## 1. Introduction

Self-organized patterns of self-propelled entities are found at every scale in biology: from *in vitro* cytoskeletal extracts [1–3], bacterial systems [4–8], and insect swarms [9] to schools of fish [10], herds of mammals [11], and human crowds [12]. Collective phenomena are also observed in artificial active systems such as synthetic swimmers [13, 14], phoretic colloids [15, 16], and active rollers [17–19]. Active systems in the absence of torques—often referred as scalar active matter—can phase separate [20–22]. Other complex, self-organized, collective motion patterns exhibit either polar or nematic (orientational) order. Given the vectorial nature of orientational order, active systems sharing such a feature are often referred to as ‘vectorial active matter’ [23]. It is widely believed that self-organized pattern formation in vectorial active matter requires particles to interact via some type of velocity alignment mechanism [24, 25], though alternative mechanisms exist [26], such as combined attraction-repulsion and active forces [27]. The best known example of an active system with polar velocity alignment is the Vicsek model [28], where complex, active polar flows spontaneously emerge. A different class of active systems is obtained by changing the velocity alignment in the Vicsek model from polar to nematic [29, 30]. Self-propelled (point-mass) particles with nematic alignment do not generate polar flows, but self-organize into high-density, nematic bands [30]. The theoretical understanding
of vectorial active matter has been strongly based on the study of off-lattice Vicsek-like models by means of numerically costly, particle-based simulations [28, 30–34] and by the derivation of hydrodynamic equations for these systems [35–41]; for comprehensive reviews, see [24, 25]. There is, however, a renewed interest in understanding the behavior of active matter on lattices given the suitability of efficient numerical schemes, and in particular due to the large repertoire of available analytical tools to study lattice models that are expected to help to improve our theoretical understanding of active matter physics [42–47].

Here, we introduce an active lattice-gas cellular automaton (LGCA) for compressible (i.e. without volume exclusion) vectorial active matter, including both polar and nematic velocity alignment, where interactions are, by construction, zero-range. We show that the introduced LGCA with polar alignment undergoes a non-equilibrium phase transition that is fundamentally different from the non-local, incompressible active LGCA, where the onset of order is characterized by the presence of travelling polar bands [48] as occurs in the off-lattice Vicsek model [31]. In contrast, in the present model (with polar alignment), the emergence of order occurs via a process by which a macroscopic fraction of the particles coalesce into the same state, and thus sharing the same position and momentum, i.e. a polar condensate. For nematic alignment, we find that the system displays a condensation, but with a fundamental difference: a macroscopic fraction of the particles in the system collapses into a nematic filament. Note that recently, condensation behavior was found in scalar active matter in the presence of an external potential by a mechanism known edge-diffusivity [49, 50]. In scalar active matter—i.e. in self-propelled disks—it has been observed that phase separation into dense and dilute phases can result from volume exclusion interactions among the particles, a phenomenon known as motility-induced phase separation (MIPS; see [51] for a review). Our study shows that in vectorial active matter, condensation behavior is based on velocity alignment and does not require an external potential or external field, even in the absence of volume exclusion interactions. Thus, the emergence of polar and nematic order in the active LGCA is fundamentally different from the onset of order in other models [28, 31, 35, 42–47]. Similarly, phase separation in the active LGCA, that takes the form of a condensation, is fundamentally different from MIPS [51]. At high interaction strengths and low effective temperatures, these models show a great variety of collective behaviors such as homogeneous, polar-ordered states, formation of spatially-extended clusters, emergence of long-range order, filament-linked clusters, phase separation into extended jammed clusters, the formation of moving bands, and rippling; whereas our active LGCA models consistently show only two phases: a disordered, homogeneous state, and an aligned, highly condensed state in either single lattice sites, or in filaments of a single lattice site in width, depending on the nature of the alignment interaction. Furthermore, the hydrodynamic equations for the introduced LGCA allows us to show that generically the discrete lattice symmetries—at least for square and hexagonal lattices—have a strong impact on the emergent large-scale properties of on-lattice active systems. In summary, our study unveils the existence of condensed phases in vectorial active matter on the lattice, without attractive or excluded volume interactions.

## 2. Model definition

We study an hexagonal LGCA, where each node possesses six velocity channels,

$$ \mathbf{c}_{\vartheta} = \left( \cos \left[ (\vartheta - 1) \frac{\pi}{3} \right], \sin \left[ (\vartheta - 1) \frac{\pi}{3} \right] \right) $$

with $\vartheta \in [1, 2, \dots, 6]$, which can be occupied by a number $s_{\vartheta}$ of cells moving in the direction specified by $\mathbf{c}_{\vartheta}$. In contrast to standard passive [52] and active [48, 53] LGCA that consider that each velocity channel of a node can either be unoccupied or occupied by at most one cell, i.e. $s_{\vartheta} \in \{0, 1\}$—similar to the Pauli exclusion principle— here we assume particles without volume exclusion by letting any velocity channel $\mathbf{c}_{\vartheta}$ of a node contain an unrestricted number of particles $s_{\vartheta} \in \mathbb{N}$: i.e. many particles can share the same state (position and momentum). At each discrete time step, as illustrated in figure 1, two processes take place:

(i) Particles modify their velocity by transitioning to a channel $\vartheta$, according to the transition probability $T_{\vartheta}$:

$$ T_{\vartheta} = \frac{1}{Z} \exp[\beta H_{\vartheta}(\mathbf{r}, k)] \qquad (1) $$

with $H_{\vartheta}(\mathbf{r}, k)$ defined as:

$$ H_{\vartheta}(\mathbf{r}, k) = \sum_{\alpha=1}^{6} s_{\alpha}(\mathbf{r}, k) (\mathbf{c}_{\vartheta} \cdot \mathbf{c}_{\alpha})^{m} \qquad (2) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** LGCA dynamics. Particles reside in velocity channels within lattice sites. During the reorientation step (R), particles within each lattice site are assigned new velocity channels according to the transition probability $T_\vartheta(s)$. Subsequently, during the migration step (M), particles are deterministically moved to neighboring lattice sites.</figcaption></figure>

where $s(\mathbf{r},k) = [s_1(\mathbf{r},k), s_2(\mathbf{r},k), \dots, s_6(\mathbf{r},k)]$ is the configuration of node $\mathbf{r}$ at time $k$, which is given by the number of particles $s_\alpha(\mathbf{r},k)$ occupying channel $\alpha$ (with $\alpha \in [1, 2, \dots, 6]$), the strength of the velocity alignment is denoted by $\beta \in \mathbb{R}_+$, and $Z$ is the node normalization constant. Finally, $m=1$ corresponds to polar and $m=2$ to nematic velocity alignment. These transition probabilities are the volume exclusion-free analogues of those derived previously [54], which were obtained as discretized, steady state distributions of flocking off-lattice Langevin equation models with pairwise velocity alignment. Note that equation (1) is a Maxwell-Boltzmann distribution, where the parameter $\beta$ is inversely proportional to the effective temperature. Indeed, as seen in [54], $\beta$ is proportional to the interaction strength, and inversely proportional to the noise amplitude of off-lattice, stochastic differential equation models. Thus, in the limit $\beta \to 0$ (corresponding to infinitely high temperatures or noise amplitudes), equation (1) reduces to a discrete, uniform distribution, with particles exhibiting no preference in their reorientation. On the other hand, in the limit $\beta \to \infty$ (corresponding to an absolute zero temperature or no noise), $T_\vartheta = 1$ for a single orientation $\vartheta$ and $T_\vartheta = 0$ for all other orientations, resulting in completely deterministic behavior, for $m=1$, and two opposite orientations with $T_\vartheta = \frac{1}{2}$, and $T_\vartheta = 0$ for all other orientations, resulting in a one-dimensional, symmetric random walk along a single axis, for $m=2$. We stress that equations (1) and (2) define a purely local velocity alignment process: the reorientation dynamics does not consider the configuration of neighboring nodes. Note that the number of particles at node $\mathbf{r}$ at time $k$ is given by $n(\mathbf{r},k) = \sum_\alpha s_\alpha(\mathbf{r},k)$. This implies that equation (2) for $m=1$ can be expressed as $H(\mathbf{r},k)_\vartheta = n(\mathbf{r},k)\mathbf{P}(\mathbf{r},k) \cdot \mathbf{c}_\vartheta$ where $\mathbf{P}(\mathbf{r},k) = [\sum_\alpha s_\alpha(\mathbf{r},k)\mathbf{c}_\alpha]/n(\mathbf{r},k)$.

(ii) Particles move spatially according to the following rule: particles in node $\mathbf{r}$ at channel $\vartheta$ are transported into the same channel into the neighboring node $\mathbf{r} + \mathbf{c}_\vartheta$; see figure 1. It should be noted that lattice-gas cellular automata are mesoscopic, rather than microscopic, models. Therefore, time and space are coarse-grained. At this level of resolution, the reorientation processes that occur continuously in microscopic models are assumed to have reached a steady state at every discrete time step [54], and thus can only move at a fixed speed of one lattice site per time step.

### 3. Emergence of polar condensates

For polar alignment ($m=1$), we observe, in strong contrast to other polar active systems whose polar ordered phase consists of either a traveling polar band or an spatially homogeneous polar phase [28, 31, 35], the emergence of polar, moving condensates; see figure 2(a). In the following, we explain the mechanism that leads to the emergence of such condensates.

We make use of a mean-field argument and define the probability of finding a particle in channel $\vartheta$ of a certain node at time $k$ as $q_\vartheta(k) = s_\vartheta(k)/n(k)$, with $n(k) = \sum_{\alpha \in \Theta} s_\alpha(k)$ being the number of particles in the node. The node configuration, after reorientation, takes the form:

$$q_{\vartheta}(k+1) = \frac{1}{Z} \exp \left[ \beta n \sum_{\alpha \in \Theta} \cos(\alpha - \vartheta) q_{\alpha}(k) \right]. \quad (3)$$

Let us assume a steady state condition: $q_\vartheta(k+1) = q_\vartheta(k)$ for all $\vartheta$. Under this assumption, the transcendental equation (3) can be easily solved; see figures 2(b) and (c). Note that it is possible to recast equation (3) as $q_\vartheta = \exp[\beta n \mathbf{P} \cdot \mathbf{c}_\vartheta]/Z$ in order to obtain an equation for the polar order parameter of the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Dynamics of the polar, BA-LGCA. (a). Simulation snapshot; node occupancy is color coded, with dark color indicating large particle number. The inset shows the state of a node. See [55] for movies. (b). Global polar order $||\mathbf{P}||$ vs. average occupation per node $n$. (c). Probability of a particle escaping a polarized condensate given its particle number. Inset: Probability of a particle staying in a polarized cluster. (d) Histogram of the number of particles in all lattice sites at two different time steps. (e). Mean cluster size (blue) and polar order (red) as function of time; the dotted line corresponds to $t^{1/2}$, the characteristic growth in classic ZRPs [56], while the dashed line corresponds to $t^{3/4}$. Average over 200 simulations with $\beta = 4.3$ and $\bar{\rho} = 12$. (f). Polar order $P_g$ and condensate size $N_c$ as a function of $\beta\bar{\rho}$ for average $\bar{\rho} = 12$. (g). Global polar order $P_g$ and local polar order $P_\ell$ as a function of $\beta\bar{\rho}$ for average $\bar{\rho} = 12$. 150 simulations were averaged. A lattice size of 60 × 60 was used for all simulations. Parameter values were recorded after 1000 time steps.</figcaption></figure>

node: $\mathbf{P} = \sum_{\alpha \in \Theta} q_{\alpha} \mathbf{c}_{\alpha}$. And thus, $\mathbf{P} = \sum_{\alpha \in \Theta} \frac{1}{Z} \exp[\beta n \mathbf{P} \cdot \mathbf{c}_{\alpha}] \mathbf{c}_{\alpha}$. This transcendental equation describes a phase transition. The disorder state, given by $||\mathbf{P}|| = 0$, and corresponding to $q_{\vartheta} = 1/6$ for all $\vartheta$, is a solution of the above transcendental equation. This solution is stable for a number of particle $n$ below the critical $n^*$ value. This is evidenced by the behavior of $||\mathbf{P}||$ in figure 2(b). Above $n^*$, $q_{\vartheta} = 1/6$ is not stable and particles accumulate preferentially in one velocity channel and polar order in the node emerges. It is very important to realize that the probability of finding all particles in a node occupying the same velocity channel, for $n > n^*$ grows exponentially with the number of particles $n$ in the node. This implies that once all particles in a node occupy the same velocity channel, the probability that single particles transition to other velocity channels decays exponentially with $n$, as does the probability of not finding the $n$ particles in the same node. This mechanism leads to a condensation of particles. To understand this, assume, without loss of generality, that all particles are initially located in channel 0; then $q_0 = e^{\beta n}/Z = (1 + e^{\beta n(\gamma_1-1)} + \dots + e^{\beta n(\gamma_5-1)})^{-1}$, with $\gamma_j = \cos(-\pi j/3)$ so $\gamma_j - 1 < 0$. The probability that particles change channel is $1 - q_0 = \frac{e^{\beta n(\gamma_j-1)}}{[1 + \sum_{j=1}^5 e^{\beta n(\gamma_j-1)}]}$, which decays exponentially fast with $n$, while the probability of finding all particles again in channel 0 after one time step behaves as $(1 + e^{\beta n(\gamma_1-1)} + \dots + e^{\beta n(\gamma_5-1)})^{-n}$. Consequently, the probability for a particle in a condensate of mass $n \gg n^*$ to leave the condensate decays exponentially with $n$; figure 2(c). This feature is shared by the standard zero-range process (ZRP) [56]. It should be mentioned that, in this case, randomness decreases exponentially with particle number $n$ due to the pairwise, additive definition of $H_\psi(\mathbf{r}, k)$ (equation (2)), which increases linearly as more aligned particles are present. Conversely, other population models may exhibit similar fluctuations originating from multiplicative noise terms [57], which we do not consider in this model. Particle condensation can be evidenced by looking at the distribution of particles per node $P(n, t)$; figure 2(d) (see [58]). The second peak that shifts over time is indicative of an aggregation dynamics. Note that here, it is the alignment mechanism that holds particles
together. The emergent condensate moves ballistically in contrast to the processes studied in [56, 58]. This means that the condensation dynamics is not solely mediated by particle exchange with the bulk, but results from the motion of the condensate: as the condensate moves cross the system collects clusters in its way and grows in size, figures 2(d) and (e). This is further supported by the fact that the mean condensation time increases with system size (see figure S8 of the supplementary material), since clusters need to travel longer distances to accumulate particles in distant lattice sites. Note that, as argued in [59], in the absence of velocity correlations among clusters (and only in this scenario), polar order becomes intimately connected with the number and distribution of moving clusters. As expected, in figure 2(f) we observe that as particles coalesce and the mean cluster size grows with time, the global order parameter $P_g$—defined as $P_g = \| \sum_r \sum_{\alpha \in \Theta} q_\alpha(\mathbf{r})c_\alpha \|$ also does it. Furthermore, note that the mean cluster size scales with time with an exponent larger than 1/2; see [56]. Moreover, figure 2(f) shows that the condensate size $N_c$, defined as $N_c = \langle \max\{s_\vartheta(\mathbf{r}) : \vartheta \in \Theta, \mathbf{r} \in \mathcal{L}\} \rangle$, and $P_g$ exhibit almost identical dependency with $\beta\bar{\rho}$, with $\bar{\rho} = \langle n(\mathbf{r}) \rangle = N/L^2$, where N is the number of particles in a system of linear size L. An additional observable obtained from simulations was the local polar order, defined as $P_\ell = \sum_r \| \sum_{\alpha \in \Theta} q_\alpha(\mathbf{r})c_\alpha \| / N_\ell$, where $N_\ell$ is the number of nodes $r \in \mathcal{L}$ with at least one particle in any channel. Thus, $P_\ell$ is the average of the polar order within individual nodes, across all occupied nodes in the lattice. It is evident that, if particles are homogeneously distributed among all channels in every non-empty node, then $P_\ell = 0$. Conversely, if for every node, all particles are contained within a single channel, then $P_\ell = 1$. Figure 2(g) shows a plot of the global order $P_g$ and the local polar order $P_\ell$, as a function of the model parameters after several time steps. After the transition, the local polar order $P_\ell$ saturates up to its maximum value, while the global polar order $P_g$ exhibits non-monotonic behavior, reaching a maximum exactly at the transition point. This means that, after the transition, most occupied nodes accumulate particles within a single node. At the transition point, few condensates form and therefore the condensate with the most particles increases the overall global polar order. After the transition, several condensates with their own orientation form at once. The difference among the condensates' orientation drives the goal polar order down. Thus, disconnected condensates are present for all parameter values after the transition.

## 4. Emergence of nematic filaments

The strict local nature imposed by zero-range interactions seems to prevent nodes to exchange orientational information in order to establish extended ordered structures. This is true for polar alignment: the polar order of two neighboring nodes is not coupled. Counterintuitively, for nematic alignment, and in contrast to polar alignment, nodes can exchange orientational information and particles self-organize into stable, elongated, nematic filaments that lead to global nematic order; figure 3(a). We stress that these nematic filaments should not be confused with, and are different from the nematic bands observed in off-lattice models [30]: (i) the width of these filaments is always of one node and does not vary with the distance to the critical point as occurs in nematic bands [30], and (ii) as consequence of this, the transversal instability observed for nematic bands, precluding global nematic order, is not present [39, 60].

How can an extended nematic filament emerge from zero-range interactions? Picture an initial configuration where a large number $n$ of particles are placed on a single channel of a lattice site. As in the polar case, $q_\vartheta = 1/6$ is not stable for $n > n^*$ and particles accumulate in two opposite velocity channels. In the next time step, particles move to neighboring sites associated to these two velocity channels. If the density in each new lattice site is high, the same process occurs in the new sites, and thus information on the initial orientational order, propagates forward and backward, as shown in the sketch of figure 3(c). If we consider that this process starts taking place in an area surrounded by a disordered population of particles, then as the nematic filament expands enlarging its scattering cross-section, these particles feed the growing filament until it eventually percolates. The described mechanism leads to global nematic order as evidenced by the global nematic order parameter $|Q| = \left| \sum_r \sum_{\alpha \in \Theta} q_\alpha(\mathbf{r}) e^{2i\alpha} \right|$. Note that the emergent filament contains a macroscopic fraction of the particles in the system, as evidenced by the behavior of the filament density $N_f$ with $\beta$, which is similar to the one of $N_c$ for polar alignment; figure 3(d). Furthermore, we observe that for fixed $\bar{\rho}$, the $N$ grows as $L^2$, and thus $N_f$ as $L$, which implies that for $L \to \infty$, this quantity diverges: the nematic filament acts as a condensate with a macroscopic particle fraction located at it moving in two opposite directions. The local nematic order, $Q_\ell = \sum_r \left| \sum_{\alpha \in \Theta} q_\alpha(\mathbf{r}) e^{2i\alpha} \right| / N_\ell$ was also observed in simulations. Here as well, we observed a monotonic increase of the local nematic order and a non-monotonic behavior of the global nematic order; figure 3(b). As before, this indicates the formation of few filaments with most particles at the transition point, driving polar order up, and several filaments with different orientations far above the transition point, maintaining a high local nematic order, but inhibiting nematic order at a global level.
<figure><img src="image_5.png" /><figcaption>Figure 3. Dynamics of the nematic, BA-LGCA. (a). Snapshot of a simulation, with high node occupation numbers displayed in dark color; inset shows the state of a node. See [55] for movies. (b). Global nematic order |Q| and local nematic order Ql as a function of βρ for average ρ = 12. (c). Schematic depicts the filament formation process. Arrows indicate the direction of particle groups, with colors, from red to magenta, decreasing particle number. (d). Filament density vs. βρ averaged over 150 simulations for fixed ρ = 0.6. a Lattice of size 60 × 60 was used for all simulations. Parameter values were recorded after 1000 time steps.</figcaption></figure>

5. Hydrodynamic equations

The onset of the (macroscopic) symmetry breaking, while the system remains near the disordered phase, can
be described in an approximate manner by calculating the expected occupation numbers of lattice velocity
channels. However, note that well-developed condensates and filaments cannot be described by smooth PDE
solutions. As a first step, let us define $f_{\vartheta}(\mathbf{r}, k) := \langle s_{\vartheta}(\mathbf{r}, k) \rangle$. By performing an in-node mean-field
approximation and employing the binomial statistics of velocity channel occupation numbers,
one obtains the lattice-Boltzmann equation (LBE): $f_{\vartheta}(\mathbf{r} + \mathbf{c}_{\vartheta}, k + 1) = \rho(\mathbf{r}, k) T_{\vartheta}(\mathbf{r}, k)$, where
$\rho(\mathbf{r}, k) := \sum_{\vartheta \in \Theta} f_{\vartheta}(\mathbf{r}, k) = \langle n(\mathbf{r}, k) \rangle$ is the density, equal to the expected value of the total number of
particles in each lattice site, and $T_{\vartheta}(\mathbf{r}, k) := \frac{1}{Z} \exp[\beta \sum_{\alpha \in \Theta} f_{\alpha}(\mathbf{r}, k) (\mathbf{c}_{\vartheta} \cdot \mathbf{c}_{\alpha})^m]$ is the mean-field alignment
probability. We perform a discrete Fourier transform on $\vartheta$ to arrive to hydrodynamic fields that depend on
position and time. For polar alignment ($m=1$), the relevant hydrodynamic fields reduce to the local density
$\rho := \sum_{\vartheta \in \Theta} f_{\vartheta}$ and local polar order (or momentum) $\mathbf{M} := \sum_{\vartheta \in \Theta} f_{\vartheta} (\cos \vartheta, \sin \vartheta)$. The temporal evolution
of these fields, in dimensionless form, leading to unit active speed, is given by:

$$
\partial_t \rho + \nabla \cdot \mathbf{M} = 0 \qquad (4a)
$$

$$
\partial_t \mathbf{M} + \frac{1}{2} [\nabla \varrho - \mathcal{D}(\nabla \varrho)] + \frac{\beta}{\rho} \left[ (\mathbf{M} \cdot \nabla) \mathbf{M} - \frac{1}{2} \nabla M^2 + (\nabla \cdot \mathbf{M}) \mathbf{M} \right] = \rho \langle \mathbf{u} \rangle - \mathbf{M}, \quad (4b)
$$

where $\langle \mathbf{u} \rangle := \sum_{\vartheta \in \Theta} T_{\vartheta} (\cos \vartheta, \sin \vartheta)$, $\mathcal{D}(\nabla \rho) := \frac{1}{\rho^2} [\mathbf{M}(\mathbf{M} \cdot \nabla \rho) + \mathbf{M} \times (\mathbf{M} \times \nabla \rho)]$. For a detailed
derivation, see [55]. The term $\langle \mathbf{u} \rangle$ can be approximated, assuming $\beta \ll 1$, by a Taylor series. If we
approximate it up to third order in $\beta$, we find $\langle \mathbf{u} \rangle = \frac{\beta}{2} \mathbf{M}(1 - \frac{\beta^2}{8} \mathbf{M}^2)$, and by replacing this expression into
equation (4b), we obtain an equation that is invariant under arbitrary rotations of $\mathbf{M}$ as occurs in Toner–Tu
theory [35]. However, if the expansions is performed up to 5th order, we uncover that a steady solution of the
form $\mathbf{M} = M(\cos \theta, \sin \theta)$, with $M$ a scalar constant and $\theta$ an angle, cannot exist for arbitrary values of $\theta$. The
only possible solutions correspond to $\theta = \frac{a\pi}{3}$, with $a \in [0, 1, \dots, 5]$, i.e. to one of the possible 6 directions of the lattice; for details see [55].

Differences with the standard off-lattice Vicsek model can be shown even using the third order expansion in $\beta$ mentioned above. There are, as expected, two spatially homogeneous solutions, corresponding to $M=0$ (disorder) and $M>0$ (order). Performing a linear stability analysis around the disordered homogeneous steady state, and assuming perturbations along a single direction, we find that the system is stable against uniform perturbations for an initial average density $\beta\bar{\rho} < 2$, showing that the product $\beta\bar{\rho}$ is the effective bifurcation parameter. The second steady state corresponds to spatially homogeneous polar order state with $M = [8(1 - \frac{2}{\beta\bar{\rho}})]^{1/2}$. The ordered homogeneous state only exists when $\beta\bar{\rho} > 2$, where the disordered state becomes unstable. For non-uniform perturbations, the spatially homogeneous ordered state, that is never observed in the LGCA, is stable only for the parameter range $2.564 \lesssim \beta\bar{\rho} < 3$. Otherwise, the real parts of the eigenvalues grow with increasing wavenumber, indicating that small wavelength perturbations dominate, as confirmed by direct integration of the hydrodynamic equations. This instability indicates a drastic density increase in localized areas of the space, which can be associated to the emergence of condensates at this level of description. Numerical evidence strongly suggests that spatially homogeneous ordered states, and ordered polar bands are not possible in the LGCA model, even when forcing the system by using these configurations as initial conditions (see Supplementary Data).

For nematic alignment ($m=2$), the relevant hydrodynamic fields include nematic order, defined for simplicity as $Q := \sum_{\vartheta \in \Theta} f_{\vartheta}(\cos 2\vartheta + i \sin 2\vartheta)$, and the hydrodynamic equations read:

$$ \partial_t \rho = \mathfrak{D}\mathfrak{D}^* \rho + \mathfrak{R}(\mathfrak{D}^2 Q^*) \quad (5a) $$

$$ \partial_t Q = \frac{1}{2} \{ \mathfrak{D}^2 \rho + \mathfrak{D}\mathfrak{D}^* Q + 2\mathfrak{D}^* [\mathfrak{R}(\mathfrak{D}Q)] \} + \frac{\beta\rho}{4} \left[ Q + \frac{\beta}{8} Q^* \left( Q^* - \frac{\beta}{4} Q^2 \right) \right] - Q, \quad (5b) $$

where $\mathfrak{D} := \partial_x + i\partial_y$, * indicates complex conjugate, $\mathfrak{R}(\cdot)$ denotes real part, and fast relaxation of $M$ was assumed. It is worth stressing, as we explain below, that according to equation (5) nematic order is constrained to exist, as occurs for $M$, on the orientations allowed by the lattice structure. Equation (5) possesses spatially homogeneous disordered and nematically ordered solutions. The disordered state is stable when $\beta\bar{\rho} < 4$, where the bifurcation depends only on the product $\beta\bar{\rho}$, as before. The order in the spatially homogeneous steady state is given by $Q = (2 \pm 2\sqrt{9 - \frac{32}{\beta\bar{\rho}}}) e^{i\frac{2k\pi}{3}}$, $k \in \{0, 1, 2\}$. Two features are worth discussing. First, the ordered state exists when $\beta\bar{\rho} > \frac{32}{9}$, i.e. it exists even when the disordered state is stable, indicating a subcritical bifurcation. Second, while for $\frac{32}{9} < \beta\bar{\rho} < 4$ three pairs of solutions with different moduli, but same argument $\arg(Q) \in \{\frac{2k\pi}{3}, k \in \mathbb{N}\}$ exist, for $\beta\bar{\rho} > 4$, six different solutions exist, each one with a different argument corresponding to one of the lattice directions $\vartheta \in \Theta$. For $\beta\bar{\rho} > 4$, solutions with $\arg(Q) = \frac{2k\pi}{3}$ are stable, while those with $\arg(Q) = \frac{(2k+1)\pi}{3}$ are unstable. This is a consequence of the constraints imposed on $Q$ by the lattice, such that nematic structures can only emerge along the three hexagonal lattice axes, as observed in the LGCA simulations. Again, simulations suggest that spatially homogeneous ordered states, as well as spatially extended bands are not stable, and consistently break into nematic filaments, or revert back to the homogeneous disordered state. In this case, however, ordered states are much less unstable than in the polar case, and need long times in order to break apart (see Supplementary Data).

## 6. Discussion

The original LGCA was introduced with the intention of reproducing basic fluid mechanics properties. However, the initial study performed on a square lattice [61] displayed a series of anomalies due to the use a finite number of velocity channels. Later on, it was found that such anomalies vanish in hexagonal lattices [52], where the phenomenology of Navier-Stokes equations is recovered, despite the use of only 6 velocity channels. Similarly, it is tempting to speculate that an active LGCA with polar velocity alignment in an hexagonal lattice (and short-range interactions) has to reproduce the large-scale properties predicted by the Toner-Tu equations [35] for a (compressible) active polar fluid. After all, it is expected that most microscopic details are, at macroscopic scales, irrelevant except for the range and symmetry of the interactions. However, our study shows that spatially homogeneous ordered phases (cf Toner-Tu phase [35]) or polar bands (see [31, 46]) cannot exist in the LGCA model defined in this work. Moreover, in this LGCA, order emerges via a fundamentally different phase transition from the ones previously reported [28, 30, 31, 35, 46–48] in active matter. In LGCA, the emergence of global order is intrinsically related to a condensation
transition—that occurs in the absence of an external potential—which shares similarities with the ZRP. For the polar LGCA, it is associated to the formation of a polar macroscopic condensate, while for the nematic LGCA to the collapse of a macroscopic number of particles in the systems into a nematic filament (that leads to global nematic order; see [30, 39]).

After their introduction as fluid-dynamical models, LGCA have been used to model biological systems, specifically for modeling collective cell migration (see [62] for a review). The original LGCA formulation considers an exclusion principle, which limits the maximum number of cells which may be present at a single spatial point at a given time. On the one hand, our model could be considered as a 2D projection of a 3D system. For example, the accumulation of particles at certain lattice sites could correspond to the formation of an elongated structure by the piling up of agents, such as during the fruiting body and slug formation in *D. discoideum* [63]. On the other hand, cells are quite compressible, and are able to accumulate in great numbers. Experimentally, it has been found that cells can greatly reduce their volume during tissue growth [64, 65]. Cells also regularly change their size to maintain epithelial homeostasis [66]. For example, the photoreceptor epithelium shows a wide variety of receptor packing densities, corresponding to photoreceptor cells of varying sizes [67]. Furthermore, invasive cancer cells can be ten times as compressible as their healthy counterparts [68], and cell compressibility has been shown to correlate with invasive potential [69]. Moreover, it has been shown that densely packed cells confined to narrow environments can still move collectively, as long as they maintain sufficiently high adhesive contacts and velocity correlation levels [70]. The present LGCA model extension, which does not consider an exclusion principle, represents a first step towards a better mathematical representation of these cellular aggregation and compression processes.

Finally, we strongly suspect that the condensation reported here is related to the recent findings that indicate that here exist fundamental differences between additive and non-additive interactions in vectorial active matter [38, 71]. In particular, we speculate that the emergent macroscopic behavior can be affected, not only by modifying the symmetry or range of the interactions (see [48] for an analogous model with extended range and an exclusion principle), as it is usually expected, but by eliminating the dependency on $n(\mathbf{r}, k)$ in equation (2), for instance by defining for $m = 1$, $H(\mathbf{r}, k)_\vartheta = \mathbf{P}(\mathbf{r}, k) \cdot \mathbf{c}_\vartheta$. This will be the focus of future investigations.

## Data availability statement

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

## Acknowledgments

F P acknowledges financial support from CY Initiative through the Ambition project ‘Collective Intelligence’ and Labex MME-DII. J M N S acknowledges support from the PAPIIT-UNAM Grant, Project IA104821.

## Author contributions

J M N S: Formal analysis, investigation, methodology, software, visualization, writing (original draft), and writing (review and editing). H H: Funding acquisition, supervision, and writing (review and editing). A V B and L B: supervision, methodology, and writing (review and editing). A D: Funding acquisition, supervision, project administration, and writing (review and editing). F P: Conceptualization, formal analysis, investigation, methodology, project administration, supervision, and writing (review and editing).

## Conflict of interest

There are no conflicts to declare.

## ORCID iDs

Josué Manik Nava-Sedeño https://orcid.org/0000-0001-9451-8676

Haralampos Hatzikirou https://orcid.org/0000-0002-1270-7885

Anja Voß-Böhme https://orcid.org/0000-0001-6510-8694

Fernando Peruani https://orcid.org/0000-0003-2985-8935

## References

[1] Schaller V, Weber C, Semmrich C, Frey E and Bausch A R 2010 Polar patterns of driven filaments *Nature* **467** 73–77
[2] Sanchez T, Chen D T N, DeCamp S J, Heymann M and Dogic Z 2012 Spontaneous motion in hierarchically assembled active matter *Nature* **491** 431–4

[3] Nagai K H, Sumino Y, Montagne R, Aranson I S and Chaté H 2015 Collective motion of self-propelled particles with memory *Phys. Rev. Lett.* **114** 168001

[4] Dombrowski C, Cisneros L, Chatkaew S, Goldstein R E and Kessler J O 2004 Self-concentration and large-scale coherence in bacterial dynamics *Phys. Rev. Lett.* **93** 098103

[5] Sokolov A, Aranson I S, Kessler J O and Goldstein R E 2007 Concentration dependence of the collective dynamics of swimming bacteria *Phys. Rev. Lett.* **98** 158102–4

[6] Peruani F, Starruß J, Jakovljevic V, Søgaard-Andersen L, Deutsch A and Bär M 2012 Collective motion and nonequilibrium cluster formation in colonies of gliding bacteria *Phys. Rev. Lett.* **108** 098102

[7] Ariel G, Rabani A, Benisty S, Partridge J D, Harshey R M and Be'er A 2015 Swarming bacteria migrate by Lévy walk *Nat. Commun.* **6** 1–6

[8] Be'er A and Ariel G 2019 A statistical physics view of swarming bacteria *Mov. Ecol.* **7** 9

[9] Attanası A *et al* 2014 Finite-size scaling as a way to probe near-criticality in natural swarms *Phys. Rev. Lett.* **113** 238102

[10] Ward A J W, Sumpter D J T, Couzin I D, Hart P J B and Krause J 2008 Quorum decision-making facilitates information transfer in fish shoals *Proc. Natl Acad. Sci. USA* **105** 6948–53

[11] Ginelli F, Peruani F, Pillot M-H, Chaté H, Theraulaz G and Bon R 2015 Intermittent collective dynamics emerge from conflicting imperatives in sheep herds *Proc. Natl Acad. Sci. USA* **112** 12729–34

[12] Moussaïd M, Helbing D and Theraulaz G 2011 How simple rules determine pedestrian behavior and crowd disasters *Proc. Natl Acad. Sci. USA* **108** 6884–8

[13] Palacci J, Sacanna S, Steinberg A P, Pine D J and Chaikin P M 2013 Living crystals of light-activated colloidal surfers *Science* **339** 936–40

[14] Theurkauff I, Cottin-Bizonne C, Palacci J, Ybert C and Bocquet L 2012 Dynamic clustering in active colloidal suspensions with chemical signaling *Phys. Rev. Lett.* **108** 268303

[15] Buttinoni I, Bialké J, Kümmel F, Löwen H, Bechinger C and Speck T 2013 Dynamical clustering and phase separation in suspensions of self-propelled colloidal particles *Phys. Rev. Lett.* **110** 238301

[16] Golestanian R, Liverpool T B and Ajdari A 2007 Designing phoretic micro- and nano-swimmers *New J. Phys.* **9** 126

[17] Bricard A, Caussin J-B, Desreumeaux N, Dauchot O and Bartolo D 2013 Emergence of macroscopic directed motion in populations of motile colloids *Nature* **503** 95–98

[18] Bricard A, Caussin J-B, Das D, Savoie C, Chikkadi V, Shitara K, Chepizhko O, Peruani F, Saintillan D and Bartolo D 2015 Emergent vortices in populations of colloidal rollers *Nat. Commun.* **6** 1–8

[19] Kaiser A, Snezhko A and Aranson I S 2017 Flocking ferromagnetic colloids *Sci. Adv.* **3** e1601469

[20] Barberis L and Peruani F 2019 Phase separation and emergence of collective motion in a one-dimensional system of active particles *J. Chem. Phys.* **150** 144905

[21] Digregorio P, Levis D, Suma A, Cugliandolo L F, Gonnella G and Pagonabarraga I 2018 Full phase diagram of active Brownian disks: from melting to motility-induced phase separation *Phys. Rev. Lett.* **121** 098003

[22] Fily Y and Marchetti M C 2012 Athermal phase separation of self-propelled particles with no alignment *Phys. Rev. Lett.* **108** 235702

[23] Gompper G *et al* 2020 The 2020 motile active matter roadmap *J. Phys.: Condens. Matter* **32** 193001

[24] Vicsek T and Zafeiris A 2012 Collective motion *Phys. Rep.* **517** 71–140

[25] Marchetti M, Joanny J, Ramaswamy S, Liverpool T, Prost J, Rao M and Simha R 2013 Hydrodynamics of soft active matter *Rev. Mod. Phys.* **85** 1143–89

[26] Barberis L and Peruani F 2016 Large-scale patterns in a minimal cognitive flocking model: incidental leaders, nematic patterns and aggregates *Phys. Rev. Lett.* **117** 248001

[27] Caprini L and Löwen H 2023 Flocking without alignment interactions in attractive active Brownian particles *Phys. Rev. Lett.* **130** 148202

[28] Vicsek T, Czirók A, Ben-Jacob E, Cohen I and Shochet O 1995 Novel type of phase transition in a system of self-driven particles *Phys. Rev. Lett.* **75** 1226–9

[29] Peruani F, Deutsch A and Bär M 2008 A mean-field theory for self-propelled particles interacting by velocity alignment mechanisms *Eur. Phys. J.: Spec. Top.* **157** 111–22

[30] Ginelli F, Peruani F, Bär M and Chaté H 2010 Large-scale collective properties of self-propelled rods *Phys. Rev. Lett.* **104** 184502

[31] Grégoire G and Chaté H 2004 Onset of collective and cohesive motion *Phys. Rev. Lett.* **92** 025702

[32] Chaté H, Ginelli F and Montagne R 2006 Simple model for active nematics: quasi-long-range order and giant fluctuations *Phys. Rev. Lett.* **96** 180602

[33] Chaté H, Ginelli F, Grégoire G, Peruani F and Raynaud F 2008 Modeling collective motion: variations on the Vicsek model *Eur. Phys. J. B* **64** 451–6

[34] Kürsten R and Ihle T 2020 Dry active matter exhibits a self-organized cross sea phase *Phys. Rev. Lett.* **125** 188003

[35] Toner J and Tu Y 1995 Long-range order in a two-dimensional dynamical XY model: how birds fly together *Phys. Rev. Lett.* **75** 4326–9

[36] Bertin E, Droz M and Grégoire G 2006 Boltzmann and hydrodynamic description for self-propelled particles *Phys. Rev. E* **74** 022101

[37] Ihle T 2011 Kinetic theory of flocking: derivation of hydrodynamic equations *Phys. Rev. E* **83** 030901

[38] Kürsten R and Ihle T 2021 Quantitative kinetic theory of flocking with three-particle closure *Phys. Rev. E* **104** 034604

[39] Peshkov A, Aranson I S, Bertin E, Chaté H and Ginelli F 2012 Nonlinear field equations for aligning self-propelled rods *Phys. Rev. Lett.* **109** 268701

[40] Peshkov A, Bertin E, Ginelli F and Chaté H 2014 Boltzmann-Ginzburg-Landau approach for continuous descriptions of generic Vicsek-like models *Eur. Phys. J.: Spec. Top.* **223** 1315–44

[41] Großmann R, Schimansky-Geier L and Romanczuk P 2013 Self-propelled particles with selective attraction-repulsion interaction: from microscopic dynamics to coarse-grained theories *New J. Phys.* **15** 085014

[42] Peruani F, Klauss T, Deutsch A and Voss-Boehme A 2011 Traffic jams, gliders and bands in the quest for collective motion of self-propelled particles *Phys. Rev. Lett.* **106** 128101

[43] Partridge B and Lee C F 2019 Critical motility-induced phase separation belongs to the Ising universality class *Phys. Rev. Lett.* **123** 068002

[44] Kourbane-Houssene M, Erignoux C, Bodineau T and Tailleur J 2018 Exact hydrodynamic description of active lattice gases *Phys. Rev. Lett.* **120** 268003
[45] Manacorda A and Puglisi A 2017 Lattice model to derive the fluctuating hydrodynamics of active particles with inertia *Phys. Rev. Lett.* **119** 208003

[46] Solon A P and Tailleur J 2015 Flocking with discrete symmetry: the two-dimensional active Ising model *Phys. Rev. E* **92** 042119

[47] Nesbitt D, Pruessner G and Lee C F 2021 Uncovering novel phase transitions in dense dry polar active fluids using a lattice Boltzmann method *New J. Phys.* **23** 043047

[48] Bussemaker H J, Deutsch A and Geigant E 1997 Mean-field analysis of a dynamical phase transition in a cellular automaton model for collective motion *Phys. Rev. Lett.* **78** 5018

[49] Golestanian R 2019 Bose-Einstein-like condensation in scalar active matter with diffusivity edge *Phys. Rev. E* **100** 010601

[50] Meng F, Matsunaga D, Mahault B and Golestanian R 2021 Magnetic microswimmers exhibit Bose-Einstein-like condensation *Phys. Rev. Lett.* **126** 078001

[51] Cates M E and Tailleur J 2015 Motility-induced phase separation *Annu. Rev. Condens. Matter Phys.* **6** 219–44

[52] Frisch U, Hasslacher B and Pomeau Y 1986 Lattice-gas automata for the Navier-Stokes equation *Phys. Rev. Lett.* **56** 1505

[53] Deutsch A and Dormann S 2005 *Cellular Automaton Modeling of Biological Pattern Formation* (Springer)

[54] Nava-Sedeño J, Hatzikirou H, Peruani F and Deutsch A 2017 Extracting cellular automaton rules from physical Langevin equation models for single and collective cell migration *J. Math. Biol.* **75** 1075–100

[55] See supplemental material for movies of the active LGCA and for further technical details on the derivation of the hydrodynamic field equations

[56] Evans M R and Hanney T 2005 Nonequilibrium statistical mechanics of the zero-range process and related models *J. Phys. A: Math. Gen.* **38** R195

[57] Scheffer M, Van Bavel B, van de Leemput I A and van Nes E H 2017 Inequality in nature and society *Proc. Natl Acad. Sci. USA* **114** 13154–7

[58] Majumdar S N, Krishnamurthy S and Barma M 1998 Nonequilibrium phase transitions in models of aggregation, adsorption and dissociation *Phys. Rev. Lett.* **81** 3691

[59] Peruani F and Bär M 2013 A kinetic model and scaling properties of non-equilibrium clustering of self-propelled particles *New J. Phys.* **15** 065009

[60] Großmann R, Peruani F and Bär M 2016 Mesoscale pattern formation of self-propelled rods with velocity reversal *Phys. Rev. E* **94** 050602

[61] Hardy J, De Pazzis O and Pomeau Y 1976 Molecular dynamics of a classical lattice gas: transport properties and time correlation functions *Phys. Rev. A* **13** 1949

[62] Deutsch A, Nava-Sedeño J M, Syga S and Hatzikirou H 2021 BIO-LGCA: a cellular automaton modelling class for analysing collective cell migration *PLoS Comput. Biol.* **17** e1009066

[63] Marée A F M and Hogeweg P 2001 How amoeboids self-organize into a fruiting body: multicellular coordination in *Dictyostelium discoideum* *Proc. Natl Acad. Sci.* **98** 3879–83

[64] Puliafito A, Hufnagel L, Neveu P, Streichan S, Sigal A, Fygenson D K and Shraiman B I 2012 Collective and single cell behavior in epithelial contact inhibition *Proc. Natl Acad. Sci.* **109** 739–44

[65] Aland S, Hatzikirou H, Lowengrub J and Voigt A 2015 A mechanistic collective cell model for epithelial colony growth and contact inhibition *Biophys. J.* **109** 1347–57

[66] Gradeci D, Bove A, Vallardi G, Lowe A R, Banerjee S and Charras G 2021 Cell-scale biophysical determinants of cell competition in epithelia *eLife* **10** e61011

[67] Jiao Y, Lau T, Hatzikirou H, Meyer-Hermann M, Corbo J C and Torquato S 2014 Avian photoreceptor patterns represent a disordered hyperuniform solution to a multiscale packing problem *Phys. Rev. E* **89** 022721

[68] Swaminathan V, Mythreye K, O'Brien E T, Berchuck A, Blobe G C and Superfine R 2011 Mechanical stiffness grades metastatic potential in patient tumor cells and in cancer cell lines *Cancer Res.* **71** 5075–80

[69] Cheburkanov V, Pendyala K, Parlani M, Lele T, Friedl P and Yakovlev V 2022 Imaging mechanical properties of cancer cells during metastasis with Brillouin microspectroscopy *Proc. SPIE* **11944** 71–79

[70] Ilina O et al 2020 Cell–cell adhesion and 3D matrix confinement determine jamming transitions in breast cancer invasion *Nat. Cell Biol.* **22** 1103–15

[71] Chepizhko O, Saintillan D and Peruani F 2021 Revisiting the emergence of order in active matter *Soft Matter* **17** 3113–20