Circling crystals in chiral active matter with self-alignment

Abstract

We study a crystal composed of active units governed by self-alignment and chirality. The first mechanism acts as an effective torque that aligns the particle orientation with its velocity, while the second drives individual particles along circular orbits. We find that even a weak degree of chirality, when coupled with self-alignment, induces collective motion of the entire crystal along circular trajectories in space. We refer to this phase as a circling crystal. When chirality outweighs self-alignment, the circular global motion is suppressed in favor of vortex-like regions of coordinated motion. This state is characterized by oscillating spatial velocity correlations, a power law decay of the energy spectrum, and oscillatory temporal correlations. Our findings can be tested experimentally in systems ranging from epithelial tissues to swarming robots, governed by chirality and self-alignment.

---

I. Introduction

In the last few years, active matter^1–3 has emerged as an exciting new field for applying and testing concepts of non-equilibrium physics. Active systems are typically composed of self-propelled units that convert energy from internal or external sources into directed motion. These units often tend to follow circular or helical orbits, systematically rotating clockwise or counterclockwise. These kinds of systems are often termed chiral^4,5 because their trajectory cannot be superimposed by their mirror images. Examples of chiral active matter are widespread in biology: they range from living microswimmers, such as bacteria^6–9 moving on a substrate to malarial parasites,¹⁰ and chiral microtubules driven by molecular motors,¹¹ sperm cells¹² and small cell clusters confined in circular islands.¹³ Chirality has also been implemented in the design of artificial systems by breaking the rotational symmetry in the body shape or in the self-propulsion mechanisms, as in the case of L-shaped colloids,^14,15 isotropic droplets swimming in circles due to chiral ordering and geometrical frustration,^16,17 or chiral active granular particles^18–27 often termed spinners.^28–31

The impact of chirality on single-particle properties, such as odd diffusion^32–35 and edge currents under external potentials,^36,37 has been extensively explored, while recent efforts have focused on chirality-induced collective phenomena.^38–43 Both theoretical analyses and numerical simulations have shown that adding chirality suppresses motility-induced phase separation,^44–46 leading to the formation of hyperuniform fluids^47–49 and the emergence of intriguing structures such as rotating vortices in the velocity field⁵⁰ or self-reverting vorticity.⁵¹ Introducing alignment interactions in chiral systems^52–56 leads to an enhanced flocking behavior and pattern formation⁵⁷ or micro-flock pattern formation⁵⁸ due to the fact that the interplay between self-rotation and alignment fosters synchronization among particles and stabilizes large-scale coherent chiral motion. A recent line of research focuses on chiral active crystals, where each fundamental unit of the solid consists of chiral active particles. In this case, chirality reduces the correlation length of the spatial velocity correlations,^59,60 typical of active matter and induces cross-spatial correlations among different velocity components, resulting in a net angular momentum.⁶⁰ Moreover, odd-elastic or chiral crystals exhibit spontaneous grain rotation, autonomous propulsion, and dynamic fragmentation.⁶¹ In addition, the chirality-induced suppression of displacement fluctuations associated with hyperuniformity generates truly long-translational order even in two dimensions.⁶²

Several experimental systems, such as cells,⁶³ epithelial tissues,^64–66 and active granular particles,^25,67–69 are characterized by a coupling between translational and rotational motion commonly referred to as self-alignment.⁷⁰ This mechanism acts as an effective torque that tends to align the particle's orientation with its velocity.^71–77 Self-alignment is one of the key mechanisms that can lead to flocking behavior^74,76 even if this is a single-particle mechanism that does not couple the orientation of different particles.

The influence of self-alignment on collective phenomena has been recently investigated in several studies. Specifically, self-alignment suppresses motility-induced phase separation in favor of a homogeneous flocking phase,⁷⁶ and it is responsible for a re-entrant glass transition⁷⁴ due to emergent migratory patterns. In addition, recent studies focus on the case of self-aligning active crystals, where each unit of the solid is subject to a self-aligning activity. While, in the absence of self-alignment, active crystals excite collective modes at every wave vector,⁷⁸ self-alignment activates a few elastic modes, as proved experimentally⁷³ and theoretically.⁷⁹ In this case, the system displays flocking behavior which can be analytically predicted and interpreted as a continuous phase transition described by an effective velocity-dependent free energy with a Landau–Ginzburg profile.⁸⁰

However, the interplay between chirality and self-alignment remains scarcely explored. It is still unclear if the competition between these two mechanisms influences the collective dynamics of active systems and generates different collective behavior.

In this study, we explore the competition between self-alignment and chirality on collective behaviors in dense active systems, with crystalline order. We discover that even a weak chirality introduces a global rotation, giving rise to a new collective phenomenon, which we term a circling crystal. In this phase, all the particles in the system are able to coordinate their motion and move collectively along circular trajectories. By contrast, strong chirality suppresses the global circling in favor of vortex structures in the velocity field. Here, the system is characterized by oscillating spatial velocity correlations with a finite swirl size and a power-law decaying energy spectrum, while the velocity auto-correlations show a time oscillating profile.

The structure of the paper is as follows: in Section II, we present the model used to perform the simulations of a solid composed of active units characterized by self-alignment and chirality. The corresponding results are discussed in Section III. Finally, in Section IV, we provide a general discussion together with possible perspectives, applications, and directions for future studies and experimental validations. The details of the model employed in the numerical simulations are presented in Appendix A, while the derivations supporting our theoretical predictions are given in Appendix B.

II. Model

We consider a two-dimensional solid consisting of N interacting chiral active Brownian particles subject to a self-alignment mechanism. Particles are active and they move with a constant velocity v₀ along their intrinsic orientation n̂ (black arrow in Fig. 1(a)), indicated by the black cap on each particle, and are chiral since their orientation rotates at constant frequency, Ω (purple circular arrow in Fig. 1(a)), also termed chirality. This rotating mechanism pushes a free particle to move along a circular trajectory whose typical radius is determined by the interplay between the chirality strength and the particle's self-propulsion, r_Ω = v₀/Ω. In this system, the crystal units are subject to a self-alignment mechanism, modeled as an effective torque, T_sa (blue circular arrow in Fig. 1(a)), that tends to align the particle orientation, n̂, with its velocity v (red arrow in Fig. 1(a)). The particle motion is described by an inertial translational dynamics and an overdamped rotational dynamics: where x_i is the position of the i particle, v_i = ẋ_i is its velocity, and θ_i is the orientation angle setting the direction of the active force. The parameters γ and γ_r represent the translational and rotational friction coefficients, respectively, while D_t and D_r are the corresponding diffusion coefficients. The terms ξ_i and η_i denote independent Gaussian white noises with zero mean and unit variance. Each particle, of mass m, is self-propelled at a constant speed along the direction n̂_i = (cos θ_i,sin θ_i), determined by its orientation angle θ_i. Particle motion is influenced by a chiral term, Ω, which induces circular orbits, and by a deterministic torque, T^sa_i = β(n̂_i × v_i), that aligns the orientation n̂_i with the velocity v_i. The self-alignment strength is set by the parameter β, which defines the typical distance, γ_r/β, that a particle travels before its orientation aligns with its velocity.

