Repulsive torques alone trigger crystallization of constant speed active particles

Abstract

We investigate the possibility for self-propelled particles to crystallize without reducing their intrinsic speed. We illuminate how, in the absence of any force, the competition between self-propulsion and repulsive torques determines the macroscopic phases of constant-speed active particles. This minimal model expands upon existing approaches for an improved understanding of crystallization of active matter.

---

Can crystals form out of constituents always in motion? Thermal systems, such as atomic and molecular matter, typically crystallize when the temperature is decreased, slowing down their constituents. This correlation between slowing down and crystallization also manifests itself in active materials,^1–4 despite their constituents being self-propelled by intrinsic forces. This connection even extends to other solidification processes in such materials. A prominent example is collection of active Brownian particles which exhibit partial solidification known as motility-induced phase separation.^5,6 As its name conveys, solidification in this case is entangled with speed reduction.⁷ Similarly, self-propelled Voronoi models^8,9 describing biological tissues exhibit jamming, another solidification transition with reduction of the cells' mobility.¹⁰ All those systems share the common feature of constituents interacting through forces, inducing the reduction of their speed. Yet, active particles are also found to display rich collective behaviours when only interacting through torques. Then, only the orientation of self-propulsion is altered and their speed is preserved.^11–15 From the seminal Vicsek's model,¹⁶ those systems have proved fruitful from a theoretical perspective^17–20 and successful to account for living²¹ and synthetic^22,23 systems. Since such active systems cannot, by design, comply with the above paradigm of solidification, their possibility to crystallize is intriguing. In this Communication, we prove that a minimal system of active particles interacting only via repulsive torques does crystallize. We rationalize its fluid and crystal phases through a detailed account of binary scattering collisions. We expect this system to complement force-based models for an ultimate comprehensive description of crystallization in out-of-equilibrium systems,^24–30 as well as to account for experimental observations.³¹

Constant-speed, torque-only, model system

We introduce a minimal model of active particles moving at constant speed v₀ and only interacting through repulsive torques. Particles, with instantaneous position r(t) and orientation p̂(t) = (cos θ, sin θ), evolve in a two-dimensional box with periodic boundaries:

ṙ_i(t) = v₀p̂_i(t), (1)

where r_ij = r_i − r_j = |r_i − r_j| exp(iψ_ij), and the noise with amplitude D_θ obeys 〈ξ_i(t)〉 = 0 and 〈ξ_i(t)ξ_j(t′)〉 = δ_ijδ(t − t′). We use the forward Euler method to carry out the numerical integration of eqn (1) and (2) with N = 270 particles. Particle i reorients to align in a time τ along r_ij away from its neighbours j located within an interaction range R_int, see Fig. 1a.

Fig. 1 Liquid-to-crystal transition follows from fast re-orientation/slow self-propulsion. (a) Schematics of the repulsive interaction between two particles. Particles (black) re-orient so as to align (orange) opposite to their center-to-center vector (blue). (b) Instantaneous configurations (left) and particles trajectories (right) increasing the repulsion over self-propulsion (top to bottom). Particles are coloured according to |y₆| evidencing the building-up of orientational order. Trajectories are all of the same duration illustrating the increasing localization of the dynamics. Dot diameter: R_int. Black rectangle indicates trajectories magnification.

Here, we focus on the competition between self-propulsion at v₀ and repulsion of strength τ⁻¹ in the limit of small noise D_θ = 0.4 × 10⁻⁴v₀/R_int, thereby keeping the intrinsic persistence length constant to a value 2.5 × 10⁴R_int. This competition is parametrized by the dimensionless parameter τ⁻¹/(v₀/R_int) which compares the time for a particle to travel a distance of one interaction range R_int/v₀ with the time for completion of re-orientation τ. In what follows, we express all quantities in units of R_int for lengths and R_int/v₀ for times, and study the system varying its two control parameters τ⁻¹ and the number density ρ.

