Slipping and fluidisation in active crystalline rotors

Abstract

In equilibrium, commensurability between sliding surfaces underpins our understanding of nanoscale friction and yielding in crystals. However, these concepts have only recently begun to be imported into the realm of active matter. Here, we develop an experimental platform and a theoretical description for active colloids confined and self-assembled into small crystals. In our experiments with Quincke Rollers, these crystallites form due to cohesive interactions between the particles. We apply a circular confinement which controls the population of the system. We focus on perfect hexagonal crystallites of 61 particles whose dimensions are compatible with the confinement. Competition between solidity and self-propulsion leads to self-shearing and complex flow-inversion behaviour, along with self-sliding states and activity-induced melting. We discover active stick-slip dynamics, which periodically switch between a commensurate static state and an incommensurate self-sliding state characterised by a train of localised defects and describe the steady-state behaviour using a discretised model of active hydrodynamics. We then investigate the intermittent stick-slip dynamics using an extension of the Frenkel–Kontorova (FK) model, a fundamental workhorse of slipping and flow in crystals appropriate to our geometry. Our findings in a colloidal model system point to a wealth of phenomena in active solids as design principles for both assembly and robotics down to the nanoscale.

---

1 Introduction

Self-assembly of active systems goes beyond the minimisation of free energy to design exotic structures and dynamics and takes many forms, from flock formation,^1–5 to clustering through motility-induced phase separation.^6,7 Perhaps the simplest self-assembled structures are crystals. When created from active components, crystalline systems exhibit significant modifications to the defect behaviour⁸ and synchronisation of collective modes.⁹ Study of active crystals opens a new perspective on living systems from bacteria¹⁰ to starfish embryos, the latter in particular exhibiting beautiful chiral oscillations.¹¹ However, such biological systems are very complex, and considerable insight can be gained from examples of these exotic phenomena in a minimally simple setting. To this end, tunable colloidal^12,13 and granular¹⁴ model systems provide a means to realise such behaviour in a system which is amenable to a comprehensive understanding.

In passive matter, crystalline solids present a well-understood microscopic framework in which yielding and slip may be understood, through defects, dislocations and gliding of crystal planes.¹⁵ By analogy therefore, active colloidal crystals provide a promising means to explore mechanisms underlying failure and fluidisation. Computer simulation predicts a rich behaviour, of defect mediated transitions,¹⁶ extreme deformations¹⁷ and strong fluctuations while the two-step melting behaviour in passive systems in 2d is preserved.¹⁸ Activity also introduces novel phonon-like collective modes¹⁹ and coarsening mechanisms.²⁰

There are rather few experiments with active colloidal crystals. Some insights have been gained from active granular systems, where the microscopic mechanism of spontaneous flow may be related to defects and dislocations.¹⁴ Recently, improved control over active colloids at high concentration has permitted some colloidal solids to be investigated.^21–23 Active colloidal crystals exhibit novel annealing²⁴ and (2d) melting phenomena.^25–27 Among the few studies to explore yielding and flow have considered rotors with strong non-reciprocal coupling through the use of spinners^28,29 where crystal structure interplays with active dynamics to form characteristic patterns.^28,30 Meanwhile, active colloidal amorphous solids have received some attention with the observation of re-entrant dynamical behaviour.³¹

Here we address the behaviour of a well-characterised experimental colloidal model system of Quincke rollers in a dc external electric field under circular confinement.^2,32,33 These particles are passive until a threshold field strength E_Q is exceeded, yet even under these conditions,³⁴ long-ranged electrohydrodynamic attractions lead to the formation of well-defined crystallites. This quasi-2d model system of active colloids enables us to explore the competition between crystalline order and self-propulsion. In this way, our system may be regarded as prototypical for studies of crystallisation in (polar) active matter.

We characterise the activity-induced rotation, melting, and complex flow behaviour due to the frustration between cohesive interactions which promote a crystalline solid and activity which causes yielding. In our model system, we focus on a perfect crystallite with N_max = 61 particles, which acts as a rotor under active motion of the particles. By rotor, here we refer to solid body undergoing rotation. Note that the rotor spins in the xy plane, perpendicular to the plane of rotation of each Quincke roller particle. Here we observe that increasing self-propulsion leads to a sequence of rigid-body rotation, faster angular velocity at the edge and, finally, an inversion to faster angular velocity at the centre. We are motivated to gain insight into this flow inversion behaviour, and therefore we introduce a version of the Toner-Tu model³⁵ which describe the hydrodynamics in active fluids. Here we discretise the model to reflect the layers in our experimental system.

We then move to a time-resolved description of the rotor and consider the intermittent self-sliding behaviour. To do so, we apply the Frenkel–Kontorova (FK) model of nanotribological phenomena^36,37 to our active clusters. In the FK model, slipping is associated with a breakdown in local crystalline order. We expect analogues of sliding friction to be generically present in active colloids that form crystallites, due to the frustration between particle motility and periodic order, and our system presents a suitable opportunity to probe this.

This paper is organised as follows. In Section 2, we describe our experimental methodology, our discretised version of the Toner-Tu equations is explained in Section 3. Our results are presented in Section 4, which we divide as follows. In Section 4.1 we outline the non-equilibrium phase behaviour that we observe. This is followed by an analysis of the steady-state slipping behaviour in Section 4.2 and the predictions of our minimal model in Section 4.3. We then probe the stick-slip behaviour in Section 4.4, before applying the Frenkel–Kontorova model in Section 4.5. We discuss and conclude in Section 5.

