A particle-field approach bridges phase separation and collective motion in active matter

Whereas self-propelled hard discs undergo motility-induced phase separation, self-propelled rods exhibit a variety of nonequilibrium phenomena, including clustering, collective motion, and spatio-temporal chaos. In this work, we present a theoretical framework representing active particles by continuum fields. This concept combines the simplicity of alignment-based models, enabling analytical studies, and realistic models that incorporate the shape of self-propelled objects explicitly. By varying particle shape from circular to ellipsoidal, we show how nonequilibrium stresses acting among self-propelled rods destabilize motility-induced phase separation and facilitate orientational ordering, thereby connecting the realms of scalar and vectorial active matter. Though the interaction potential is strictly apolar, both, polar and nematic order may emerge and even coexist. Accordingly, the symmetry of ordered states is a dynamical property in active matter. The presented framework may represent various systems including bacterial colonies, cytoskeletal extracts, or shaken granular media.

Smooth particle representation and overlap calculation. A single rod at position r k with orientation ϕ k is considered to be an elliptical object which is represented by a smooth field ψ(r; r k , ϕ k ). This field can be rewritten as ψ(r; r k ,ϕ k ) = e −(r−r k )· 1+Q[ϕ k ] 4l 2 using the matrix Q[ϕ] = cos(2ϕ) sin(2ϕ) sin(2ϕ) − cos(2ϕ) (2) that encodes the orientation of a rod in two dimensions. Its eigenvectors e [ϕ] = (cos ϕ, sin ϕ) T and e ⊥ [ϕ] = (− sin ϕ, cos ϕ) T determine the major and minor axis, respectively, as a function of the polar angle ϕ. As in the main text, the dimensions of the major and minor axis are parametrized by l and l ⊥ . The Gaussianity of the field ψ facilitates the analytical calculation of the overlap I(r k , r j ; ϕ k , ϕ j ) = d 2 r ψ(r; r k ,ϕ k ) ψ(r; r j ,ϕ j ) of two rods. Due to the convolution structure of this integral, it may be solved via Fourier transformation methods yielding with the position-independent prefactor I 0 = I 0 (ϕ k , ϕ j ) = πl l ⊥ 1 − ε 2 1 − ε 2 cos 2 (ϕ k − ϕ j ) (5) which reflects the overlap of two rods whose centers of mass are identical. The mean nematic matrixQ(ϕ k , ϕ j ) = Interaction energy. The central question from the modeling point of view is how to derive the interaction from the overlap I. The basic idea is that spatially extended objects should repel if they come close to each other due to steric effects. We recall that the interaction energy was assumed to be binary, i.e. the total interaction energy is the sum of binary contributions: The interaction energy u 2 should be an increasing function of the overlap I. A natural assumption could therefore be the direct proportionality of the overlap to the interaction energy: u 2 ∝ I. Accordingly, pairs of particles which are far apart do not contribute to the interaction energy, whereas there is an energy cost associated to their mutual overlap as they approach each other. This model works well as long as overlaps remain small. However, spurious torques emerge in this model if large overlaps can occur, e.g. at very high particle densities or for high activity. To illustrate this issue, let us consider two rods whose centers of mass are close-by ( Supplementary Fig. 1); in particular, we compare the overlap for two particles whose body axis are aligned with a configuration where their orientations form a cross. It is self-evident that the cross configuration would correspond to the configuration of lower potential energy as the overlap is smaller than in the aligned situation. Physically, it is, however, questionable why two rods which slide past each other should align in a perpendicular fashion. Therefore, we argue that the ansatz u 2 ∝ I is questionable for driven nonequilibrium systems if significant overlaps are likely to occur. At the end of this section, we give the precise mathematical expression for these spurious torques for completeness. In order to avoid the emergence of spurious torques for close-by particles, we rethink the problem again. According to the Buckingham π theorem [1], the most general dependence of the binary interaction energy on the overlap is u 2 = κF[ξ], a large overlap in aligned configuration b reduced overlap in perpendicular arragement Supplementary Figure 1 | Sketch of two rods whose centers of mass are close-by. a The overlap (shown in color) of two nematically aligned rods is large. In contrast, the overlap is minimized by perpendicular alignment of the long body axis as illustrated in panel b. Therefore, spurious torques favoring perpendicular alignment would be implied if the interaction energy was chosen directly proportional to the overlap itself as this configuration corresponds to a state of minimal energy.
where κ is a constant with the units of energy and F[ξ] is a nonlinear function of the dimensionless variable ξ = I/a 0 , where I is the overlap and a 0 is a characteristic surface. We often use F[ξ] = ξ γ , with γ > 0 since we expect u 2 to be an increasing function of the overlap. Other functional forms may be used without affecting the observed large-scale behavior of the system qualitatively. It is also possible to implement a dependence F[ξ] that diverges for ξ > ξ * to ensure a maximal compressibility of the particles. Moreover, we want Onsager's argument [2] to be valid for a system of elongated passive particles represented by the proposed phase field, i.e. we expect the particles to get nematically aligned at the mean-field level locally. It turns out that the choice a 0 = I 0 is consistent with Onsager's argument. In summary, the interaction energy u 2 is expressed as Accordingly, the binary interaction energy reduces to a constant if two particles are close by (r j r k ) and, thus, the emergence of spurious torques is excluded.
Derivation of force and torque. In order to derive the force and torque, the interaction energy should be differentiated with respect to the position and orientation of a focal particle: F k = −∇ k U and M k = −∂ ϕ k U. As the interaction is binary by construction, we can write force and torque as a sum via where In order to keep the notation compact, we abbreviate the argument of the function F by Applying the chain rule, we thus obtain the intermediate result Performing the remaining differentiation eventually yields the expressions for force and torque. For the force, the differentiation is trivial: Accordingly, the pairwise force exerted by one rod on another is thus determined by where the kernel was abbreviated for convenience. For the torque, the differentiation of the nominator in Supplementary Eq. (11b) yields the collision avoidance term, whereas the nematic alignment follows from the differentiation of the denominator: Note that in the limit ε = 0, the torque vanishes identically and the force reduces to an isotropic central body force.
Spurious torques in the model u 2 =κI. We eventually give the expressions for the torque for a model where the binary interaction energy is directly proportional to the overlap, i.e.
The derivation is technically identical to the one of the previous paragraph. The binary torque in this model, which is obtained from m 2 (r, ϕ, ϕ ) = −∂ ϕ u 2 (r; ϕ, ϕ ), reads with the prefactorS These equations are the analogues to Supplementary Eqs. (14)-(15). The structural similarity of the corresponding terms is obvious. However, the prefactor in front of the nematic alignment term can turn negative, signaling the emergence of perpendicular alignment for This inequality is never true for large inter-particle distances. However, if the particles come close, the nematic alignment term changes sign. This reflects the emergence of spurious torques which were qualitatively discussed above along with Supplementary Fig. 1. As long as large overlaps (small inter-particle distances) cannot occur due to high repulsion or weak self-propulsion, spurious torques are irrelevant. However, qualitative differences are expected if particles can cross each other, e.g. because there is a high level of fluctuations, a high level of activity or the particle density is very large. If, on the other hand, the repulsion strength is very high, those configuration will never occur and both models become equivalent.