Fig. 1b and Supplementary Videos (ESI†) illustrate that varying τ⁻¹ at constant density ρ = 1.03 alters the nature of the macroscopic state of the system: it spontaneously transitions from a liquid state to a hexagonal crystal state. Indeed as τ⁻¹ increases, the system structures itself and local orientational order builds up, as shown in Fig. 1b. The local bond-orientational order parameter defined as ψⁱ₆ = 〈e^i6ψ_ij〉_{j∈[scr N, script letter N]i} goes from a broad distribution peaked around |ψ₆| ≃ 0.4 at τ⁻¹ = 0.4 to a narrowly peaked distribution at |ψ₆| ≃ 1 at τ⁻¹ = 200.‡ Together with this structuring, the dynamics of the constant-speed particles localize, as illustrated by the equal-duration trajectories shown in Fig. 1b.

In this Communication, we explore first the properties of these liquid and crystal states as τ⁻¹ increases, and then the mechanisms that shape the full two-dimensional phase diagram shown in Fig. 6.

Fluid phase and binary collisions

We first characterize the properties of the fluid phase occuring at low τ⁻¹ and constant ρ = 1.03 and elucidate their microscopic origin. We quantify the dynamics of the particles by measuring their mean-squared displacement 〈Δr²〉 (t) = 〈(r_i(t + t₀) − r_i(t₀))²〉_i,t0, shown in Fig. 2a. For all τ⁻¹ < 28, the dynamics is diffusive at long times: 〈Δr²〉 = 4Dt^α with α = 1, see Fig. 2a inset. The dynamics, however, slow down as τ⁻¹ increases, which is captured by the decrease of the diffusion coefficient D on Fig. 2b. Note that any value of D is much lower than the translational diffusion caused by the orientationnal noise of eqn (2). This intrinsic diffusivity D_noise≡v₀²/(2D_θ) = 1.25 × 10⁴ does not explain the liquid properties which can follow only from the interactions between particles. Interestingly, the diffusion coefficient decays algebraically as D ∼ (τ⁻¹)^−0.6, see Fig. 2b, thereby suggesting a unique origin for this decrease.

Fig. 2 Dynamics in the fluid phase. (a) Mean-squared displacement, relative to the center of mass, as a function of time for various τ⁻¹ and ρ = 1.03 showing the diffusive dynamics characterized by the exponent α measured at long times (pink-shaded region), see inset. The mean-squared displacements are measured at steady state from systems initialized in random configurations. Thin blue line: crystal phase at τ⁻¹ = 80. (b) The diffusion coefficient decreases with τ⁻¹ (red dots measured at long time) quantifying the slowing down of the system dynamics. Black dots: diffusion coefficient assuming only uncorrelated binary scattering collisions, showing good agreement up to τ⁻¹ ≃ 5. Dotted line: D = D_noise. Dashed line indicates algebraic decay with −0.6 slope.

A naive but simple hypothesis is that particles experience uncorrelated binary scattering events alone. In this case, the diffusion coefficient would be given by the variance δθ² of the scattering angle distribution as D = v₀²/(2v₀R_intρδθ²).¹² To assess the validity of this description, we investigate scattering events numerically. We integrate the trajectories of two particles upon collisions, with initial distance R_int and orientations θⁱⁿ₁ and θⁱⁿ₂ both uniformly distributed in [−π, +π]. Fig. 3a illustrates such collisions for three different values of τ⁻¹ and two sets of initial orientations. Fig. 3b shows that increasing τ⁻¹ broaden the distribution of scattering angles δθ = θ^out_i − θⁱⁿ_i, with θ^out_i the orientation after collision. The diffusion coefficients computed from the variance of these distributions collapse on the fluid coefficients up to τ⁻¹ ≃ 5, as shown in Fig. 2b. This agreement proves that the properties of the fluids are quantitatively captured by uncorrelated binary collisions, a surprising result given their relatively high density ρ = 1.03.