Fig. 1 Schematic of dynamics of chiral self-aligning active solids. (a) Schematic illustration of an active solid composed of active particles characterized by self-alignment and chirality. Here, the particle orientation, indicated by the black cap on the particle, is represented by n̂, while v denotes its velocity. The rotation direction, i.e. the particle chirality, is indicated through a purple circular arrow, while the self-alignment torque T_sa by a blue circular arrow. (b)–(e) Trajectories of a target particle within the active solid as the strengths of self-alignment and chirality are varied. The color along the trajectory encodes time, with darker points representing earlier positions and lighter points indicating later positions. The points along each trajectory are connected by gray lines to highlight the particle's path. The simulation results presented here were obtained for reduced self-alignment values of B = 0.1 and B = 0.5, and reduced chirality values of ω = 0.5 and ω = 40. The other dimensionless parameters of the simulations are: Pe = 10, M = 10⁻⁴, , and Φ = Nπσ²/4L² = 1.1.

The force F_i accounts for purely repulsive interactions between particles and is derived from a Weeks–Chandler–Andersen (WCA) potential: F_i = −∇_iU_tot, with and U(r) = 4ε[(σ/r)¹² − (σ/r)⁶] if r < 2^1/6σ and zero otherwise. Here, ε sets the energy scale, and σ denotes the particle diameter. We consider a regime of high packing fraction ϕ = Nσ²π/4L² ≈ 1.1, where L is the side of the simulation box. In two dimensions, the high value of ϕ forces the particles to occupy the vertices of a triangular lattice (Fig. 1(a)) arranging in an almost perfect crystal-like configuration. As a consequence, they cannot exchange positions with their neighbors. Each individual unit is confined in a cage formed by the nearest particles, and the only way a particle can move beyond the typical interparticle distance is through collective motion. The high density considered in our work allows us to reproduce conditions relevant to many biological systems, such as epithelial tissues⁸¹ or chiral cellular aggregates¹³ while avoiding effects like density fluctuations^82,83 or density inhomogeneities.

In our simulations, we set D_t = 0, as in active matter systems this term is typically much smaller than the effective diffusion arising from activity. Chirality introduces a characteristic time for a particle to complete a circular orbit in the absence of other effects, given by 1/Ω. Self-alignment defines another characteristic time, γ_r/(βv₀), i.e., the time required for the particle's orientation to align with its velocity. These times compete with other intrinsic time scales of the system, namely the persistence time of the particle, τ = 1/D_r, and the translational inertial time, τ_d = m/γ.

Simulations are performed in a box of side L with periodic boundary conditions by rescaling lengths by the particle diameter σ and time by τ. With these choices, the dynamics are governed by several dimensionless parameters: the Péclet number, Pe = v₀/(D_rσ), which quantifies the persistence length of a particle relative to its size; the reduced mass, M = D_rm/γ, which determines the relevance of inertia and is chosen to be small, M = 10⁻⁴; the reduced chirality, ω = Ωτ, which compares the typical time for a particle to complete a circular orbit with the persistence time; the reduced self-alignment strength, B = βσ/γ_r, which compares the self-aligning length to the particle size; and, finally, the reduced interaction strength, . Details of the numerical simulations are discussed in Appendix A.

III. Results

In order to investigate the interplay between chirality and self-alignment in an active crystal, we perform simulations by fixing the Péclet number at Pe = 10 and systematically varying the strength of the reduced self-alignment B and the reduced chirality ω (see Fig. 2(a)). In ref. 80, a non-chiral active solid characterized by self-alignment was already studied and revealed a transition from a disordered to a flocking state as the alignment strength increases. Here, we aim to investigate how the inclusion of a chiral term modifies the dynamics of the crystal, ultimately leading to the complete suppression of the flocking state observed for sufficiently strong chirality.

Fig. 2 Phase diagram. (a) Phase diagram in the plane of reduced self-alignment strength B and reduced chirality ω. Colors denote different phases which are visualized by the four snapshots (b)–(e), where vectors, colored using the same color as the corresponding phase in the phase diagram, represent the normalized particle velocities. (b) Circling crystal (CC) phase (light brown), where all particle velocities are aligned, and all the system performs synchronized circular orbits in space. (c) Vortex phase (VP) (brown), where the system develops vortex structures and the global circling is suppressed. (d) Disordered system (DS) (black), where particle velocities form finite-size domains where particles are aligned. (e) Chiral disordered system (CDS) (grey), where the particles’ velocity is aligned in finite-size domains, but the particles show oscillatory dynamics. The orange region that separates CC and VP phases represents metastable configurations, showing circling crystal or vortex structures depending on the initial conditions. Colored rectangles in the phase diagram outline the values of parameters analyzed in successive figures. The remaining dimensionless parameters of the simulations are Pe = 10, M = 10⁻⁴, , and ϕ = Nπσ²/4L² = 1.1.

A. Single-particle dynamics in a self-aligning chiral active solid

We start by monitoring the dynamics of a target particle inside the crystal (see the Supplementary Video). The particle trajectories, which in Fig. 1(b)–(e) are represented by colored dots that become lighter with time and are connected by gray lines, strongly depend on the relative strengths of reduced chirality ω and reduced self-alignment B. For weak B (i.e., before the flocking transition⁸⁰) and small ω, the particle irregularly moves with a displacement that exceeds the particle diameter σ (Fig. 1(d)). This observation suggests that the system does not flock but the particles move in finite-size domains as in ref. 84, implying the absence of a global synchronization among the particles of the crystal. In this regime, the increase in chirality induces more localized trajectories (Fig. 1(e)), with small position fluctuations even smaller than the particle diameter σ. At low chirality and strong self-alignment (Fig. 1(b)), particles are able to globally coordinate their motion, giving rise to well-defined circular orbits with a radius that scales as the ratio v₀/Ω, much larger than the particle diameter σ. These extended orbits emerge from the collective motion of the system, in which all particles perform a circular coordinate motion, suggesting the occurrence of a new flocking-like state with rotating features. Finally, upon increasing chirality in this regime of strong self-alignment, the radius of the circular orbits decreases (Fig. 1(c)), becoming even smaller than the particle diameter σ. The reduced size of the orbits effectively suppresses the fluctuations of the particles' displacement within their lattice cages. However, local alignment seems to be preserved, as reflected by the nearly perfect circular trajectories followed by each particle.

This single-particle trajectory study implies that introducing a sufficiently strong self-alignment enables the particles to coordinate their motion, either globally or locally, depending on the intensity of the chiral term, and effectively reduces noise effects. It is reflected in the almost perfectly circular trajectories shown in Fig. 1(b) and (c).

B. Phase diagram

