Introduction

Active matter has gained much attention as a field rich with unique and potentially useful non-equilibrium phenomena1. These systems consist of self-driven units that have relatively simple individual dynamics but collectively exhibit macroscopic coherent motions2. While flocks of birds and schools of fish are the commonly given examples of such emergent collective behavior, similar phenomena have been observed on much smaller scales with active microagents such as bacteria3,4,5,6,7,8,9,10, chemically activated motile colloids11,12,13,14,15, microtubule–kinesin bundles (active nematics)16, living liquid crystals (suspensions of motile bacteria in liquid crystals)17,18, and field-driven colloids19,20,21,22,23,24,25,26,27.

Most studies of active suspensions have focused on translational particles such as bacteria and their inanimate mimics28,29,30,31,32,33,34,35. However, there has been rising interest in systems containing self-rotating particles36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51. An isolated micro-rotor has no means of self-motility. But when a rotor is coupled with other rotors, motility becomes possible. Theoretical studies showed that coupling mediated by contact in dry systems36 or hydrodynamic interactions37,52 can induce propulsion and self-organization. While biological examples of rotors are limited, e.g., the bacterium Thiovulum majus53, there has been a growing number of synthetic designs of rotors fueled by chemical reaction54 or external forcing such as magnetic fields39,44,46,51,55, electric fields56,57,58, oscillating platform59, ultrasound60, and applied air flow61,62.

Simulations based on agent-based models, without hydrodynamics, for gear-shaped spinners36,63 and rotating dimers64 showed that a mixture of clockwise (CW) and counterclockwise (CCW) rotors segregate into same-spin phases. This behavior was actualized in an experiment65 using a large number of small robots spinning on a table. The collective dynamics of fluid-suspended rotors, however, can be quite different from that observed in the dry systems due to hydrodynamic interactions. For example, unlike the frozen state of gear-shaped spinners at low density36, fluid-embedded spinners exhibit a gas-like phase, where the particles move randomly in the stirred fluid38. Furthermore, while both the dry and the fluid systems form similar large-scale dynamic patterns, such as lanes and vortices, their quantitative characteristics are different. Much insight has been gained by studying these systems using discrete-spinner models, however, this approach becomes computationally expensive when considering large systems with a high spinner density. In such cases, it becomes favorable to utilize continuum models that replace the discrete individual particles with a continuous distribution. Such continuum approach allows for the computationally feasible study of the collective dynamics of large systems27,66,67,68,69.

The existing continuum models of self-rotating systems are phenomenological and typically based on the simplified hydrodynamic theory that includes either chiral components in the stress tensor or additional chiral forcing terms in the Navier–Stokes equation45,46,70. These models capture phenomena such as the surface waves in active spinner systems46, active coarsening, and the emergence of self-propelled vortex doublets70, etc. However, the phenomenological implementation of active rotation makes it challenging to assess if the observed dynamics are physical or dependent on the postulated continuum formulation. It also becomes difficult to directly compare experiment and theory.

Here we derive a continuum model for a monolayer of spherical spinners with CW- and CCW spins suspended in liquid in a three-dimensional domain. The model is motivated by experimental realizations such as spinners trapped at a fluid interface39,71 or on a solid substrate46. We present numerical simulations assuming rotation due to an external DC (direct current) electric field, known as Quincke rotation. We explore the mixtures of CW and CCW spinners and study same-spin phase separation. We find that the dynamics is sensitive to fluid inertia: while in the inertialess system after a long transient turbulent-like motion the spinners settle and form steady traffic-like lanes, at small but finite Reynolds number the turbulent-like motion persists. For even larger Reynolds numbers, a population initially composed of an equal amount of CW/CCW spinners exhibits a chirality-breaking transition and collapses into a single rotation sense state. While our results are obtained for Quincke spinners, a similar continuum approach can be applied to model a monolayer of magnetic spinners energized by an applied AC magnetic field39,71. Our results advance a fundamental understanding of active systems, where large-scale chiral states emerge as a result of spontaneous symmetry breaking, as for example spinners, rotating bacteria, microtubules assays, etc.

Results and discussion

Summary of the continuum model

In the “Methods” section, we derive a continuum model for a two-dimensional system of N active spinners of mass m and radius a suspended in a fluid with viscosity μ and density ρ confined to a plane of area A, see Fig. 1a. The spinner concentration, c = N/A, is assumed to be uniform. The system is governed by the following equations for the spinner angular velocity ω(r, t), fluid velocity u(r, t), and pressure p(r, t):

$${I}_{p}\left(\frac{\partial \omega }{\partial t}+{\bf{u}}\cdot \nabla \omega \right)=D{I}_{p}{\nabla }^{2}\omega +8\pi \mu {a}^{3}\left({(\nabla \times {\bf{u}})}_{z}-2\omega \right)+2\tau $$
(1)
$$\rho \left(\frac{\partial {\bf{u}}}{\partial t}+({\bf{u}}\cdot \nabla ){\bf{u}}\right)=-\nabla p+\mu {\nabla }^{2}{\bf{u}}+4\pi \mu {a}^{2}c\nabla \times (\omega \hat{{\bf{z}}})$$
(2)
$$\nabla \cdot {\bf{u}}=0$$
(3)

