Lattice-dependent orientational order in active crystals

Abstract

Via mechanisms not accessible at equilibrium, self-propelled particles can form phases with positional order, such as crystals, and with orientational order, such as polar flocks. However, the interplay between these two types of order remains relatively unexplored. Here, we address this point by studying crystals of active particles that turn either towards or away from each other, which can be experimentally realised with phoretic or Janus colloids or with elastically-coupled walker robots. We show that, depending on how these interactions vary with interparticle distance, the particles align along directions determined by the underlying crystalline lattice. To explain the results, we map the orientational dynamics of the active crystal onto a lattice of spins that interact via (anti-)ferromagnetic alignment with each other plus nematic alignment with the lattice directions. Our findings indicate that orientational and positional order can be strongly coupled in active crystals, thus suggesting strategies to control orientational order by engineering the underlying crystalline lattice.

---

In active matter, microscopic constituents inject mechanical energy, thus driving the system out of equilibrium. As a result, active particles can self-organize in ways not accessible at equilibrium. In particular, the field has focused on how positional and orientational order can emerge.^1–7

Orientational order, such as the polar order found in flocks, can arise from direct alignment interactions between the orientations of self-propelled particles, as originally demonstrated in the Vicsek model.⁸ More recent work showed that flocking can also emerge when active particles attract each other,⁹ align their orientation with their velocity^7,10,11 or, alternatively, when particles turn away from one another.^12–14

Regarding positional order, self-propelled particles have been found to crystallise via either motility-induced phase separation,^15–20 attractive interactions,^9,21–25 or simply at densities approaching close packing.^26–32 Recent work also showed that, in confinement, self-propelled particles can form Wigner crystals that emerge through repulsive interactions, which keep the particles at a distance.^{13,14,33–36} Particles in active crystals were also found to orient and move collectively as a flock, thus displaying not just positional but also orientational order.^{9,13,14,29,37–41} Beyond such flocking crystals, the interplay between positional and orientational order in active matter has been recently explored in the XY model with vision–cone interactions,^42–45 in crystals of self-aligning walker robots,^39,40 and in crystallites of Quincke rollers.²⁵

Here, we address this question by studying crystals of self-propelled particles that turn either towards or away from each other. These interactions, which emerge for example in metal-dielectric Janus colloids,^13,46 couple the polarity of one particle to the orientation of the bond with a neighboring one. Hence, such polarity-bond interactions produce a crosstalk between positional and orientational order. We show that, on a lattice, polarity-bond interactions yield either effective alignment or anti-alignment between particle polarities, like in the XY model. In addition, they also produce nematic alignment of the particle polarities and the lattice axes. We first study the interplay between these two effects for particles on a chain. We find that the particles can achieve either local ferro- or antiferromagnetic order, either along or perpendicular to the chain. We then consider a square lattice and find that the particles can orient locally along the lattice axes and/or form domains of polar order, depending on the distance dependence of the underlying interactions. On the triangular lattice, the polarity-bond interactions can be frustrated. Overall, our findings show that, through polarity-bond interactions, the orientational order of active crystals can depend strongly on the lattice structure. Thus, our work suggests strategies to obtain desired states of orientational order in active crystals by engineering specific particle interactions and crystalline lattices.

Active crystals with polarity-bond interactions

We consider active particles on a fixed crystalline lattice. Neighboring lattice sites i and j are separated by the vector r⁽⁰⁾_ij = a(cosϕ_ij, sinϕ_ij), where a is the lattice constant and ϕ_ij define the lattice angles (Fig. 1a). The particles are bound to lattice sites by elastic forces −kΔr_i, with elastic constant k and displacement Δr_i (Fig. 1a). These elastic forces correspond to the harmonic approximation of any force that confines the particles to their lattice sites. In addition, the particles self-propel at speed v₀ along their orientation n̂_i = (cos θ_i, sin θ_i).

Fig. 1 Active crystal with polarity-bond interactions. (a) Schematic of an active crystal made of self-propelled particles bound to lattice sites. The green arrows indicate distance vectors. The purple arrows indicate the self-propulsion force ξ_tv₀n̂_i and the elastic restoring force −kΔr_i, whose balance sets the equilibrium displacement of a particle from its lattice site. (b) Polarity-bond interactions are torques, defined in eqn (1), whereby particles turn either towards or away from others. (c) Examples of distance dependences f(r) of the interaction torques for different systems, with their corresponding dimensionless parameter Ω ≡ af′(a).

The particles interact through turn-towards or turn-away torques given by

Γ_ji = Γ₀f(|r_ij|)n̂_i × r̂_ij, (1)