The analysis of the single-particle trajectories suggests the occurrence of collective behavior. Therefore, we systematically build a phase diagram in the plane of reduced chirality ω and reduced self-alignment B, which allows us to identify four distinct phases (Fig. 2(a)). For weak values of B and ω, the system remains in a disordered system (DS) phase, where particles are unable to globally coordinate their motion and chirality is too weak to significantly affect the system dynamics: the system shows finite-size domains where particles locally align their velocity (Fig. 2(d)). This local velocity alignment is a consequence of the high-density and active force persistence and results in short-range velocity correlations, consistently with previous findings without self-alignment and chirality.⁸⁴ By keeping self-alignment weak and increasing chirality, the system enters the chiral disordered system (CDS) phase. In this regime, particles still exhibit local alignment in the velocity field due to the interplay between active forces and the high packing fraction of the system (Fig. 2(e)). The main difference between the disordered phase (DS) and the CDS is that increasing the strength of the chiral term modifies the system's dynamics by reducing the size of velocity-aligned domains and introducing an oscillatory behavior in the temporal dynamics. Further details will be discussed in the following sections.

On the other hand, by increasing the value of the self-alignment strength B, the system enters a flocking-like phase typical of crystals composed of self-aligning units (see ref. 80), where the particles' velocities are globally aligned. In this regime, while chirality does not alter the transition line from the disordered to the ordered phase, it alters the flock dynamics. Indeed, here, active units do not follow an almost straight trajectory as in the standard case of flocking behavior but display collective circular trajectories, as suggested by tracking a single-particle trajectory inside the crystal (Fig. 1(b)). In this regime, all the particles' velocities are aligned (see Fig. 2(b)) as an effect of self-alignment, but the whole group of particles collectively follow circular orbits due to chirality. Consequently, we term this collective behavior, the circling crystal (CC) phase. The characteristic radius of these circular orbits is determined by the interplay between the active force and chirality and scales as v₀/Ω. For a fixed self-propulsion of the active units, increasing the value of Ω leads to a decrease in the typical radius of the orbits, as illustrated by the single-particle trajectories inside the active solid shown in Fig. 2(b) and (c).

To theoretically predict the circling crystal phase, we use a method similar to the one used in ref. 80 by including the additional effect of chirality. As derived in Appendix B, the particle dynamics (eqn (1)) can be mapped onto a coarse-grained dynamics in the continuum limit, with the following form

where η(r,t) is a vector representing Gaussian white noise satisfying 〈η(r,t)η(r′,t′)〉 = δ(t − t′)δ(r − r′) and z corresponds to the unit vector normal to the plane of motion. The term denotes the Landau–Ginzburg free energy, which reads

The first term on the right-hand side of eqn (3) represents the kinetic contribution, which tends to restore random configurations. The remaining two terms correspond to a single-particle free energy with a Mexican-hat shape, whose form depends on the sign of the mass term, i.e. the term proportional to |v(r,t)|². Compared to the results of ref. 80, chirality generates an additional force Ωz × v(r,t) which cannot be expressed as the functional derivative of an effective free energy. The structure of the chirality-induced term is reminiscent of a Lorentz force due to an external magnetic field and therefore it is intuitively responsible for a circular current, applying at the single-particle level.

To shed light on the effect of the chiral term on the system dynamics, we analyze the complex velocity ψ(r,t) = v_x(r,t) + iv_y(r,t), which is then convenient to express in polar coordinates ψ(r,t) = v(r,t)e^iφ(r,t), where v = v(r,t) is the velocity modulus and φ = φ(r,t) the polar angle. With this choice the dynamics (2) can be rewritten as

The time evolution of the velocity modulus, v, is not affected by the chirality and remains consistent with the Mexican-hat free-energy profile described by eqn (3). However, we find that the chiral term uniquely influences the time evolution of the velocity orientation φ through a constant term. Therefore by considering a reference frame that rotates with angular velocity Ω, i.e. implementing the change of variables φ′ = φ − Ωt in eqn (4b), chirality disappears from the dynamics. This is equivalent to looking for a solution with a time-shifted phase, which can hold only for weak chirality. In this regime, assuming that spatial gradients of the velocity modulus are negligible, the constantly rotated angle φ′ is described by a diffusive dynamics

As a consequence, in the rotating reference frame, the angle field gives rise to the Goldstone modes associated with the spontaneous breaking of rotational symmetry. This theory explains the circling crystal phase because the solution for the polar angle φ evolves diffusely but is subject to a constant temporal increase with angular velocity Ω. Since the only effect due to chirality is imposing a constant rotation to the velocity field, our theory predicts that the transition line from the disordered to the circling crystal is not affected by chirality and occurs when the mass term changes sign, i.e., when β exceeds the critical value β_c, given by

This theoretical prediction is consistent with our numerical findings (Fig. 2(a)) and with the flocking transition observed in the absence of chirality.⁸⁰

By maintaining a high value of self-alignment B and selecting strong chirality values ω, the system enters the vortex phase (VP). In this regime, we observe that chirality competes with self-alignment, thereby suppressing the flocking behavior. For ω ≳ 10, the orbital radius r_Ω is comparable to the one of the particle, r_Ω ≲ σ/2. This indicates that the circular motion of the particles is almost localized within the individual lattice cage (see Fig. 1(c)), which corresponds to low effective temperature of the system. As a consequence, global alignment – which generally emerges when self-alignment allows random persistent fluctuations of the particles to break the global rotational symmetry – is suppressed. Instead, local ordered domains emerge, where within each of the particle velocities are aligned and exhibit coordinated circular motion (see Fig. 2(c)). To ensure that the VP is not an artificial finite-size effect, we have systematically performed simulations with different particle numbers, verifying that the vortex spacing observed in the velocity field is unaffected by the size of the simulation box. We remark that, between the circling crystal (CC) phase and the vortex phase (VP), there exists an intermediate region, highlighted in orange in the phase diagram, where the system may show a stable global rotating phase or vortex-like structures, depending on the randomly chosen initial conditions.

In Fig. 2(a), three different colors (blue, red, and green) are used to highlight specific lines of the phase diagram corresponding to different values of the reduced self-alignment (B = 0.05, 0.2, and 0.5). Along these lines, we observe transitions between all the four phases discussed in this work. The same color scheme will be consistently employed in the subsequent figures.

C. Spatial velocity correlation functions and energy spectra

In order to characterize the phases previously described, we analyzed both the spatial velocity correlation function and the kinetic energy spectrum (Fig. 3). The spatial velocity correlation function is defined aswith r_ij = |r_i − r_j|. Here, the symbol 〈⋯〉 corresponds to the time average. In the case of a crystal composed by passive units, C(r) does not have a spatial dependence and is approximately given by a Dirac delta function centered at r = 0. In the absence of active forces, the particles' velocities are not correlated since the system is governed by the Maxwell–Boltzmann distribution. In an active solid without both chirality and self-alignment, C(r) typically exhibits a spatial decay with an exponential profile and a correlation length increasing with the persistence time τ of active force.⁸⁴ When chirality is included in the dynamics of individual particles, the characteristic length of the correlated domains decreases as chirality increases consistently with theoretical and numerical results of ref. 59 and 60 in the absence of self-alignment. On the other hand, for strong enough self-alignment the system undergoes a flocking transition. In this case, the velocity correlation function has an almost constant profile.⁸⁰