where \({I}_{p}=\frac{2}{5}m{a}^{2}\) and \(D=\frac{{k}_{B}T}{6\pi \mu a}\). The thermal energy of the system is given by kBT. Details of the model derivation are provided in the “Methods” section.

Fig. 1: Schematics of the discrete system.
figure 1

a Illustration of spinners confined to a plane. Blue and yellow color corresponds to clockwise (CW) and counterclockwise (CCW) rotation. The axis of rotation for each spinner is shown as gray vertical line oriented in the z-direction. b The Quincke effect is an electrohydrodynamic instability that gives rise to a torque on a dielectric particle in a uniform DC (direct current) electric field94 if the free-charge polarization is antiparallel to the applied field. Above a critical field strength E > EQ rotation about an axis in the plane perpendicular to the electric field (Ω ⋅ E = 0, where Ω is the individual spinner angular velocity) is generated by the misaligned induced dipole of the particle, P. The direction of rotation can be either CW or CCW, depending on the initial perturbation in the dipole direction.

Full size image

Equation (1) is the spin angular momentum balance of the spinners; the three right-hand side terms account for the thermal noise, the torque on the spinners due to fluid drag, and the externally applied torque τ. Equation (2) is the linear momentum conservation for the fluid (the Navier–Stokes equations) with the addition of the last term on the right-hand side, which drives the fluid flow due to the rotation of the spinners. Equation (3) is the incompressibility condition since we assume constant fluid density ρ.

It is worth noting that Eqs. (1)–(3) are generic for any spherical spinners. The specifics of the particular physical system are introduced in the form of the externally applied torque τ, which is what causes the particles to spin. In this paper we consider the case of Quincke spinners, studied experimentally in ref. 72, τ = c−1(P×E)z, where P is the polarization vector of the individual spinner and E is the externally applied constant (DC) electric field, see Fig. 1b. The governing equations for the polarization vector P are provided in the “Methods” section, see Eqs. (28)–(29).

The complete model is provided in non-dimensional form in the Methods section, see Eqs. (31)–(36). All equations are nondimensionalized using the spinner radius a as the characteristic length scale and the Maxwell–Wagner polarization time tQ, which sets the magnitude of the rotation rate, as the characteristic time scale. The scaling of the physical variables and the definitions of the model parameters are summarized in Table 1. The scaling analysis shows that the most important parameter governing the behavior of the system is the Reynolds number for the fluid, \({\rm{Re}}=\rho {a}^{2}/\mu {t}_{Q}\). For \({\rm{Re}}\ll 1\), the effect of inertial terms in Eq. (2) of is small; fluid inertia becomes important for \({\rm{Re}}\gg 1\). In the rest of the paper, we explore the system behavior as Re increases. Unless otherwise noted, all simulations are performed with α = 0.25, γ = 1.1, κ = 0.1, D* = 1, \({D}_{P}^{* }=0.1\), box size L = 480a, N = 1024, and time step dt = 0.00025tQ.

Table 1 The dimensionless parameters and their definitions in terms of the properties of the spinners and the environment.
Full size table

Spinners stir large-scale fluid flow

Figure 2a–c and Supplementary Movie 1 show that the spinner fluid undergoes a phase separation into clusters of CCW (yellow) and CW spinners (blue), in agreement with the simulations using a discrete spinner model38. These clusters themselves are rotating, growing, separating, and reconnecting as also observed in the active coarsening continuum model in ref. 70.

Fig. 2: Predicted collective behavior at intermediate Reynolds number \({\rm{Re}}=0.1\).
figure 2

a Spinner angular velocity, ω, showing clusters of clockwise spinners (ω < 0) and counter-clockwise spinners (ω > 0). b Squared magnitude of the fluid velocity, ∣u∣2. c Plot of the fluid vorticity (∇ × u)z. Black arrows show the direction of fluid flow. d Sketch illustrating that in the interior of a region of co-rotating spinners flow is negligible due to the canceling effect of the same sense of rotation of neighboring spinners. At the interface between regions of opposite rotation there is a net flow (black arrows) generated by the counter-rotating spinners. Yellow spheres rotate counterclockwise and blue spheres rotate clockwise, gray arrow indicates the direction of rotation.

Full size image

The flow is localized the interface between CW and CCW clusters, see Fig. 2b and Supplementary Movie 2. Physically, the fluid flow is canceled out in the region between the same spin spinners. Hence, even though the spinners might rotate very fast, the flow might be very weak, especially in the interior of the same spin regions; accordingly, while inertia might be very important in the angular conservation equation, it might be less so in the linear momentum conservation (Navier–Stokes equation for the fluid flow). At the boundary between clusters, opposite spin spinners pump fluid in the same direction causing bands of flow or edge currents, see an illustration of this in Fig. 2d.

Emergence of lanes depends on Reynolds number

An important question concerns the long-term behavior of a system composed of spinners that are initially randomly 50 CW : 50 CCW. Previous theoretical studies38,41 of discrete micro-spinners have shown that the system exhibits separation into same-spin clusters that eventually evolve into emergent patterns, such as lanes or vortices. This work however was restricted to Stokes flow (\({\rm{Re}}=0\)), while recent colloidal experiments39,44 suggest a non-negligible Reynolds number.