which arise in self-aligning active particles,⁷ as well as from electrostatic interactions in Janus particles with a metallic (dark) and a dielectric (light) hemisphere^13,46 (Fig. 1b). The torque Γ_ji exerted by particle j on particle i, with amplitude Γ₀ and a general distance-dependence given by f(r) > 0 (Fig. 1c), tends to turn particle i either towards (Γ₀ > 0) or away from (Γ₀ < 0) the distance vector r_ij = r_j − r_i = |r_ij|r̂_ij connecting it with particle j. We define Γ₀ as the torque amplitude at a distance given by the lattice constant, such that f(a) = 1. All together, the particles follow the overdamped Langevin equationswhere ξ_t and ξ_r are the translational and rotational friction coefficients, and D_t and D_r are the translational and rotational diffusivities associated with the corresponding Gaussian white noises η^t_i and η^r_i. Here, we indicated the torque as a scalar quantity as it only has a component along the ẑ axis.

Self-propulsion displaces particles away from the lattice sites. Particles reach a displacement Δr_i = ln̂_i, with displacement length l = ξ_tv₀/k, in a time scale τ_e = ξ_t/k set by the elastic restoring force (Fig. 1a). As in recent work,⁴⁰ we assume that this elastic relaxation time is much smaller than the time scale of the angle evolution: τ_e ≪ τ_θ = ξ_r/Γ₀. Under this approximation, particle positions adiabatically follow the slower orientation dynamics. Ignoring translational noise, which is negligible in front of rotational noise for Janus particles,^13,46 particle positions are given by

r_i(t) = r⁽⁰⁾_i + Δr_i(t) = r⁽⁰⁾_i + ln̂_i(t), (4)

where r⁽⁰⁾_i is the position of the lattice site of particle i.

Active crystals as spin lattices

Under the approximation of fast elastic relaxation, particle positions can be eliminated in favor of the orientations; hence, the active crystal reduces to a spin lattice. To this end, we insert the positions of eqn (4) in eqn (1) and obtainHere, we made the torque dimensionless by rescaling time as t̃ = D_rt. As a result, the torque strength is parametrized by the (signed) dimensionless torque amplitude . In eqn (5), the original polarity-bond interaction n̂_i × r̂_ij between the particles decomposes into two effects: (i) turning either towards or away from the neighbouring lattice site, n̂_i × r⁽⁰⁾_ij, and (ii) either alignment or anti-alignment with the neighbour's orientation, n̂_i × n̂_j.

Assuming nearest-neighbor interactions, and that the displacement l is much smaller than the lattice constant a, we expand the radial dependence in powers of l/a as (see Section S1 of the SI)

Here, we defined the dimensionless distance-dependence parameter Ω ≡ af′(a), which quantifies how the interaction torque depends on distance. It is negative (positive) for torques that decay (grow) with distance (Fig. 1c). Introducing eqn (6) in eqn (5), we obtainto first order in l/a. The first two terms represent the orienting towards neighbouring lattice sites, at the zeroth and first order of the l/a expansion. The third term describes the effective neighbour alignment or antialignment.

We now sum over nearest neighbours to obtain the torque on particle i:

Here, the first term of eqn (7) has cancelled because the lattice vectors r⁽⁰⁾_ij add up to zero for a regular lattice. Thus, all the remaining contributions are of order l/a. Hence, this geometric factor sets the magnitude of the torques on a lattice together with the dimensionless torque magnitude . Moreover, we expressed all the contributions in terms of particle orientations θ_i and lattice angles ϕ_ij (Fig. 1a, Section S2 of the SI). The first term in eqn (8) corresponds to an alignment or anti-alignment torque ∝ sin(θ_j − θ_i) like that of the XY model with ferro- or anti-ferromagnetic interactions (Fig. 2a). The second term produces nematic alignment of a particle i with the lattice axes, given by ϕ_ij, which effectively behave as an external nematic field acting on the spins (Fig. 2b). Finally, the third term produces alignment or anti-alignment of particle i with the mirror image of the neighboring particle j; the lattice axis connecting them, encoded in the angle ϕ_ij, acts as the mirror plane (Fig. 2c).

Fig. 2 Polarity-bond interactions on a lattice. On a lattice (dashed axis), and given a reference particle j (gray), the original interaction torques yield three contributions (eqn (8)): alignment or anti-alignment with the neighbour's orientation ∝ sin(θ_j − θ_i) (a), with the lattice axes ∝ sin2(ϕ_ij − θ_i) (b), and with the neighbour's mirror image ∝ sin(2ϕ_ij − θ_j − θ_i) (c). Green (red) particles indicate the orientations favoured when the prefactor of the corresponding term in eqn (8) is positive (negative), which depends on the sign of the interaction torques, Γ₀/|Γ₀|, and the distance-dependence parameter Ω.

