Introduction

Active matter consists of agents, which include one or more building blocks that can convert energy into forces or torques or are externally driven, leading to an inherent motion1,2. The interactions among the individual active units can largely vary making the different systems exhibit an extensive number of structures, together with the emergence of a broad range of collective motions. These phenomena can be found across a wide range of length scales, e.g., in flocks of bird3, traffic dynamics4, swarming in bacterial colonies5,6, or cluster, swarm and lane formation in self-phoretic colloids7,8. Understanding the non-equilibrium physics of active matter is of major importance in unravelling the processes of life, since all living biological systems are driven far from equilibrium.

Active micrometre-scaled particles powered by externally imposed electric-magnetic manipulation serve as promising candidates for designing agents capable of applications such as nanomedical drug delivery9,10. A prominent example is given by particles carrying a magnetic moment exposed to externally applied magnetic fields9,11,12,13,14,15,16, which can be used to steer tracer particles through a solvent9, or be brought into a rotating state by imposed torques11,12,13,14,15,16,17 becoming then an archetypical synthetic chiral active fluid. Biological examples of chiral active fluids are membranes with rotating macromolecules18, or Volvox colonies which have been shown to hydrodynamically interact and organize in bound states of co-rotation19. Chiral particles are receiving much attention recently20,21 since they display a plethora of novel behaviours, which are far from fully explored. For example the formation of crystals at high rotor densities22, where dislocations move with whorl-like patterns. Most chiral systems show intrinsic density inhomogeneities, which can be related to the presence of an additional term in the system stresses that is odd under mirror or time-reversal symmetry23,24,25, which indicates the presence of the so-called odd viscosity. In soft matter chiral active systems, odd viscosity effects have so far only sparsely been realized experimentally26.

Active turbulence refers to the spontaneous collective complex spatiotemporal motion emerging in systems with active components27. It has been observed in systems such as bacterial suspensions28,29, swarming sperm30, or active nematics31,32. In classical turbulence, kinetic energy is externally inserted into the system at length scales in which viscous dissipation is negligible, and inertial effects transport the energy to smaller length scales where it is viscously dissipated. This is the so-called energy cascade33. In active turbulence, energy input is linked to the motion of the microscopic components and inertia is expected not to play an important role34,35,36. Self-organization phenomena of the active components and related instabilities develop correlated flows at various length scales resulting in a descriptively similar behaviour than inertial turbulence, where energy transfer across scales is not required, but possible. For both, inertial and active turbulence, self-similar dynamics are observed with the formation of vortices over a range of length scales, together with power-law dependencies of the energy spectrum.

Systems consisting of rotating or spinning particles have shown to either exhibit active turbulence, or odd viscosity26,37, but not both. Whether both features may occur simultaneously, and if they are even related is still an open question, answered affirmatively in this paper. Here, we report an experimental and numerical study of ensembles of rotating colloidal rods interacting dominantly via hydrodynamic interactions, where both the odd viscosity and turbulence can be simultaneously observed and quantified. We introduce a chiral active system composed of spinning silica rod-like colloids with a magnetic tip at one end, with the direction of the magnetic moment perpendicular to the rod axis (see the “Methods” section). The sedimented colloids follow an externally applied magnetic field almost instantaneously, and in a particular frequency range, they precess perpendicular to the container basis, resulting in an ensemble of synchronously rotating cylinders with parallel symmetry axes. In parallel, we conduct large-scale simulations with a model of discs rotating at a fixed angular velocity Ω, immersed in an explicit hydrodynamic solvent, and interacting with other discs by steric interactions38,39,40,41 (see the “Methods” section). This assumes that the relevant interactions between the rotating rods take place only perpendicularly to their cylinder axes. The fluid in contact with the colloid’s surface co-rotates with it, generating long-range hydrodynamic interactions together with two additional stresses between rotors: rotational and odd. The rotation stress exerted between pairs of rotors results in their propulsion, which leads to ensemble dynamics and rich cooperative effects, such as dynamic vortex formation of different sizes with an energy spectrum scaling similar to 2D inertial turbulence, all function of the packing fraction. Odd viscosity is a specific property of chiral rotor fluids, and therefore absent in usual fluids. The odd stresses appear orthogonally to the standard shear stresses leading to an inwards accumulation in the vortices, which are dynamically formed due to the rotation stresses. The weak compressibility of our system permits the emergence of clear density-vorticity correlations which allows us to reliably measure the odd viscosity of the system in bulk. We expect this to become a standard method to quantify odd viscosity in many other chiral active systems. The persistent formation of vortices with sizes ranging from a few particle diameters to almost the whole system size results in a decay of the energy spectrum with a power-law of approximately a decade in experiments, and significantly larger in simulations, enabling us to demonstrate the self-similar dynamics in this case of exceedingly small Reynolds number.

In the remainder of the paper, we (i) introduce the dynamics of two interacting rotors, (ii) study the translational propulsion of rotors in dilute to dense ensembles, (iii) characterize the vorticity and density fields, which provides a measurement of the system’s odd viscosity, and (iv) evaluate the active turbulence via the energy spectra.

Results

Hydrodynamics of rotor pairs

The characterization of the dynamics of two rotors allows a first quantification of the relevant interactions in the system. Trajectories of two interacting rotors are shown in Fig. 1a–d for a duration of Ωt = 30, which corresponds to t = 3 s in experiments. When two rotors approach, they perform an orbital rotation around their joint centre of mass. This orbital rotation is fast at short separations, as can be seen in the trajectories shown in Fig. 1a and c, and slows down with increasing distance (Fig. 1b and d), becoming Brownian and independent of each other for even larger separations. This behaviour can be understood as a result of the advection mutually generated by the hydrodynamic flows, which become weaker for increasing distances. The simulated hydrodynamic flow fields of two rotors placed at a fixed distance are shown in Fig. 1e. The angular velocity of the pair Ωpair around their centre of mass is quantified from the measured trajectories and presented in Fig. 1f, where both show a similar power-law decay as a function of the separation between the rotors. For an infinitely long isolated cylindrical rotor, the radial velocity profile can be estimated as uφ(r) = (σ/2)2Ω/r, and for two, the position-dependent flow field can be calculated as a superposition of the individual flow fields, yielding to a prediction of the angular velocity Ωpair/Ω = σ2/(2r2), where the boundary conditions on the surface of the rotors are disregarded. This quantitative prediction is roughly a factor of two higher than the experimental results, while it is in almost perfect agreement with the simulation measurements as can be seen in Fig. 1f: the analytical expression for the flow follows from the infinite cylinder approximation, which is valid for the 2D simulations.