Using the continuum model, if we consider \({\rm{Re}}=0\) and an initial 50 CW : 50 CCW random configuration, lanes always eventually form, see Fig. 3a and Supplementary Movie 3, regardless of the initial configuration. Therefore, the continuum model is in agreement with the discrete spinner model. The formation of lanes can be qualitatively interpreted as a result of the fact that the Stokes flow corresponds to least dissipation. Segregation into the same spin regions is thus favored since the flow vanishes in the interior of the same spin region and the only remaining flow is confined to the interface between opposite spin regions. The flow (and thus the dissipation) is reduced if the interface between these same spin domains is minimized, which in the case of the periodic boundary conditions considered in our simulations, is a straight line. Conservation of mass further requires two interfaces with opposite flow directions resulting in a lane configuration.

Fig. 3: Predicted spinner angular velocity at increasing Reynolds numbers \({\rm{Re}}\).
figure 3

Plots of spinner angular velocity ω for the long-term behavior for a \({\rm{Re}}=0\), b \({\rm{Re}}=0.01\), and c \({\rm{Re}}=0.1\). Yellow regions have counterclockwise (CCW) spinners and blue regions have clockwise (CW) spinners. The black arrows show the direction of fluid flow.

Full size image

A small, non-zero Reynolds number, such as \({\rm{Re}}=0.01\), leads to distorted lanes, see Fig. 3b and Supplementary Movie 4. These quasi-stable lanes have turbulent-like flow near the interfaces between the two opposite-spin clusters. The lanes can also spontaneously dissolve into to a transient turbulent-like state and then reappear.

For small but finite Reynolds number, e.g., \({\rm{Re}}=0.1\) and greater, the system only exhibits a turbulent-like motion with no emergent collective large-scale structure, as shown in Fig. 3c and Supplementary Movie 1. If the initial state is lanes, instead of a random configuration, the lanes become unstable and dissolve into turbulent-like motion, see Supplementary Movie 5.

For an even larger Reynolds number, such as \({\rm{Re}}=1\), the ratio of CW to CCW spinners is no longer conserved and the system exhibits a chirality-breaking transition. As the inertial terms of the Navier–Stokes equation become more significant, the fluid vorticity is not solely determined by the spinner angular velocity distribution. Instead, the inertial-induced component of the fluid vorticity can cause the spinners to flip spin orientation. For this case, the fluid flow initially exhibits turbulent-like motion, but eventually begins to favor a single spin orientation, see Supplementary Movie 6. Clearly this behavior is not present in dry system of spinners.

We can quantitatively distinguish between the lane and the turbulent-like cases by studying the probability distribution function (PDF) of the fluid velocity, u. As shown in Fig. 4a, the PDF for lanes exhibits a clear bi-modal distribution due to the two sides of the lane. For the turbulent-like case, Fig. 4b, the PDF is Gaussian unlike the non-Gaussian behavior observed in experiments with colloidal spinners71,73. Overly populated tails in our simulation were found to be due to underresolved simulations. Velocity and vorticity correlation functions for turbulent-like regimes are displayed in Fig. 4c, d.

Fig. 4: Probability distribution functions (PDF) for the fluid velocity and correlation functions for the fluid velocity and vorticity.
figure 4

Reynolds numbers are: a \({\rm{Re}}=0\) and b \({\rm{Re}}=0.1\). Spatial correlation of the c fluid velocity 〈u(r) ⋅ u(r + r0)〉 and d fluid vorticity \(\langle \left(\nabla \times {\bf{u}}({\bf{r}})\right)\cdot \left(\nabla \times {\bf{u}}({\bf{r}}+{{\bf{r}}}_{0})\right)\rangle \). Here r = (x, y) is the spatial position, and r0 is the displacement vector. The vertical axis for (c) and (d) are normalized so the curves have intercepts equal to one. The gray curve shows the lane formation case, \({\rm{Re}}=0\), and the blue curve shows the turbulent-like case, \({\rm{Re}}=0.1\). In the case of lanes, the PDF (a) exhibits a clear bi-modal distribution as expected with formation of lanes, while for sustained turbulent-like motion, the PDF (b) exhibits a Gaussian distribution. For reference, the dashed curve provides the Gaussian distribution with the same mean and standard deviation.

Full size image

Finally, it is worth mentioning that the discrete spinner model38,41 also observed the formation of vortices in the case of larger spinner concentrations, which is not captured by our continuum model. This deviation is likely because steric repulsion is neglected in the continuum model but becomes non-negligible at higher spinner concentrations. Steric repulsion can be included in the continuum model by introducing a particle stress tensor27.

Turbulent-like flow has a power-law energy spectrum