2. Experimental methods

We employ a colloidal suspension of sterically stablised poly(methyl methacrylate) colloids, of diameter σ = 2.92 μm and polydispersity 3.1% as determined with SEM. The particles are suspended in a low conductivity solution of AOT surfactant at a concentration of 0.15 mM in hexadecane. The suspension is then injected into a sample cell made of ITO-coated glass slides (Solems, ITOSOL12), separated by a layer of adhesive tape (100 μm in thickness). Sedimentation occurs due to the density mismatch with the solvent, and the colloids form a quasi two-dimensional layer, as represented in Fig. 1(a). Quincke rotation is achieved with the application of a dc field E using a potential amplifier (TREK 609E-6) connected to the sample cell.

Fig. 1 Confinement induced self-assembly of Quincke rotors. (a). The confinement of Quincke rollers into circular regions is achieved by the application of a localised electric field E. (b). Different populations N, indicated by the labels, are found across the confining regions (E ≈ E_Q). Scale bar represents 50 μm. (c). Order parameters as a function of N/N_max as the system becomes active, i.e. at E ≈ E_Q. Polarisation Π_φ and bond orientational order, ψ₆, parameters are shown (see Section 4.1 for definitions of order parameters). (d). Schematic dynamical phase diagram summarising structural and dynamical states. (e) and (f). Phase diagrams indicating (e) the hexagonal ψ₆ and (f) polar Π_φ order for passive and active rollers under strong confinement. For electric fields above E_Q (solid horizontal line), clusters become active and give rise to apolar (disordered) gas and polar (ordered) fluid phases. Colourbars indicate the magnitude of the order parameters. Symbols are: (black hexagons) passive aggregates, (black squares) active apolar (Π_φ ≤ 0.5) and (black circles) polar states (Π_φ > 0.5). White circles indicate polar fluids (ψ₆ ≤ 0.5), whereas solid black circles indicate active crystallites (ψ₆ > 0.5). Colours in data points for N_max correspond to data in Fig. 2. Representative coordinates are shown in the bottom row for the state points indicated.

Confinement is introduced by patterning the top electrode of the sample cell. Circular regions of radius R_c ≈ 5σ are produced employing conventional photolithography methods. The thickness of the photoresist layer is ≈1σ, sufficient for charge screening on the top electrode. A non-zero electric current thus results only within the circular regions upon the application of the field E. This yields an electro-osmotic flow which strongly confines rollers within the regions of interest. The origin of this flow is the transport of electric charges in the direction of the circular regions, which produces an inward flow at the bottom electrode.³⁸ In the absence of the electric field, colloids are prone to escape the circular regions due to thermal motion.

At low amplitudes of the field, e.g. E < E_Q, where E_Q ≈ 0.8 V μm⁻¹ for an isolated particle, a population of N colloids is dragged towards the confining regions and self-assemble into clusters. Fig. 1(b) shows nine confining regions, each with a different population N. The maximum population observed corresponds to N_max = 61 which is a perfect hexagonal arrangement with four layers around a central particle. The Quincke rollers are imaged using brightfield microscopy (Leica DMI 3000B) with a 10× magnification, and recorded at high speed (900 fps, Basler ACE) in regions of 128 × 128 px. The motion of individual rollers is reconstructed using a Python version of a common tracking algorithm.³⁹

3. Model of discrete hydrodynamics in circular confinement

Here we introduce our discretised description of the rotor, which we approximate as a circle. This minimal model has two key ingredients and is constructed as follows. (i) a constant preferred self-propulsion velocity for an individual particle layer in isolation (captured by preferred speed u₀ and a “stiffness” ξ that speeds up or slows down a particle towards the preferred speed) and (ii) a sliding friction between adjacent layers (captured by [small eta, Greek, tilde]), generated by particle–particle interactions, which favour rigid-body rotation of the cluster as a whole. In addition, the particles at the boundary have velocity u_R ≠ u₀. This different velocity at the boundary captures the experimental observation that although each interior layer has two neighbouring layers, the outermost layer has only one neighbouring layer and therefore is affected less by its neighbours. The use of an effective friction reflects our Toner-Tu type treatment of the system as being one-component, as we neglect the multicomponent nature of our experimental system, which consists of colloids, solvent and ions.

Consider the steady state of linearised Toner-Tu type equations of homogeneous dry active matter for the local self-propulsion velocity v(r), which in this setting is

[small eta, Greek, tilde]∇²v(r) − ξ(v(r) − u₀) = 0, (1)

where the preferred interior velocity is azimuthal, u₀ = u₀[small phi, Greek, circumflex]. For a circular domain, in plane polars we have r = (r,ϕ) or r = rr̂ + ϕ[small phi, Greek, circumflex] and v = v_r(r)r̂ + v_ϕ(r)[small phi, Greek, circumflex] with 0 ≤ r ≤ R and ϕ ∈ [0,2π]. Here ϕ is the azimuthal angle. We consider the situation where v_r = 0 and radially symmetric, such that v_ϕ(r) satisfies

[small eta, Greek, tilde](∂²_rv_ϕ + r⁻¹∂_rv_ϕ − r⁻²v_ϕ) − ξ(v_ϕ − u₀) = 0 (2)

where the boundary conditions are v_ϕ(0) = 0,v_ϕ(R) = u_R.

