Emergence of Lanes and Turbulent-like Motion in Active Spinner Fluid

Assemblies of self-rotating particles are gaining interest as a novel realization of active matter with unique collective behaviors such as edge currents and non-trivial dynamic states. Here, we develop a continuum model derived from coarse-grained equations of motions for a system of discrete spinners. We apply the model to explore the mixtures of spinners and same-spin phase separation. We ﬁnd that the dynamics is strikingly sensitive to ﬂuid inertia: In the inertialess system, after transient turbulent-like motion the spinners segregate and form steady traﬃc lanes. Contrary, at small but ﬁnite Reynolds number, the turbulent-like motion is sustained and the spinner population exhibit a chirality breaking transition: only population with a certain sense of rotation survives. The results shed light on the dynamic behavior of non-equilibrium materials exempliﬁed by active spinners.

Simulations based on agent-based models, without hydrodynamics, for gear-shaped spinners [36,53] and rotating dimers [54] showed that mixture of clockwise and counterclockwise rotors segregate into same-spin phases. This behavior was actualized in an experiment [55] us-ing 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 [37][38][39]43]: large-scale dynamic patterns, such as lanes and vortices. While much progress has been gained by studying these systems using these discrete-spinner models, this approach becomes computationally expensive when considering large systems with a large spinner density. In such cases, it becomes favorable to utilize continuum models that replace the discrete individual particles with a continuous distribution for the entire suspension of particles. Such continuum approach allows for larger systems to become computationally feasible as well as easier to study the collective dynamics of the system [27,[56][57][58][59]. The existing continuum approaches for self-rotating systems are phenomenological [40,53,54,60]. This makes it challenging to assess if the observed dynamics are physical or dependent on the chosen continuum formulation. It also becomes difficult to directly compare experimental results against the phenomenological model. Here we consider a monolayer of spherical spinners with clock-and counterclockwise spins suspended in liquid in a three-dimensional domain. The model is motivated by experimental realizations such as spinners trapped at a fluid interface [43] or on a solid substrate [40]. Assuming the spinners remain confined to a plane and using the mechanics of individual spinners, we derive a continuum model. We then present numerical simulations assuming rotation due to an external DC electric field, known as Quincke rotation. We have explored the mixtures of clockwise (CW) and counterclockwise (CCW) spinners and studied the same-spin phase separation. We have found that the dynamics is strikingly sensitive to fluid inertia: while in the inertialess system af- Blue and yellow color corresonds to clockwise (CW) and counterclockwise (CCW) rotation. The axis of rotation for each spinner is shown as gray vertical line oriented in the zdirection. (b) The Quincke effect is an electrohydrodynamic instability which gives rise to a torque on a dielectric particle in a uniform DC electric field [61] if the free-charge polarization is antiparallel to the applied field. Above a critical field strength E > EQ rotation in the plane perpendicular to the electric field (Ω · E = 0) is induced 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.
ter long transient turbulent-like motion the spinners settles and form steady traffic-like lanes, at small but finite Reynolds number the turbulent-like motion is sustained indefinitely long. For even larger Reynolds number, initially equal CW/CCW spinner population exhibits a chirality breaking transition and collapses to a single sense of rotation state. While our results are obtained for socalled Quincke spinners, similar models can be applied to magnetic spinners as well. Furthermore, our results shed light on fundamental properties of generic active systems(e.g., spinners, rotating bacteria, microtubules assays, etc.) where large-scale chiral states emerge as a result of spontaneous symmetry breaking.

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 the xy-plane, see Fig. 1a. We then assume a uniform spinner concentration, c = N/A, where A is the area of the two-dimensional domain. The system is then modeled by the following equations for the spinner angular velocity ω(r, t), fluid velocity u(r, t), and pressure p(r, t): where I p = 2 5 ma 2 and D = k B T 6πµa . The thermal energy of the system is given by k B T .
Equation (1) is the spin angular momentum balance of the spinners; the right hand terms account, respectively, for thermal noise, the torque due to fluid drag on the spinners, and the externally applied torque dependent on the physical system τ . Eq.(2) is the linear momentum equation of the fluid or Navier-Stokes with the addition of the final term, which drives the fluid flow due to the rotation of the spinners. Eq.(3) is the incompressibility condition since we assume constant fluid density ρ.
It is worth noting Eq.(1)-(3) are general for any spherical spinners. The specifics of the individual physical system are introduced in the form of the externally applied torque τ , which is what causes the spinners to spin. For the purposes of this paper, we will consider the case of Quincke spinners studied experimentally in [62], τ = c −1 (P × E) z , where P is the polarization vector of the spinners 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 Eq.(27)-Eq. (28). The full model we will now be considering is provided in non-dimensional form in the Method section, see Eq.(29)- (34).
All equations are nondimensionalized using the spinner radius a as the characteristic length scale and the Maxwell-Wagner polarization time t Q , which sets the magnitude of the rotation rate, as the characteristic time scale, see Methods section. The scaling analysis shows that the most important parameter governing the behavior of the system is the Reynolds number for the fluid, Re = ρa 2 /µt Q . For Re 1, the effect of inertial terms in Eq.(2) of is small; fluid inertia becomes important for Re 1. In the rest of the paper, we explore the system behavior as Re increases. Figure 2a and Supplementary Video 1 show that the spinner fluid exhibits phase separation into clusters of CCW (yellow) and CW spinners (blue), also seen in the discrete spinner model [38]. These clusters themselves are rotating, growing, separating, and reconnecting as also observed in the 'active coarsening' case in Ref [60]. The flow is localized the interface between CW and CCW clusters, see Fig. 2b and Supplementary Video 2. Physically, the fluid flow is canceled out in the region between same spin spinners. Hence, even though the spinners might rotate very fast, the flow might be very weak, especially in the interior of 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). A 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. 2c.