In many quasi-two-dimensional active systems exhibiting turbulent-like behavior, such as magnetic rollers39, or swimming bacteria74, the energy is injected at the microscopic scale. For example, rotation of bacterial flagella results in the onset of large-scale collective behavior and so-called "active turbulence”6, characterized by the inverse cascade75. It is tempting to assume that in the Quincke system, the inverse cascade should be present as well due to the microscopic rotation of individual spinners. However, the situation is less clear in our model because the microscopic energy injection scale is coarse-grained in the continuum approximation and replaced by the effective driving torque τ in Eq. (1). This torque may slowly vary in space and does not produce the inverse cascade.

As the simulations show turbulent-like behavior, it is of interest to compute the energy spectrum, E(k), and the energy flux in k-space, ΠE. The energy spectrum E(k) is obtained by taking the Fourier transform of the spatial auto-correlation of u:

$$E(k)=\frac{k}{2\pi }{\left\langle {\int}_{A}{e}^{-i{\bf{k}}\cdot {{\bf{r}}}_{0}}\left\langle {\bf{u}}({\bf{r}})\cdot {\bf{u}}({\bf{r}}+{{\bf{r}}}_{0})\right\rangle d{x}_{0}d{y}_{0}\right\rangle }_{| {\bf{k}}| = k}$$
(4)

such that \({{\bf{r}}}_{0}={x}_{0}\hat{{\bf{x}}}+{y}_{0}\hat{{\bf{y}}}\) and 〈. . . 〉∣k∣=k is the average over all wavenumbers k such that ∣k∣ = k. Correspondingly, the energy flux is obtained as follows

$${{{\Pi }}}_{E}(k)=\langle {{\bf{u}}}^{\,{<}\,k}[(| {\bf{u}}\cdot \nabla ){\bf{u}}]\rangle $$
(5)

Here u<k is the low-pass filtered velocity field with the wave numbers outside ∣k∣ < k set to zero:

$${{\bf{u}}}^{\,{<}\,k}({\bf{r}},t)=\mathop{\sum}\limits_{| {\bf{k}}| <{\bf{k}}}\tilde{{\bf{u}}}({\bf{k}},t){e}^{i{\bf{k}}\cdot {\bf{r}}}$$
(6)

and \(\tilde{{\bf{u}}}({\bf{k}},t)\) is the Fourier transform of u(r, t)76.

A hallmark of the inverse cascade in 2D turbulence is the energy scaling E(k) ~ k−5/3 and negative value of the energy flux ΠE(k) for small wavenumbers k. As shown in Fig. 5, the energy cascade scales as k−5/3. However, the energy flux is negative only in a very narrow range, ΠE(k) < 0 for k < k0 ≈ 0.016, which is not consistent with the energy cascade observed in 2D classical turbulence76,77,78 and colloidal experiments44,71. We suspect that this behavior is due to the relatively small system size in the computations and coarse-gaining of the microscopic energy infection scale. In the cases where lanes form at long times, the fluid flow initially exhibits transient turbulence-like motion with the same k−5/3 scaling, prior to settling into lanes.

Fig. 5: Energy spectrum E(k).
figure 5

The spectrum is computed by Eq. (4) for a flow at the intermediate Reynolds number \({\rm{Re}}=0.1\). The spectrum shows an energy cascade and scaling, k−5/3. Inset: dependence of the energy flux ΠE(k) vs k. Computations performed with L = 300a to increase resolution, L is the integration domain size, a is the particle radius.

Full size image

The turbulence inertial range is characterized by two lengthscales, the Taylor microscale76, and the integral scale79. The Taylor microscale falls in between the large-scale eddies and the small-scale eddies. The integral scale LI can be extracted from the velocity autocorrelation function in Fig. 4c. It gives LI ≈ 100 for particle Reynolds number Re = 0.1. The corresponding r.m.s. (root mean square) velocity urms can be determined from Fig. 4b and is about 5, resulting in the integral scale Reynolds number of the order of 50. The Taylor microscale λ is several times smaller, and is determined from the fluctuations of velocity and velocity gradients, \({({\partial }_{{\bf{r}}}\langle {\bf{u}}\rangle )}^{2}={{\bf{u}}}_{{\rm{rms}}}^{2}/{\lambda }^{2}\). It roughly can be estimated from the velocity gradients (vorticity) autocorrelation function, Fig. 4d, which gives λ ≈ 15.

Initially, the spinners begin to phase separate into larger and larger same spin clusters. However, the energy spectrum exhibits a peak that occurs at the length scale of the domain, ~L. This demonstrates the cluster size becomes limited by the size of the domain. The spatial correlation functions for the fluid flow73 and fluid vorticity are another useful metric for studying turbulent flows. As shown in Fig. 4, the turbulent flow case exhibits no characteristic length scale as seen in meso-scale turbulence. This behavior is consistent with colloidal spinner experiments71,73.

Activity enhances transport

The capability of the spinner system for active transport can be quantified by tracking the position of point particles si that are being passively advected by the fluid flow, ∂si/∂t = u(si, t) and computing the mean-square-displacement (MSD), \(\left\langle | {{\bf{s}}}_{i}(t+{t}_{0})-{{\bf{s}}}_{i}({t}_{0}){| }^{2}\right\rangle \). In addition, the mean square velocity (MSV) for the fluid flow, \(\left\langle | {\bf{u}}{| }^{2}\right\rangle \), is another useful metric for the studying the active transport of the system.