This model combines an effective particle–particle friction, [small eta, Greek, tilde] which penalises relative motion between nearest neighbour particles, and self-propulsion with average speed u₀. For the circular domain, we assume that the motion is only in the azimuthal direction (i.e. the radial component of the velocity of each particle is zero) and that this velocity is the same in each layer which we label by the number n ∈ [1, 4]. Taking the particle diameter σ, we can denote the radial position of the centre of each particle in layer n by r_n = σ(n − 1). Given an azimuthal velocity v_n ≡ v_ϕ(r_n), we can define the discrete radial difference operator,

Applying this twice we obtainBalancing the local relative friction with self-propulsion leads to the equation,

η([Δ²_rv]_n + r_n⁻¹[Δ_rv]_n − r_n⁻²v_n) − (v_n − u₀) = 0, (5)

where η = [small eta, Greek, tilde]/ξ measures the local friction. The boundary conditions are v₀ = 0 and v₅ = u_R. Due to its position at the boundary with lower friction, we expect that the critical field is lower for the particles at the boundary than those in the bulk. In the bulk, at high area fraction, each particle will have around six neighbours, which impede the flow of ions that underlies Quincke rotation. Hence as the electric field is increased, the particles at the boundary will begin to move before the particles in the bulk, i.e. they will have a lower effective critical field. Above this boundary critical field but below the bulk critical field, u_R > u₀ = 0. If the bulk critical field is E_Q, then at E = E_Q, the boundary layer is already moving with a finite speed, hence above E_Q, we can assume that u₀ = A₀ΔE^α where α = 1/2² for ΔE = (E − E_Q) and . We note that because α < 1, (and assuming A_i, , B_i are of similar size) that the bulk layer speed can increase faster than the boundary layer speed as ΔE ≪ 1 increases from zero. This leads to the inverted angular velocity profiles observed.

To see this, we note that the dynamics of eqn (5) can be expressed by the matrix equation

where S_mn is a tridiagonal matrix and . Note that R = 5σ. We fix the boundary conditions v₀ = 0, v₅ = u_R by setting the values of the 1st and 6th row of the matrix. The matrix is
where is chosen to implement the boundary conditions. Then we can invert the matrix to solve for v_n. Once we have v_n we can obtain the local angular velocity, [small omega, Greek, vector] of each layer: ω_n = v_n/r_n is the angular velocity of layer n, where [[small omega, Greek, vector]]_n = ω_n. Parameters here are σ = 1, η = 10, A₀ = 1, A_R = 0.7, B_R = 0.5, so that

4. Results

4.1. Phase behaviour in the population-activity plane

The structural and dynamic behaviour in the population-activity (N–E) plane is summarised in Fig. 1b–f and Movie S1. To characterise the srtucture and dynamics of the active system, we shall use two order parameters to characterise the system, one of which is structural (ψ₆) and the other is dynamical (Π_φ).

To elucidate the local structural (hexagonal) order, we use the particle-averaged, time-averaged, bond orientational order parameter

whereThe quantity z_j is the co-ordination number of particle j from a Voronoi tessellation, and θ_k^j is the angle between the bond from i to k and a reference axis. The parameter ψ₆ runs between perfect ordering (ψ₆ = 1) and complete disorder (ψ₆ = 0).

The collective dynamics are characterised by the time-averaged polar order parameter,

which quantifies the degree of alignment between rollers. Here û_i is the local orientation of roller i, and ê_φ is a unit vector along the azimuthal direction. Π_φ = 1 indicates perfect azimuthal alignment, and Π_φ = 0 represents random particle orientations.

Fig. 1c shows the dynamical order parameter Π_φ and structural order parameter ψ₆ as a function of N at the onset of activity E ≈ E_Q. A transition from the isotropic gas with little alignment in the dynamics (and a low ψ₆) to a strongly aligned polar state is seen from the polarisation order parameter Π_φ at N/N_max ≈ 0.4, even for populations with relatively little structural order. (Recall that N_max = 61 for a perfect hexagonal crystal of four layers around a central particle, see Fig. 1e.) The structural order (ψ₆) does increase at slightly higher densities N/N_max ≈ 0.5 but the intermediate regime may be considered to be a polar fluid state. These two measures indicate the emergence of active crystallites at high densities.

In Fig. 1e, we consider the hexagonal order parameter ψ₆ in the (N, E) plane. For the passive case, E < E_Q, we see moderate values of ψ₆ even at low populations, N/N_max < 0.5. This corresponds to inactive clusters formed due to electrohydrodynamic attractions.³⁴ In this range of field strength ψ₆ increases with population, and the clusters become increasingly crystalline for N approaching N_max. Switching on activity by increasing the field strength to E > E_Q leads to a drop in ψ₆ as hexagonal ordering competes with activity.³⁴ For the highest field strengths of E/E_Q = 2.3, we observe the formation of a layered fluid structure in which particles are organised in concentric layers⁴⁰ (see right side of Fig. 1f). At the onset of Quincke rotation (E = E_Q), the role of density becomes apparent. With population , activity causes hexagonal ordering ψ₆ to decrease, which has also been seen for bulk systems,³⁴ but for larger populations, ψ₆ remains at values similar to those of passive crystallites.