SUPPLEMENTARY NOTE 2: ADDITIONAL LANGEVIN SIMULATIONS
Phase diagrams -identification of motility-induced phase separation. Here, we provide additional quantitative results of numerical Langevin simulations that support the main findings reported in the main text. At first, we present phase diagrams for different aspect ratios to identify regions in parameter space where motility-induced phase separation emerges. We focus on the plane of density ρ 0 and activity/persistence, quantified here by the Péclet number defined as v 0 /( D ϕ ), with the characteristic length = √ l l ; note, however, that we fix l l ⊥ = 1 in all simulations. Different measures to identify MIPS were proposed in literature, such as the structure of the local density distribution (unimodal vs. bimodal) or cluster size distributions. In order to identify MIPS beyond the visual inspection of snapshots of simulations, we measured two order parameters that should capture central properties of MIPS. Motility-induced phase separation, as it is commonly described for self-propelled discs, is a phenomenon where a system spontaneously phase separates into a dense phase with local hexatic order and a low density, disordered gas phase. In short, two structural properties are relevant: (i) density fluctuations and (ii) the level of hexatic order. For the measurement of density fluctuations, the system is subdivided into N b boxes of edge length l b = 10 and the number of particles n j is counted in each box j. The expected number of particles in a box is ρ 0 l 2 b for the particle density ρ 0 . We calculate the standard deviation of the local particle number std(n) and construct an order parameter φ n by dividing by the standard deviation for a Poissonian point pattern that is given by N p(1 − p) with p = l 2 b /L 2 . Accordingly, this parameter equals one for a Poissonian point pattern and it is significantly larger than one for a strongly inhomogeneous distribution of particles; values smaller than one indicate regular arrangements of particles (e.g. a lattice). The level of density fluctuations in the density-activity plane is shown in Supplementary Fig. 2 for three different aspect ratios.
The measurement of the level of local hexatic order is based on a Voronoi tesselation [3]. The Voronoi neighbors Ω j for each particle j are identified first. Then, the complex order parameter is calculated for each particle, where α kj is the polar angle of the displacement vector r k − r j . The level of local hexatic order is given by |ψ j |. The breakdown of MIPS for anisotropic particles as shown in Fig. 2 of the main text is accompanied by the breakdown of local hexatic order: Supplementary Fig. 3 shows the same states like Fig. 2a,b, but the color coding corresponds to the local hexatic order parameter. In addition, we measured the distribution of the local particle number that reveals a bimodal structure for self-propelled spheres, whereas it is unimodal for rods. The average level of local hexatic order is quantified by averaging ψ j over particles all particles: These considerations eventually allow to construct an order parameter for the identification of motility-induced phase separation, that is characterized by the simultaneous presence of high density fluctuations and the emergence of local hexatic order, via the product of φ n and φ h6 . Corresponding phase diagrams are shown in Supplementary Fig. 4 for different aspect ratios in the density-activity plane. These quantitative phase diagrams reveal that motility-induced phase separation does only emerge for self-propelled spheres, whereas the order parameter is significantly reduced for self-propelled rods. A mechanistic explanation why MIPS breaks down is discussed in the main text along with analytical arguments.  Supplementary  Fig. 2). Color-coded is the order parameter φn × φ h 6 , that is the product of the order parameter for density fluctuations and the average level of local hexatic order. The characteristic parabolic parameter regime where MIPS emerges is observed for selfpropelled discs (panel a), whereas it is not observed for self-propelled rods (panels b and c).
Particle transport. To underline the different structural properties of the states that are illustrated in Fig. 2, we quantified the transport properties of particles by measuring the mean-squared displacement as a function of time ∆t = t − t 0 , graphically represented in Supplementary Fig. 5. As initial conditions, we chose the states that are shown in Fig. 2.
As commented in the context of Eq. (29) in the main text, one could, in principle, alternatively use the expressions for effective transport coefficients to quantify numerically the transport properties of the system, e.g. the collective diffusion coefficient Γ of the density field [4]. Note, however, that the expression for the collective diffusion coefficient was derived in the context of the linear stability analysis of the spatially homogeneous, isotropic state. This requires weak density fluctuations as well as the absence of both, local and global order since the local polar order parameter was adiabatically eliminated in this derivation.
The resulting expression for Γ is therefore, strictly speaking, only applicable in the disordered phase; it can, moreover, be used to detect the linear threshold to the emergence of MIPS [4]. However, as soon as local order emerges, the expression for Γ cannot be interpreted as the collective diffusion coefficient any longer, as the assumptions which led to its derivation are not longer justifiable. Thus, the applicability of this measure is limited to a small range of aspect ratios only (disc-like particles). For this reason, we decided to quantify the transport properties of the system via the mean-squared displacement at the single-particle level. For disordered states, where the collective diffusion coefficient can be measured numerically such as in the case of self-propelled discs, its measurement and the transport characteristics derived from the single-particle mean-squared displacement are expected to yield identical predictions.
The numerical study reveals a key difference between self-propelled discs and self-propelled rods: whereas the aggregates that are formed via MIPS rarely move -they are globally disordered and move diffusively with a diffusion coefficient that decreases as the size of the aggregate grows -polar clusters formed by self-propelled rods are highly dynamic and move ballistically. This is because polar ordering is inherently related to convective mass transport. In fact, polar clusters are becoming more persistent, i.e. ballistic, as they grow in size as rotational diffusion is suppressed for particles inside clusters. Technically speaking, MIPS is based on diffusion mediated transport -drops grow in size by exchanging particles/clusters via a diffusion process -while, in sharp contrast, the cluster dynamics of elongated active particles results from the interaction of polar cluster moving ballistically, with a persistence that increases with cluster size [5,6]. upplementary Figure 5 | Mean-squared displacement as a measure for collective transport properties. The simulation parameters correspond to the snapshots shown in Fig. 2 of the main text. Notably, the effective speed and persistence of the collective dynamics depends on the aspect ratio φ of particles. In particular, elongated self-propelled rods move highly persistent and ballistically.
Positional ordering and substructure of polar bands. In the main text, we highlighted that the polar band formed by self-propelled rods for high aspect ratios has a nontrivial substructure (cf. Fig. 2e): particles that are aligned tend to move next to each other, thereby forming smectic layers. To explain this type of ordering, we point out the fixed points of Eq. (3) that contains both, an alignment-like torque as well as a torque that couples the orientation of a rod to the relative position of its interaction partner. The latter interaction mechanism favors configurations where two rods arrange such that their body axis and the vector that connects the two centers of mass are perpendicular to each other. This torque is responsible for the formation of a non-trivial substructure which, in contrast, cannot be observed in Vicsek-like alignment models for self-propelled rods [7] where this type of coupling is absent.
We constructed two order parameters to quantify the formation of these substructures in more detail. Supplementary  Fig. 6 shows the same snapshot as Fig. 2e of the main text, but the color-coding corresponds to different order parameters that quantify the local structural properties as follows. Both order parameters are based on a Voronoi tesselation, i.e. local quantities are associated to particles by summing over Voronoi neighbors. Our analysis reveals that the band is composed of patches within which they form quadratic, anisotropic lattices that are revealed by the order parameter This definition is structurally related to the local hexatic order parameter discussed above. It is shown in Supplementary  Fig. 6a. In line with our argument above that the torqueφ i ∝ sin[2(ϕ i − α ij )] favors perpendicular alignment of the relative distance vector connecting the centers of mass of two rods and the direction of motion of a rod, one may define an effective interaction energy proportional to cos 2 (ϕ i − α ij ). We associate to each particle the local average of this quantity, as shown in Supplementary Fig. 6b that reveals a high level of positional ordering in almost all parts of the band.