The sign of each of these terms depends on the value of Ω, which is determined by the distance dependence f(r) of the interaction torques (Fig. 1c). For metal-dielectric Janus colloids,^13,46 their electrostatic dipole–dipole interactions give f(r) = a⁴/r⁴, which gives Ω(a) = −4. For particles reorienting in the chemical concentration field produced by others,¹⁴ we have f(r) = a²/r², which gives Ω = −2. Similarly, systems where torques arise from short-ranged repulsive interactions will have Ω < 0. Other possible cases are topological interactions, which are distance-independent, and hence have Ω = 0. Such topological interactions could either be programmed in robots or arise naturally in animals that turn towards or away from their nearest neighbors regardless of their distance. Yet another option is torques due to elastic forces,¹¹ for which f(r) = (r − [small script l])/(a − [small script l]), and hence Ω = 1/(1 − [small script l]/a) can be either positive or negative depending on the ratio between the spring's rest length [small script l] and the lattice constant a. Elastic forces were proposed to model the soft interactions between cells,^10,33,47 and they were realised in crystals made of hexbugs connected by springs.³⁹ Overall, different systems realise different values of the distance-dependence parameter Ω (Fig. 1c). Hence, below we explore its role and we find that it controls the type and strength of orientational order in our active crystals.

Interestingly, the torque in eqn (8) can be derived from an effective energy H, such that the dynamics of the particle orientations θ_i read

andThis effective energy has contributions due to an XY-type alignment H^XY_ij = −cos(θ_j − θ_i), lattice alignment H^LA_ij = [cos2(ϕ_ij − θ_i) + cos2(ϕ_ij − θ_j)]/2, and mirror alignment H^MA_ij = −cos(2ϕ_ij − θ_i − θ_j) (Fig. 2). Whereas the distance-dependence parameter Ω controls the sign and relative strength of these different contributions as discussed above, the turn-towards or turn-away character of the interaction torques controls the global sign of the effective energy function. Therefore, switching between turn-towards and turn-away torques^13,46 causes a complete inversion of the energy landscape, whereby stable equilibrium points become unstable and vice versa. Such a switch is known as a landscape-inversion phase transition,⁴⁸ which is of mixed order⁴⁹ and displays unique phase-ordering processes.⁵⁰

One-dimensional chain

To study the emerging orientational order in active crystals with polarity-bond interactions, we start by considering a one-dimensional chain with periodic boundary conditions (Fig. 3a). We perform Brownian dynamics simulations of eqn (3) and (8) with N = 10⁵ particles using the Euler method with a time step . From the simulations, we characterize the emergence of nematic order as a function of the dimensionless parameters and Ω of the torque interactions (Fig. 3b). In Section S3 of the SI we show that N = 10⁵ is sufficiently large to avoid finite-size effects, and that the time evolution of the nematic order is not sensitive to different realisations of the random orientations in the initial condition.

Fig. 3 States of particles on a chain. (a) Schematic of a chain of active particles. (b) Nematic order S as a function of the dimensionless parameters of the torque interactions: the distance-dependence parameter Ω and the torque amplitude . (c) and (d) Polar and nematic order parameters as a function of Ω for turn-away (c, ) and turn-towards (d, ) torques. Colour shadings indicate different states, shown in snapshots in (e)–(k).

On a chain, each particle has two neighbours with lattice angles ϕ_ij = 0,π. In this case, the effective energy eqn (10) reduces to that of an anisotropic XY model for spins n̂_i = (cos θ_i, sin θ_i) in a nematic field which aligns them either parallel or perpendicular to the chain axis (see Section S4 of the SI):

Here, the superscripts x and y indicate spatial components. The distance-dependence parameter Ω controls both the anisotropy of the interactions, reflected in the first two terms, as well as the alignment with the chain axis, encoded in the last term.

Despite the presence of effective alignment interactions, the chain does not exhibit global polar order; the polar order parameter vanishes (dashed lines in Fig. 3c and d). The situation is reminiscent of the XY model in 1d, for which the Hohenberg–Mermin–Wagner theorem forbids the breaking of the continuous rotation symmetry, and hence the emergence of long-range polar order.^51–53 Here, however, the theorem does not apply because the anisotropy of the interactions as well as the lattice alignment in eqn (11) already break the continuous rotational symmetry.

To rationalise the absence of polar order in our model, we adapt Peierls' argument for the lack of ferromagnetic order in the 1d Ising model.^53–55 We can extend this argument to our case because the effective energy has the discrete symmetry θ_i → θ_i + π. Hence, we consider an excitation in the form of a domain of π-flipped spins, such that the system configuration looks like …↗↙…↙↗…. Because of the symmetry of the effective energy function, the domain bulk costs no extra energy; the only energy penalty comes from the domain walls, whose relative contribution decreases with increasing system size N. However, the number of ways to place the domain walls, which determines the system's entropy, increases with system size. In the thermodynamic limit N → ∞, and for any non-zero temperature (here noise strength D_r > 0), this entropic contribution wins and prevents the emergence of polar order. This argument does not rule out the existence of local polar order, as seen in Fig. 3k. On large scales, however, no polar order persists.