Fig. 3 Correlation functions and energy spectra. (a), (c), and (e) Spatial velocity correlations as a function of the distance r for reduced self-alignment values of B = 0.05, 0.2, and 0.5, respectively. (b), (d), and (f) Kinetic energy spectra of the system, E(k), defined in eqn (8), as a function of k for B = 0.05, 0.2, and 0.5. The simulations at B = 0.5 are performed in a larger box (180 × 180) in order to avoid finite-size effects. The other dimensionless parameters of the simulations are: Pe = 10, M = 10⁻⁴, , and ϕ = Nπσ²/4L² = 1.1.

To focus on the spread of information across scales, we study the kinetic energy spectrum, defined as

here, k = |k| corresponds to the modulus of the wave vector and coincides with the Fourier transform of the velocity field. The maximum and minimum accessible values of k are defined, respectively, by the inverse of the particle diameter and by the inverse of the box size, both scaled by 2π. In general, when the energy spectra exhibit a power-law decay in the range of wavenumbers k between the value corresponding to the spectral peak and the smallest k defined by the inter-particle distance, this indicates the emergence of a scale-free regime in real space. The presence of a peak at a well-defined wavenumber k signals the existence of structures with a characteristic size in the velocity field. When the peak is not visible, it occurs at smaller values of k that are not accessible due to the finite system size, which imposes a lower limit on the range of k values studied. The increase of the energy spectrum always observed at large k values originates from the intrinsically disordered packing at short length scales. For a non-chiral active crystal without self-alignment, the energy spectrum profile scales as E(k) ∼ k⁻¹ for large k, consistently with an Ornstein–Zernike profile.^85,86 In ref. 59, where chirality is introduced at the level of single-particle dynamics, the energy spectra reveal the emergence of correlated velocity fields for low noise and weak chirality.

In our study, we combine both chirality and self-alignment in a crystal composed by active units. For B = 0.05 (Fig. 3(a) and (b), blue panel), the spatial velocity correlation functions for ω = 0.5, 50, and 100 (Fig. 3(a)) exhibit an exponential decay, with a characteristic length that decreases as chirality increases. This reflects the fact that, as chirality increases, particles tend to move along orbits of smaller radius, resulting in a system characterized by short range correlations in the velocity field. The corresponding energy spectra (Fig. 3(b)) show a slower decay at larger k for ω = 0.5 (DS phase), characterized by a power-law behavior of the form E(k) ∼ k⁻¹, in accordance with an Ornstein–Zernike profile, as discussed in recent works on chiral and non-chiral active crystals.^59,60,86 As chirality increases (CDS phase), the spectral peak shifts to larger k, reflecting the emergence of smaller structures in the velocity field (Fig. 2(e)), while the power-law decay becomes steeper and, for extreme chirality strength, scales as E(k) ∼ k⁻². Increasing chirality leads to smaller particle domains where velocities are correlated, but also to a progressive loss of local fluctuations and instabilities, which is reflected in a faster decay of the energy spectrum.

At B = 0.2 (Fig. 3(c) and (d), red panel) and ω = 0.5, the system enters the circling crystal (CC) phase: all particles are aligned and move collectively along circular orbits (see Fig. 2(b)), resulting in a nearly constant spatial velocity correlation function close to one as shown in ref. 80. The corresponding energy spectrum displays a delta-like peak at k = 0, which cannot be properly extrapolated from the data since the minimum accessible k is determined by the inverse of the simulation box size, scaled by 2π. Increasing ω to 50 or 100 drives the system back into the chiral disordered (CDS) phase, where correlations decay exponentially, and the energy spectrum follows a power-law behavior of the form E(k) ∼ k⁻².

For B = 0.5 (Fig. 3(e) and (f), green panel), the CC phase persists at ω = 0.5 with similar correlation and spectral features. At higher chirality (ω = 50, 100), the system transitions to the vortex phase (VP), characterized by a correlation function that crosses zero and becomes negative at scales set by the domain size (Fig. 2(c)), and energy spectra with pronounced peaks at finite k. Both the zero-crossing in the spatial velocity correlation function and the peak positions in the spectrum indicate that the characteristic size of coherent velocity structures decreases with increasing chirality. In this case, the energy spectrum at large k scales as a power law, E(k) ∼ k⁻⁴, indicating an almost complete suppression of local instabilities in favor of well-defined vortex structures in the velocity field.

The fast power-law decay observed in the energy spectra for all regimes characterized by high chirality reflects the fact that this term introduces a characteristic correlation length in the velocity field, defining the spatial scale on which most of the energy is concentrated.

To summarize, in the DS and CDS phases the velocity correlation functions decay rapidly, especially as the strength of the chiral term increases. In the CC phase, the alignment of the velocities is reflected in the almost constant behavior of the function. The same function, which oscillates in the VP phase, reflects finite-size domains in which the particle velocities are aligned. Finally, all four phases are also characterized by different exponents in the decay of the energy spectra, reflecting on which spatial scale the energy of the system is concentrated.

D. Velocity autocorrelation function

To characterize the dynamics of the phases described in the phase diagram, we have also analyzed the velocity autocorrelation function, defined as:where 〈⋯〉 corresponds to the time and particle average. Starting from the lowest self-alignment value (Fig. 4(a)–(c), blue panel), C_t for ω = 0.5 (corresponding to the DS phase in Fig. 4(a)) exhibits a rapid exponential decay. This behavior indicates that the particle dynamics in this regime are dominated by the rotational noise, which tends to randomize the particle orientation over the characteristic time scale τ, in agreement with previous studies on non-chiral active solids. In contrast, for ω = 50, 100 (Fig. 4(b) and (c)) we are in the CDS phase and the velocity autocorrelation functions display clear oscillatory profiles, with a periodicity set by the chiral term Ω, and an evident damping behavior as also shown in ref. 59. These chirality-induced oscillations are consistent with the oscillating mean-square displacement encountered even in a free chiral active particle^87–89 and can be related to the occurrence of the odd diffusive behavior typical of chiral systems.^32,33,90

Fig. 4 Velocity auto-correlation functions. (a)–(c) Velocity auto-correlation functions C_t(t) for B = 0.05 and ω = 0.5, 50, 100. (d)–(f) C_t(t) for B = 0.2 and ω = 0.5, 50, 100. (g)–(i) C_t(t) for B = 0.5 and ω = 0.5, 50, 100. The simulations at B = 0.5 are performed in a larger box (180 × 180) in order to avoid finite-size effects. The other dimensionless parameters of the simulations are: Pe = 10, M = 10⁻⁴, , and ϕ = Nπσ²/4L² = 1.1.