Fig. 1f depicts the experimental phase diagram obtained from the polar order parameter, Π_φ, in the same (N, E) plane. The passive and active phases are separated by the onset of activity at field E = E_Q. Below E_Q, the polar order parameter is close to zero as the motion of the passive colloids is largely uncorrelated. However, when the system is active, we find a transition between an apolar gas (low Π_φ) at low population to a polar fluid (high Π_φ) as N increases. The boundary between these states appears to be insensitive to field strength.

Comparing Fig. 1e and f, we find the active apolar gas, which has a low value for both structural and dynamical order parameters. This then transitions to two ordered regimes reflecting the two order parameters: a polar fluid at moderate density (N/N_max ≳ 0.4) and an active crystalline phase at even higher densities in which both parameters Π_φ and ψ₆ approach 1. At the highest densities, we observe that sufficiently strong activity is still able to melt the crystalline order, a phenomenon we return to below. We note that using structural order parameters such as ψ₆ to classify a state point as a crystal does not by itself demonstrate that it has a long relaxation time, i.e. that it is solid.²³ We summarise the results we obtain for the full system in the schematic structural-dynamical phase diagram in Fig. 1d. This illustrates the states we observe, of passive ordered states for E < E_Q. In the active case at low density we find an active gas, which gives way to a polar fluid and finally ordered, active crystallites are found at high area fraction.

4.2. Steady-state slipping: activity vs. the boundary

To investigate the interplay between crystalline order and dynamics, we focus on a perfect crystal, which for this confinement radius corresponds to N = N_max = 61, and the activity-dependent rotational dynamics of both hexagonal and layered-fluid structures. Note that here the cohesive interactions lead to the hexagonal crystallites. The effect of the confinement (whose radius is a little larger than the crystallite) is mainly to control the number of particles, rather than to act as a confining wall which might perturb the structure.⁴¹

In its quiescent state for , this population forms a nearly perfect hexagonal shape, see Fig. 1. At high field strengths, i.e., E ≈ 2.3E_Q, the structure is fluid-like with four concentric layers around a central stationary particle. The active state exhibits field-dependent structures with different dynamics, see Fig. 2a. For both hexagonal and layered-fluid cases, the different layers are well resolved by the radial density profiles ρ(r)/[small rho, Greek, tilde] in Fig. 2b where [small rho, Greek, tilde] is the mean density. One signature of hexagonal structure is a split peak which merges into a single peak as the system becomes more fluidised at a high field strength.

Fig. 2 Perfect crystal rotor. N = N_max = 61. (a). Time-averaged angular velocity as a function of E in layers n = 2 (black squares) and n = 4 (white circles), with n counting outwards from the central particle (labelled n = 0). At intermediate field strength E ≈ 1.5E_Q the outer layer rotates faster than the inner layer, but this situation is inverted at high E > 1.6E_Q, where the inner layers rotate faster than the outer layer. Colours correspond to the state points shown in Fig. 1(e and f). (b) Scaled radial density profiles characterising structures obtained at different fields E, indicated in the legend. Numbers over the peaks indicate the layer number n. Curves are offset vertically for clarity. (c) Angular velocity profiles as function of the layer number n, where n = 0 is the central particle. Symbols are experimental data, solid lines are fits obtained using eqn (5) in the SI, and dotted lines are step guides. Data are offset for clarity. (d) Distortion order parameter, 〈Θ〉, as a function of layer n. (e) Hexagonal order parameter, 〈ψ₆〉, as a function of layer n. For all panels, colours correspond to the electric field values indicated by the colourbar in (b).

Upon increasing the electric field strength E, we observe a range of dynamical scenarios, including rigid-body rotation and periodic rearrangements, as illustrated by the angular velocity, ω, profiles in Fig. 2c. In particular, in order of ascending field strength: (i), for E < E_Q, no rotation is measured (grey triangles); (ii), at the onset of activity, we observe rigid-body rotation (yellow pentagons, Movie S2); (iii), the cluster rotates fastest at the outside (magenta circles, Movie S3); (iv), we find rigid-body rotation (purple triangles); (v), the cluster rotates fastest in the interior (blue hexagons and teal squares, Movies S4 and S5); therefore, we find an inversion in the flow profile as a function of field strength. This complex dynamical behaviour in Fig. 2c contrasts with the radial density, Fig. 2b, which shows little change until the very highest field strength. The bond-orientational parameter ψ₆ is a more sensitive measure of structure (Fig. 2e), and shows additional detail such as a drop in hexagonal order when the outer layers rotate faster (magenta data). However, the lack of sensitivity of these time-averaged structural measures motivate a time-resolved approach, which we pursue below.

4.3. Comparison with the minimal model

Before exploring the time-resolved behaviour, we consider the discretised Toner-Tu model introduced in Section 3 to describe the steady-state behaviour.³⁵ We recall that our model has two main ingredients: (i) a constant preferred self-propulsion velocity for an individual particle layer in isolation (captured by a preferred speed u₀ and a “stiffness” ξ that brings a particle towards the preferred speed) and (ii) a sliding friction between adjacent layers (captured by [small eta, Greek, tilde]), generated by particle–particle interactions, which favour rigid-body rotation of the cluster as a whole. In addition, the particles at the boundary have velocity u_R ≠ u₀.

The results from the model are shown as the solid lines in Fig. 2c. In the passive case E < E_Q, by construction the angular velocity ω = 0 (grey data in Fig. 2c). At the onset of activity (E = E_Q, yellow data), although there is a slight upturn in the model for the outer layer (n = 4), overall there is good quantitative agreement between the model and the experiment. Conversely at E = 1.4E_Q (red data), the model overestimates the angular velocity at small radii (n = 1, 2).