Yet, our active chains are not always disordered. They can display global nematic order (Fig. 3b). We quantify it through the scalar nematic order parameter , which is the largest eigenvalue of the nematic order-parameter tensor Q_αβ = 〈2n^α_i(t)n^β_i(t) − δ_αβ〉_i,t, where α and β are indices for spatial components. In our system, nematic order arises from the lattice-alignment contribution in the last term of eqn (11), which acts as an external nematic field with strength controlled by the distance-dependence parameter Ω. For Ω = 1, the lattice-alignment contribution vanishes. In this case, the effective energy eqn (11) corresponds to that of the XY model, for which the Hohenberg–Mermin–Wagner theorem forbids global order. Accordingly, we obtain states with no global nematic order (black horizontal stripe in Fig. 3b) but with local order, either ferromagnetic or anti-ferromagnetic (Fig. 3f and j).

For other values of Ω, there can be global nematic order (Fig. 3b). For small torque amplitudes Γ₀, fluctuations allow the system to sample different configurations. For large torque amplitudes Γ₀, the interactions favour specific configurations (Fig. 3e–k), which we describe and label with arrow symbols below.

For turn-away interactions (Γ₀ < 0), we find two states (Fig. 3c): anti-aligned perpendicular to the chain (Fig. 3e, ↑↓) and anti-aligned along the chain (Fig. 3g, →←), in addition to the state with only local anti-ferromagnetic order for Ω = 1 (Fig. 3f, ↗↙). Respectively, for turn-towards interactions (Γ₀ > 0), we find three states (Fig. 3d): anti-aligned along the chain (Fig. 3h, →←), aligned perpendicular to the chain (Fig. 3i, ↑↑), and aligned along the chain (Fig. 3k, ←←), in addition to the state with only local ferromagnetic order for Ω = 1 (Fig. 3j, ↗↗). We note that any of the aligned states described here displays only local polar order. In the following, we explain the emergence of these states by analyzing the equilibrium configurations of two spins.

Two-particle configurations

To understand the states on the chain, we study the equilibrium configurations of two particles on a lattice, described as coupled spins θ₁,θ₂ governed by the effective energy in eqn (11). Each point in the θ₁,θ₂-plane corresponds to a spin configuration as shown in Fig. 4a. We obtain their effective energy from eqn (11) and show it in Fig. 4b. For turn-away (turn-towards) interactions, the ground state is given by the minimum (maximum) of .

Fig. 4 Equilibrium states for two particles. (a) Arrow representations of the two-particle states in the θ₁θ₂ plane. (b) Effective energy H for different values of Ω. Light (dark) areas are energetically favourable for turn-towards (turn-away) torques, corresponding to . The turquoise (orange) contours enclose regions with 90% of the probability p(θ₁,θ₂) ∝ e^−H for Γ₀ = 10a/l (Γ₀ = −10a/l). (c) Most probable two-particle states for , which explain the observations made from the simulations in Fig. 3.

Beyond the ground state, since the angle dynamics in eqn (9) is equivalent to a system of interacting Brownian particles, the probability density follows the Boltzmann distribution p(θ₁,θ₂) ∝ e^{−H(θ₁,θ₂)}. The turquoise (orange) contour lines in Fig. 4b enclose regions with 90% of the probability for turn-towards (turn-away) torques at . Using Fig. 4a as a reference, we identify the states corresponding to these high-probability regions. The results, shown in Fig. 4c, match with the states found in our simulations (Fig. 3). Thus, the equilibrium behavior of two particles explains the variety of states found for the many-body system.

Note that, for Ω = 1, the effective energy reduces to that of the XY model, which is rotationally invariant, and hence the minimum becomes degenerate. Accordingly, the probability for turn-towards (turn-away) torques concentrates around the ferromagnetic θ₂ = θ₁ (anti-ferromagnetic θ₂ = θ₁ + π) ground state, without any preferential alignment with the chain axis (Fig. 4b, Ω = 1). For turn-away interactions, the ground states are non-degenerate for Ω ≠ 1: they are the ↑↓ configuration for Ω < 1 and →← for Ω > 1. For turn-towards interactions, for Ω < 0, the ground state is also the →← configuration. For Ω ≥ 0, the ground state is degenerate, given by any ferromagnetic configuration θ₂ = θ₁ (Fig. 4b, turquoise diagonal stripes). However, this degeneracy is broken once fluctuations are taken into account, as they allow the particles to explore the shape of the effective energy around the minimum. Analyzing the probability p(θ₁,θ₂) ∝ e^{−H(θ₁,θ₂)} reveals the most likely configurations: ↑↑ for 0 < Ω < 1, and ←← for Ω > 1, as shown in Fig. 4c, which match those in Fig. 3.