Fig. 3 Binary scattering features. (a) Trajectories of pairs of particles upon binary collisions for two sets of initial orientations {θⁱⁿ₁, θⁱⁿ₂} (depicted on the schematics) and three values of τ⁻¹. (b) The probability density function of the scattering angles [scr P, script letter P](δθ) broadens with τ⁻¹. (c) The probability density function of the penetration depth [scr P, script letter P](d) reflects the strengthening of excluded-volume upon binary collisions. The inset shows the decrease of the mean value of the distribution with τ⁻¹.

Above τ⁻¹ > 5, however, the diffusion coefficient deviates from this simplified dynamics. The observation of binary scattering trajectories of Fig. 3a gives a hint into the mechanism at play: as τ⁻¹ increases, not only scattering is enhanced but particles also stay further apart upon colliding. This observation is quantified in Fig. 3c through the distributions of the minimal distance during collision . These distributions reveal that, around τ⁻¹ ≃ 5, an excluded area around the particle position emerges: they virtually never come at contact, [scr P, script letter P](d = 0) ≃ 0. Importantly, this change in binary collision behaviour survives at the scale of the many-body system where depleted area forms and strengthen around each particles, as captured by the pair-correlation function g(r) in Fig. 4. Driven by this exclusion effect, the fluid evolves from an uncorrelated gas to a more and more structured liquid as τ⁻¹ increases.

Fig. 4 Structuring of the fluids. Radial pair-correlation functions signal an evolution from a gas to a liquid state of the fluid upon increasing τ⁻¹ coinciding with the enhancement of exclusion effect during binary collisions. g(r) is the orthoradial average of . Thin blue line: crystal phase at τ⁻¹ = 80. ρ = 1.03.

Altogether, the slowing down of the dynamics and the structuring of the system result from the specifics of collisions between pairs of constant-speed particles repelling via torques. They put the system on the path of crystallization.

Crystal phase and effective mean-field cages

Crystallization indeed occurs at high τ⁻¹, as captured by both structural and dynamical features. Bond-orientational order (Fig. 1) and positional order (Fig. 4) extend over the whole system size. The dynamics of the particles localize as shown by the saturation of the mean-squared displacement at τ⁻¹ = 80 in Fig. 2a.

The spontaneous formation of a stable yet dynamical crystal is quite unexpected. Indeed, a constant-speed particle interacting via repulsive torques with random lattices of quenched obstacles does never localize.^12,13 We note two main differences between the many-body system described by eqn (1) and (2) and these isolated dynamics: particles are interacting with (i) constantly moving particles, which, (ii) assemble into a crystalline lattice instead of being randomly positioned. To elucidate which of this two differences accounts for the successful localization of the particles dynamics, we now introduce a mean-field description of the system in its crystalline state. We investigate the dynamics of a constant-speed particle initially positioned within a mean-field unit cell consisting of a regular hexagon of quenched particles. This configuration is depicted in Fig. 5a where typical trajectories are displayed for various values of τ⁻¹. At low τ⁻¹ particles escape the cell whereas at high τ⁻¹ the cell act as an efficient trap. To assess the validity of this mean-field description, we quantify the localization of the dynamics through the typical cage sizes , measured both from 10³ mean-field trajectories and from many-body systems initialized into crystalline configurations, see Fig. 5b. Overall, the mean-field description predicts well the variations of the cage sizes in the crystal; with fluctuations in the positions of particles further accounting for the slightly larger sizes in the latter case. Besides, the threshold for efficient caging τ_c⁻¹ = 16 agrees with the value yielding the loosest, possibly metastable, crystals, τ⁻¹ = 36, see Fig. 5b and c.

Fig. 5 Efficient caging in the crystals. (a) Six quenched particles (black dots) on a regular hexagon constitute mean-field cages. The lattice spacing respects the crystal density. Pink circles are of radius R_int and denote the range at which the active particles, initially at the center, interact with the quenched ones. The typical trajectories (red to yellow) of constant-speed particles show two behaviours: at low τ⁻¹ particles escape whereas at high τ⁻¹ they are trapped. (b) The cage size decreases with τ⁻¹. Black dots are the sizes of mean-field cages. Purple dots are effective cage sizes in the crystals. (c) Mean-squared displacements, relative to the center of mass, saturate at long-time translating the localization of particles' dynamics. The mean-squared displacements are measured at steady state from systems at ρ = 1.03 initialized in crystalline configurations.