Collective Behavior
An important question concerns the long-term behavior of the system assuming the spinners are initially chosen to be randomly 50 CW : 50 CCW. Previous theoretical studies [38,39] of discrete micro-spinners have shown that the system exhibits separation into same-spin clusters and then eventually evolve into emergent patterns, such as lanes or vortices. This work however was restricted to Stokes flow (Re = 0), while recent colloidal experiments [43,47] suggest a non-negligible Reynolds number.
Using the continuum model, if we consider Re = 0 and an initial 50 CW : 50 CCW random configuration, lanes always eventually form, see Fig. 3a and Supplementary Video 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 by the fact that Stokes flow corresponds to the least dissipation. Segregation into same spin regions is thus favored since the flow vanishes in the interior of same spin region and the only flow remaining is confined to the interface between opposite spin regions. The flow (and thus the dissipation) is further 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 direction resulting in a lane configuration. A small, non-zero Reynolds number, such as Re = 0.01, leads to distorted lanes, see Fig. 3b and Supplementary Video 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.
Strikingly, for small but finite Reynolds number, e.g. Re = 0.1 and greater, the system just exhibits a turbulent-like motion with no emergent collective largescale behavior, as shown in Fig. 3c and Supplementary Video 1. If the initial state is lanes, instead of a random configuration, the lanes become unstable, for Re = 0.1, and dissolve into turbulent-like motion, see Supplementary Video 5. For an even larger Reynolds number, such as 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, these fluid flows initially exhibit turbulent-like motion, but eventually begin to favor a single spin orientation, see Supplementary Video 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 flow, u. As shown in Fig. 4, the PDF for lanes exhibits a clear bi-modal distribution due to the two sides of the lane. For the turbulent-like case, the PDF is Gaussian unlike the non-Gaussian behavior observed in experiments with colloidal spinners [63,64]. Overly populated tails in our simulation were found to be due to underresolved simulations.
Finally, it is worth mentioning that the discrete spinner model [38,39] 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 tensor [27].

Turbulent-like Flow
In many quasi-two-dimensional active systems exhibiting turbulent-like behavior, such as magnetic rollers [43], or swimming bacteria [65], 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 cascade [66]. 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 in the by taking the Fourier transform of the spatial auto-correlation of u: such that r 0 = x 0x +y 0ŷ and ... |k|=k is the average over all wavenumbers k such that |k| = k. Correspondingly, the energy flux is obtained as follows Here u <k is the low-pass filtered velocity field with the wave numbers outside |k| < k are set to zero: andũ(k, t) is the Fourier transform of u(r, t) [67]. 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) from small wavenumbers k. As shown in Fig. 5, the energy cascade scales as k −5/3 . However, the energy flux is negative only in very narrow range, Π E (k) < 0 for k < k 0 ≈ 0.016, which is not consistent with the energy cascade observed in 2D classical turbulence [67][68][69] and colloidal experiments [47,64]. We anticipate that this behavior is due to relatively small system size, and coarse-gaining of the microscopic energy infection scale. In the cases that exhibit lane formation at long times, the fluid flows initially exhibit transient turbulence-like motion with the same k −5/3 scaling, prior to settling into lanes.
The the turbulence inertial range is characterized by two lengthscales, the Taylor microscale [67] and the integral scale [70]. The Taylor microscale falls in between the large-scale eddies and the small-scale eddies. The integral scale L I can be extracted from the velocity autocorrelation function in Fig. 4c. It gives L I ≈ 100 for particle Reynolds number Re= 0.1. The corresponding r.m.s. (root mean square) velocity u rms 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, (∂ x u ) 2 = u 2 rms /λ 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 en-ergy spectrum exhibits a peak which 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 flow [63] and fluid vorticity are another useful metric for studying turbulent flows. As shown in Fig. 4, the turbulent flow case exhibits no characteristic lengthscale as seen in mesoscale turbulence. This behavior is consistent with colloidal spinner experiments [63,64].