Lattice-dependent order in two dimensions

To illustrate that the connection between lattice structure and orientational order extends to two dimensions, we now consider a square and a triangular lattice, respectively with lattice angles , where n = 0,1,2,…. For both lattices, the lattice-alignment term of the effective energy in eqn (10) vanishes, because particles tend to (anti-)align with two perpendicular or three symmetric axes. This allows us to study the competition between the neighbour-alignment and mirror-alignment terms in eqn (10), tuned by the distance-dependence parameter Ω.

Square lattice

On the square lattice, there is a mapping between any given configuration with turn-away interactions (Γ₀ < 0) and another one with turn-towards interactions (Γ₀ > 0), as shown in Section S5 of the SI. Therefore, we focus on the turn-towards case with and perform simulations of 500 × 500 spins with a time step dt = 0.0005 for a time t = 5000 from an initial condition with random orientations.

To explore the role of the distance-dependence parameter, we consider the values Ω = −0.5,0,0.5, for which Fig. 5a and b show snapshots at large and small scales. The amplitude of neighbour and mirror alignment is proportional to |Ω + 1| and |Ω − 1|, respectively. For Ω = 0.5, neighbor XY alignment is stronger. Consequently, the system forms polar domains and topological defects, similar to the XY model. However, the weak contribution of mirror alignment creates a preference to orient along the lattice, as reflected in the orientational distribution function shown in the inset. For Ω = 0, the neighbor alignment and mirror alignment contributions have equal strengths. In this case, particles orient along one lattice axis, forming a state with nematic order consisting of oppositely-pointing stripes of different widths. For Ω = −0.5, mirror alignment dominates, and the particles are anti-aligned along the direction of orientation and aligned perpendicular to it. This arrangement results in regular stripes with alternating orientation.

Fig. 5 Ordering on a square lattice. (a) and (b) Snapshots of 500 × 500 spins (a) and close-up views of 10 × 10 spins (b) for turn-towards interactions with and different distance-dependence parameters Ω = −0.5, 0, 0.5. Color represents particle orientation. The insets in (a) show the angular distribution. (c) Effective energy H between a particle oriented along the x axis and two neighbors: one along the x axis with orientation θ_→, and one along the y axis with orientation θ_↑ (see schematic). The turquoise contours enclose regions with 90% of the probability p(θ_→,θ_↑) ∝ e^−H.

To understand these patterns, we consider a particle with a fixed orientation along the x direction (θ = 0) and we study the effective interaction energy eqn (10) for varying orientations of the neighbour in the direction of orientation, θ_→, and of the neighbour perpendicular to it, θ_↑. For Ω ∈ (−1,1) the neigbour alignment term tends to align both θ_→ and θ_↑ with the reference particle, while the mirror alignment term tends to anti-align θ_→ and align θ_↑ with θ = 0. Fig. 5c shows the interaction energy H, with contour lines enclosing 90% of the probability p(θ_→,θ_↑) ∝ e^{−H(θ_→,θ_↑)}. For Ω = −0.5, mirror alignment prevails, creating anti-alignment along the orientation direction and alignment perpendicular to it. This is consistent with the aligned stripes of alternating direction seen on large scales (Fig. 5a, left). For Ω = 0.5, neighbour alignment is stronger, resulting in aligned regions (Fig. 5a, right). For Ω = 0, both interaction terms tend to align θ_↑ with the reference particle. In contrast, the alignment and misalignment effects on θ_→ cancel, such that the interaction does not set the orientation θ_→. This is consistent with our observation of stripes that are strongly correlated perpendicular to the particle orientation, but that randomly alternate in the direction of orientation. In all cases, the configurations predicted from this 3-particle picture based on the interaction energy agree with the simulation results.

Triangular lattice

We now consider a triangular lattice and perform simulations for both turn-away and turn-towards torques of 500 × 500 spins with a time step dt = 0.0005 for a time t = 1000 starting from an initial condition with random orientations (Fig. 6).

Fig. 6 Ordering and frustration on a triangular lattice. (a), (d), (e) and (h) Close-up snapshots of simulations of 500 × 500 particles on a triangular lattice. Particle orientations are shown in color as in Fig. 5. The insets show the distribution of particle orientations. (a) and (d) The case with Ω = 1 corresponds to the XY model, with either antiferromagnetic (Γ₀ < 0) or ferromagnetic (Γ₀ > 0) interactions. (b) and (c) Schematics of frustrated interactions. (b) For turn-away torques, satisfying the antiferromagnetic interactions of a central particle (green bonds) results in unfavourable interactions between its neighbours (red bonds). (c) The system then reaches a compromise state. (e) and (h) The case with Ω = −1 corresponds to only mirror-alignment interactions, and hence the particles align relative to the lattice. (f) For turn-away torques, satisfying the mirror-alignment interactions of a central particle (green bonds) results in unfavourable interactions between neighbours (red bonds). (g) The system then finds a frustrated compromise state. (h) For Ω = −1, frustration is also present for turn-towards torques.

