Inertial clustering and emergent phase separation of spherical spinners

We study the hydrodynamics of spherical spinners suspended in a Newtonian fluid at inertial regime. We observe a spontaneous condensation of the spinners into particle rich regions, at low but finite particle Reynolds numbers and volume fractions. The particle clusters have a coherent internal dynamics. The spinners form colloidal vortices surrounded by the fluid depleted of the particles. The formation of vortices is observed both in periodic simulation box and when the spinners are confined between two flat walls. The stabilisation of the observed states relies only on hydrodynamic interactions between the spinners and requires a finite amount of inertia. The observations pave the way for the realisation of 3-dimensional spinner materials, where coherent structures and collective dynamics arise only from the rotational motion of the constituents.

Another possibility is provided by systems, where the dynamics of the individual building blocks is purely rotational. Nature's examples of this include the dancing of Volvox [25] and the formation of vortex arrays [26] of spinning bacterial cells. Previous studies of artificial colloidal spinner materials [27][28][29][30][31][32][33] have concentrated either to very low Reynolds numbers or to 2 dimensions (2D). At vanishing Reynolds numbers and at high area fractions simulations have predicted a phase separation of binary mixtures [28][29][30], the emergence of edge currents [31] as well as the stabilisation of 2D crystals [32], while the presence of an odd viscosity has been predicted for chiral active fluids [34]. Experiments of circular disks spinning on a gas-liquid interface at finite Reynolds numbers, demonstrated a dynamic ordering arising from the competition between magnetic attraction and hydrodynamic repulsion [5,6]. While attractive interactions has been observed for inertial spinners in a dense passive media [35,36] and ordered structures have been predicted for spinning disks at finite Reynolds numbers in 2D [33]. The effects of inertia in 3-dimensional (3D) spinner solutions are currently unknown.
We consider a simple 3D spinner system, consisting of spherical particles suspended in a Newtonian fluid and subjected to torques. Our results demonstrate that the initially uniform particle density is unstable: a spontaneous condensation of particle rich and particle poor domains from the ini- * juho.lintuvuori@u-bordeaux.fr tially uniform distribution of the spinners is observed. When both the particle volume fraction and Reynolds number are low but finite, the spinners spontaneously organise into colloidal vortices surrounded by a pure fluid. We demonstrate that the observed emergent phase separation and the collective motion originates solely from the hydrodynamic interactions between the spinners and requires a finite amount of inertia.

II. METHODS
We use a lattice Boltzmann method (LBM) [37], to study the dynamics of suspensions of rotationally driven spherical particles. The LBM was used to solve the quasiincompressible Navier-Stokes equation for the fluid flow [37]. The no-slip boundary condition on the particle surface is realised by bounce back on links methods [38,39] which can be modified to take into account the movement of the particles [40]. We considered a density matched solution ρ = ρ fluid = ρ particle . We set the LBM lattice spacing ∆x = 1, time unit ∆t = 1 and density ρ = 1, as customary in LBM. The spinners were modelled as spherical particles, with radius R = 6∆x (R = 2.1∆x in Fig. 4), with a very short range repulsion between them to avoid particle-particle overlaps [32,41]. A constant torque T (Fig. 1a) is applied on each particle. This leads to a spinning motion around a unique axis (X) which gives the (particle) vorticity direction (Fig. 1a). For a small particle Reynolds number (Re) an isolated spinner reaches a steady state rotational frequency ω 0 = T /8π µR 3 , where µ is the dynamic viscosity of the fluid. When Re is increased a deviation from ω 0 is expected due to inertial effects with the observed ω < ω 0 [42][43][44]. For the Re 20 considered in this work the deviation is reasonably small (Fig. 2), thus ω 0 was used to calculate the particle Reynolds numbers used in the text. The particles are placed in a three dimensional periodic rectangular box with a square cross section in Y Z plane and squashed along the vorticity direction X (Fig. 1b) (Two parallel solid walls are added in Y Z plane in Fig. 7a, b and c).
The dynamical state of the system is characterised by the particle Reynolds number Re = ρω 0 R 2 µ and the particle volume fraction φ = N 4/3πR 3 L X L Y L Z × 100%, where L X|Y |Z are the simulation box lengths.
Assuming a particle with a radius 100µm and using the kinematic viscosity of water 10 −6 m 2 /s, a Re = 1 would require a spinning frequency ω 0 = 100Hz. Using these values we can map a single simulation time unit to ∆t ≈ 5 × 10 −5 s. A typical simulation run consisting of 10 6 LB steps, corresponds to 50s in real time.