Together, these results indicate that trapping by quenched crystalline lattices suffices to explain the localization of particles dynamics in the crystal phase. In particular, no emergent correlation between particles' velocities is required for crystallization.

Mechanisms of crystallization

The above analyses put two mechanisms forward for the spontaneous crystallization: the increase of excluded-volume and the efficient caging within hexagonal cells. Which of this two mechanisms sets the onset of crystallization?

To elucidate the role of each mechanism, we systematically construct the two-dimensional phase diagram varying the density ρ and the repulsion strength τ⁻¹, see Fig. 6. It shows that crystals emerge in a confined region of the phase diagram.

Fig. 6 Necessary conditions for stable crystals. Phase diagram of the system in the (ρ, τ⁻¹)-plane. The dots reflect the state of the system according to the bond-orientational order parameter 〈|ψ₆|〉, measured at steady state from systems initialized in random configurations. The crystal is confined by both criteria on excluded-volume (pink line) and on efficient caging (purple line). Dashed lines at ρ = 0.92 and ρ = 1.15 indicate respectively the crystallization threshold and the close-packed limit for hard-disks.

First, we rely on the measure of the mean minimal distance upon binary collision (see Fig. 3c) to estimate excluded volume effects. Precisely, we map our system of soft particles onto hard disks of effective radius R_eff(τ⁻¹) = 〈d〉(τ⁻¹) and define the effective density ρ_eff(τ⁻¹) = ρ_HD/R_eff², where ρ_HD ≈ 0.92 is the number density of crystallization in hard-disk systems.³² Plotting the curve ρ = ρ_eff in Fig. 6 defines a region where crystallization would occur were the excluded-volume effect to trigger the transition. At high τ⁻¹ ≥ 10², and upon increasing the density, crystallization occurs at ρ_eff: this perfect agreement with the hard disks description suggests that entropy dominates in this region.

In stark contrast, this hard-disk ansatz does not hold at smaller τ⁻¹. This discrepancy is elucidated by accounting for efficient caging instead. Repeating the mean-field analysis for various densities we determine the threshold for caging τ_c⁻¹(ρ), which varies little with ρ, see Fig. 6. This criterion based on interactions with nearest neighbours accurately defines the onset of crystallization in this region.

Overall, the above analysis evidences that the spontaneous crystallization of the system requires that both mechanisms, excluded-volume and efficient caging, simultaneously operate. These two simple lines of arguments, based on binary exclusion and mean-field caging, perform remarkably well to locate the crystallization region in the (ρ, τ⁻¹)-plane.

Considering how the crystal accommodates densities higher than closed-packed hard disks value shows how non-trivial this agreement is. At those densities, the crystal does not reduce its lattice spacing. Instead some lattice positions get occupied by pairs of particles, leading to a small reduction of 〈|ψ₆|〉, see Fig. 6. This genuine manifestation of the softness of the interactions between particles stresses the emergent character of the system crystallization.

Conclusion

Active particles moving at constant speed and interacting only through repulsive torques transition from fluids to crystals at high density and high repulsion strength (or, equivalently, low speeds). The properties of the fluid phase follow from the specifics of the binary collisions: scattering sets the diffusion up until the structuring of the fluid originating from exclusion during collisions. While the individual particles always keep their constant speed, they crystallize through the interplay of two mechanisms: sufficient excluded volume and efficient caging by the lattice neighbours.

By demonstrating the crystallization of a minimal system of constant-speed active particles and gaining insights into the microscopic origins of its various states, we provide an approach complementary to repulsive-force systems, such as active Brownian particles, for investigating out-of-equilibrium melting. Refining the transition region and looking for possible hexatic³² or other²⁴ intermediate phases are sensible next steps in this direction but require accessing larger system sizes and thus appropriate computational resources. Besides this fundamental perspective, this simple model can serve as a starting point to build and decipher eventual experimental realizations.