We start with the well-known case of the XY model, which we retrieve by setting Ω = 1 (see eqn (10)). In this case, the energy is rotation-invariant and the particles do not align with the lattice, as shown by the orientational distributions in the insets in Fig. 6a and d. For turn-towards torques (Γ₀ > 0), the particles experience ferromagnetic XY interactions, and hence they develop local polar order (Fig. 6d). For turn-away torques (Γ₀ < 0), the XY interactions are antiferromagnetic. In a triangular lattice, not all particle pairs can be simultaneously antiparallel (Fig. 6b), which is known as geometric frustration. As a result, the system reaches states like the one shown in Fig. 6a, which emerge as a compromise between achieving some anti-alignment between particles while avoiding alignment, as sketched in Fig. 6c.

For Ω ≠ 1, the effective energy eqn (10) breaks rotational invariance, and the particles orient relative to the lattice (Fig. 6e and h). To showcase the effects of the lattice, we focus on Ω = −1, for which only the mirror-alignment term in eqn (10) is present. For turn-away torques, Fig. 6f shows that satisfying the mirror-alignment interactions for the central particle, again, results in unfavourable interactions between the neighbours. To avoid them, the system reaches a compromise state consisting of alternating aligned stripes, shown in Fig. 6e and sketched in Fig. 6g. For turn-towards torques, the interactions are also frustrated (Fig. 6h). Thus, mirror-alignment interactions are frustrated by the triangular lattice for both signs of the interaction. For Ω ≠ ±1, an interplay between neighbour and mirror alignment results in generally anisotropic and frustrated states. Particles tend to align relative to the lattice but, due to frustration, the resulting states can no longer be predicted by minimising the interaction energy of two particles as in previous sections.

Overall, by extending our analysis of the one-dimensional chain, these results show that lattice-dependent orientational order can also arise in two-dimensional lattices. In addition, we found that the polarity-bond interactions between our particles can be frustrated in the triangular lattice. An interesting question for future work is whether our system can exhibit long-range polar order or not, as discussed for non-reciprocal XY models.^42–45 Another direction is to consider three-dimensional lattices.

Discussion and outlook

In summary, we studied crystals of active particles that turn either towards or away from one another. Because these interactions, which we call polarity-bond interactions, couple the orientation of a particle with the position of another, they establish a link between positional and orientational order. We showed that, when particle positions equilibrate fast compared to their orientations, the orientations can be described as spins that evolve according to an energy. In this energy, the original polarity-bond interactions give rise to both conventional aligning terms like those of the XY model but also unconventional terms that couple the particle orientations to the lattice directions. This energy allowed us to predict the variety of states that we found in direct Brownian dynamics simulations. Thus, our work contributes to ongoing efforts to establish a Hamiltonian description for systems with non-reciprocal interactions which, like our turn-towards or turn-away torques (eqn (1)), do not obey Newton's law of action and reaction.⁵⁶

Recent work on active solids showed that the interplay between positional and orientational dynamics gives rise to activity-driven oscillations termed collective actuation.^7,39,57 Here, we explored a different regime by focusing on the limit in which particle positions equilibrate fast compared to their orientations.⁴⁰ In this regime, our results show that active crystals can display several states with orientational order, with particles aligned in a variety of ways with respect to the lattice directions. The precise state that is favoured depends on whether the interaction torques are turn-towards or turn-away, as well as how they vary with distance.

Thus, our findings reveal that polarity-bond interactions enable one to control the orientational order of active crystals through the lattice structure. Experimentally, such control could be achieved in systems of either metal-dielectric Janus colloids,^13,46,58 which interact electrostatically through turn-towards or turn-away torques, or macroscopic robots, which can be programmed to do so. Under confinement, active Janus colloids form crystals at high densities due to their repulsive interactions.¹³ These repulsive interactions, when approximated for small displacements of the particles around their lattice sites, would give rise to the elastic forces considered in our model. Alternatively, the particles can be placed in engineered lattices made either with grooved substrates^59–61 or with periodic optical potentials generated with interfering lasers.^62–64 In such lattices, both the structure and the lattice constant can be controlled. In our model, these changes would affect the lattice angles and the value of the distance-dependence parameter Ω, which would then impact the orientational order of the active crystal.