Fig. 1: Hydrodynamic pair interactions.
figure 1

a–d Trajectories of two interacting colloidal rotors at various separations, with lines colour coding time evolution, and scale bars 2 μm in experiments and 2σ in simulations (see Supplementary Movie 1); a and b experimental trajectories for a duration of Ωt = 30, where t = 3 s, c and d simulated trajectories with Ωt = 5, with background points illustrating solvent positions. e Flow field of two colloids placed at a fixed distance, as calculated from hydrodynamic simulations. All simulation data is time-averaged over tav steps and a number nav of realizations. For averages of the fluid, there are δav = 10 MPC collisions in between each measuring step, and here tav = 3 × 105 and nav = 54. f Co-rotation angular velocity Ωpair as a function of the separation of the two rotors: red, experimental results; blue, simulations. For averages of the colloids δav = 103, and here with tav = 3 × 104 and nav = 240. The dashed line is the theoretical prediction Ωpair/Ω = σ2/(2r2). g Forces perpendicular to the line connecting two rotors in units of Fs = 4π2ησΩ. The dashed-dotted line is the analytical prediction from the Stokes equation in ref. 42, blue squares correspond to the explicitly measured forces in simulations, and circles are estimates obtained from the angular velocity in f via F⊥ = (kBT/D)πσΩ, with up to tav = 6 × 106 and nav = 108 for the largest separations.

We further evaluate the hydrodynamic force acting on a rotor, F⊥, known as the thrust force, which is responsible for the pair orbital rotation due to the presence of a neighbour, as illustrated in Fig. 1e. The results are normalized with Fs which assumes the Stokes drag of the colloid, i.e., ζv, with the solvent friction ζ = kBT/D, obtained from the diffusion coefficient D of isolated rotors, the ambient thermal energy kBT, and the velocity imposed at the colloid surface vs = πσΩ, such that Fs = 4π2ησΩ. We first assume that F⊥ exactly balances the Stokes hydrodynamic drag, which can be estimated by using the velocity from the measured pair angular frequency in Fig. 1f, v⊥ ≃ πΩpairr, and the solvent friction, such that F⊥ ≃ ζv⊥. Results are shown in Fig. 1g, where a similar trend for both, experiments and simulations can be observed. Furthermore, in simulations, the thrust force F⊥ can be explicitly measured by placing pairs of rotors at fixed positions as a function of their centres separation r. These measurements are in excellent agreement with the analytical expression for infinitely long cylinders in a fluid of vanishing \({{{{{{{\rm{Re}}}}}}}}\)42, as shown in Fig. 1g. The decrease of F⊥ with r indicates that the thrust force is significant over a long-range since it decays to the strength of thermal fluctuations Fth = kBT/σ only around r ≈ 10σ, since Fth/Fs ≈ 0.01 in our case.

Some small differences between experiments and simulations emerge from two effects present in experiments but not considered in simulations. In experiments sedimented magnetic rods rotate in a container of height much larger than the rod length, such that the induced flow partially escapes into the upper fluid layers above the rods. Furthermore, the friction between the solvent and the substrate could also lead to a diminution in the effective flow experienced by nearby rotors. Both effects are more pronounced at larger inter-particle distances. The overall agreement between experiments and simulations is very satisfactory, especially at short distances, where pair interactions are substantially stronger than the thermal noise.

Conversion of rotation into propulsion

The motion of more than two interacting rotors is now investigated as a function of the packing fraction ϕ. Rotor configurations and typical trajectories are shown in Fig. 2a–c for experiments and in Fig. 2d–i for simulation results. Note that accurate and automatic tracking of a large population of fast-moving colloidal particles in a dense suspension represents a significant challenge due to the limited image acquisition rate, dynamic heterogeneity and blinking of particles43. Therefore, we typically track 15–30 rotors manually for each experiment and overlay them onto raw images, with only a quarter of the images shown for clarity (see Fig. 2a–c and Supplementary Movies 2–4). At low ϕ, rotors are far from each other, displaying merely Brownian motion, as can be seen in Fig. 2a and d. At medium densities, rotors are more likely to interact in pairs, triplets, or larger ensembles, and the trajectories become ballistic and curved, reminiscent of active Brownian particles (Fig. 2b and e). Two neighbouring colloids rotate around each other in the same angular direction as the intrinsic rotation. When additional rotors get in their proximity, they get incorporated into the same motion nucleating a rotating group, or vortex. These vortices grow in size until they start to interact with neighbouring vortices. The ensemble exhibits therefore rich vortex dynamics where the rotating groups coalesce or break up dynamically, and individual rotors vividly interchange between different rotating groups. Although all colloids are spinning in sync with the external magnetic field, the resulting vortices rotate at different speeds, and they can eventually rotate in opposite directions when they need to adjust to the boundary conditions imposed by neighbouring vortices. The curvature of the trajectories shown in Fig. 2b and e depends on the configuration of nearby rotors, i.e., on the size of the rotating group the respective rotor belongs to. Small groups lead to strongly curved trajectories, whereas larger groups lead to almost straight trajectories. Curiously, a few experimental and simulated trajectories show a looping behaviour, resulting from the change of orbiting motion of one rotor from one to a different partner, triplet, or small group. For larger densities, the phenomenology is basically the same, but trajectories are on average longer and some of them might even exhibit segments bending into different directions, which occurs when rotors change from one vortex to another. Figure 2c and f show the system at the highest experimental trackable density (ϕ=0.14), and Fig. 2g–i, correspond to full simulation domains where the ensemble dynamics of multi-scale vortex formation is more obvious.

Fig. 2: Collective dynamics of hydrodynamically interacting rotors.
figure 2