Active Transport
The capability of this system for active transport can be measured by tracking the position of point particles s i that are being passively advected by the fluid flow, ∂s i /∂t = u(s i , t) and computing the mean-squaredisplacement (MSD), |s i (t + t 0 ) − s i (t 0 )| 2 . Additionally, the mean square velocity (MSV) for the fluid flow, |u| 2 , 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 Re values. Counterintuitively, the MSD and MSV increases as Re decreases. For increasing 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.
The curves for different Re values begin to qualitatively diverge in the long term. The lane formation cases, Re = 0 and Re = 0.01, show sustained ballistic motion in MSD (∼ t 2 ) and a relatively constant MSV, as all the motion is confined to the interface of the lanes. Yet for sustained turbulent-like behavior, Re = 0.1 and 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 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.

CONCLUSION
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 lane formations and turbulent-like behavior for different parameter values. Our work sheds light on active spinner materials and makes valuable predictions for experiments on the effects of fluid inertia and chirality symmetry breaking.
For the case of Stokes flow, the collective formation of the spinners into lanes is preferred as it minimizes the interfaces between counter-spin clusters and thus minimizes the energy of the system. This behavior is in agreement with the preexisting discrete spinner model. Yet for small, non-zero Reynolds flow, we observe increasing turbulent behavior with increasing Re: turbulent lanes for Re ∼ 0.01 and sustained turbulent behavior for Re ≥ 0.1. Furthermore, for non-zero Reynolds numbers we observed 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. Furthermore, we have found the mean square displacement (MSD), mean square velocity (MSV), the energy spectrum statistics for the fluid in order to quantitatively analyze the cases of lane formation and sustained turbulent-like behavior.
This continuum model can be applied to other experimental spinner system as discussed previously. While this work used periodic boundaries, future work could easily apply this model with both fixed, impenetrable boundaries and soft, deformable boundaries, as well as arrays of artificial obstacles. Another unexplored possible avenue is the case of non-uniform spinner densities. As the spinner densities, in experiments, are often nonconstant in space and time, it would be of interest to see how this would affect the collective and turbulent behavior seen in this paper. This derivation could also be easily generalized for a three dimensional suspension of neutrally-buoyant micro-rotors.

Derivation of Continuum Model
Let us consider a distribution of N spherical spinners of radius a and mass m confined to a 2D plane of area A. The translational velocity of a single spinner located at position r is then The translational motion of the spinner is the result of being passively advected by the surrounding fluid flow, u, and Brownian motion of the spinner. 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 given as D = k B T /6πµa, where k B T is the thermal energy and µ is the fluid viscosity. We are also assuming a lower spinner concentration such that the effect of steric repulsion is negligible.
It is useful to 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 Using Eq. (7)- (14) and the Fokker-Planck equation [71], it can be shown that the governing equation for Ψ is given by the following Smoluchowski equation assuming an incompressible fluid (∇ · u). From this point on the derivation of our continuum model resembles Ref [57]. It is useful to define a function . Ω such that We will therefore define the local spinner density as c(r, t) = 1 Ω , the local angular velocity per unit area as ω(r, t) = Ω Ω , and the local external torque per unit area as τ (r, t) = τ e Ω . Integrating Eq.(9) leads to We can ignore the boundary term in Eq.(11) as spinners with |Ω| − → ∞ are non-physical and inherently have zero probability, Ψ = 0.
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 A c dA = 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 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 leads to I pΩ (r, t) = ζ p 1 2 (∇ × u) z (r, t) − Ω(r, t) + τ e (Ω, t) (14) where I p = 2 5 ma 2 is moment of inertia and ζ p = 8πµa 3 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 Ω/D r 1, where D r = k B T /(8πµa 3 ) Using Eq. (14), Eq.(13) simplifies to Furthermore if we assume a constant spinner density c and define ω = ω/c, τ = τ /c, Eq.(15) can be simplified to The above derivation is based on the assumption that a spinner is not self-propelled. In other words, u(r , t) is dependent only on spinners located at r = r .