III. RESULTS
When starting from a random initial positions ( Fig. 1b) we observe that the hydrodynamic coupling between the spinners leads to the formation of small clusters and eventually to a phase separation (Fig. 1c). The particles spontaneously organise themselves into rotating clusters surrounded by clear fluid when φ ∼ 1% and Re ∼ 1 (see e.g. Fig. 1c for φ ≈ 1.1% and Re ≈ 1.2; as well as Movie 1 in [45]).
In order to characterise the formation dynamics of the vortices, simulations with a periodic box of 8R × 54R × 54R were implemented (Fig. 3). Here, the formation of only a single particle vortex is observed (Fig. 3a). In the steady state, the vortices span the periodic boundary along the vorticity direction, and the particles translate around the centre of the cluster with a tangential velocity V θ (r) (Fig. 3a). The spinners entrain the fluid leading to the formation of a liquid vortex with the same handedness as the particle one (blue arrows in Fig. 1c). Inside the vortex, a solid body rotation V θ (r) ∼ r is observed for both the particles and the surrounding fluid (symbols and solid lines in Fig. 3b, respectively). At the edge of the vortex the particle velocities drop slightly, while the fluid velocities start to decay (Fig. 3b). The particles exhibit random motion along the vorticity direction (X) showing diffusion-like dynamics ( Fig. 3d and   from the simulations (black symbols) as a function of the theoretical predictions calculated using Stokes' limit ω 0 = T /8π µR 3 . At values Re 3 the measured values start to deviate from the Stokes limit (green dotted line), following a prediction T /8π µωR 3 = 1 + Re 2 /1200 − 15647Re 4 /20744640000 (blue dash line) as expected [42][43][44]. Throughout this work the Reynolds number based on the Stokes limit Re = ρω 0 R 2 /µ is used for clarity. The simulation box is 27R × 27R × 27R in this validation.
The normalised rotational frequency (the slope of Fig. 3b) stays constant when the global volume fraction φ is increased (black and green symbols in Fig. 3b), but decreases when Re is increased (blue symbols in Fig. 3b). This implies that the energy conversion between particle rotation and translational motion is more efficient at lower Reynolds numbers. This is also manifested by the local volume fraction of the cluster (Fig. 3a): increasing the global volume fraction φ leads to a larger cluster at a constant Re (top three panels in Fig. 3a), while increasing Re with a constant φ results to less dense clusters (bottom two panels in Fig. 3a; see also Fig. 4 for detailed density maps).
The drive T gives a rotational energy input into the system, which is damped by a hydrodynamic friction. A total energy of the particles is defined as the sum of the rotational and translational energies: E tot = E rot + E tra . These are given by The m and I are the mass and moment of inertia of a spinner. The instantaneous velocity V i and the rotational frequency ω i for a particle i are measured from the simulations. When the vortex formation is observed, the total energy is dominated by the translational motion of the particles. In a steady state a typical value E tra /E tot ∼ 0.8 is observed (Fig. 3f). In the absence of the vortices, the total energy E tot is dominated by the rotation of the spinners (Fig. 3f).
Using the ratio E tra /E tot as an order parameter, we estimate the range of the vortex formation in the thin samples, when either Re is varied for a constant volume fraction φ or as a function of the φ for a constant Re. When the volume fraction is fixed at φ ≈ 2.2%, the vortex formation is observed for 0.3 Re 10 (Fig. 3f). Conversely, for a constant Reynolds number Re ≈ 1.2, a vortex is formed for φ 8% (Fig. 3g). For very low volume fractions φ 2%, the translational energy rests finite, typically E trans /E rot ∼ 0.5 (Fig. 3g), and the formation of small clusters is observed.

A. Inertial effects
The spinning particles of radius R are coupled to the surrounding fluid velocity field u by a no-slip boundary condition at the particle surface. When the flow is slow (Re ≈ 0) an isolated particle creates a rotating flow field with only an azimuthal component [27,46] where ψ and θ are the polar and azimuthal angles in spherical polar coordinates. This leads to a constant rotational motion of a particle pair around a central point [27,46]. As it does not include attraction nor repulsion, it is not expected to lead clustering. Indeed, no clustering is reported in simulations of spinners at low Re and φ [27,28] in agreement with our simulations (Fig. 3f).
Expanding to small but finite Reynolds numbers gives rise FIG. 4. The local particle volume fraction maps and the averaged local volume fraction ρ c (r) as a function of the distance r from the centre of the vortex, in the plane perpendicular to the particle vorticity direction (X axis), for various Reynolds numbers Re and global particle volume fractions φ . For a steady state in the phase separation regime, we observe a rotating cluster with a local particle volume fraction ρ C ∼ 0.1 surrounded by a pure fluid depleted of the particles. Both the radial and azimuthal particle density are isotropic and reasonably constant inside the vortex with a small peak at the periphery.

R r
The radial component (eq. 2) creates an advection towards the particle in polar regions |ψ| 55 • , while at the equator the fluid is advected away from the particle (Fig. 5a). The effects of the secondary flow has been attributed to the repulsion between two spheres spinning side-by-side [35,36,46] as well as to an attraction of a single spinner towards a flat wall along the vorticity direction [43]. The stabilisation of the clusters, requires the formation of structures via hydrodynamic interactions at Re ∼ 1. As a smallest possible building block, we considered a spinner pair (Fig. 5c-f). We measured the vertical ∆h and radial ∆r separations between the two particles (Fig. 5b). The trajectories show a general trend akin to the single particle flow fields: an attraction along the vorticity direction and repulsion at the equatorial region ( Fig. 5c and d). However, for a Re ≈ 1.2 the appearance of a limit cycle in the ∆h-∆r space is observed ( Fig. 5c; Movie 3 in [45])). The inertial forces stabilise the particles into circular orbits around each other in the plane perpendicular to the vorticity (Fig. 5e)  creased, the inertial repulsion between the spinners increases and the limit cycle becomes less pronounced and eventually disappears: no hydrodynamic bound state was observed for Re > 10 (Fig. 5d, f and g).
In the thin periodic simulation boxes, the formation of columnar vortices along the vorticity direction was observed (Fig. 1) and the particles have 3D trajectories inside the cluster ( Fig. 3c and d). These suggest complex 3D flow effects and the importance of the interactions along the vorticity direction. To highlight the importance of the 3rd dimension, we implemented a simulation where the particles were initialised randomly in a monolayer perpendicular to the vorticity axis (Fig. 6). In this case no vortex condensation was observed: the randomly distributed particles stay as a monolayer and initially explore the 2D space due to the mutual advection from the azimuthal flow fields. The radial repulsion leads to the formation of a stable hexagonal spinner crystal with no translational motion E tra ≈ 0 (see time ω 0 t = 7700 in Fig. 6a and  b). At ω 0 t ≈ 8000 a small perturbation of 2 particles along the vorticity axis (X) was introduced. This renders the monolayer unstable and the onset of vertical particle motion arising from the interactions along the vorticity direction is observed. The spinner trajectories become 3-dimensional and the formation of small hydrodynamically bound clusters is observed. Eventually, the hydrodynamic interactions lead to the formation of a stable spinner vortex (see ω 0 t > 8000 in Fig. 6).
It should be noted that a stabilisation of rotating clusters has been predicted to occur through hydrodynamic attraction between spinning disks in strictly 2D [33]. However, when the third dimension is opened, the interactions become repulsive [33]. Our results indicate that spheres spinning in 3D are markedly different. The single particle flow field gives rise to a hydrodynamic attraction along the vorticity direction and to a repulsion at the equatorial plane. The pair dynamics show more complex inertial effects. A hydrodynamic bound state is obtained for Re ∼ 1 in the rotating (∆r, ∆h)-frame ( Fig. 5c and g). The particles move around a central point with a velocity corresponding to a translational Reynolds number Re T ∼ 0.1. Similar states are also observed at early times in the bulk simulations (see e.g. Movie 1 in [45]). The stabilisation of the large clusters likely relies on more complex inertial many body interactions in the presence of a large scale vortical motion. In the stable vortex state, both the rotational and translational particle Reynolds numbers are Re ∼ 1 (See Appendix A). The spinners do not have clear pair dynamics, but instead the particles have helical trajectories on the plane perpendicular to the vorticity axis, while diffusive behaviour is observed along the vorticity direction (Fig. 3c, d and e).