Author contributions

M. L. B. and A. M. designed and performed the research, and wrote the manuscript.

Conflicts of interest

There are no conflicts to declare.

References

1. A. M. Menzel and H. Löwen, Phys. Rev. Lett., 2013, 110, 055702.
2. C. A. Weber, C. Bock and E. Frey, Phys. Rev. Lett., 2014, 112, 168301.
3. G. Briand and O. Dauchot, Phys. Rev. Lett., 2016, 117, 098004.
4. S. E. Moran, I. R. Bruss, P. Shoenhofer and S. C. Glotzer, Soft Matter, 2022, 18, 1044–1053.
5. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
6. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 055701.
7. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 235702.
8. A. Pasupalak, L. Yan-Wei, R. Ni and M. P. Ciamarra, Soft Matter, 2020, 16, 3914–3920.
9. S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek and E. Bertin, Nat. Commun., 2020, 11, 1–9.
10. D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, Phys. Rev. X, 2016, 6, 021011.
11. H. Chaté, F. Ginelli, G. Grégoire and F. Raynaud, Phys. Rev. E, 2008, 77, 046113.
12. O. Chepizhko and F. Peruani, Phys. Rev. Lett., 2013, 111, 160604.
13. A. Morin, D. L. Cardozo, V. Chikkadi and D. Bartolo, Phys. Rev. E, 2017, 96, 042611.
14. N. Bain and D. Bartolo, Nat. Commun., 2017, 8, 1–6.
15. Y. Zhao, T. Ihle, Z. Han, C. Huepe and P. Romanczuk, Phys. Rev. E, 2021, 104, 044605.
16. T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen and O. Shochet, Phys. Rev. Lett., 1995, 75, 1226.
17. J. Toner and Y. Tu, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 58, 4828.
18. A. Peshkov, E. Bertin, F. Ginelli and H. Chaté, Eur. Phys. J.: Spec. Top., 2014, 223, 1315–1344.
19. A. P. Solon, J.-B. Caussin, D. Bartolo, H. Chaté and J. Tailleur, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 062111.
20. R. Kürsten and T. Ihle, Phys. Rev. Lett., 2020, 125, 188003.
21. W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale and A. M. Walczak, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 4786–4791.
22. A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot and D. Bartolo, Nature, 2013, 503, 95–98.
23. A. Morin and D. Bartolo, Phys. Rev. X, 2018, 8, 021037.
24. J. Bialké, T. Speck and H. Löwen, Phys. Rev. Lett., 2012, 108, 168301.
25. J. U. Klamser, S. C. Kapfer and W. Krauth, Nat. Commun., 2018, 9, 1–8.
26. G. Briand, M. Schindler and O. Dauchot, Phys. Rev. Lett., 2018, 120, 208001.
27. P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella and I. Pagonabarraga, J. Phys.: Conf. Ser., 2019, 012073.
28. E. S. Bililign, F. Balboa Usabiaga, Y. A. Ganan, A. Poncet, V. Soni, S. Magkiriadou, M. J. Shelley, D. Bartolo and W. Irvine, Nat. Phys., 2021, 1–7.
29. C. J. Reichhardt and C. Reichhardt, Nat. Phys., 2021, 1–2.
30. P. Digregorio, D. Levis, L. F. Cugliandolo, G. Gonnella and I. Pagonabarraga, Soft Matter, 2022, 18, 566–591.
31. J. Zhang, R. Alert, J. Yan, N. S. Wingreen and S. Granick, Nat. Phys., 2021, 17, 961–967.
32. E. P. Bernard and W. Krauth, Phys. Rev. Lett., 2011, 107, 155704.

---

Footnotes

† Electronic supplementary information (ESI) available: Videos. See DOI: 10.1039/d2sm00256f

‡ The average is performed over the neighbourhood [scr N, script letter N]_i given by Delaunay triangulation.

---