At higher field strength, the model predicts a flow inversion, where the inside of the rotor has a higher angular velocity than the exterior. This is apparent for E = 1.8E_Q (purple data), where we see good agreement between the model and the experiment. Further increase of the field to 2.0E_Q (blue data) yields a slightly worse prediction of the experiment. At the highest field strength (2.3E_Q, teal data), the model significantly overestimates the angular velocity for layers n = 2 to 4.

The key insight (see Section 3) is the fact that the threshold E field for self-propulsion is different at the boundary and the bulk and above the threshold the dependence on the self-propulsion speed on applied field is nonlinear (a power law).² The boundary begins to move at lower applied field than the interior behaviour but once the interior begins to move, at higher applied field, its velocity increases faster than the boundary with increasing applied field, which can lead to the flow inversion observed. Overall, our experimental discovery of flow field inversion is reproduced in a model with only two ingredients: a preferred active particle velocity, and friction between adjacent particle layers.

4.4. Stick-slip behaviour

The steady-state description that we have discussed so far assumes velocities that are constant in time. In reality, particle motion is intermittent, characterised by periods without relative motion between adjacent layers (“sticking”) interspersed with slipping events between adjacent layers. Slipping is associated with the creation and annihilation of localised defects in the hexagonal order. To identify these defects, we introduce the local bond-angle distortion parameter,computed at time t, for every simplex obtained from a Delaunay triangulation, where . Here, Θ = 0 represents equal bond angles in a perfectly ordered structure, whereas Θ = 1 indicates a severe distortion to . In Fig. 2d, we show the time-averaged radial profile of angular distortion Θ for a range of field strengths. As anticipated, Θ is largest under conditions that generate significant slipping between layers, both in the case that the outer layers move faster than the inner layers (E ≈ 1.4E_Q) and after the flow inversion when the inner layers move faster than the outer layers (E ≈ 2.3E_Q). At low field (E ≈ E_Q), Θ is low as the system exhibits rigid-body rotation, and it is again low prior to the flow inversion for (E ≈ 1.8E_Q).

The non-monotonic behaviour observed in Θ is reflected in the layer-resolved bond orientational order parameter ψ₆, shown in Fig. 2e. At low field strength, the system rotates rigidly, preserving hexagonal ordering and high ψ₆ in all layers. At the onset of slipping (E ≈ 1.4E_Q), hexagonal order is disrupted and ψ₆ is reduced. Hexagonality recovers at the higher field strengths , as the system quickly relaxes to a locally hexagonal structure following each slip event (see Movie S3). At the highest field strength (E ≈ 2.3E_Q), there is a significant drop in ψ₆, indicating fluidisation.

We investigate slipping in more detail in Fig. 3. To do so, we plot spatially resolved maps of Θ and ψ₆ as functions of time for a field strength of 1.4E_Q in Fig. 3a and b (see also Movie S3). During a slipping event, the two outer layers exhibit deviations from local hexagonal order at the vertices of the hexagon. This is reminiscent of slipping in driven assemblies of passive discs.⁴⁰ We see that locally high values of the distortion parameter Θ are correlated with low values of ψ₆. These deviations from perfect hexagonal order are confined to the outer two layers. The solid-body nature of the interior protects the hexagonal order from any bond distortion. We illustrate the intermittent dynamics of ψ₆ showing sharp transitions between high and low values that are highly correlated between the two layers in Fig. 3c. (We explore these correlations further in Fig. 5 below.)

Fig. 3 Slipping layers and order evolution. (a) and (b) Time sequence for a single layer-slipping event, showing the change of the distortion parameter Θ (a) and the hexagonal order parameter ψ₆ (b). (c) Evolution of ψ₆ in layers n = 2, 3, and 4. Incommensurate and commensurate configurations are indicated by the loss and recovery of ψ₆. Horizontal dashed lines indicate ψ₆ = 0.5. Shaded regions in c correspond to the slipping in sequence in (a) and (b).

4.5. Slipping mechanism: Frenkel–Kontorova model in periodic geometry

Slip in crystalline materials, in the form of layers sliding past each other, has long been understood in the context of the Frenkel–Kontorova model.¹⁵ This model describes frictional dynamics in the transition between commensurate and incommensurate states when thermal fluctuations are neglected.^37,42 While our system is of course thermal, insight can nevertheless be gained through the Frenkel–Kontorova model. Fig. 4a–c illustrates the one-dimensional Frenkel–Kontorova model, in which a harmonic chain of particles with interactions modelled by springs of rest length a is situated on a periodic potential corresponding to a substrate with period b as shown in Fig. 4. For a sliding chain in one dimension, the mechanism by which two lattices slip past each other involves the propagation of topological solitons known as kinks or anti-kinks. A kink corresponds to a localised compression of the chain that promotes a substrate spacing shared by two particles, as illustrated in Fig. 4b. Every kink propagation corresponds to the displacement of the local chain compression. On the other hand, anti-kinks are local extensions of the particle chain (Fig. 4c), which significantly reduce the friction between the particles.