B. Role of boundaries:
To highlight the robustness of our predictions and exclude periodic effects along the (attractive) vorticity direction as a driving force behind the condensation, we carried out simulations using a large rectangular box with two flat walls in the Y Z plane perpendicular to the spinning direction X, as well as with periodic boundary conditions (Fig. 7). Similarly to the thin samples, a spontaneous condensation of particles rich and particle poor regions is observed for Re ∼ 1 and φ ∼ 1% (Fig. 7).
When the spinners are confined (Fig. 7a), a re-entrant dynamics where a large scale vortical motion appears periodically is observed (Fig. 7a, b and c; see also Movie 4 in [45]). In the case of high Reynolds number, the system rests isotropic (see orange line in Fig. 7a for Re ≈ 19 and φ ≈ 2% sample). Spinning spheres have shown to be attracted to the no-slip surfaces perpendicular to the spinning direction [43]. This can have an effect on the stability and the formation dynamics of the vortices. The ratio E tra /E tot , calculated for the particles far away from the confining surface (wall-particle distance > 8R) (Fig. 7a) is lower than in the case of periodic box (Fig. 7d) and shows strong oscillations corresponding to the spontaneous appearance and disappearence of large scale vortices (blue and green lines in the Fig. 7a). These suggest that our predictions should be readily observable by experiments of a density matched colloidal or granular spinner suspensions confined in a Hele-Shaw type geometries.