As shown in Fig. 6a, the MSD initially exhibits a transient ballistic regime for all shown \({\rm{Re}}\) values. Counterintuitively, the MSD and MSV decrease as \({\rm{Re}}\) increases, Fig. 6b. For increasing \({\rm{Re}}\), the spinner-driven flow begins to compete with the effects from the inertial terms. Also, the force term in Eq. (2), accounting for the spinners disturbing the fluid, is proportional to viscosity. Therefore, the disturbance in the fluid generated by the spinners does not get damped down with higher fluid viscosity.

Fig. 6: Time evolution of the mean-square-displacement (MSD) and mean square velocity (MSV).
figure 6

a MSD and b MSV. Each curve representing a different value of Reynolds number \({\rm{Re}}\): \({\rm{Re}}=0\) (gray), \({\rm{Re}}=0.01\) (red), \({\rm{Re}}=0.1\) (blue), \({\rm{Re}}=1\) (yellow). Reference scalings are shown with black dashed lines.

Full size image

The curves for different \({\rm{Re}}\) values begin to qualitatively diverge in the long term. The lane formation cases, \({\rm{Re}}=0\) and \({\rm{Re}}=0.01\), show sustained ballistic motion in MSD (~t2) and a relatively constant MSV, as all the motion is confined to the interface of the lanes. Yet for sustained turbulent-like behavior, \({\rm{Re}}=0.1\) and \({\rm{Re}}=1\), the MSD shifts away from a ballistic regime closer towards a diffusive regime (~t). However, a reliable evaluation of the diffusion coefficient would require a much longer integration time. For Re = 1, the MSV also exhibits an eventual decrease. This is due to the chirality breaking of the spinner orientation.

Conclusions

To summarize, we develop a continuum coarse-grained model for a suspension of discrete microspinners. Assuming Quincke rotation, we then present numerical simulations that demonstrate the emergence of lanes of same-spin fluid and turbulent-like behavior depending on the Reynolds number. Our work sheds light on active spinner materials and makes testable predictions for experiments on the effects of fluid inertia and chirality symmetry breaking.

For the case of the Stokes flow, the formation of lanes of spinners is favored as it minimizes the interfaces between counter-spin clusters, thus minimizing the energy of the system. This behavior is in agreement with the existing discrete spinner models36,38. Yet, even small inertia causes turbulent-like behavior, which becomes more pronounced with increasing Re: turbulent lanes for \({\rm{Re}} \sim 0.01\) and sustained turbulent-like flow for \({\rm{Re}}\ge 0.1\). We also quantitatively characterize lane formation and the turbulent-like flow through the mean square displacement, MSV, and the energy spectrum statistics. Furthermore, for non-zero Reynolds numbers we observe that initially 50:50 CW/CCW population of spinners exhibit a chirality symmetry breaking transition and consequent condensation to a single sense of rotation state.

While this work used periodic boundaries, future work could apply this model with both fixed, impenetrable boundaries and soft, deformable boundaries, as well as arrays of artificial obstacles. Another possible avenue is the case of non-uniform spinner densities. As the spinner densities, in experiments, are often non-constant in space and time, it would be of interest to see how this would affect the collective and turbulent behavior seen in this study. The model can be modified to analyze other experimental systems, e.g., magnetic spinners at water-air interface25,39 and could also be easily generalized for a three-dimensional suspension of neutrally buoyant micro-rotors.

Methods

Derivation of the continuum model

Let us consider N spherical spinners of radius a and mass m confined to a 2D plane of area A. The translational motion of a spinner is the result of being passively advected by the surrounding fluid flow, u, and Brownian diffusion. The translational velocity of an individual spinner located at position r is then:

$${{\bf{v}}}_{p}({\bf{r}},t)={\bf{u}}({\bf{r}},t)+\sqrt{2D}{\boldsymbol{\xi }}(t).$$
(7)

The thermal noise of the system is modeled by a random variable ξ(t) that is Wiener process governed by the standard normal distribution. In the low-Reynolds number limit, the diffusion constant is D = kBT/6πμa, where kBT is the thermal energy and μ is the fluid viscosity. We also assume low spinner concentration such that the effect of steric repulsion is negligible.

We introduce a probability density, Ψ(r, Ω, t), of finding a spinner at position r and with angular velocity Ω. Since there are N total spinners, the normalization for Ψ is given as

$${\iint }_{A}{{\Psi }}({\bf{r}},{{\Omega }},t)d{{\Omega }}dA=N.$$
(8)

Using Eqs. (7)–(14) and the Fokker–Planck equation80, it can be shown that the governing equation for Ψ is given by the following Smoluchowski equation

$$\frac{\partial {{\Psi }}}{\partial t}+{\bf{u}}\cdot \nabla {{\Psi }}-D{\nabla }^{2}{{\Psi }}+\frac{\partial }{\partial {{\Omega }}}\left(\dot{{{\Omega }}}{{\Psi }}\right)=0,$$
(9)

assuming an incompressible fluid (∇ ⋅ u).

From this point on the derivation of our continuum model follows ref. 67. It is useful to define a function \(\left\langle \right..{\rangle }_{{{\Omega }}}\) such that