From a theoretical standpoint, our findings introduce the notion of lattice-dependent orientational order, which describes states in which rotational symmetry is broken through a coupling to the lattice structure. By revealing that the lattice structure can impact the orientational order in active crystals, our work complements previous studies of active solids, which mainly focused on how activity distorts or even melts their crystalline structure.^29,30,65 Our work also complements recent studies on the impact of spatial anisotropy, such as the one imposed by a lattice, on flocks.^44,45,66 More generally, our findings call for further developments of general continuum theories of active solids:^25,67–70 A challenge for future work is to generalise them to incorporate information about the lattice structure which, as we have found, can affect orientational order.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The code produced for this paper is available at this link https://github.com/tillwelker/Lattice-dependent-orientational-order-in-active-crystals.

Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm00627a

Acknowledgements

We thank Marín Bukov for discussions, and for pointing out the mapping of the active chain onto the anisotropic XY model. T.W. thanks Patrick Pietzonka, Holger Stark, and Sarah A.M. Loos for their ideas, inspiration, and support. Open Access funding provided by the Max Planck Society.

References

1. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189.
2. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
3. M. C. Marchetti, Y. Fily, S. Henkes, A. Patch and D. Yllanes, Curr. Opin. Colloid Interface Sci., 2016, 21, 34–43.
4. É. Fodor and M. C. Marchetti, Physica A, 2018, 504, 106–120.
5. M. Bär, R. Großmann, S. Heidenreich and F. Peruani, Annu. Rev. Condens. Matter Phys., 2020, 11, 441–466.
6. H. Chaté, Annu. Rev. Condens. Matter Phys., 2020, 11, 189–212.
7. P. Baconnier, O. Dauchot, V. Démery, G. Düring, S. Henkes, C. Huepe and A. Shee, Rev. Mod. Phys., 2025, 97, 015007.
8. T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen and O. Shochet, Phys. Rev. Lett., 1995, 75, 1226–1229.
9. L. Caprini and H. Löwen, Phys. Rev. Lett., 2023, 130, 148202.
10. B. Szabó, G. Szöllösi, B. Gönci, Z. Jurányi, D. Selmeczi and T. Vicsek, Phys. Rev. E, 2006, 74, 061908.
11. E. Ferrante, A. E. Turgut, M. Dorigo and C. Huepe, Phys. Rev. Lett., 2013, 111, 268302.
12. M. Kneževic, T. Welker and H. Stark, Sci. Rep., 2022, 12, 19437.
13. S. Das, M. Ciarchi, Z. Zhou, J. Yan, J. Zhang and R. Alert, Phys. Rev. X, 2024, 14, 031008.
14. A. G. Subramaniam, S. Adhikary and R. Singh, arXiv, 2025, preprint, arXiv:2504.07050 DOI:10.48550/arXiv.2504.07050.
15. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 055701.
16. J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936–940.
17. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
18. P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella and I. Pagonabarraga, Phys. Rev. Lett., 2018, 121, 098003.
19. M. N. van der Linden, L. C. Alexander, D. G. A. L. Aarts and O. Dauchot, Phys. Rev. Lett., 2019, 123, 098001.
20. A. K. Omar, K. Klymko, T. GrandPre and P. L. Geissler, Phys. Rev. Lett., 2021, 126, 188002.
21. B. M. Mognetti, A. Šaric, S. Angioletti-Uberti, A. Cacciuto, C. Valeriani and D. Frenkel, Phys. Rev. Lett., 2013, 111, 245702.
22. R. Singh and R. Adhikari, Phys. Rev. Lett., 2016, 117, 228002.
23. S. Thutupalli, D. Geyer, R. Singh, R. Adhikari and H. A. Stone, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 5403–5408.
24. A. Mauleon-Amieva, M. Mosayebi, J. E. Hallett, F. Turci, T. B. Liverpool, J. S. Van Duijneveldt and C. P. Royall, Phys. Rev. E, 2020, 102, 032609.
25. S. J. Kole, X. Chao, A. Mauleon-Amieva, R. Hanai, C. P. Royall and T. B. Liverpool, arXiv, 2025, preprint, arXiv:2501.15996 DOI:10.48550/arXiv.2501.15996.
26. J. Bialké, T. Speck and H. Löwen, Phys. Rev. Lett., 2012, 108, 168301.
27. C. A. Weber, C. Bock and E. Frey, Phys. Rev. Lett., 2014, 112, 168301.
28. G. Briand and O. Dauchot, Phys. Rev. Lett., 2016, 117, 098004.
29. G. Briand, M. Schindler and O. Dauchot, Phys. Rev. Lett., 2018, 120, 208001.
30. X.-q Shi, F. Cheng and H. Chaté, Phys. Rev. Lett., 2023, 131, 108301.
31. M. F. Zhang, B. Y. Fan, C. Y. Zhang, K. Chen, W. D. Tian and T. H. Zhang, Soft Matter, 2025, 21, 927–934.
32. N. Sakaï, K. Skipper, F. J. Moore, J. Russo and C. P. Royall, Soft Matter, 2025, 21, 5204–5213.
33. B. Smeets, R. Alert, J. Pešek, I. Pagonabarraga, H. Ramon and R. Vincent, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 14621–14626.
34. J.-B. Delfau, C. López and E. Hernández-García, New J. Phys., 2017, 19, 095001.
35. M. Le Blay and A. Morin, Soft Matter, 2022, 18, 3120–3124.
36. Q. Yang, M. Jiang, F. Picano and L. Zhu, Nat. Commun., 2024, 15, 2874.
37. G. Grégoire, H. Chaté and Y. Tu, Phys. D Nonlinear Phenom., 2003, 181, 157–170.
38. A. M. Menzel and H. Löwen, Phys. Rev. Lett., 2013, 110, 055702.
39. 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.
40. C. Hernández-López, P. Baconnier, C. Coulais, O. Dauchot and G. Düring, Phys. Rev. Lett., 2024, 132, 238303.
41. M. Musacchio, A. P. Antonov, H. Löwen and L. Caprini, arXiv, 2025, preprint, arXiv:2506.12967 DOI:10.48550/arXiv.2506.12967.
42. S. A. M. Loos, S. H. L. Klapp and T. Martynec, Phys. Rev. Lett., 2023, 130, 198301.
43. G. Bandini, D. Venturelli, S. A. Loos, A. Jelic and A. Gambassi, J. Stat. Mech. Theory Exp., 2025, 053205.
44. P. Popli, A. Maitra and S. Ramaswamy, arXiv, 2025, preprint, arXiv:2503.06480 DOI:10.48550/arXiv.2503.06480.
45. D. Dopierala, H. Chaté, X.-Q. Shi and A. Solon, arXiv, 2025, preprint, arXiv:2503.14466 DOI:10.48550/arXiv.2503.14466.
46. J. Zhang, R. Alert, J. Yan, N. S. Wingreen and S. Granick, Nat. Phys., 2021, 17, 961–967.
47. S. Henkes, Y. Fily and M. C. Marchetti, Phys. Rev. E, 2011, 84, 040301.
48. R. Alert, J. Casademunt and P. Tierno, Phys. Rev. Lett., 2014, 113, 198301.
49. R. Alert, P. Tierno and J. Casademunt, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 12906–12909.
50. R. Alert, P. Tierno and J. Casademunt, Nat. Commun., 2016, 7, 13067.
51. N. D. Mermin and H. Wagner, Phys. Rev. Lett., 1966, 17, 1133–1136.
52. P. M. Chaikin and T. C. Lubensky, Principles of condensed matter physics, Cambridge University Press, 1995.
53. N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group, Addison-Wesley, 1992.
54. R. Peierls, Math. Proc. Cambridge Philos. Soc., 1936, 32, 477–481.
55. K. Huang, Statistical Mechanics, John Wiley & Sons, 2nd edn, 1987.
56. Y.-B. Shi, R. Moessner, R. Alert and M. Bukov, arXiv, 2025, preprint, arXiv:2505.05246 DOI:10.48550/arXiv.2505.05246.
57. H. Xu, Y. Huang, R. Zhang and Y. Wu, Nat. Phys., 2023, 19, 46–51.
58. J. Yan, M. Han, J. Zhang, C. Xu, E. Luijten and S. Granick, Nat. Mater., 2016, 15, 1095–1099.
59. A. van Blaaderen, R. Ruel and P. Wiltzius, Nature, 1997, 385, 321.
60. K.-H. Lin, J. Crocker, V. Prasad, A. Schofield, D. Weitz, T. Lubensky and A. Yodh, Phys. Rev. Lett., 2000, 85, 1770–1773.
61. A. Ortiz-Ambriz and P. Tierno, Nat. Commun., 2016, 7, 10575.
62. M. M. Burns, J.-M. Fournier and J. A. Golovchenko, Science, 1990, 249, 749–754.
63. M. Brunner and C. Bechinger, Phys. Rev. Lett., 2002, 88, 248302.
64. K. Mangold, P. Leiderer and C. Bechinger, Phys. Rev. Lett., 2003, 90, 158302.
65. L. Ophaus, E. Knobloch, S. V. Gurevich and U. Thiele, Phys. Rev. E, 2021, 103, 032601.
66. A. Solon, H. Chaté, J. Toner and J. Tailleur, Phys. Rev. Lett., 2022, 128, 208004.
67. A. Maitra and S. Ramaswamy, Phys. Rev. Lett., 2019, 123, 238001.
68. C. Scheibner, A. Souslov, D. Banerjee, P. Surówka, W. T. M. Irvine and V. Vitelli, Nat. Phys., 2020, 16, 475–480.
69. Y. Shen, J. O'Byrne, A. Schoenit, A. Maitra, R.-M. Mège, R. Voituriez and B. Ladoux, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2421327122.
70. Y. E. Keta and S. Henkes, Soft Matter, 2025, 21, 5710–5719.

---