IV. CONCLUSIONS
Using large scale numerical simulations we have studied a simple spinner system consisting of spherical particles spinning in a Newtonian fluid at weakly inertial regime in 3 di- The simulations were initialised with the particles in a monolayer at t = 0. In this case, the particles do not leave the monolayer and only experience in-plane hydrodynamic repulsion, leading to the formation of a 2-dimensional hexagonal structure at ω 0 t = 7700. A small perturbation along the X for 2 particles was introduced at tω 0 ≈ 8000. This leads to an onset of vertical particle motion, where the spinners start to explore the full 3-dimensional space (snapshots at ω 0 t = 8800 and ω 0 t = 9900). The particle 3-dimensional hydrodynamic interactions lead to initially the formation of small and short lived clusters (snapshots at ω 0 t = 11000 to ω 0 t = 16500). Eventually, the formation of a stable rotating cluster is observed at ω 0 t = 24750. (The simulations were carried with Re ≈ 1.2 and φ ≈ 2.2% using a periodic simulation box 8R × 54R × 54R). mensions. Our results demonstrate an unexpected instability where the rotational dynamics of the colloidal spinners leads to a large spatial density and velocity variations of the initially uniform suspension. At low but finite Reynolds numbers and volume fractions we observe a spontaneous phase separation between particle rich and particle poor areas. The formation of colloidal vortices rotating the same direction as the spinners is observed. The clustering is due to complex 3D flow fields and it requires a finite amount of fluid inertia. At early times, the simulations demonstrate the formation of small hydrodynamically bound clusters (early times at Movie 1 in [45]; see also Fig. 5c and Movie 3 for a bound state of an isolated spinner pair). More complex many body interactions then lead to the onset of collective motion and eventually to the formation of particle vortices.
The predictions should be readily observable in any density matched spinner suspension, provided that hydrodynamic interactions are dominant. The simulations of the confined sam-  (Fig. 7) suggest that interesting experimental possibilities are provided by reasonably large PMMA particles in transformer oil rotating due to Quincke effect [47] or millimeter sized particles with an embedded magnet in a rotating magnetic field [48], confined in a microfluidic chambers.