a–c Fluorescent optical images of rotors at ϕ = 0.008, 0.031 and 0.144 with an overlay of a few typical trajectories. Here, only a quarter of the full probe chamber is shown (see Supplementary Movies 2–4) and \({{{{{{{\rm{Re}}}}}}}}\simeq 1{0}^{-5}\). d–f Simulated rotors configurations at ϕ = 0.007, 0.04 and 0.14 with all trajectories in grey and few highlighted ones for comparison with experimental images. Trajectories are all drawn for a time Ωt = 6. Only a part of the full simulation domain is shown and \({{{{{{{\rm{Re}}}}}}}}\simeq 0.09\). g–i Full simulation domains with trajectories, showing the formation of multiscale vortices. Dashed inset frames indicate the domains in (d–f). j Mean-square displacements for experiments (orange triangles) and simulations (purple bullets) at different packing fractions in a square simulation box of length L = 100σ with periodic boundary conditions. The lines correspond to the respective least-square fits to the active Brownian particle mean-square displacement. Simulation data was obtained with tav = 7 × 105 and nav = 1, and experiments with tav = 150. k, l Active velocity va, normalized by the rotor surface velocity vs = πσΩ, and rotational diffusion time τr with values and errors corresponding to the fits in (j). Dotted (black) lines correspond to the analytical prediction \({v}_{{{{{{{{\rm{a}}}}}}}}}\simeq \pi {{\Omega }}\sigma {(1-\phi /{\phi }_{{{{{{{{\rm{cp}}}}}}}}})}^{2{\phi }_{{{{{{{{\rm{cp}}}}}}}}}}\sqrt{\phi /\pi }\); dash-dotted lines take into account an additional multiplicative fit factor of 3.3.

Each individual rotor is propelled through the system with an instantaneous linear and rotational velocity, which depends on the local rotor configuration. On average the motion can be mapped to that of active Brownian particles, and the averaged mean squared displacement (MSD) shows the three expected regimes: purely thermal diffusion at very short times, actuated at intermediate times, and enhanced diffusion at long times. Rotational diffusion is related here to the change of direction of the motion instead of to the change of orientation of the rotor axis. Experimental and simulation results are shown in Fig. 2j for several available packing fractions. The measured MSD for each fixed ϕ are then fitted to that of an active Brownian particle44 which provides well-defined values for the actuated velocity va and the rotational time τr both in experiments and in simulations, as summarized in Fig. 2k, l. At low ϕ, the rotors barely interact with others, and in the limit of ϕ → 0, there is no propulsion on average. In fact, at very low densities, the experimental MSD barely shows active behaviour, rendering a reasonable fit impossible. This is remarkably different to ordinary active Brownian particles, which at low ϕ exhibit an active gas phase45 due to their inherent propulsion. With increasing density, the interactions between rotors become more frequent such that the active velocity first grows with ϕ. Upon further increase of the rotor density, the motion is increasingly restricted by steric interactions, and the effective fluid viscosity experienced by the rotors also grows, which eventually completely impedes their motion. These trends qualitatively explain the maximum of va at intermediate values of ϕ seen in Fig. 2k for both simulations and experiments. In order to provide a quantitative estimate of this dependence46, we assume first that the thrust velocity v⊥ of a rotor is simply due to pair interactions \({v}_{\perp }\simeq \pi {{{\Omega }}}_{{{{{{{{\rm{pair}}}}}}}}}\bar{r}\), as shown in Fig. 1f, and secondly, that the average distance between the rotors is that of a homogeneous system \(\bar{r}\simeq (\sigma /2)\sqrt{\pi /\phi }\), such that \({v}_{\perp }\simeq (\pi \sigma {{\Omega }})\sqrt{\phi /\pi }\). Simultaneously, the rotors dissipate momentum via mutual interactions and the rotor density increases the effective fluid viscosity experienced by the rotors47 which results in a decrease of the velocity at high density, as expected also for active Brownian particles48. The increase in viscosity for a 2D system of passive colloidal particles in linear order49 is η(ϕ)/η0 = 1 + 2ϕ, although in order to account for the abrupt increase when approaching the close packing density ϕcp, it is more appropriate to consider the phenomenological equation \(\eta (\phi )/{\eta }_{0}\simeq {(1-\phi /{\phi }_{{{{{{{{\rm{cp}}}}}}}}})}^{-2{\phi }_{{{{{{{{\rm{cp}}}}}}}}}}\)50. The drag force considered for the thrust velocity F⊥ ∝ η0v⊥ can then be considered to balance with the drag of the effective velocity in a dense system Fa ∝ η(ϕ)va which provides a full estimate of the effective active velocity va ≃ η0/η(ϕ) v⊥ in the full density range, as shown in Fig. 2k. This estimate agrees qualitatively very well with the simulation results, and perhaps somewhat surprisingly even quantitatively with the experimental measurements.

The measured τr obtained from the MSD in Fig. 2j and shown in Fig. 2l corresponds to a rotational diffusion time only in the dilute limit. At larger densities, this time is related to rotational diffusion but also to an intrinsic rotation or change of direction, due to the orbital motion in the vortices. To estimate τr, we consider the average time a rotor orbits another rotor before moving into a different orbit. In the regime of small densities, this time decreases with density; at intermediate densities, the presence of vortices of multiple sizes leads on average to τr remaining constant, while for very high densities, where the rotors are basically not moving, τr rapidly increases. In the limit of small densities, the time can then be approximated as \({\tau }_{{{{{{{{\rm{r}}}}}}}}}\simeq \bar{r}/{v}_{\perp }\simeq 1/(2{{\Omega }}\phi )\). This estimate works for both experiments and simulations, even at medium densities, as shown in Fig. 2l.

Vorticity and density fields: odd viscosity

To obtain the vorticity experimentally, we track all the particles (up to 104) between two frames with the assistance of ImageJ, to obtain the velocity field, and average over 6 realizations. Similar averaging is also performed in simulations. With the ensemble configurations and the related velocity fields, the corresponding coarse-grained vorticity ω = ∂xvy−∂yvx and local density ϕloc fields can be computed, as shown in Fig. 3. Areas with positive vorticity correspond to underlying vortices rotating in the same direction as the imposed magnetic field, areas with negative vorticity appear for vortices rotating in the opposite direction, which essentially fill the space in between the positive vorticity areas. Although qualitatively similar, the vorticity field measured in experiments is weaker than in simulations, possibly due to the friction with the substrate.

Fig. 3: Vorticity and density fields.
figure 3