For higher self-alignment values, B = 0.2 (Fig. 4(d)–(f), red panel) and B = 0.5 (Fig. 4(g)–(i), green panel), all velocity autocorrelation profiles, even the ones at small chirality, exhibit oscillatory behavior consistent with the chirality-induced periodicity. Furthermore, increasing the self-alignment strength reduces the damping of these oscillations, until, for the strongest self-alignment, the curves no longer display any decay. Indeed, strong self-alignment allows the particles to coordinate their motion in an extremely persistent manner, preventing the loss of coordination that is normally induced by rotational noise. The introduction of the alignment term effectively forces particle trajectories to no longer be dominated by thermal noise, but instead to coordinate their motion over periods of time much longer then the characteristic time scales of the system, globally in the circling crystal phase or locally in the vortex phase, allowing the emergence of undamped periodic motion whose period is determined by the chirality of individual particles.

In the end, the chiral term is able to affect the system dynamics in all phases except the DS phase, in which both chirality and self-alignment are still dominated by the thermal noise. Moreover, if the self-alignment term is sufficiently strong, the chiral particle motion becomes persistent in time, as observed in the CC and VP phases.

IV. Conclusions

In this paper, we have studied the interplay between self-alignment and chirality in active crystals. The first mechanism drives the particle orientation to align with the velocity, while the second one drives the particles along circular orbits. Consistent with the previous results, for weak chirality and strong self-alignment, a circling crystal phase was found where the individual particles synchronize and all move along the same circular orbit. When chirality is strong, the circular collective motion is suppressed and the system develops vortex-like structures in the velocity field. In this state, the system displays an oscillatory behavior in the spatial velocity correlations and time-oscillating velocity autocorrelations. It is important to note that the four phases presented in this work differ in the properties of the velocity field and in the presence or absence of a global motion of all the particles forming the active crystal. Owing to the high packing fraction, each particle is confined within a lattice cage with no possibility of escape, thereby preserving the hexagonal crystal order independently of the values of self-alignment or chirality. Further studies could aim to develop a field theory capable of explaining and predicting the suppression of this transition numerically observed, for instance by generalizing the theory developed in ref. 91 to self-aligning chiral systems.

The results presented in this paper suggest that a natural extension of this work would be to study the effects of self-alignment and chirality in systems of active Brownian particles at lower densities. Recent studies on active turbulence^92,93 indicate that the low-density counterpart of the system presented here may exhibit active turbulent-like behavior for strong chirality. Active turbulence^94,95 differs from the usual inertial turbulence observed in fluids at high Reynolds numbers. In the latter case, the interplay between inertial and viscous forces destabilizes the laminar flow, leading to the formation of vortices. The situation is different in active turbulence: here, no external driving is needed, as the internal propulsion of the active units is sufficient to induce complex dynamics, giving rise to self-driven flows that can even become chaotic.

Our findings are amenable to experimental validation in both biological systems and in active granular media.^{18,31,64,65,73,96–100} In the biological context, it is well established that many cellular assemblies exhibit self-alignment mechanisms,⁶⁴ while others display intrinsic chirality.¹³ Exploring systems that combine these two features could provide new insights into their collective behaviors. Regarding granular systems, several experimental realizations already exhibit self-alignment. Examples include Hexbugs,^19,23,73,101 commercial toys characterized by a certain degree of self-alignment, or 3D-printed granular walkers,⁷⁵ where the strength of self-alignment can be tuned by adjusting the mass distribution and the position of the center of friction. Including chirality in these systems can be easily achieved by simply breaking the rotational symmetry of the particle body, as demonstrated in light-driven robots propelled by internal vibrations,²⁰ intelligent robots equipped with internal motors generating horizontal oscillations²⁶ or air-fluidized disks consisting of a set of blades with a constant tilt.^30,31 The possibility of combining these two ingredients, self-alignment and chirality, in dense systems could enable the experimental observation of the phases described in this work and provide new insights into the interplay between collective motion, chirality, and structural ordering in active matter. The experimental realization of these phases may pose technical challenges due to boundary effects, particularly for the circling crystal phase, in which the entire system exhibits circular orbits in space. In the presence of hard or soft walls, particles near the boundaries lack sufficient space to perform such orbits; as a consequence, bulk vortices and edge currents may emerge. Overall, our results suggest that a careful control of microscopic activity parameters offers a promising route to stabilize novel collective states in active crystalline systems.

Author contributions

M. M. performed the numerical study; L. C. derived the theory; and M. M and L. C. wrote the first draft of the paper. M. M., A. A., H. L., and L. C. equally contributed to the data interpretation, presentation, and writing of the manuscript.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5sm01135c.

Data for this article, including simulation data, are available at ‘Zenodo’ at https://doi.org/10.5281/zenodo.17592187.

Appendices

Appendix A: details of the numerical simulations

The dynamics described in the main text using eqn (1a) and (1b) are numerically implemented by using the Euler integration scheme with a time step of Δt = 5 × 10⁻⁶τ. Rescaling time by the persistence time of the particle's trajectory, τ, and positions by the particle's diameter, σ, the integration scheme is as follows:

In the equations above, the prime symbol denotes dimensionless variables used in our simulations. In this formulation, the particle orientation is given by n̂_i = (cos θ_i,sin θ_i), and Y_i are Gaussian random variables with zero mean and unit variance. The force arising from interactions can be written as F_i = −∇_iU_tot, where , with U chosen as a WCA potential, as described in the main text.

The particle dynamics are controlled by the dimensionless parameters discussed in the main text that we summarize again here: the Péclet number, Pe = v₀/(D_rσ), that is fixed to Pe = 10 in our numerical study; the reduced mass, M = D_rm/γ, which compares the inertial time m with the persistence time τ = 1/D_r. This parameter corresponds to the inverse of the coefficient multiplying the second and last terms on the right-hand side of eqn (A1b). This is fixed to M = 10⁻⁴ to explore the low-inertia regime; and the dimensionless interaction strength, , which can be extracted from the third term of the same equation.

Chirality and self-alignment introduce additional dimensionless parameters: the reduced chirality, ω = Ωτ, varied between 10⁻¹ ≤ ω ≤ 10², and the reduced self-alignment strength, B = βσ/γ_r, varied in the range 5 × 10⁻² ≤ B ≤ 10⁰. Simulations are performed for a total time of 2 × 10²τ in a square box of size L = 125 with periodic boundary conditions. Simulations are performed with a number of particles N = 2.2 × 10⁴, yielding a packing fraction of Φ = Nπσ²/4L² = 1.1.

Each point in the phase diagram corresponds to a single simulation, which is sufficiently long (t = 2 × 10²τ) to ensure that the system has reached the steady state. For the points close to the phase boundaries, we have performed multiple simulations in order to determine the boundaries more clearly.

Appendix B: derivation of the theoretical prediction of the circling crystal

a. Algebraic derivation of the chiral self-aligning active dynamics. The dynamics for chiral active Brownian particles interacting with self-alignment in two dimensions can be expressed in Cartesian components as followswhere v_i = ẋ_i. Here, we have neglected the translational noise to be consistent to the model used for the simulations.

As a first step, we rewrite the self-alignment interactions in the following form (as already done in ref. 80)

β(n̂_i × v_i) × n̂_i = β(v_i − n̂_i[v_i·n̂_i]), (B2)