Fig. 4 Frenkel–Kontorova behaviour in Quincke rotors. (a) Schematic representation of the Frenkel–Kontorova model. A slipping layer is modelled by a chain, comprising particles connected by springs of length a. The lattice periodicity is b, and the interaction between the static lattice and the sliding layer is given by a commensurate sinusoidal potential. (b) Kink and (c) anti-kink representation. A kink forms when two particles are contained in the same well, producing a local compression of the chain. The opposite case of an anti-kink is a local expansion of the chain, which leads to a vacancy. For our confined geometry, the slipping layers are treated as chains forming (d) commensurate and (e) incommensurate configurations. (f) The evolution of the rescaled interaction potential βU, (g) chain extension δa follows the oscillations of ψ₆ in Fig. 3c (see Fig. 5). In (d)–(f), the blue lines are the outermost layer (n = 4) and the dashed crimson lines are the adjacent internal layer, n = 3. Shaded regions in panels (f)–(h) correspond to sequences shown in Fig. 3a and b. Here the field strength is 1.4E_Q.

Fig. 5 Correlations between structural measures and energy. (a) and (b) Correlation between bond-orientational order ψ₆ and rescaled interaction potential βU/U_max in layers 4 (a) and 3 (b). (c) and (d) Correlation between bond-orientational order ψ₆ and stretching a in layers 4 (c) and 3 (d). Here the field strength is 1.4E_Q.

Unlike the classical Frenkel–Kontorova model with an infinite chain, our confined system consists of concentric layers as shown in Fig. 4d and e. These illustrations show commensurate and incommensurate configurations during stick-slip dynamics and demonstrate that our system allows only for anti-kinks due to local extensions of the chain. The hard-core interactions between particles prohibit compression and the formation of Frenkel–Kontorova kinks in these close-packed crystallites. However, in our system, anti-kinks correspond to stretching the system against the cohesive force between the particles, i.e., the electrohydroynamic interaction which binds the crystallite together.

To interpret the slipping in our system in terms of the Frenkel–Kontorova model, we estimate the energetics of the process. While the interactions between Quincke rollers are highly complex, involving hydrodynamic coupling and electrodynamic contributions,^2,43 we have found that in fact the structure and dynamical behavior of our system is quite well described by treating the electrohydrodynamic interactions with a long-ranged Yukawa-like attraction.³⁴ As a first approximation, here we follow the same approach and model the direct interactions between the particles with an attractive Yukawa potential. We stress that this is a major simplification and that the parameters we use should be taken as approximate. Here we use a well depth −βu_yuk(σ) = 10, where β = (k_BT)⁻¹ is the inverse of the thermal energy and a screening length equal to the particle diameter σ, which is somewhat shorter than that used previously,³⁴ taken to represent stronger screening of electrohydrodynamic interactions at the high area fractions we consider here. We use this to determine the potential energy of the assembly, which we plot as a function of time in Fig. 4f. The time-evolution is highly correlated with the structural order parameter ψ₆ shown in Fig. 3c, which we explore further below.

To measure the commensurability, we estimate the “chain extension”, or stretching as the nearest-neighbour separation in each slipping layer, with reference to the layer towards the centre. That is, the layer towards the centre is treated as a “substrate” in the context of the Frenkel–Kontorova model. We term the stretching δa, such that the distance between two stretched particles is δa + a with a the separation of two which are not stretched as shown in Fig. 4a and, for our geometry, in Fig. 4d and e. Note that, in a perfect hexagon, all particles touch their neighbours, so a = σ. In Fig. 4g, we see that the stretching δa is anti-correlated with the potential energy. To consider the formation of anti-kinks, i.e., chain extensions, we determine the quantity which is of order unity at the onset of slipping in the Frenkel–Kontorova model.¹⁵ To evaluate this quantity, we estimate the difference in potential energy between the slipping and non-slipping states ΔβU and the stiffness κ for our system.

We find ΔβU by subtracting the potential energy per particle, βU, of the minima (corresponding to the commensurate states) from that of the maxima, corresponding to the incommensurate states, see Fig. 4d–f. We determine κ by fitting a parabola to βU(r) for a bulk system. (see Fig. A1 in the Appendix, SI).

We evaluate for layers n = 2 and 3 as 17 (cf., measured value from Fig. 4g is δa = 0.094) and between the outer layers n = 3 and 4 as 18 (cf., measured value from Fig. 4g is δa = 0.079). While in the classic Frenkel–Kontorova model, slipping is expected when here we find a rather larger value. It is quite possible that we have over-estimated ΔβU or indeed that we have underestimated κ. While of course we can use other parameters which may provide an improved agreement, we note that an accurate interaction potential has yet to be determined for this system, although work in this direction is in progress. We believe that this topic could well be investigated in more detail in the future with a more accurate analysis than we have been able to perform here. We further observe that whether or not it is even possible to apply the FK model to active matter would be worth investigating, perhaps using computer simulation with precisely-controlled interactions between the particles.

Relating to the slipping phenomena that we observe is coupling between structural measures (ψ₆) and the stretching parameter δa and the energy βU. In Fig. 5a and b, we investigate correlations between the energy and the structure. These appear to be quite well correlated, although βU/βU_max is a nonlinear function of ψ₆, having a strong response at high ψ₆. In Fig. 5(c and d), we also consider the relationship between ψ₆ and δa. Here the scatter is larger, and the relationship between ψ₆ and δa is more linear.

5. Discussion and conclusion

We have investigated the behaviour of a model system of Quincke rollers confined in circular regions, which fix the population of each realisation. At low particle populations, this active system forms an apolar gas with a transition to a polar fluid at a higher density. Further increase of particle population leads to a hexagonally ordered packing which rotates.