Coarse-grained vorticity and local packing fraction fields of the rotor dynamics on a 2D square grid of binning length l0 = 10σ for simulation and experiment. Both, vorticity and density fields for simulations are generated from the same system configuration, such that the correlation can be directly observed. Similar to the experimental data. Superimposed streamlines and rotor positions are displayed. The packing fraction used in experiments, and simulations is ϕ = 0.144, ϕ = 0.14, respectively. The simulations have been performed in a square simulation box of length L = 300σ with periodic boundaries.

The corresponding density fields in Fig. 3 also show clear inhomogeneities. Positive vorticity areas tend to be more populated than negative vorticity areas, both for experiments and simulations. This accumulation indicates the presence of a radial pressure on the rotating vortex originating from a non-vanishing odd vicosity24,51. As a first quantification of this effect, the probability density distribution is calculated separately for areas of positive and negative vorticity, displaying that both distributions are clearly displaced relative to each other, as shown in Fig. 4a. For positive vorticity, the maximum of the distribution occurs for densities larger than the average density, and conversely for negative vorticity. This means that areas with ω > 0 tend to attract particles, and that they are depleted from ω < 0 areas. The separation of the maximum distribution for a given density is larger for the lower average density in the investigated cases, which also shows a broader distribution, and the trend can be clearly seen both in experiments and simulations.

Fig. 4: Vorticity–density correlations.
figure 4

a Normalized probability distribution of the local packing fraction ϕloc of the areas with positive vorticity (solid lines and solid symbols) and negative vorticity (dashed lines and open symbols). Lines correspond to simulation results and symbols to experimental ones. Results for lower volume fraction are shown in magenta (experiments ϕ = 0.075 and simulations ϕ = 0.08). Higher volume fraction results are shown in cyan (experiments ϕ = 0.144 and simulations ϕ = 0.14). b Normalized variation of the local rotor mass density ρ as a function of the local vorticity ω. The lines indicate least-squares fits according to Eq. (1) with which νodd can be determined. Averages made with tav = 1.2 × 104 and nav = 1, and experiments with tav = 1 and nav = 5, nav = 6 for ϕ = 0.075 and ϕ = 0.144, respectively. The inset shows the estimates for the normalized odd viscosity as a function of the density for experiments (red symbols) and simulations (blue symbols), error bars are the standard deviation of the least-squares fit.

In order to provide a quantitative characterization of the odd viscosity, at a coarse-grained level, the rotors can be described using the continuum theory for chiral active fluids. Then, for an incompressible fluid, the vorticity spreads diffusively, and the related stresses are compensated by the pressure. However, if the system permits weak density inhomogeneities, the stresses due to a non-vanishing odd viscosity point to the centre of circulation, which translates (see the “Methods” section) into an increase of the density linearly proportional to the given vorticity as24

$$\frac{{{\Delta }}\rho }{\rho }={\nu }^{{{{{{{{\rm{odd}}}}}}}}}\frac{\omega }{{c}^{2}}\,.$$
(1)

This expression is also valid when the vorticity ω is locally varying, such that Δρ = ρ(r)−〈ρ(r)〉 is the local density change with respect to the average density in the system, with c the propagation velocity of a colloidal density inhomogeneity, and νodd the kinematic odd-viscosity, i.e., the momentum diffusivity due to the presence of the odd stresses. Therefore, a circular flow of vorticity ω > 0 is experiencing stress forces pointing to the centre of circulation, leading to rotor accumulation in the centre of the circular flow24. Similarly, a circular flow ω < 0 leads to rotor depletion in the centre of circulation. The advantage of Eq. (1) is that it can be directly employed to quantify the system’s odd viscosity. From data such as in Fig. 3, a histogram relating local vorticities and local densities can be evaluated, both in experiment and in simulation. The results presented in Fig. 4b show a linear dependence between local density changes and vorticity in the simulations as well as in the experiments. Fitting the lines in Fig. 4b with Eq. (1) allows the accurate quantification of νodd/c2 at different densities in the system ϕ, with resulting values shown in the inset of Fig. 4b. The experimentally measured odd viscosity is in reasonable agreement with the numerical results, which reflects that in both systems the two-dimensional hydrodynamic rotation is responsible for the propagation of the stresses. In the investigated density range, νodd/c2 decreases with density, due to the decrease of the compressibility of the system with increasing density. This decrease of νodd/c2 can also be observed in the larger separation of the maxima for lower densities in Fig. 4a. In order to provide a quantification of the odd viscosity, we need an evaluation of c, the propagation velocity of a colloidal density inhomogeneity. A lower limit for this value is given by the colloids’ diffusion which opposes the density inhomogeneities with corresponding diffusion coefficient D, as c ≃ D/σ ≃ 0.2 μm/s. In that case, and for the employed system parameters at a density of ϕ = 0.075, the lowest limit for the odd viscosity is νodd ≳ 1.5 × 10−2 μm2/s. Note that viscous stresses are transported much faster than odd stresses, and the density inhomogeneities resulting from the odd stresses can only be observed in long-lived vortex flow (see Supplementary Movie 5). To the best of our knowledge, the odd viscosity has not been quantified in bulk but was only quantified before in a few very different systems, such as the so-called edge-pumping effect26 in experiments, via the power spectra of the surface waves, measuring the deformation of a flexible boundary, or as for a system of granular rotors, in simulations, measuring the normal stresses37. The direct relation between vorticity and density inhomogeneities as expressed in Eq. (1) can be though related to a much larger range of systems20,21.

Active turbulence

The colloid trajectories in Fig. 2g–i and the vorticity-density fields in Fig. 3 show the simultaneous presence of eddies of different sizes, which is clearly reminiscent of turbulence52. To provide a quantification of the turbulent dynamics, we investigate the rotor velocity field v(r) and its corresponding energy spectra Eq as a function of the wave vector q, defined in Eq.(11) (see Fig. 5).

Fig. 5: Energy spectra of the colloid dynamics.
figure 5