$$\left\langle \right.X{\rangle }_{{{\Omega }}}({\bf{r}},t)=\int_{-\infty }^{\infty }X({\bf{r}},{{\Omega }},t){{\Psi }}({\bf{r}},{{\Omega }},t)d{{\Omega }}.$$
(10)

We will therefore define the local spinner density as c(r, t) = 〈1〉Ω, the local angular velocity per unit area as \(\overline{\omega }({\bf{r}},t)={\langle {{\Omega }}\rangle }_{{{\Omega }}}\), and the local external torque per unit area as \(\overline{\tau }({\bf{r}},t)={\langle {\tau }_{e}\rangle }_{{{\Omega }}}\).

Integrating Eq. (9) leads to

$$\frac{\partial c}{\partial t}+{\bf{u}}\cdot \nabla c-D{\nabla }^{2}c+\left(\dot{{{\Omega }}}{{\Psi }}\right){| }_{-\infty }^{\infty }=0.$$
(11)

We can ignore the boundary term in Eq. (11) as spinners with ∣Ω∣ → ∞ are non-physical and inherently have zero probability, Ψ = 0.

$$\frac{\partial c}{\partial t}+{\bf{u}}\cdot \nabla c=D{\nabla }^{2}c.$$
(12)

Therefore Eq. (12) gives the evolution equation for the spinner density of the system. It is worth noting that, with appropriate boundary conditions, c is conserved and ∫AcdA = N according to Eq. (8).

Spin angular momentum equation

If we multiply Eq. (9) by Ω and integrate with respect to Ω, it is straightforward to show that

$$\frac{\partial \overline{\omega }}{\partial t}+{\bf{u}}\cdot \nabla \overline{\omega }=D{\nabla }^{2}\overline{\omega }+2{\langle\dot{{{\Omega }}}\rangle}_{{{\Omega }}}-\left({{\Omega }}\dot{{{\Omega }}}{{\Psi }}\right){| }_{-\infty }^{\infty }.$$
(13)

We can again ignore the boundary term as spinners with ∣Ω∣ → ∞ are non-physical and have zero probability.

The conservation of spin angular momentum of a single spinner is

$${I}_{p}\dot{{{\Omega }}}({\bf{r}},t)={\zeta }_{p}\left(\frac{1}{2}{(\nabla \times {\bf{u}})}_{z}({\bf{r}},t)-{{\Omega }}({\bf{r}},t)\right)+{\tau }_{e}({{\Omega }},t)$$
(14)

where \({I}_{p}=\frac{2}{5}m{a}^{2}\) is moment of inertia and ζp = 8πμa3 is the friction coefficient for a sphere and the externally applied torque on the spinner is τe(Ω, t). The effect of thermal noise is ignored in Eq. (14) because the spinners rotation rate is very high and thus the rotational Peclet number Ω/Dr ≫ 1, where Dr = kBT/(8πμa3).

Using Eq. (14), Eq. (13) simplifies to

$$\frac{\partial \overline{\omega }}{\partial t}+{\bf{u}}\cdot \nabla \overline{\omega }=D{\nabla }^{2}\overline{\omega }+\frac{{\zeta }_{p}}{{I}_{p}}\left(c{(\nabla \times {\bf{u}})}_{z}-2\overline{\omega }\right)+2{I}_{p}^{-1}\overline{\tau }.$$
(15)

Furthermore, if we assume a constant spinner density c and define \(\omega =\overline{\omega }/c\), \(\tau =\overline{\tau }/c\), Eq. (15) can be simplified to

$$\frac{\partial \omega }{\partial t}+{\bf{u}}\cdot \nabla \omega =D{\nabla }^{2}\omega +\frac{{\zeta }_{p}}{{I}_{p}}\left({(\nabla \times {\bf{u}})}_{z}-2\omega \right)+2{I}_{p}^{-1}\tau .$$
(16)

The above derivation is based on the assumption that a spinner is not self-propelled. In other words, \({\bf{u}}({\bf{r}}^{\prime} ,t)\) is dependent only on spinners located at \({\bf{r}}\ne {\bf{r}}^{\prime} \).

Fluid flow

Assuming Stokes flow, the fluid flow generated by a single spinner can be approximated as a point rotlet81, with torque \(8\pi \mu {a}^{3}{{\Omega }}\hat{{\bf{z}}}\). The resulting fluid flow is given by the solution to the following forced Stokes equation:

$$-\nabla {p}_{s}+\mu {\nabla }^{2}{{\bf{u}}}_{s}=-\nabla \times \left({\zeta }_{p}{{\Omega }}\delta ({\bf{r}}-{{\bf{r}}}_{0})\delta (z)\hat{{\bf{z}}}\right).$$
(17)

For distribution of N spinners with center of mass position r0i and angular velocity Ωi, the fluid flow can be written as

