Repulsive torques alone trigger crystallization of constant speed active particles

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 determine the macroscopic phases of constant-speed active particles. This minimal model expands upon existing approaches for an improved understanding of crystallization of active matter.

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][25][26][27][28][29][30] , as well as to account for experimental observations 31 .

Constant-speed, torque-only, model system
We introduce a minimal model of active particles moving at constant speed v 0 and only interacting through repulsive torques. Particles, with instantaneous position r(t) and orientationp(t) = (cos θ , sin θ ), evolve in a two-dimensional box with periodic boundaries:ṙ where r i j = r i − r j = r i − r j exp(iψ i j ), and the noise with amplitude D θ obeys ξ i (t) = 0 and ξ i (t)ξ j (t ) = δ i j δ (t −t ). We use the forward Euler method to carry out the numerical integration of Eqs. 1 and 2 with N = 270 particles. Particle i reorients to align in a time τ along r i j away from its neighbours j located within an interaction range R int , see Fig. 1a.
Here, we focus on the competition between self-propulsion at v 0 and repulsion of strength τ −1 in the limit of small noise D θ = 0.4 × 10 −4 v 0 /R int , thereby keeping the intrinsic persistence length constant to a value 2.5 × 10 4 R int . This competition is parametrized by the dimensionless parameter τ −1 /(v 0 /R int ) which compares the time for a particle to travel a distance of one interaction range R int /v 0 with the time for completion of reorientation τ. In what follows, we express all quantities in units of R int for lengths and R int /v 0 for times, and study the system varying its two control parameters τ −1 and the number density ρ. Figure 1b and Supplementary Videos † illustrate that varying τ −1 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 τ −1 increases, the system structures itself and local orientational order builds up, as shown in Fig. 1b. The local bond-orientational order parameter defined as ψ i 6 = e i6ψ i j j∈N i goes from a broad distribution peaked around |ψ 6 | 0.4 at τ −1 = 0.4 to a narrowly peaked distribution at |ψ 6 | 1 at τ −1 = 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 τ −1 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 τ −1 and constant ρ = 1.03 and elucidate their microscopic origin. We quantify the dynamics of the particles by measuring their mean-squared displacement ∆r 2 (t) = (r i (t + t 0 ) − r i (t 0 )) 2 i,t 0 , shown in Fig. 2a. For all τ −1 < 28, the dynamics is diffusive at long times: ∆r 2 = 4Dt α with α = 1, see Fig. 2a inset. The dynamics, however, slow down as τ −1 increases, which is captured by the decrease of the diffusion coefficient D on Fig. 2b. § The average is performed over the neighbourhood N i given by Delaunay triangulation.  Note that any value of D is much lower than the translational diffusion caused by the orientationnal noise of Eq. 2. This intrinsic diffusivity D noise ≡ v 2 0 /(2D θ ) = 1.25 × 10 4 does not explain the liquid properties which can follow only from the interactions between particles. Interestingly, the diffusion coefficient decays algebraically as D ∼ (τ −1 ) −0.6 , see Fig. 2b, thereby suggesting a unique origin for this decrease.
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 δ θ 2 of the scattering angle distribution as D = v 2 0 /(2v 0 R int ρδ θ 2 ) 12 . 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 θ in 1 and θ in 2 both uniformly distributed in [−π, +π]. Figure 3a illustrates such collisions for three different values of τ −1 and two sets of initial orientations. Figure 3b shows that increasing τ −1 broaden the distribution of scattering angles δ θ = θ out i − θ in 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 τ −1 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.
Above τ −1 > 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 τ −1 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 d = min t (r 12 (t)) {θ in 1 ,θ in 2 } . These distributions reveal that, around τ −1 5, an excluded area around the particle position emerges: they virtually never come at contact, 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 paircorrelation 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 τ −1 increases. 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 τ −1 , 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 τ −1 = 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 Eqs. 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 τ −1 . At low τ −1 particles escape the cell whereas at high τ −1 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 ∆r 2 ∞ = lim t→∞ ∆r 2 , measured both from 10 3 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 τ −1 c = 16 agrees with the value yielding the loosest, possibly metastable, crystals, τ −1 = 36, see Figs. 5b-c. Together, these results indicate that trapping by quenched crys-talline 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 τ −1 , see Fig. 6. It shows that crystals emerge in a confined region of the phase diagram. First, we rely on the measure of the mean penetration depth 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 (τ −1 ) = d (τ −1 ) and define the effective density ρ eff (τ −1 ) = ρ HD /R 2 eff , where ρ HD ≈ 0.92 is the number density of crystallization in hard-disk systems 32 . 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 τ −1 ≥ 10 2 , 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 τ −1 . This discrepancy is elucidated by accounting for efficient caging instead. Repeating the mean-field analysis for various densities we determine the threshold for caging τ −1 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 (ρ, τ −1 )-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 |ψ 6 | , 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 32 or other 24 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.