Energy spectra obtained from the power-spectral-density of the rotor dynamics as in Eq. (10) divided by the mean-square velocity. Blue bullets depict simulation data whereas red triangles denote experimental data. Dotted and dashed experimental and simulation lines correspond to comparable densities. The experimental data is multiplied by a constant factor to match the high-q simulation values. The minimum and maximum wave numbers are \({q}_{\min }\sigma /(2\pi )=3\times 1{0}^{-3}\) and \({q}_{\max }\sigma /(2\pi )=5\times 1{0}^{-1}\). Statistics here are similar to Fig. 4. Inset: Averaged flow field in a simplified vortex, rotors arranged in a circle of radius Rv, with a fixed separation \(\bar{r}\) between nearest neighbours, with tav = 105 and nav = 80.

Energy is injected into the system on the rotor size scale which determines a limiting maximum wave vector \({q}_{\max }\simeq 2\pi /\sigma\). Since the colloids only propel when nearby rotors break their local symmetry, energy propagates because of the hydrodynamic forces and is then transported across length scales following self-similar dynamics. On the other hand, the vortex size has a maximum limit given by the system size, such that the minimum possible wave vector is \({q}_{\min }\simeq 2\pi /(L/2)\). In the limit ϕ → 0, the amount of large-scale cooperative motion is almost negligible and energy is predominantly stored in small vortices such that the energy spectrum slightly increases towards large q-values. Consequently, at low densities, energy is not transferred to large-scale motions. For intermediate ϕ, vortices of all sizes are taking part in the dynamics. The energy spectrum then follows a power-law decay q−5/3. Approaching \({q}_{\max }\), the energy distribution deviates from the power-law and approaches a constant, which is related to the particle size. Simulations show the presence of very big vortices in Fig. 3, and the power law decay in Fig. 5 is valid almost up to \({q}_{\min }\). Experiments in Fig. 3 show that the maximum size of vortices is clearly smaller than the system size, and the power-law decay is valid in a smaller range, namely until q1σ/(2π) ≳ 0.02. This is most likely due to the dissipation of energy through the substrate friction, similar to the truncation of the inverse energy cascade witnessed in 2D classical turbulence. Here the truncation occurs at \({l}_{{{{\Gamma }}}_{{{{{{{{\rm{subs}}}}}}}}}}\simeq \rho {v}_{0}(\phi )/{{{\Gamma }}}_{{{{{{{{\rm{subs}}}}}}}}}\), with a linear substrate friction density Γsubs53. From the value q1, the maximum experimental vortex size can be estimated to be 50σ, which is consistent with the vortices in Fig. 3, and implies Γsubs = 7.63 kg m−2 s−1. This approach might therefore also be a way to quantify the substrate friction, for which alternatively specific experiments could be also conducted26. Finally, at very high densities, steric interactions suppress rotor motion, and the free formation of vortices, such that it becomes increasingly difficult to measure the energy spectra. In contrast to bacterial turbulence, there is no dominant vortex scale that is introduced by the interplay of hydrodynamics, alignment interactions, activity, and rotational noise34,35,54.

To get a more intuitive insight into the self-similar behaviour of the system, we consider a simplified picture of a vortex formed by a few rotors with fixed positions in a circle of radius Rv, as shown in the inset of Fig. 5. If we approximate the drag on each rotor induced by the flow of its neighbours, the only contribution remaining is, due to symmetry, tangential to the circular trajectory. This is the reason for the emergence of the circular arrangements, independent of their size, which become then unstable given the presence of collisions, compressibility effects, and thermal fluctuations.

Discussion

We have shown that rotating micrometre-sized particles are an interesting model system of chiral active matter, where the emergence of turbulence and odd viscosity can be simultaneously observed, as demonstrated here in experiments and simulations. While different types of rotating colloids are known to convert rotational into translational energy in symmetry-breaking situations, such as in the presence of confinement26,40, investigations of this effect in bulk were fragmentary, and the precise measurement of odd viscosity effects in low-Reynolds-number soft-matter systems was elusive until now11,12,46,55. Results of simulations and experiments are to a large degree in agreement, showing very similar behaviour and dependence on system variables. Individual rotors behave similarly to active Brownian particles with their propulsion and rotational diffusion dependent on the configuration of neighbouring rotors. Translational velocity and rotational characteristic time can be measured and satisfactorily compared with an analytical prediction, for all the range of available concentrations. We furthermore present an effective method that allows the quantitative measurement of the bulk odd viscosity, shown here for experiments and simulations, and of use for a wide range of systems, such as roller liquids16,55, chiral granular gases56, or colloidal Janus asymmetric rotors57,58.

Active turbulence in rotor materials is due to long-ranged hydrodynamic interactions, which here have been shown to lead to active turbulent dynamics without the emergence of a dominant vortex scale, and corresponding energy spectra with a −5/3 power-law, with the same exponent as predicted by Kolmogorov for inertial turbulent systems33. The same dependence has also been found in non-chiral active spinner systems consisting of a mixture of clockwise and counter-clockwise rotating particles at moderate Reynolds numbers, for a system of particles approximately three orders of magnitude larger than ours12, and also for a continuum theoretical approach at zero Reynolds number36. A recent classification of active turbulence27 distinguishes the phenomenology of systems with either polar or nematic order. Active turbulence in polar systems shows typically a dominant vortex length and the scaling regimes of the energy power spectra are nonuniversal. Meanwhile, active turbulence in nematic systems is characterized by universal exponents, together with the continuous creation and annihilation of topological defects. Spinners might therefore constitute a separate category of active turbulence since orientational order is not intrinsic to the system, there are no topological defects, and the scaling exponents seem universal.

The simultaneous emergence of odd viscosity and active turbulence is an inherent property of the system, provided that the rotors setup is chiral. The hydrodynamic stresses among colloids include rotational and odd contributions absent in usual fluids59. Rotational viscosity couples the internal rotation to the translation of the colloids leading to rich vortical dynamics, while odd viscosity implies an orthogonal coupling of shear stresses inducing an effective pressure pointing into or out of the centre of emergent vortices. Moreover, both mechanisms mutually support each other, since the odd viscosity-induced compactification at the centre of the turbulence-induced vortices stabilizes them, enhancing then the turbulent effect. Both rotational and odd viscosity have the same origin, appear together, and are here measured from the same system setup such that the simultaneous emergence of both phenomena is foreseen to be found in other chiral active systems where the active and odd stresses are mediated by long-ranged hydrodynamics. However, this could be hampered if compressibility is negligible, as is the case when colloids have strong cohesive, repulsive, or simply steric interactions, as those in systems with very high densities, or in compact clusters. We also expect the concepts here discussed to be useful for various related hydrodynamic rotor systems of biological relevance, such as rotating membrane macromolecules, algae, or sperm, and also for the design of microrobots, and microrobots assemblies18,19,60,61.