$${\bf{u}}=\mathop{\sum }\limits_{i=0}^{N}{{\bf{u}}}_{s}({\bf{r}}-{{\bf{r}}}_{0i},z,{{{\Omega }}}_{i}).$$
(18)
$${\bf{u}}({\bf{r}},z,t)={\int}_{A}\mathop{\int}\nolimits_{-\infty }^{\infty }{{\bf{u}}}_{s}({\bf{r}}-{{\bf{r}}}_{0},z){{\Psi }}({{\bf{r}}}_{0},{{\Omega }},t)d{{\Omega }}\,d{{\bf{r}}}_{0}.$$
(19)

The pressure distribution, p(r, z, t) can be computed similarly. By multiplying Eq. (17) by Ψ(r0) and integrating, the Stokes equation for the spinner distribution can be shown to be

$$0=-\nabla p+\mu {\nabla }^{2}{\bf{u}}+4\pi \mu {a}^{2}\nabla \times (\overline{\omega }\hat{{\bf{z}}}).$$
(20)

For flows with non-negligible inertia, Eq. (20) can be approximated for low Reynolds flow as

$$\rho \left(\frac{\partial {\bf{u}}}{\partial t}+({\bf{u}}\cdot \nabla ){\bf{u}}\right)=-\nabla p+\mu {\nabla }^{2}{\bf{u}}+4\pi \mu {a}^{2}\nabla \times (\overline{\omega }\hat{{\bf{z}}}).$$
(21)

Since the majority of experiments and previous theoretical work involve an effectively two-dimensional system25,39,60,61, we will therefore consider a two-dimensional fluid flow in the plane of the spinners. For the case of two-dimensional flow, a frictional drag term, −βu, can be added to the right side of Eq. (21). For small β values, this term does not qualitatively change the observed behaviors. The effect of larger β is more extensively studied in ref. 70. The form of Eqs. (1)–(3) is consistent with the continuum model derived for ferrofluids82 and the model used to analytically study the hydrodynamics for a system of micro-rotors83,84.

Spinner’s torque

A common spinner system consists of ferromagnetic colloids in either a rotating or oscillating magnetic field42. The external torque is then given by

$$\tau ={c}^{-1}{\mu }_{0}{({\bf{M}}\times {\bf{H}})}_{z},$$
(22)

where μ0 permeability of vacuum, magnetization is M, and the magnetic field is H. An oscillating magnetic field allows for spontaneous-symmetry breaking in the direction of rotation. Therefore, spinners can have either a CW or CCW spin44,85. However, a rotating magnetic field predetermines the spinner’s direction of rotation, causing the spinners to be either all CW or all CCW46.

We, instead, consider a system of Quincke spinners driven by a DC electric field. Quincke rotation has the advantage of being well studied86,87,88,89 and it allows for simultaneous populations of CW and CCW spinners.

Quincke rotation

For a fluid with conductivity and permittivity (σ1,ϵ1) and a suspended spherical particle with conductivity and permittivity (σ2,ϵ2), Quincke rotation occurs when

$$\frac{{\epsilon }_{2}}{{\sigma }_{2}}> \frac{{\epsilon }_{1}}{{\sigma }_{1}}$$
(23)

and the external electric field ∣E∣ is above a threshold given by

$${E}_{c}={\left[\frac{{\epsilon }_{1}{t}_{Q}}{2\mu }\left(\frac{{\epsilon }_{2}-{\epsilon }_{1}}{2{\epsilon }_{1}+{\epsilon }_{2}}-\frac{{\sigma }_{2}-{\sigma }_{1}}{2{\sigma }_{1}+{\sigma }_{2}}\right)\right]}^{-1/2}\ ,$$
(24)

where

$${t}_{Q}=\frac{({\epsilon }_{2}+2{\epsilon }_{1})}{({\sigma }_{2}+2{\sigma }_{1})}$$
(25)

is the Maxwell–Wagner relaxation time. The characteristic angular velocity for a Quincke spinner is

$${{{\Omega }}}_{Q}=\frac{1}{{t}_{Q}}\sqrt{{\left(\frac{| {\bf{E}}| }{{E}_{c}}\right)}^{2}-1}.$$
(26)

Neglecting the higher-order electromagnetic interactions89 between spinners, we instead focus on the spinners’ hydrodynamic interactions. If R is the distance between two rotors, the hydrodynamic forces scale as ~R−2 while the electromagnetic forces86 scale as ~R−4.

The external torque from the electric field is then given as

$$\tau ={c}^{-1}{({\bf{P}}\times {\bf{E}})}_{z}$$
(27)

with electric polarization density P(r, t) (dipole moment per unit area) governed by

$$\frac{\partial {\bf{P}}}{\partial t}+({\bf{u}}\cdot \nabla ){\bf{P}}={D}_{P}{\nabla }^{2}{\bf{P}}+{t}_{Q}^{-1}({{\bf{P}}}_{eq}-{\bf{P}})+{\boldsymbol{\omega }}\times {\bf{P}}$$
(28)

where the equilibrium polarization density is Peq is

$${{\bf{P}}}_{eq}=\frac{-c\epsilon | {\bf{E}}| }{1+{t}_{Q}^{2}{{{\Omega }}}_{Q}^{2}}\left(\hat{{\bf{x}}}+\frac{\omega }{| \omega | }{t}_{Q}{{{\Omega }}}_{Q}\,\hat{{\bf{y}}}\right).$$
(29)