For the magic number N_max = 61, corresponding to a perfect hexagonal packing, these particles exhibit a complex sequence of steady states as the activity is increased with a flow inversion and fluidisation at high activity. At the level of individual layers, a phenomenological model based on self-propulsion and friction captures this behaviour. As activity is increased, the quiescent state gives way first to a state with rigid rotation, and then subsequently to a state with sliding layers. At lower fields, the outermost layer moves fastest but at higher fields the inner regions rotate faster than the boundary, i.e. flow inversion. The model assumes a lower threshold to Quincke rotation (ie the onset of activity) at the boundary than the interior of the rotor. So at low activity the boundary moves faster. However, as the field is increased, the difference in linear active velocity between the interior and the boundary drops. The angular velocity is then reduced for the boundary relative to the interior and flow inversion results.

At the particle-resolved level, our experiments reveal intermittent stick-slip dynamics between the layers. We use measurements of the interactions between the particles³⁴ to compare with predictions for continuous slipping from the Frenkel–Kontorova model. The FK model predicts that when the square root of the ratio of the potential energy related to the displacement of the slipping layer from its optimum position to the stretching energy is around unity (i.e. ), however here we find rather higher values. Application of the Frenkel–Kontorova model to active matter would clearly be worthy of further investigation.

A related passive system of a hard-sphere colloidal “corral” being driven from the boundary, like a tiny rheometer, has been investigated previously.^40,41,44 Both the passive and active systems exhibit a fluidisation transition from a hexagonal configuration to a layered fluid. However, unlike the active system considered here, the passive rotors show neither flow-field inversion nor intermittent slipping. This highlights how the new phenomenology characteristic of active solids results in our current experiments from a combination of interparticle cohesion and activity.

Activity allows for solids in which complex dynamics emerge from simple ingredients. In our Quincke roller system, the formation of coherent active rotation due to confinement suggests a route towards the extraction of useful work from synthetic active matter, for example with the use of a gearwheel and axle as investigated previously.⁴⁰ We have demonstrated how cohesion and activity can be used to design soft nanomachines in which flow fields and intermittent dynamics can be manipulated at the microscopic scale.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article, including representative microscopy images compiled into movies are available at Zenodo at https://doi.org/10.5281/zenodo.14781732.

Description of Supplementary Movies: Movie 1 – Experimental set-up. Movie shows multiple wells of diameter R = 30 μm containing different populations of Quincke rollers. Electric strength is E = 1.4E_Q. Scale bar represents 100 μm. Movie is reproduced at 30 fps. Movie 2 – Rigid body rotation behaviour at E = E_Q. Movie shows a roller population of N = 61 showing rigid body rotation. (a) Experimental acquisition, and (b) shows guides to the eye for individual rollers located at different layers and over symmetry lines. Scale bar in (a) represents 10 μm. Movie is reproduced at 6 fps. Movie 3 – Slipping behaviour at E = 1.4E_Q. Same roller population showing steady-state slipping of the outermost layers. (a) Experimental acquisition. (b) Guides to the eye for individual rollers in different layers. (c) Shows local hexagonal order ψ₆ and (d) is the local distortion Θ. Scale bar in (a) represents 10 μm. Movie is reproduced at 6 fps. Movie 4 – Onset of the flow inversion at E = 2.0E_Q. Movie shows the overtaking of internal layers past over the outermost layers. (a) Experimental acquisition. (b) Guides to the eye for individual rollers. (c) Local hexagonal order ψ₆. Scale bar in (a) represents 10 μm. Movie is reproduced at 6 fps. Movie 5 – Fluid behaviour at E = 2.3E_Q. Movie shows the complete fluidisation of the structure, with every layer rotating with an independent angular velocity ω. The flow inversion becomes more evident at this value of the field strength, where the internal layers rotate faster. (a) Experimental acquisition. (b) Guides to the eye for individual rollers at different layers. (c) Local hexagonal order ψ₆. Scale bar in (a) represents 10 μm. Movie is reproduced at 6 fps.

Acknowledgements

The authors would like to thank Leticia Cugliandolo, Olivier Dauchot, Petia Vlahovska and Srikanth Sastry for helpful discussions. CPR would like to acknowledge finding from ANR, DiViNeW. AS, CPR and TBL would like to thank the Kavli Center for Theoretical Physics for hospitality and support during the Active Solids: From Metamaterials to Biological Tissue program. This work was supported by EPSRC grants EP/R014604/1 and EP/T031077/1.

Notes and references