Methods

Experimental setup

Silica rods with a magnetic head (Fe3O4) are used. The ferromagnetic material is grown on the Janus rods to impose a permanent magnetic dipole moment perpendicular to the rod’s long axis62. Directional growth of silica from nanoparticle encapsulated microemulsion droplets63,64 is employed, followed by seeded growth of silica layers, which is a modification of the synthesis protocol of bare silica rods65, the Stober process. A scanning electron microscopy (SEM) image of these match-stick-like magnetic silica Janus rods is shown in Fig. 6a, where the inset shows a transmission electron microscopy (TEM) image highlighting the doping of magnetic nanoparticles at the head. This resulted in slightly tapered colloidal rods of 3.5 ± 0.3 μm in length, and a head and a tail of diameters of 0.77 ± 0.08 and 0.61 ± 0.05 μm, respectively, measured from SEM images of ~70 particles. The ζ-potential of these particles are ~−65 mV. Each rod possesses a permanent magnetic dipole moment at the end approximately perpendicular to its long axis. These particles were then suspended in deionized water (Millipore, 18.2 MΩ) and loaded into a custom sample chamber built by glueing a Teflon cylinder (internal dimension: 1 cm; outer dimension: 2 cm; height: 1 cm) onto a piece of the coverslip. The chamber was cleaned with isopropyl alcohol and DI water thoroughly before being dried with nitrogen gas. The sample was allowed to rest on a microscope stage for 10 min until all the particles sediment to the bottom. A rotating magnetic field of 150 Gauss at constant angular velocity was applied by a pair of Helmholtz coils. All experiments were conducted at room temperature on an inverted light microscope (Olympus IX73) equipped with a ×60 oil-immersion lens (NA of 1.42) and the images were captured by a Ximea colour camera (MQ042Cg-CM). The centroids of rotors were determined by a standard Matlab routine43.

Fig. 6: Experimental system and dynamics of single rotors.
figure 6

a An SEM image of ferromagnetic Janus rods; the inset is a transmission electron microscopy image of the magnetic head of a rod. b Experimental setup with two pairs of orthogonally placed Helmholtz coils for generating a rotating magnetic field. c State diagram of a single rod under a rotating magnetic field as a function of the frequency and the field strength. Insets show images of rods corresponding to the three dynamic states observed, horizontal rod rotating synchronously (light green squares), standing rod rotating synchronously (red circles), and horizontal rod rotating asynchronously (dark green triangles) with applied dynamic magnetic field. Most of the experiments in the later work are conducted at a dynamic magnetic field of 10 Hz and 150 Gauss (red bullet); Scalebars are 2 μm in a, c, d, 10 cm in b, and 200 nm in the inset of (a). d Experimentally measured mean squared displacement of a single rotor (red circles) and its fit to 〈r2〉 = 4Dτ, with D = 0.171 ± 0.001 μm2/s; Inset shows typical trajectories of rotors. e Averaged fluid velocity uφ(r) in an axial direction as a function of distance from the colloid centre. Symbols are simulation results and the solid line corresponds to the analytical prediction, uφ(r) = σ2πΩ/(2r), with tav = 106 and nav = 4. Inset: Simulated 2D fluid averaged velocity field around a single rotor.

The magnetic field is generated by two pairs of orthogonally placed Helmholtz coils hosted on a home-built microscope stage, Fig. 6b. Under a slowly rotating magnetic field (B = 150 Gauss), the rod lays flat on the substrate due to gravity, with the long axis aligned perpendicular to the applied field, Fig. 6c. As a rotating magnetic field is applied between a certain range of frequency, typically 2−20 Hz, the rod stands up against gravity and rotates synchronously with the applied field, Fig. 6c, which we therefore term a rotor. In this study, we primarily focus on a rotating frequency of 10 Hz. We first examine the translational motion of a single rotor or a dilute suspension of rotors when they are far apart by measuring the mean squared displacements (MSD, 〈r2〉), which grows linearly with time, Fig. 6d, allowing the determination of the translational diffusion coefficient D. Typical trajectories are shown in the inset of Fig. 6d.

Simulation method

The employed numerical method, multiparticle collision dynamics, is a mesoscopic simulation technique to simulate fluids that does not rely on the microscopic degrees of freedom and thus is not in need for calculating all the interactions between the solvent particles. The method includes hydrodynamic interactions and thermal fluctuations. Provided that suitable parameters are employed, the correct low-Reynolds number behaviour is reproduced66. The algorithm basically consists of two alternating steps. In the streaming step, the positions of the fluid particles are ballistically updated, i.e., ri(t + h) = ri(t) + vi(t)h. The second step is the collision step, in which the fluid particles are sorted into square collision boxes of length a = σ/6 and exchange momentum with all particles in a given collision box according to a certain protocol. The collision routine we employ has already been introduced in ref. 67. It builds on the basic collision routine in which each fluid particle’s relative velocity is rotated by an angle of ±π/2 with respect to the centre of mass velocity in a collision box, with equal probability. When studying rotating objects, and to avoid the occurrence of unphysical torques in the fluid, angular momentum conservation is necessary. We employ a variant of the collision routine that conserves linear and angular momentum68 but not energy. The colloid rotation constitutes a persistent input of energy in the system, which is compensated by considering a thermostat in each time step to all collision boxes individually69, becoming also a guarantee of a constant system temperature kBT. For the fluid, we employ an average number of particles per collision cell n = 10 and take the collision time, i.e., the time between two collisions, as \(h=0.02a\sqrt{m/({k}_{{{{{{{{\rm{B}}}}}}}}}T)}\), yielding a viscosity of \(\eta =17.9\sqrt{m{k}_{{{{{{{{\rm{B}}}}}}}}}T}/a\), according to ref. 68.

The rotors are modelled as impenetrable moving and rotating no-slip boundaries that exchange linear and angular momentum with the fluid in the streaming step and in the collision step, by introducing virtual particles39,40. From the mean squared displacement, the diffusion coefficient in the dilute regime can be determined to be \(D=3.73\times 1{0}^{-4}{\sigma }^{2}/(a\sqrt{m/({k}_{{{{{{{{\rm{B}}}}}}}}}T)})\). In the concentrated regime, rotors sterically interact with each other via a purely repulsive Lennard–Jones potential,