where we have used that n_i² = 1. At first, we define the acceleration s_i = v̇_i and then apply the time derivative to eqn (B1a). By using eqn (B1b) to calculate , we obtain

In order to express the dynamics as a function of position, velocity, and acceleration only, we eliminate n_i by using eqn (B1a):

The dynamics (B4) are equivalent to the equations of motion (B1a) and (B1b) numerically implemented and result from a change of variables from (x_i,v_i,n̂_i) to (x_i,v_i,s_i). We remark that chirality solely appears in the last line of eqn (B4), while the remaining terms are simply due to the interplay between self-alignment activity and potential forces.

b. Vanishing the inertia limit. In the limit of small inertia, the leading contributions to the dynamics (B4) are the terms proportional to γ/m, and neglecting higher-order contributions we obtain

In the overdamped limit m/γ → 0, one can eliminate the variable s_i, neglecting the time evolution of the acceleration, i.e. ṡ_i = 0. In this way, we obtain an equation of motion for v_i, which reads

The first term in eqn (B6) is a local alignment term and by calling N_i the instantaneous number of particles interacting with the particle i, we have

In this interaction, the particle's velocity v_i is attracted to the velocity of particle j depending on the second derivative of the potential. The other terms in eqn (B6) are a force F_i normalized by the persistence time τ, the chiral term and a noise one. In addition, self-alignment induces other single-particle non-linear terms that depend on v_i and F_i.

c. Lattice approximation. In the crystalline configuration, particles occupy the sites of an ordered triangular lattice and exhibit small fluctuations around their equilibrium positions. To proceed with the derivation, we adopt the lattice approximation, effectively freezing the particle positions at their lattice sites. Under this assumption, the net force on each particle vanishes, F_i = 0. In addition, the second derivative of the potential becomes constant and can be expressed as

Here, N_i = 6 for a two-dimensional triangular lattice, and the constant matrix , which depends on the distance between particles i and j, is obtained by computing the second spatial derivatives of the potential U:

where r^α_ij = r^α_i − r^α_j, with α = x,y. In this expression, each prime on the potential U means a spatial derivative. In addition, the following relation holds:

∇^α_i∇^β_jU = −∇^α_i∇^β_iU. (B10)

The Cartesian components of r^α_ij/|r_ij| can be written as cos(δ_j) and sin(δ_j), where δ_j is the angle between the vector r_ij and the horizontal axis. Owing to the regular hexagonal structure of the cluster, this angle can be expressed as δ_j = δ₀ + jπ/3, with j = 0, 1, …, 5, where δ₀ denotes the orientation of the hexagon with respect to the reference frame, which we set to zero without loss of generality. In this way, we finally obtain the elements of the matrix

where we have suppressed the explicit dependence on j. Finally, we remark that the sum over the six neighbors of the out-of-diagonal elements of vanisheswhile the same protocol for the diagonal elements gives a constant with the following expression

In this way, we have

where is the identity matrix. At the end, the particle dynamics in the lattice approximation can be expressed aswhere the constant is defined by eqn (B15) and we have neglected the directional contribution of the lattice. In this way, the first term on the right-hand side of eqn (B17) corresponds to a discrete Laplacian.

d. Continuum limit. Here, we consider the continuum limit by replacing v_i → v(r,t), where r is a continuous coordinate on the lattice. In this way, the discrete Laplacian in eqn (B17) is replaced bywhere the particle diameter σ also corresponds to the lattice constant. The dynamics (B4) in the continuum limit becomes a dynamics for the velocity field v(r):

The same dynamics can be written also in terms of the functional derivative of a free-energy functional

where the term η(r,t) is a vector representing white noise such that 〈η(r,t)η(r′,t′)〉 = δ(t − t′)δ(r − r′), and Ωz × v(r,t) is the term that introduces the chirality in the dynamics and is the Landau–Ginzburg free energy, which reads:

The first term on the right-hand side of eqn (B21) represents the kinetic contribution, which tends to restore random configurations. The two other terms correspond instead to a single-particle free-energy with a Mexican-hat shape, whose form depends on the sign of the mass term.

e. Circling crystal prediction. At this point, in order to investigate the effect of the chirality, we can express the velocity field as a complex number, ψ(r,t) = v_x(r,t) + iv_y(r,t), and decompose the dynamics into its cartesian components as follows:

Recombining the two equations, we obtain

Recalling the complex expression for ψ(r,t) = v_x(r,t) + iv_y(r,t) and neglecting the noise term, which is not essential for the following results, we can write the following equation:

where we have used |v(r,t)|² = |ψ(r,t)|². Expressing the complex velocity in its exponential form, ψ(r,t) = v(r,t)e^iφ(r,t), in terms of the velocity modulus v(r,t) and the polar angle φ(r,t), and substituting this into eqn (B25), we can obtain evolution equations for v(r,t) and φ(r,t) by calculating the real and imaginary parts in eqn (B25), obtainingwhere we have used:

[small psi, Greek, dot above](r,t) = v̇(r,t)e^iφ(r,t) + i[small phi, Greek, dot above](r,t)v(r,t). (B28)

Finally, since v > 0 and that it does not vary in space, we obtain:

Integrating this last expression over time, we obtain

From this result, it is clear that an additional contribution appears in the time evolution of the velocity orientation: the chiral term Ωt, which induces a constant rotation of the particle velocity. At this point, it is convenient to define a new angle φ′ = φ − Ωt, corresponding to a reference frame rotating with angular velocity Ω. In this new referent frame the dynamics of the system remain the same described in ref. 80 and the transition from a disordered to an ordered state occurs when the mass term in eqn (B21) changes sign, i.e., when β exceeds the critical value β_c, given by

Appendix C: supplementary video details

The supplemental video provided in the SI illustrates the temporal evolution of a chiral active crystal subject to self-alignment, as described using eqn (1). In the video, the velocity of each particle is represented by a black arrow, while the trajectory and behavior of a selected target particle are highlighted in red. Four different time evolutions of the system are shown in parallel, corresponding to distinct values of the self-alignment strength and chirality. The stars displayed in each panel indicate the phases under consideration, and their colors follow the same convention as in Fig. 2. Specifically, the top and bottom panels correspond to reduced self-alignment parameters of B = 0.5 and B = 0.1, respectively, while the left and right panels correspond to chirality values of ω = 0.5 and ω = 40.

For B = 0.1 and ω = 0.5 (black star) or ω = 40 (grey star), the system exhibits the disordered system (DS) or chiral disordered system (CDS) phases, respectively, where particle velocities display only local alignment. In the CDS phase, the target particle undergoes a periodic rotation induced by the chiral term, as revealed by the trajectory of the tagged particle. For B = 0.5 and ω = 0.5 (beige star), the system enters the circling crystal (CC) phase, characterized by global alignment of particle velocities and coherent circular motion of the entire crystal, as evidenced by the trajectory of the target particle. Finally, for B = 0.5 and ω = 40 (brown star), the system reaches the vortex phase (VP), where particle velocities remain aligned and coordinated within finite-size domains, and the velocity field exhibits vortex-like structures. The tagged particle still performs a periodic rotation, but due to strong chirality its motion becomes almost localized within the crystal cage.