Fluid Flow
Assuming Stokes flow, the fluid flow generated by a single spinner can be approximated as a point rotlet [72], with torque 8πµa 3 Ωẑ. The resulting fluid flow is given by the solution to the following forced Stokes equation: For a distribution of N spinners with center of mass position r 0i and angular velocity Ω i , the fluid flow can be written as The pressure distribution, p(r, z, t) can be computed similarly. By multiplying Eq.(17) by Ψ(r 0 ) and integrating, the Stokes equation for the spinner distribution can be shown to be For flows with non-negligible inertia, Eq.(20) can be approximated for low Reynolds flow as ρ ∂u ∂t +(u·∇)u = −∇p+µ∇ 2 u+4πµa 2 ∇×(ωẑ). (21) Since the majority of the previous experiments and theoretical work involves an effectively two-dimensional system, we will therefore consider a two-dimensional fluid flow in the plane of the spinners. For the case of twodimensional 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 [60]. The form of Eq.(1)- (3) is consistent with the continuum model derived for ferrofluids [73] and the model used to analytically study the hydrodynamics for a system of micro-rotors [74,75].
Spinner's Torque A common spinner system consists of ferromagnetic colloids in either a rotating or oscillating magnetic field [42]. The external torque is then given by with permeability of vacuum µ 0 , magnetization M, and magnetic field H. An oscillating magnetic field allows for spontaneous-symmetry breaking in the direction of rotation. Therefore, spinners can have either a clockwise or counterclockwise spin [47,76]. However, a rotating magnetic field predetermines the spinner's direction of rotation, causing the spinners to be either all clockwise or all counter-clockwise [40]. We, instead, consider a system of Quincke spinners driven by a direct-current (DC) electric field. Quincke rotation has the advantage of being well studied [77][78][79][80] and it allows for simultaneous populations of clockwise and counter-clockwise spinners.

Quincke Effect
For a fluid with conductivity, and permittivity (σ 1 , 1 ) and suspended spherical particles with conductivity and permittivity (σ 2 , 2 ), Quincke rotation occurs when and the external electric field |E| is above the critical field threshold, where t Q = ( 2 +2 1 )/(σ 2 +2σ 1 ) is the Maxwell-Wagner relaxation time. Then, the characteristic angular velocity for Quincke spinners is given as Neglecting the higher order electromagnetic interactions [80] 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 forces [77] scale as ∼ R −4 . The external torque from the electric field is then given as with electric polarization density P(r, t), the dipole moment per unit area, governed by ∂P ∂t + (u · ∇)P = D P ∇ 2 P + t −1 Q (P eq − P) + ω × P (27) such that the the equilibrium polarization density is P eq is and = 4π 1 a 3 Eq.(27)-Eq. (28) are provided in Ref [74] with the inclusion an artificial diffusion term for numerical stability.

Quincke Model
We will assume a constant, uniform electric field E parallel to xy plane. Without loss of generality, we can define E = Ex.
It is then of interest to nondimensionalize the system of equations in order to reduce the number of parameters. Eq.(1)-Eq.(3) and Eq.(25)-Eq. (28) can be nondimensionalized as Eq.(29)-Eq.(34) using the scalings and dimensionless parameters shown in Table 1. X * is the corresponding dimensionless value of a variable X.
The parameters γ and α simply effect how fast the spinners rotate and corresponding effect on the fluid flow. Furthermore, D * determines the characteristic length for the interface between two opposite-spinning clusters. Therefore, the only parameters that significantly change the qualitative behavior of the system are Re and κ.
The spontaneous symmetry breaking between CW and CCW is introduced in the initialization of ω * . Each grid point of ω * is randomly assigned a small positive value (CW) or a small negative value (CCW).

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 pseudospectral method [81]. All Fourier and inverse transforms are computed using the Fast Fourier Transform (FFT). Eq.(29)-Eq.(32) are then time-stepped using first order exponential time differences [82]. For Stokes flow, the fluid velocity is solved using the stream-function formulation [83]. For non-Stokes flow, the fluid velocity is timestepped using the method of Chorin Projection [84]. Details about the numerical implementation can be found in the Supplementary material.

DATA AVAILABILITY
The data that support the findings of this study are available in the main text and supplementary information. Additional information is available from the corresponding author upon request.

CODE AVAILABILITY
The code is available from the corresponding author upon reasonable request.