$$\begin{array}{l}{{{{{{{\mathcal{U}}}}}}}}({l}_{ij})=\left\{\begin{array}{ll}4\varepsilon \left[{\left(\frac{a}{{r}_{ij}-\sigma }\right)}^{12}-{\left(\frac{a}{{r}_{ij}-\sigma }\right)}^{6}\right]+\varepsilon \,,\quad &{{{{{{{\rm{for}}}}}}}}\,{r}_{ij}\le \sigma +{2}^{1/6}a\\ 0\,, \hfill \quad &{{{{{{{\rm{else}}}}}}}}\end{array}\right.\end{array}$$
(2)

The rotors are therefore simulated as 2D impenetrable discs of diameter σ, and always roughly one collision box of fluid is between two rotors to ensure proper hydrodynamic coupling. Thus, after the momentum exchange between colloid and fluid, the positions of the rotors are updated according to a molecular dynamics scheme. The mass density of the colloids’ material and the fluid mass density are taken to be the same. The angular velocity Ω of the colloids is fixed to \({{{\Omega }}}_{0}=0.01857/(a\sqrt{m/({k}_{{{{{{{{\rm{B}}}}}}}}}T)})\), if not otherwise stated. The interactions of the fluid particles with the colloid surface generates a co-rotation of the fluid with a velocity decaying in the azimuthal direction, as shown in the inset of Fig. 6e. The radial velocity profile quantitatively agrees with the measurements of the simulated flow, as shown in Fig. 6e. For the calculation of the packing fraction, we consider the hard cores of the rotors, that cannot be occupied by fluid, i.e., an ensemble of N rotors in a square simulation box of length L has a packing fraction of ϕ = Nπ(σ/2)2/L2. Simulations are performed in a square simulation box with periodic boundary conditions, with a box length L = 300σ unless otherwise stated.

Simulation code was developed in CUDA C/C++, NVIDIA A100 GPUs in the JUWELS supercomputer and is typically used to simulate 30 × 106 solvent particles and up to 70,000 rotors.

Dipolar magnetic interactions

Provided the use of magnetic colloids, we experimentally quantify pair interactions parallel to the line of centres by estimating the effective pair potential as a function of radial distance70, \(U(r)=-{k}_{{{{{{{{\rm{B}}}}}}}}}T\log g(r)\). Results are shown in Fig. 7 for a dilute suspension of rotors for three values of the magnetic rotating frequency. At the lower frequency, the rotors only display repulsive interactions, while for increasing frequencies, a weak short-range attraction on the order of 0.2kBT emerges. This proves that dipolar magnetic interactions between rotors are negligible in experiments since they should be noticeable and decrease in intensity with frequency. The measurements of F∥ in ad hoc simulations, corresponding to Fig. 1g, are shown in the inset of Fig. 7 to be attractive, but much smaller than the thermal noise, than F⊥, and of shorter range. These measurements confirm that the dynamic behaviour of the ferromagnetic rotors is dominated by hydrodynamic interactions and that these are properly accounted for in the simulations.

Fig. 7: Effective pair potential.
figure 7

Estimated from the radial distribution function as \(U(r)=-{k}_{{{{{{{{\rm{B}}}}}}}}}T\log g(r)\) at density ϕ = 0.007 for simulations (in blue) and ϕ = 0.0002 for experiments (in red) at three input spin frequency values, showing the absence of significant dipolar magnetic interactions. Statistics similar to Fig. 1, inset same as Fig. 2j. Inset: Forces parallel to the line connecting the two rotors as measured in simulations, with error bars indicating the standard deviation of 7 measurements.

Dimensionless numbers

The comparison between experiments and simulations is done via dimensionless quantities. Although single rotors cannot be considered self-propelled particles, the magnetic activation can be quantified with the velocity at the colloid surface, vs = πσΩ, such that we consider the system Péclet number as Pe = vsσ/D and the Reynolds number \({{{{{{{\rm{Re}}}}}}}}={v}_{{{{{{{{\rm{s}}}}}}}}}\sigma /\nu\) with D the colloid diffusion coefficient, and ν the fluid kinematic viscosity. With specified default values, in experiments we are working with \({{{{{{{\rm{Re}}}}}}}}\approx 1{0}^{-5}\) and Pe ≃ 38, for which the rod-head diameter has been considered. Meanwhile in simulations, \({{{{{{{\rm{Re}}}}}}}}\simeq 0.09\) and Pe ≃ 20, respectively. Clearly, the input parameters do not perfectly match, but these dimensionless numbers ensure that both simulations and experiments are performed in the regime of low Reynolds and large Péclet numbers, where the same physical behaviour is to be expected. Note that in classical turbulence typical macroscopic relevant system lengths are considered in order to determine the Reynolds number. Meanwhile, here we consider the microscopic particle diameter σ as the relevant size, which is the typical choice for active systems. Even in the case that we would consider the average particle velocity (see Fig. 2k), and the size of the largest vortex, the relevant Reynolds number would just be 5–20 times larger, which is still exceedingly small.

Vortex dynamics in Stokes flow with odd viscosity

The Stokes equation of a chiral active fluid describes the time evolution of the flow velocity u, in terms of its density ρ, pressure p, kinematic viscosity ν, vorticity ω = εαβ∂αuβ (where the Einstein notation is considered, and εαβ is the Levi–Civita sysmbol), and importantly also by the odd kinematic viscosity νodd, which is proportional to the field of intrinsic rotation \(\widetilde{{{\Omega }}}={{\Omega }}\langle {\sum }_{i}\delta ({{{{{{{\boldsymbol{r}}}}}}}}-{{{{{{{{\boldsymbol{r}}}}}}}}}_{i})\rangle\), as24,51

$$\rho {\partial }_{t}{u}_{\alpha }=-{\partial }_{\alpha }\left(p-\rho {\nu }^{{{{{{{{\rm{odd}}}}}}}}}\omega \right)+\rho \nu {\partial }_{\beta }{\partial }_{\beta }{u}_{\alpha }\,,$$
(3)

Taking the curl of Eq. (3) leads to the vorticity diffusion equation

$${\partial }_{t}\omega =\nu {\partial }_{\beta }{\partial }_{\beta }\omega \,.$$
(4)

which can be solved with Fourier transform methods. The initial condition of a line or puntual vortex ω(r, t = 0) = Λδ(r) which later diffuses due to viscosity can be considered, in order to account for the internal vortex dynamics of the system together with Eq. (4). The solution is then

$$\omega ({{{{{{{\boldsymbol{r}}}}}}}},t)=\frac{{{\Lambda }}}{4\pi \nu t}{{\rm {e}}}^{-\frac{{{{{{{{{\boldsymbol{r}}}}}}}}}^{2}}{4\nu t}}\,.$$
(5)

Integrating this expression results in an expression for the velocity profile in polar coordinates when incompressibility is ensured, ∂αuα = 0, leading to

$${u}_{\varphi }(r,t)=\frac{{{\Lambda }}}{2\pi r}\left(1-{{\rm {e}}}^{-\frac{{r}^{2}}{4\nu t}}\right)$$
(6)

and ur = 0. These expressions of the flow field together with Eq. (3) yields to the following relation in the radial direction

$$0={\partial }_{r}\left(p-\rho {\nu }^{{{{{{{{\rm{odd}}}}}}}}}\omega \right)\,,$$
(7)

and thus p−p0 = ρνoddω, with p0 ≡ p(r → ∞). This means, that the pressure compensates the antisymmetric stress stemming from odd viscosity in order to satisfy incompressibility. If small changes in density due to finite compressibility are now considered52, and assuming

$${{\Delta }}p={c}^{2}{{\Delta }}\rho$$
(8)

with Δp = p−p0 and Δρ = ρ−ρ0, we obtain the expression for the density accumulation in Eq. (1), also employed for the measurements in Fig. 4b. In this case, c is the propagation velocity of a colloidal density inhomogeneity and might thus be linked to the diffusive spreading proportional to the colloids’ diffusion coefficient D.

Note that Eq. (3) might include two additional terms, one accounting for rotational friction among the rotors and the other for the friction between the rotor fluid and the substrate. The first of which is proportional to the rotational kinetic viscosity νR and takes the form \({\varepsilon }_{\alpha \beta }{\partial }_{\beta }\rho {\nu }_{{{{{{{{\rm{R}}}}}}}}}(2\widetilde{{{\Omega }}}-\omega )\), and is the coarse-grained version of the hydrodynamic coupling of intrinsic rotation and translational degrees of freedom. The second term is proportional to the substrate friction coefficient Γsubs and takes the form −Γsubsuα. However, both terms do not alter the result in Eq. (1), if incompressibility to arrive at Eq. (7), and \(\widetilde{{{\Omega }}}={{{{{{{\rm{const.}}}}}}}}\), i.e., a homogeneous rotor density, is assumed.

Turbulence analysis

To analyse the turbulent dynamics of the system, we investigate the velocity field v(r) and its corresponding energy spectra Eq. The formal definition of Eq is

$${E}_{{{{{{{{\rm{kin}}}}}}}}}=\frac{1}{2}\langle {v}^{2}\rangle =\int\nolimits_{0}^{\infty }{{{{{{{\rm{d}}}}}}}}q\,{E}_{q}\,,$$
(9)

with reciprocal space vector q and q ≡ ∣q∣. Using the Wiener–Khintchin theorem to express the correlation function 〈v2〉 in terms of reciprocal variables and assuming isotropy, we can write

$${E}_{{{{{{{{\rm{kin}}}}}}}}}=\frac{1}{4\pi }\int\nolimits_{0}^{\infty }{{{{{{{\rm{d}}}}}}}}q\,q\,{\langle \hat{{{{{{{{\boldsymbol{v}}}}}}}}}\cdot {\hat{{{{{{{{\boldsymbol{v}}}}}}}}}}^{* }\rangle }_{{{{{{{{\boldsymbol{q}}}}}}}}}\,,$$
(10)

where \({\langle \hat{{{{{{{{\boldsymbol{v}}}}}}}}}\cdot {\hat{{{{{{{{\boldsymbol{v}}}}}}}}}}^{* }\rangle }_{{{{{{{{\boldsymbol{q}}}}}}}}}\) is the two-dimensional Fourier transform of the velocity correlation function. The integral in Eq. (10) goes over values of \({\langle \hat{{{{{{{{\boldsymbol{v}}}}}}}}}\cdot {\hat{{{{{{{{\boldsymbol{v}}}}}}}}}}^{* }\rangle }_{{{{{{{{\boldsymbol{q}}}}}}}}}\) in radially symmetric shells in q-space. In a discretised version of Eq. (10), the integral is then evaluated as a sum over the discretised values in equal-q shells, i.e.,

$${E}_{q}=\frac{1}{8\pi {{\Delta }}q}\mathop{\sum}\limits_{q-{{\Delta }}q < k < q+{{\Delta }}q}{\langle \hat{{{{{{{{\boldsymbol{v}}}}}}}}}\cdot {\hat{{{{{{{{\boldsymbol{v}}}}}}}}}}^{* }\rangle }_{{{{{{{{\boldsymbol{k}}}}}}}}}\,.$$
(11)

Calculation of coarse-grained values

Experiments and simulations provide configurations at different times where the rotor’s positions are well-defined. Additionally, in simulations we obtain the instantaneous rotor velocities. In order to obtain the necessary density, velocity, and vorticity fields required in our study a coarse-grained procedure is applied. Density ρ(r) is obtained by averaging the rotor’s positions in a grid with a bin size that might vary, but typically (10σ)2. The bin size should be at least a few colloid diameters in order to identify coarse grain effects, but not much larger, since the structure of small vortices would already be averaged out. Coarse-grained velocities are obtained with two configurations at close times, v(r, t) = (r(t + Δt)−r(t))/Δt, with the time interval has been chosen such that ΩΔt ≤ 2, where in simulations higher resolution is possible, which is necessary at high densities due to frequent inter-rotor collisions. The velocity, and thus vorticity fields are then obtained by averaging the rotor velocities in each bin. For the calculation of the energy spectra, we employ a bin size of σ2 in order to obtain the highest resolution in q-space.