Acknowledgements

HL acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the SPP 2265, under grant number LO 418/25. LC acknowledges financial support from the University of Rome La Sapienza, under the project Ateneo 2024 (RM124190C54BE48D).

References

1. M. Marchetti, J. Joanny, S. Ramaswamy, T. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189.
2. J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
3. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
4. H. Löwen, Eur. Phys. J.: Spec. Top., 2016, 225, 2319–2331.
5. B. Liebchen and D. Levis, Europhys. Lett., 2022, 139, 67001.
6. W. R. DiLuzio, L. Turner, M. Mayer, P. Garstecki, D. B. Weibel, H. C. Berg and G. M. Whitesides, Nature, 2005, 435, 1271–1274.
7. H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 268102.
8. A. P. Petroff, X.-L. Wu and A. Libchaber, Phys. Rev. Lett., 2015, 114, 158102.
9. E. Pérez Ipiña, S. Otte, R. Pontier-Bres, D. Czerucka and F. Peruani, Nat. Phys., 2019, 15, 610–615.
10. P. Patra, K. Beyer, A. Jaiswal, A. Battista, K. Rohr, F. Frischknecht and U. S. Schwarz, Nat. Phys., 2022, 18, 586–594.
11. F. Afroze, D. Inoue, T. I. Farhana, T. Hiraiwa, R. Akiyama, A. M. R. Kabir, K. Sada and A. Kakugo, Biochem. Biophys. Res. Commun., 2021, 563, 73–78.
12. I. H. Riedel, K. Kruse and J. Howard, Science, 2005, 309, 300–303.
13. B. Wang and G. Xu, Biophys. J., 2022, 121, 1931–1939.
14. F. Kümmel, B. ten Hagen, R. Wittkowski, I. Buttinoni, R. Eichhorn, G. Volpe, H. Löwen and C. Bechinger, Phys. Rev. Lett., 2013, 110, 198302.
15. Y. Shelke, N. R. Srinivasan, S. P. Thampi and E. Mani, Langmuir, 2019, 35, 4718–4725.
16. D. Seč, T. Porenta, M. Ravnik and S. Žumer, Soft Matter, 2012, 8, 11982–11988.
17. L. N. Carenza, G. Gonnella, D. Marenduzzo and G. Negro, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 22065–22070.
18. C. Scholz, M. Engel and T. Pöschel, Nat. Commun., 2018, 9, 931.
19. O. Chor, A. Sohachi, R. Goerlich, E. Rosen, S. Rahav and Y. Roichman, Phys. Rev. Res., 2023, 5, 043193.
20. F. Siebers, A. Jayaram, P. Blümler and T. Speck, Sci. Adv., 2023, 9, eadf5443.
21. C. W. Chan, D. Wu, K. Qiao, K. L. Fong, Z. Yang, Y. Han and R. Zhang, Nat. Commun., 2024, 15, 1406.
22. L. Caprini, I. Abdoli, U. M. B. Marconi and H. Löwen, Newton, 2025, 1, 100253.
23. D. Boriskovsky, R. Goerlich, B. Lindner and Y. Roichman, Soft Matter, 2026, 22, 297–305.
24. R. Kant, A. Maitra, A. K. Sood and S. Ramaswamy, arXiv, 2025, preprint, arXiv:2509.00729 DOI:10.48550/arXiv.2509.00729.
25. J. P. Carrillo-Mora, A. Garcés and D. Levis, Phys. Rev. E, 2025, 112, 065417.
26. M. Noirhomme, I. Mammadli, N. Vanesse, J. Pande, A.-Š. Smith and N. Vandewalle, Phys. Rev. Appl., 2025, 23, 064008.
27. F. Novkoski, M. Mélard, M. Delens, F. Wéry, M. Noirhomme, J. Pande, A.-Š. Maier, A.-Š. Smith and N. Vandewalle, arXiv, 2025, preprint, arXiv:2509.07437 DOI:10.48550/arXiv.2509.07437.
28. M. Workamp, G. Ramirez, K. E. Daniels and J. A. Dijksman, Soft Matter, 2018, 14, 5572–5580.
29. C. Scholz, A. Ldov, T. Pöschel, M. Engel and H. Löwen, Sci. Adv., 2021, 7, eabf8998.
30. M. Ángel López-Castaño, Á. Rodríguez-Rivas and F. V. Reyes, Front. Phys., 2022, 10, 972051.
31. M. A. López-Castaño, A. M. Seco, A. M. Seco, Á. Rodríguez-Rivas and F. V. Reyes, Phys. Rev. Res., 2022, 4, 033230.
32. C. Hargus, J. M. Epstein and K. K. Mandadapu, Phys. Rev. Lett., 2021, 127, 178001.
33. E. Kalz, H. D. Vuijk, I. Abdoli, J.-U. Sommer, H. Löwen and A. Sharma, Phys. Rev. Lett., 2022, 129, 090601.
34. E. Kalz, H. D. Vuijk, J.-U. Sommer, R. Metzler and A. Sharma, Phys. Rev. Lett., 2024, 132, 057102.
35. P. L. Muzzeddu, E. Kalz, A. Gambassi, A. Sharma and R. Metzler, New J. Phys., 2025, 27, 033025.
36. L. Caprini, H. Löwen and U. M. B. Marconi, Soft Matter, 2023, 19, 6234–6246.
37. Y. Wang and J. Berx, Soft Matter, 2025, 21, 6697–6706.
38. G.-J. Liao and S. H. L. Klapp, Soft Matter, 2021, 17, 6833–6847.
39. K. L. Kreienkamp and S. H. L. Klapp, New J. Phys., 2022, 24, 123009.
40. T. Hiraiwa, R. Akiyama, D. Inoue, A. M. R. Kabir and A. Kakugo, Phys. Chem. Chem. Phys., 2022, 24, 28782–28787.
41. S. Ceron, K. O'Keeffe and K. Petersen, Nat. Commun., 2023, 14, 940.
42. T. Lei, C. Zhao, R. Yan and N. Zhao, Soft Matter, 2023, 19, 1312–1329.
43. J. Mecke, J. O. Nketsiah, R. Li and Y. Gao, Natl. Sci. Open, 2024, 3, 20230086.
44. K. S. H. L. Liao Guo-Jun, Soft Matter, 2018, 14, 7873–7882.
45. J. J. Bickmann Jens, B. Stephan and W. Raphael, J. Chem. Phys., 2022, 156, 194904.
46. L. Langford and A. K. Omar, Phys. Rev. Lett., 2025, 134, 068301.
47. C. M. P. Lei Qun-Li and N. Ran, Sci. Adv., 2019, 5, eaau7423.
48. B. Zhang and A. Snezhko, Phys. Rev. Lett., 2022, 128, 218002.
49. Y. Kuroda and K. Miyazaki, J. Stat. Mech.:Theory Exp., 2023, 2023, 103203.
50. S. A. Zhang Bo and S. Alexey, Nat. Commun., 2020, 11, 4401.
51. L. Caprini, B. Liebchen and H. Löwen, Commun. Phys., 2024, 7, 153.
52. N. Kruk, J. A. Carrillo and H. Koeppl, Phys. Rev. E, 2020, 102, 022604.
53. Z. Bing-Quan and S. Zhi-Gang, Phys. A, 2022, 598, 127373.
54. A. Negi, K. Beppu and Y. T. Maeda, Phys. Rev. Res., 2023, 5, 023196.
55. Y. Wang, B. Ventéjou, H. Chaté and X.-Q. Shi, Phys. Rev. Lett., 2024, 133, 258302.
56. D. Jian, J. Su, J. Wang and J. Wang, Commun. Phys., 2025, 8, 386.
57. B. Liebchen and D. Levis, Phys. Rev. Lett., 2017, 119, 058002.
58. D. Levis and B. Liebchen, J. Phys.: Condens. Matter, 2018, 30, 084001.
59. A. Shee, S. Henkes and C. Huepe, Soft Matter, 2024, 20, 7865–7879.
60. U. M. B. Marconi and L. Caprini, Soft Matter, 2025, 21, 2586–2606.
61. Z.-F. Huang, M. te Vrugt, R. Wittkowski and H. Löwen, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2511350122.
62. Y. Kuroda, T. Kawasaki and K. Miyazaki, Phys. Rev. Res., 2025, 7, L012048.
63. B. Szabo, G. Szöllösi, B. Gönci, Z. Jurányi, D. Selmeczi and T. Vicsek, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 061908.
64. S. Henkes, Y. Fily and M. C. Marchetti, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2011, 84, 040301.
65. D. L. Barton, S. Henkes, C. J. Weijer and R. Sknepnek, PLoS Comput. Biol., 2017, 13, e1005569.
66. F. Giavazzi, M. Paoluzzi, M. Macchi, D. Bi, G. Scita, M. L. Manning, R. Cerbino and M. C. Marchetti, Soft Matter, 2018, 14, 3471–3477.
67. L. Giomi, N. Hawley-Weld and L. Mahadevan, Proc. R. Soc. A, 2013, 469, 20120637.
68. O. Dauchot and V. Démery, Phys. Rev. Lett., 2019, 122, 068002.
69. C. Tapia-Ignacio, L. L. Gutierrez-Martinez and M. Sandoval, J. Stat. Mech.:Theory Exp., 2021, 2021, 053404.
70. P. Baconnier, O. Dauchot, V. Démery, G. Düring, S. Henkes, C. Huepe and A. Shee, Rev. Mod. Phys., 2025, 97, 015007.
71. C. Malinverno, S. Corallino, F. Giavazzi, M. Bergert, Q. Li, M. Leoni, A. Disanza, E. Frittoli, A. Oldani and E. Martini, et al., Nat. Mater., 2017, 16, 587–596.
72. E. F. Teixeira, H. C. M. Fernandes and L. G. Brunnet, Soft Matter, 2021, 17, 5991–6000.
73. P. Baconnier, D. Shohat, C. H. López, C. Coulais, V. Démery, G. Düring and O. Dauchot, Nat. Phys., 2022, 18, 1234–1239.
74. M. Paoluzzi, D. Levis and I. Pagonabarraga, Commun. Phys., 2024, 7, 57.
75. M. Casiulis, E. Arbel, C. van Waes, Y. Lahini, S. Martiniani, N. Oppenheimer and M. Y. Ben Zion, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2502211122.
76. M. Musacchio, A. P. Antonov, H. Löwen and L. Caprini, J. Chem. Phys., 2025, 162, 244902.
77. Y. Kinoshita, N. Uchida and A. M. Menzel, J. Chem. Phys., 2025, 162, 054906.
78. L. Caprini, U. Marini Bettolo Marconi, A. Puglisi and H. Löwen, J. Chem. Phys., 2023, 159, 041102.
79. C. Hernández-López, P. Baconnier, C. Coulais, O. Dauchot and G. Düring, Phys. Rev. Lett., 2024, 132, 238303.
80. M. Musacchio, A. P. Antonov, H. Löwen and L. Caprini, arXiv, 2025, preprint, arXiv:2506.12967 DOI:10.48550/arXiv.2506.12967.
81. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. P. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212–216.
82. Y. Kuroda, H. Matsuyama, T. Kawasaki and K. Miyazaki, Phys. Rev. Res., 2023, 5, 013077.
83. D. Nishiguchi and M. Sano, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2015, 92, 052309.
84. L. Caprini, U. M. B. Marconi, C. Maggi, M. Paoluzzi and A. Puglisi, Phys. Rev. Res., 2020, 2, 023321.
85. L. Caprini, U. Marini Bettolo Marconi and A. Puglisi, Phys. Rev. Lett., 2020, 124, 078001.
86. S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek and E. Bertin, Nat. Commun., 2020, 11, 1–12.
87. T. Sven van and H. Löwen, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2008, 78, 020101.
88. F. J. Sevilla, Phys. Rev. E, 2016, 94, 062120.
89. W. Chan Chung, D. Wu and R. Zhang, Nat. Commun., 2024, 15, 1043.
90. L. Caprini, U. M. B. Marconi, B. Liebchen and H. Löwen, arXiv, 2025, preprint, arXiv:2509.05053 DOI:10.48550/arXiv.2509.05053.
91. E. Kalz, A. Sharma and R. Metzler, J. Phys. A: Math. Theor., 2024, 57, 265002.
92. Y.-E. Keta, J. U. Klamser, R. L. Jack and L. Berthier, Phys. Rev. Lett., 2024, 132, 218301.
93. J. U. Klamser and L. Berthier, arXiv, 2025, preprint, arXiv:2510.04575.
94. V. Bratanov, F. Jenko and E. Frey, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 15048–15053.
95. R. Alert, J. Casademunt and J.-F. Joanny, Annu. Rev. Condens. Matter Phys., 2022, 13, 143–170.
96. N. Kumar, H. Soni, S. Ramaswamy and A. Sood, Nat. Commun., 2014, 5, 4688.
97. N. Koumakis, A. Gnoli, C. Maggi, A. Puglisi and R. Di Leonardo, New J. Phys., 2016, 18, 113046.
98. L. Walsh, C. G. Wagner, S. Schlossberg, C. Olson, A. Baskaran and N. Menon, Soft Matter, 2017, 13, 8964–8968.
99. O. Chor, A. Sohachi, R. Goerlich, E. Rosen, S. Rahav and Y. Roichman, Phys. Rev. Res., 2023, 5, 043193.
100. A. P. Antonov, M. Musacchio, H. Löwen and L. Caprini, Nat. Commun., 2025, 16, 7235.
101. M. Musacchio, M. Felber, M. Paoluzzi, A. Gnoli, A. Puglisi and L. Angelani, arXiv, 2025, preprint, arXiv:2511.21472.

---