1. M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. Aditi Simha, Rev. Mod. Phys., 2013, 85, 1143.
2. A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot and D. Bartolo, Nature, 2013, 503, 95–98.
3. A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini and M. Viale, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 11865–11870.
4. M. J. Bowick, N. Fakhri, M. C. Marchetti and S. Ramaswamy, Phys. Rev. X, 2022, 12, 010501.
5. N. A. M. Araújo, L. M. C. Janssen, T. Barois, G. Boffetta, I. Cohen, A. Corbetta, O. Dauchot, M. Dijkstra, W. M. Durham, A. Dussutour, S. Garnier, H. Gelderblom, R. Golestanian, L. Isa, G. H. Koenderink, H. Löwen, R. Metzler, M. Polin, C. P. Royall, A. Šari, A. Sengupta, C. Sykes, V. Trianni, I. Tuval, N. Vogel, J. M. Yeomans, I. Zuriguel, A. Marin and G. Volpe, Soft Matter, 2023, 19, 1695–1704.
6. Y. Fily, A. Baskaran and M. C. Marchetti, Soft Matter, 2012, 8, 3002–3009.
7. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
8. L. Braverman, C. Scheibner, B. van Saders and V. Vitelli, Phys. Rev. Lett., 2021, 127, 268001.
9. P. Baconnier, D. Sohbat, C. Hernández López, C. Coulais, V. Démery, G. During and O. Dauchot, Nat. Phys., 2022, 18, 1234–1239.
10. A. P. Petroff, X.-L. Wu and A. Libchaber, Phys. Rev. Lett., 2015, 114, 158102.
11. T. H. Tan, A. Mietke, J. Li, Y. Chen, H. Higinbotham, P. J. Foster, S. Gokhale, J. Dunkel and N. Fakhri, Nature, 2022, 607, 287–293.
12. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Abd Volpe, Rev. Mod. Phys., 2016, 88, 045006.
13. K. J. M. Bishop, S. B. Biswal and B. Bharti, Annu. Rev. Chem. Biomol. Eng., 2023, 14, 1–30.
14. G. Briand, M. Schindler and O. Dauchot, Phys. Rev. Lett., 2018, 120, 208001.
15. P. M. Chaikin and T. C. LuLubensky, Principles of Condensed Matter Physics, Cambridge University Press, 1995.
16. C. A. Weber, C. Bock and E. Frey, Phys. Rev. Lett., 2014, 112, 168301.
17. X.-Q. Shi, F. Cheng and H. Chaté, Phys. Rev. Lett., 2023, 131, 10830.
18. J. U. Klamser, S. C. Kapfer and W. Krauth, Nat. Commun., 2018, 9, 5045.
19. L. Caprini, U. Marini Bettolo Marconi and H. Löwen, Phys Rev. E, 2023, 108, 044603.
20. C. B. Caporusso, L. F. Cugliandol, P. Digregorio, G. Gonnella, D. Levis and A. Suma, Phys. Rev. Lett., 2023, 131, 068201.
21. N. Klongvessa, F. Ginot, C. Ybert, C. Cottin-Bizonne and M. Leocmach, Phys. Rev. Lett., 2019, 123, 248004.
22. D. Geyer, D. Martin, J. Tailleur and D. Bartolo, Phys. Rev. X, 2019, 9, 031043.
23. N. Sakaï, K. Skipper, F. J. Moore, J. Russo and C. P. Royall, Soft Matter, 2025, 21, 5204–5213.
24. F. Meng, D. Matsunaga, R. Golestanian and P. Tierno, Nat. Commun., 2019, 10, 2444.
25. H. Massana-Cid, C. Maggi, N. Gnan, G. Frangipane and R. Di Leonardo, Nat. Commun., 2024, 15, 6574.
26. P. Digregorio, D. Levis, L. F. Cugliandolo, G. Gonnella and M. I. Pagonabarraga, Soft Matter, 2022, 18, 566.
27. P. Digregorio, D. Levis, L. F. Cugliandolo, G. Gonnella and I. Pagonabarra, Soft Matter, 2022, 18, 566–591.
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., 2022, 18, 212–218.
29. A. Mauleon-Amieva, M. P. Allen, T. B. Liverpool and C. P. Royall, Sci. Adv., 2023, 9, eadf5144.
30. 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.
31. A. G. Bayram, F. J. Schwarzendahl, H. Löwen and L. Biancofiore, Soft Matter, 2023, 19, 4571–4578.
32. H. Karani, G. E. Pradillo and P. M. Vlahovska, Phys. Rev. Lett., 2019, 123, 208002.
33. R. Reyes Garza, N. Kyriakopoulos, Z. M. Cenev, C. Rigoni and J. V. I. Timonen, Sci. Adv., 2023, 9, eadh2522.
34. 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.
35. J. Toner and Y. Tu, Phys. Rev. Lett., 1995, 75, 4326–4329.
36. Physics of Sliding Friction, ed. B. N. J. Persson and E. Tosatti, Springer, 1996, vol. 311.
37. T. Bohlein, J. Mikhael and C. Bechinger, Nat. Mater., 2012, 11, 126–130.
38. W. D. Ristenpart, P. Jiang, M. A. Slowik, C. Punckt, D. A. Saville and I. A. Aksay, Langmuir, 2008, 24, 12172–12180.
39. J. Crocker and D. Grier, J. Colloid Interface Sci., 1996, 179, 298–310.
40. I. Williams, E. C. Oğuz, T. Speck, P. Bartlett, H. Löwen and C. P. Royall, Nat. Phys., 2016, 12, 98–103.
41. I. Williams, E. C. Oğuz, P. Bartlett, H. Löwen and C. P. Royall, Nat. Commun., 2013, 4, 1–6.
42. A. Ward, F. Hilitski, W. Schwenger, D. Welch, A. W. W. Lau, V. Vitelli, L. Mahadevan and Z. Dogic, Nat. Mater., 2015, 14, 583–588.
43. S. Imamura, K. Sawaki, J. J. Molina, M. S. Turner and R. Yamamoto, Adv. Theory Simul., 2023, 6, 2200683.
44. I. Williams, E. C. Oğuz, H. Löwen, W. C. K. Poon and C. P. Royall, J. Chem. Phys., 2022, 156, 184902.

---