SUPPLEMENTARY NOTE 3: KINETIC THEORY
Fokker-Planck dynamics of the one-particle density. In this appendix, we summarize some technical details on the kinetic theory discussed in the main text. The central observable linking particle-based descriptions and a field theoretical treatment is the one-particle density distribution which determines the density of particles at a particular reference point in phase space {r, ϕ}. Its dynamics can be derived from the N -particle Fokker-Planck equation [8] corresponding to the particle-based rod dynamics. In general, however, the dynamics of the one-particle density is not closed but it depends on the pair correlation function g 2 via the force and torque functionals: M (r, ϕ, t) = d 2 r dϕ m 2 (r − r , ϕ, ϕ )P (r , ϕ , t) g 2 (r, r ; ϕ, ϕ , t).
Thus, an appropriate closure scheme has to be found in order to address physical questions within kinetic theory. The explanation of the phenomena addressed in the main text -active phase separation and its breakdown as well as the emergence of polar order from nematic interactions -cannot be explained on the mean-field level (g 2 ≈ 1) but requires the consideration of inter-particle correlations as they are caused by the kinetics of collisions.
The spherical limit approximation -from particles to fields. In order to obtain the field theoretical equivalents of particle-based forces and torques, the kernels f 2 (r, ϕ, ϕ ) and m 2 (r, ϕ, ϕ ), cf. Supplementary Eqs. (13) and (15), have to be inserted into the corresponding interaction integrals above. In the spherical limit approximation, these expressions are substantially simplified: Several further simplifications are helpful. First, the one-particle density in the interaction integrals is expanded into a Taylor series in spatial gradients. Furthermore, an ansatz for g 2 has to be provided. We neglected its explicit temporal dependence and, further, reduced the number of variables of the pair-distribution using global symmetries, namely spatial homogeneity and rotational invariance: In the main text, we reported that the most relevant contribution to the two-particle distribution function is the kinetic enhancement of the particle density in front of a focal particle with respect to its direction of motion due to self-propulsion: as particles move actively in a semi-dilute environment, they tend to collide with others, and consequently the probability to find a particle in front is higher than in the back. We take this effects into account considering the corresponding, first Fourier component of g 2 : Using these approximations, we obtain the substantially simplified force and torque functionals with the density field ρ(r, t) = π −π dϕ P (r, ϕ, t) (32) and the transport coefficients The coefficients ζ i are positive if the density is enhanced in front of a particle with respect to its direction of motion as argued in the main text.
Derivation of Eqs. (15,16). In the main text, the transport coefficients σ p and σ n which determine the symmetry breaking of the polar and nematic order parameter at the local level, defined bẏ are discussed [Eqs. (15,16) in the main text]. Details on their derivation are summarized in this appendix. The derivation follows essentially the steps outlined in the corresponding Methods section: the temporal dynamics of the one-particle density distribution P (r, ϕ, t) given by Supplementary Eq. (26) is successively inserted into the definition of the corresponding order parameters: In two dimensions, there is, however, a more elegant and direct way to calculate both transport coefficients at once via Fourier transform with respect to the angular variable ϕ. For this purpose, we first introduce the angular Fourier modes of the one-particle density distribution as follows: f n (r, t) = 2π 0 dϕ e inϕ P (r, ϕ, t), P (r, ϕ, t) = 1 2π n e −inϕ f n (r, t).
Order parameters and angular Fourier modes are directly related via .
Since we are solely interested in the transport coefficients at the local level, gradient terms are neglected in the following discussion. Accordingly, only the terms involving torque and angular diffusion in the where denotes the angular Fourier modes of the torque functional [Supplementary Eq. (27b)] at the coarse-grained level. As a next step, we use symmetry properties of the model to simplify the torque functional M (r, ϕ, t) = d 2 r dϕ P (r , ϕ , t)m 2 (r − r , ϕ, ϕ )g 2 (r, r ; ϕ, ϕ , t) and the subsequent calculation of Fourier modes. At first, the one particle distribution function P (r , ϕ, t) can be expanded in Taylor series in r around the point r because we focus on the local level and, moreover, interactions are local as they are based on the contact of two particles: M (r, ϕ, t) ≈ d 2 r dϕ P (r, ϕ , t)m 2 (r − r , ϕ, ϕ ) g 2 (r, r ; ϕ, ϕ , t) + O(∇) (41a) = 1 2π n f n (r, t) d 2 r dϕ e −inϕ m 2 (r − r , ϕ, ϕ )g 2 (r, r ; ϕ, ϕ , t) + O(∇).
In parallel to the discussion in the main text, we introduce an effective mean-field torquem 2 = m 2 g 2 which combines the effects of the actual interaction m 2 and the collision kinetics, described by the pair correlation function g 2 . We refer to it as an effective mean-field model because structurally identical equations were obtained ifm 2 was the actual microscopic interaction in a system that is studied within mean-field approximation, where correlations are neglected (g 2 ≈ 1). Note, however, that in the present context correlations are included into the theory as we do not specify or presuppose a particular structure of the pair-correlation function. By virtue of rotational and translational invariance of the model, the number of variables ofm 2 can be reduced viã m 2 =m 2 (|r − r | , α − ϕ, ϕ − ϕ), where α = arg(r − r). In words, this reflects that the interaction does only depend on the relative position and relative orientation of one particle with respect to another one. This can be used to simplify the torque functional as follows: M (r, ϕ, t) = 1 2π n f n (r, t) As a last step, we use the fact that the interaction is anti-symmetric with respect to reflections at the axis that is parallel to the orientation of a focal particle:m 2 (r, −α, −ϕ ) = −m 2 (r, α, ϕ ). Accordingly, only the sin(nϕ)-terms will contribute to the integration: The dynamics is not closed but represents an infinite hierarchy of Fourier modes. In general, a closure relation is needed which depends on the structure of state or pattern which shall be described. The closure scheme will, however, only affect nonlinear terms. The relevant terms, which determine σ p and σ n that were discussed in the main text, readily follow from the summand m = n: These expressions are the basis for the discussion in the main text.