and

$$\epsilon =4\pi {\epsilon }_{1}{a}^{3}\left(\frac{{\epsilon }_{2}-{\epsilon }_{1}}{2{\epsilon }_{1}+{\epsilon }_{2}}\,-\,\frac{{\sigma }_{2}-{\sigma }_{1}}{2{\sigma }_{1}+{\sigma }_{2}}\right)$$
(30)

Eq. (28)–Eq. (29) are provided in ref. 83 with the inclusion an artificial diffusion term for numerical stability.

Continuum model of a Quincke spinner fluid

We will assume a constant, uniform electric field E parallel to x–y plane. Without loss of generality, we can define \({\bf{E}}=E\hat{{\bf{x}}}\).

We nondimensionalize the system of equations in order to reduce the number of parameters. Equations (1)–(3) and Eqs. (26)–(29) can be nondimensionalized as Eqs. (31)–(36) using the scalings and dimensionless parameters shown in Table 1. X* is the corresponding dimensionless value of a variable X.

$$\frac{\partial {\omega }^{* }}{\partial {t}^{* }}+({{\bf{u}}}^{* }\cdot {\nabla }^{* }){\omega }^{* }= \, {D}^{* }{\nabla }^{* 2}{\omega }^{* }+\kappa \left({({\nabla }^{* }\times {{\bf{u}}}^{* })}_{z}-2{\omega }^{* }\right)\\ \, +2\kappa \gamma {({{\bf{P}}}^{* }\times \hat{{\bf{x}}})}_{z}$$
(31)
$$\,{\text{Re}}\,\left(\frac{\partial {{\bf{u}}}^{* }}{\partial {t}^{* }}+({{\bf{u}}}^{* }\cdot {\nabla }^{* }){{\bf{u}}}^{* }\right)=-{\nabla }^{* }{p}^{* }+{\nabla }^{* 2}{{\bf{u}}}^{* }+4\alpha {\nabla }^{* }\times ({\omega }^{* }\hat{{\bf{z}}})$$
(32)
$${\nabla }^{* }\cdot {{\bf{u}}}^{* }=0$$
(33)
$$\frac{\partial {{\bf{P}}}^{* }}{\partial {t}^{* }}+({{\bf{u}}}^{* }\cdot {\nabla }^{* }){{\bf{P}}}^{* }={D}_{P}^{* }{\nabla }^{* 2}{{\bf{P}}}^{* }+{{\bf{P}}}_{eq}^{* }-{{\bf{P}}}^{* }+{\omega }^{* }\hat{{\bf{z}}}\times {{\bf{P}}}^{* }$$
(34)
$${{\bf{P}}}_{eq}^{* }=-{\gamma }^{-1}\left(\hat{{\bf{x}}}+\frac{{\omega }^{* }}{| {\omega }^{* }| }{{{\Omega }}}_{Q}^{* }\,\hat{{\bf{y}}}\right)$$
(35)
$${{{\Omega }}}_{Q}^{* }=\sqrt{{\gamma }^{2}-1}$$
(36)

The system involves seven dimensionless parameters: the dimensionless spinner diffusion coefficient (D*), the particle Reynolds number (κ−1), the strength of the electric field relative to critical field (γ), the Reynolds number (Re) for the fluid, and the percentage of area covered by spinners (α), and a small artificial diffusion constant (\({D}_{P}^{* }\)) for the polarization evolution equation. Since Quincke rotation only occurs for E > Ec, we will operate for γ > 1. For the dimensionless form of the equations see Supplementary Note 1.

The parameters γ and α control how fast the spinners rotate and the corresponding effect on the fluid flow. D* determines the characteristic length for the interface between two opposite-spinning clusters. Therefore, the only parameters that could significantly change the qualitative behavior of the system are Re and κ.

In the initialization of ω*, each grid point is randomly assigned a small positive value (CW) or a small negative value (CCW).

Scaling and dimensionless parameters

We scale the length by the particle radius a, time by the tQ is the Maxwell–Wagner relaxation time tQ. The electric field E is normalized by the critical electric field Ec. Here N is the number of spinners, A is the area covered by spinners, and μ is the fluid dynamic viscosity. Correspondingly, the polarization P is scaled by Nϵ/A, where N is the number of spinners, A is the area covered by spinners, and ϵ is given by Eq. (30). Finally, the pressure p is normalized by μ/tQ, where μ is the fluid dynamic viscosity. Definitions of the model dimensionless parameters are given in Table 1.

Numerical methods

The reported results are numerical simulations of the continuum model using a square domain of length L, double periodic boundary conditions. The code is highly parallelized using CUDA to run on a NVIDIA graphics card. Spatial derivatives are computed using a pseudo-spectral method90. All Fourier and inverse transforms are computed using the fast Fourier transform (FFT). Equations (31)–(34) are then time-stepped using first-order exponential time differences91. For Stokes flow, the fluid velocity is solved using the stream-function formulation92. For non-Stokes flow, the fluid velocity is time-stepped using the method of Chorin Projection93. Details about the numerical implementation can be found in Supplementary Note 2.