Phase separation of rotor mixtures without domain coarsening driven by two-dimensional turbulence

Unlike in thermodynamic systems, phase separation can occur without a thermodynamic driving force in active systems. How phase separation of purely hydrodynamic origin proceeds is an intriguing physical question. To this end, we study the phase separation of a binary mixture of oppositely rotating disks in a two-dimensional (2D) viscous fluid at an athermal condition by hydrodynamic simulations, focusing on the inertia effect. At symmetric and off-symmetric compositions, phase separation forms the oppositely flowing bands and a circular rotating droplet in the disordered matrix phase. In both cases, phase separation creates the largest structure directly from a chaotic state without gradual domain coarsening, unlike in the thermodynamic and corresponding dry rotor mixtures. We show that this unusual behaviour results from the nonlinear convective acceleration, i.e., the inverse cascade phenomena characteristic of 2D turbulence. Our finding reveals nontrivial nonlinear hydrodynamic effects on the self-organisation of active/driven particles in a fluid. Understanding the precise role of hydrodynamic interactions in the self-organization of circle-swimming active matter is an exciting avenue. Here, in a mixture of counter-rotating disks, the authors show that phase separation formed the largest size structure directly, without domain coarsening, and is driven by the inverse cascade phenomena characteristic of 2D turbulence.

O ver the past decade, there has been growing interest in active and driven matter, which not only is an interesting nonequilibrium dynamic system but also serves as a bridge between biological and physical worlds 1,2 . Active matter is an assembly of active agents, each of which consumes energy to move and thus is intrinsically in an out-of-equilibrium state. They often exhibit systematic self-propulsion 3 . Examples of such systems in nature include, in the increasing order of size, bacteria 4 , cells 5 , zooplanktons 6 , insects 7 , birds, and fish shoals 1 . In many of these examples, hydrodynamic interactions play a critical role in collective behaviour [8][9][10] . However, studying complex systems such as biological systems suffers from the drawback of having no individual control over their elements and thus often requires developing coarse-grained models theoretically. It is helpful to investigate model systems with simple activity for elucidating the physical essence of such unique out-of-equilibrium hydrodynamic states of active/driven matter.
Active motion can be classified into two fundamental types: translation and rotation. For translational motion, it has recently been found that active particles can phase separate even without attractive interactions, which is now widely known as "motilityinduced phase separation" (see, e.g., Refs. [11][12][13][14][15][16]. Very recently, a new type of torque-based phase separation of polar self-propelled particles has also been reported 17,18 . For rotational motion, phase separation was also found in mixtures of counter-rotating particles in dry 19 and wet [20][21][22][23] environments. Such self-rotating particles can also be seen in many biological systems, including spinning entities like bacterium Thiovulum majus 24 , dancing Volvox algae 25 , sperm cell clusters 26 and in synthetic experimental systems [27][28][29][30][31][32] . They produce interesting collective effects regulated by purely hydrodynamic means. It was also shown that a circular swimmer, in the limit of small radii of the trajectories, exhibits similar behaviour to rotors 33 . It is desirable to understand the roles of hydrodynamics and inertia in such a collective coherent motion of self-rotating active matter 8,9,[21][22][23]25,26,29,[34][35][36][37][38][39][40][41][42][43][44][45] . Two-body hydrodynamic interactions between rotating disks at low Reynolds numbers have been relatively well understood (see, e.g., Goto and Tanaka 21 for rotating disks and Fily et al. 20 for rotating spheres). However, many-body interactions are highly complex, especially with nonlinear convective acceleration 21 . Lately, the importance of inertia effects in the self-organisation of active or driven particles has also been recognised and studied intensively [46][47][48][49][50][51] . Furthermore, the turbulent nature of pattern evolution of active and driven particles in a fluid has also attracted much attention 27,29,31,32,[52][53][54][55][56][57][58][59] , known as "active turbulence". For example, it was shown by Linkmann et al. 58,59 that there is a transition between spatiotemporally chaotic and condensate states for an effective one-fluid model of active turbulence. Mesoscopic turbulence is also not necessarily identical. For bacterial turbulence, for example, the roles of excluded volume and rotlet dipole have been investigated recently 60 . From the perspective of phase separation of active matter, how the inverse energy cascade characteristic of two-dimensional (2D) turbulence leads to the formation of a condensate state coupled with phase separation of self-rotating particles is an intriguing open question.
To this end, we study the pattern evolution of mixtures of clockwise (CW) and counter-clockwise (CCW) rotating disks immersed in a 2D fluid, focusing on the turbulent nature. Here, we note that we consider a noise-free, athermal system. We investigate how the system self-organises under nonlinear manybody hydrodynamic interactions as a function of the disk composition (i.e., composition symmetry) and rotation speed (i.e., the Reynolds number). This problem of the self-organisation of selfrotating particles in a 2D fluid is related to two fundamental physical phenomena: 2D turbulence and phase separation.
From the viewpoint of 2D turbulence [61][62][63][64][65] , we discuss the similarity and difference between our rotating disk model and Onsager's point vortex model 66,67 . Onsager considered volumeless point vortices in an inviscid unbounded 2D fluid to understand 2D turbulence. This model has gained a renewed interest because quantum vortices formed in a Bose-Einstein condensate [68][69][70] can be regarded as a realisation of Onsager's point vortex model since vortices are discrete with circulations constrained to Γ = ± h/m (h: Planck's constant; m: the mass of a superfluid particle). On the other hand, our disks have a finite size and thus, the vortex should be regarded as a Rankine vortex instead of a point vortex. Furthermore, the surrounding fluid is viscous in our case, so the viscous dissipation cannot be neglected. It is also worth stressing that we rotate disks by imposing a constant torque individually, and fluctuations of the velocity field are self-generated by hydrodynamic couplings between rotors. Thus, our system's time evolution can be regarded as the selforganisation of Rankine vortices in a viscous fluid. We also emphasise that this process is accompanied by the spatial inhomogenization of the energy injection and its chirality.
On the other hand, from the viewpoint of phase separation, it was shown numerically 19 and experimentally 71 that oppositelyrotating disks phase separate even without a thermodynamic driving force, indicating its purely dynamic origin. Thus, an interesting question is (i) how phase separation of a system where only interparticle interactions are of hydrodynamic (or dynamic) origin is different from the thermodynamic phase separation 72 and the corresponding driven system without hydrodynamic interactions 19,71 and (ii) how the nonlinear effect (convective acceleration) affects phase-separation behaviour.
In the following, we show by numerical simulations that the phase-separated dynamical state of an oppositely rotating disk mixture has a rotating disordered or crystalline droplet morphology for an off-symmetric composition (OSM), and an oppositely flowing band structure for a nearly symmetric composition (SM). Strikingly, these large-scale patterns are formed directly from the initial homogeneous chaotic state without continuous domain coarsening. Comparing results with and without the nonlinear convective acceleration term shows that such peculiar phase-separation behaviours are exclusively due to nonlinear turbulent effects. These behaviours cannot be inferred straightforwardly from pair-level interactions between rotating disks due to the intrinsically nonlinear, many-body nature of hydrodynamic interactions. We discuss how nonlinear hydrodynamic interactions select such unique phase-separation paths and dynamical states based on the knowledge of 2D turbulence.

Results and discussion
Numerical methods. Many-body problems in colloidal systems induced by hydrodynamic interactions are often complex to solve theoretically. Thus, numerical simulations are the most promising approach to the problem. For this purpose, we can use the socalled squirmer models 8,9,39,40 or simpler methods that apply torque to particles directly 20,21,23,27,35,43,[73][74][75] . The former is for modelling self-propelled particles, whereas the latter is for modelling externally actuated particles, leading to different dynamic behaviours 76,77 . For mixtures of counter-rotating particles, the latter method has mainly been used. We also employ this method due to its technical merit. Strictly speaking, particles driven by externally applied forces should be called driven matter (or passive matter) rather than active matter, and their behaviours can be very different (see, e.g., Fily et al. 20 on the difference in the Stokes regime). Nevertheless, a driven system is expected to help understand how long-range, many-body hydrodynamic interactions induce particle self-organisation.
There are methods like Brownian dynamics, lattice Boltzmann, and Stokesian dynamics, but they are subjected to theoretical complexities (e.g., moving boundary problems), violation of the incompressible condition, and/or low Re approximation. Here, we adopt the Fluid Particle Dynamics (FPD) method 21,[78][79][80][81] , which is a model developed on the basis of physically transparent assumptions and is theoretically simple. In this method, solid colloidal particles are approximated by undeformable viscous liquid particles whose inside viscosity η c is much higher than solvent viscosity η ℓ . This allows us to get rid of the solid-fluid boundary condition and simulate the dynamical behaviours of colloids in a solvent by directly solving the Navier-Stokes equation (see Methods). This method can also incorporate nonlinear hydrodynamic effects at high Re, which is critical in this study. Note that our method includes both particle and fluid inertia (see the Methods). Since our method directly solves the flow field on the lattice under a periodic boundary condition, it can treat long-range many-body hydrodynamic interactions without complex treatments like the Ewald sum. Here we note that a similar numerical simulation method directly solving the Navier-Stokes equation was also developed by Höfler and Schwarzer 82 . For the details of our method, its validity, and analyses, see Methods and Refs. 21,[78][79][80][81] .
The set-up of phase-separation simulations. Using the FPD method, we study phase-separation kinetics of a mixture of counter-clockwise (CCW) and clockwise (CW) rotating particles driven by the same torque strength under a doubly periodic boundary condition. The disks interact with hard-core-like steric repulsions (see Methods for the detail) and hydrodynamic interactions. All results in the main text are for the total particle area fraction of 0.157, which is chosen because it was previously verified for the corresponding one-component system that a reentrant phase behaviour is observed at this area fraction 21 . For our study, we vary the fraction of CCW particles, Φ, between 0 and 0.5. To characterise the angular speed of the particle, we employ the particle-scale rotational Reynolds number, which is widely used in the field of active matter: Re ¼ ρa 2 Ω=η ' , where ρ is the density, a is the disk radius, Ω is the disk's angular frequency, and η ℓ is the liquid viscosity. Since Ω is a reaction to the applied torque and may depend on the environment of a rotor, we define Re theoretically as Re ¼ jαj=0:0064 (|α|: the rotational force density) (see Methods on the details). The range of the particle-scale Reynolds number we study is O(1 − 10), which is a typical Re range for experiments with magnetic rotors 27 and camphor 83 . This range of the particle-scale Reynolds number, Re, looks small, but the macroscopic-scale Reynolds number can be much larger; for example, for Re=10.8, if we estimate macroscopic-scale Re by using the characteristic domain size and velocity of the final states in our simulations as the relevant length and velocity, it is estimated as about 1000 (see below).
Phase separation occurs under many-body hydrodynamic interactions with nonlinearity. Nonlinearity gives rise to complex interparticle interactions, which can be seen from even pair-level interactions being not that simple (see Supplementary Fig. 1 for the pair interactions between rotors). The hydrodynamic nonlinearity generates the so-called Magnus force: The Magnus force acting on a 2D disk rotating with angular velocity Ω (|Ω| = Ω) in infinite 2D fluid flowing with a uniform velocity v is given by f M = 2πρa 2 v ⊗ Ω 84 , which depends on Re (Ω). Thus, a rotating disk with translational velocity is inevitably subjected to the Magnus force.
State diagram of phase separation. First, we show a state diagram as a function of Φ and Re (Ω) in Fig. 1a. Figure 1b shows the final steady-state structures observed in our simulations at various state points (Φ, Re). We can see two phase-separation morphologies, droplet and stripe patterns, primarily depending on Φ, i.e., the chiral symmetry, above the threshold Re (~1.4). For low Φ and Re ≥ 2.4, the minority CCW rotors form a rotating droplet cluster considerably denser than the matrix phase. As the mixture's symmetry increases to Φ = 0.4~0.5 (SM) for Re ≥ 2.4, the system forms two phase-separated bands, whose interfaces flow oppositely with undulation fluctuations. The two bands have the same particle density; thus, the particle density remains uniform throughout the system. See Supplementary Notes 2A-C and Supplementary Figs. 2-3 on the essential characteristics of phase separation. Unlike thermodynamic phase separation, we emphasise that the two phases are always composed of almost pure components.
We can see phase separation for Re ≥ 2.4, but not at Re=0.3 (see Fig. 1b). Thus, the border between the mixed and demixed states may be located around Re = 1. However, the precise location of the mixing-demixing border and the mixing-demixing hysteresis behaviour is obscured by the limitation of our simulation time (see Supplementary Note 3 and Supplementary  . We can robustly reach the final steady state within the reasonable computational time for Re ≥ 2:4, but cannot for Re ≤ 1:4 due to extremely slow dynamics. Finally, we note that our system's total particle area fraction (0.157) is lower than 0.2. According to Yeo et al. 22 , no phase separation should be observed in the linear Stokes regime at this volume fraction. Thus, the above demixing behaviours are expected to originate from nonlinear hydrodynamic effects, as confirmed by different methods later.  Figure 2a and b show the phase-separation processes for Φ = 0.3 and 0.5, respectively. We cannot see any clear indication of phase separation in the particle distribution until t = 10000 for both OSM and SM, although the flow field gradually self-organises even before this time but in a chaotic manner.
For OSM (Φ = 0.1~0.3), as Re increases, the homogeneous state with composition fluctuations first transforms to a phaseseparated state, where a rotating fluid droplet of CCW (red) particles exists in the matrix of CW (blue) particles with little difference in the particle volume fraction between them (see Supplementary Movie 1). After some incubation time (between t = 10,000 and 20,000), the rotating disordered fluid droplet further transforms through condensation into a rotating droplet with a much higher density (see Supplementary Note 2C and Supplementary Fig. 3d), generating the hexatic order in the droplet due to packing effects (see Fig. 2a and Supplementary Movie 2). Thus, this second transition involves the minority droplet phase's condensation by squeezing out the fluid between the rotors in the disordered fluid droplet initially formed. The particle density and the degree of order of the rotating droplet further increases between t = 20000 and t = 50000. In this process, the liquid droplet with a rough (i.e., broad) temporally fluctuating interface transforms into the denser crystalline droplet with a much smoother (i.e., sharper) circular interface. This incubation time to form the rotating crystal state becomes shorter and shorter with increasing Re. We can also see that the sharpening of the interface is accompanied by a steep increase in the droplet interface velocity (see the bottom row of Fig. 2a).
For SM, as Re increases, the homogeneous state with composition fluctuations transforms into a phase-separated state (see Supplementary Movie 3). In Fig. 2b, phase separation into two (CCW and CW) bands can be seen around t = 10000, and then the velocity field pattern further develops with time. We also note that the banded structure is formed along either x-or ydirection under the boundary condition. The band's interfaces fluctuate with time, and their undulation fluctuations are always accompanied by a giant vortex. The formation of the two phaseseparated bands is accompanied by a steep increase in the flow velocities in the interfaces between the two bands (see the bottom row of Fig. 2b).  The flow fields shown in Fig. 2 for t ≤ 5000 are chaotic in nature, reminiscent of 2D turbulence (see also Supplementary Movies 4 and 5). Although the particle-scale Re is relatively small, the macroscopic-scale Reynolds number characterised by the velocity and the domain size (the droplet size for OSM and the bandwidth for SM) steeply increases with time: For example, at Re = 10.8 (see Fig. 2), the macroscopic-scale Re in the final state reaches about 900 and 1100 for OSM and SM, respectively, which are large enough to induce the turbulent behaviours.
Comparing the composition field in Fig. 2 and the flow fields in Fig. 3, we can immediately notice that the fluid velocity magnitude v = |v| is maximum at the domain interface for both OSM and SM. The interface becomes very sharp only when a compact droplet with internal hexagonal order is formed, i.e., only for high-Re OSM. This is because a crystal-like droplet is dense and rotates as a whole in the matrix. We can see that the Magnus force, f M = 2πρa 2 v ⊗ Ω, is strongest at the interface, reflecting that the translational particle velocity v is maximum there (see Fig. 3b).
Similarly, the strongest velocity field in SM is generated at the interfaces flowing oppositely between the two phases. This behaviour can be naturally explained by the fact that rotors with opposite chiralities tend to move in the same direction (see Supplementary Fig. 1b). In this case, the interface is intrinsically unstable and dynamically fluctuates for the following reasons. Since the rotational directions of disks in the two coexisting bands are opposite, the Magnus forces acting on particles on the two sides of the interface act in opposite directions (see Fig. 3d). These opposing forces are balanced on an average throughout the interface but not locally, causing the interface fluctuations. A remaining problem for future investigation is to elucidate what determines the interface thickness.
We can also see well-developed vortices linked to the undulation of flow bands (see Fig. 3c), which is robustly observed in every run of simulations of SM with band formation. It is a natural consequence of bands being sandwiched by the interfaces rapidly flowing in opposite directions. This situation is reminiscent of Jupiter's zonal alternating flow bands and Great Red Spot 85 : The vortex (i.e., Great Red Spot) is trapped between bands flowing in opposite directions and is thought to be unable to move perpendicular to the zonal flow (north-south direction) and thus stable. There should undoubtedly be differences 86 , but interestingly, some similarities may exist.
Absence of gradual domain coarsening. In ordinary phase separation 72 and the corresponding dry rotor mixture 19,71 , domains increase their size gradually with time, known as domain coarsening. However, the phase-separation behaviour of counterrotating disk mixtures described above fundamentally differs from these two cases. The phase-separated structures are formed without gradual domain coarsening, whose process is accompanied by the self-organisation of flow towards the largest-scale coherent structures after the initial fluctuation period (see Fig. 2). We note that this feature is independent of the system size since the inverse cascade robustly leads to the maximum-size condensate directly without accompanying domain coarsening, as will be discussed later.
To characterise the structural evolution in the phase-separation process, we show the temporal evolution of the radial distribution functions g(r) for particle distributions without distinguishing the rotational direction (see Fig. 4a and b) and the structure factors  For OSM, g(r) is initially feature-less but shows the sudden emergence of the oscillatory decay in g(r) around t = 15,000, reflecting the crystal-like order formation in the droplet (see Fig. 4a). S(k) also exhibits the increase in the low-k part around t = 15,000, reflecting the sudden formation of a dense droplet (see Fig. 4c). The absence of the peak at a finite k and its temporal shift toward lower k indicate the absence of domain coarsening.
For SM, we see little temporal change in g(r) (see Fig. 4b) since (i) particles are almost homogeneously distributed in space even for the final steady-state (see Fig. 2b) and (ii) we do not distinguish the rotational directions in the calculation of g(r). On the other hand, S(k) shows the temporal increase in the low-k region (Fig. 4d) but without a peak at a finite k, reflecting the system-size phase separation into two bands. These results clearly indicate no gradual domain coarsening during phase separation in the presence of the nonlinear term, unlike ordinary phase separation.
We note that these pattern formations accompany the selforganisation of flow, which can be confirmed by the temporal evolution of the structure factors of velocity magnitude S v (k) (compare Fig. 4c and d with Fig. 5a and b). The S v (k) peaks at a finite k for SM, reflecting the velocity magnitude modulation along the y axis after the band formation, whose characteristic k is twice the smallest k of the system (see Fig. 2b).
Origin of phase separation without domain coarsening. Now, we show that the peculiar phase-separation behaviour without domain coarsening described above is due to the inverse energy cascade characteristic of 2D turbulence [61][62][63][64] in which the kinetic energy is transferred from the small energy-injection scale to larger scales.
We can see the importance of the turbulent nature by the emergence of the universal -5/3 scaling law in the turbulent energy spectrum E(k) (see Methods for the definition) in the final steady states at low Re (Re=0.3), where phase separation does not occur, as shown in Fig. 6a and b. We stress that this scaling is observed only when a condensate state is not formed. When a system-spanning condensate is formed by phase separation, we can no longer see such a scaling due to the formation of largescale coherent flow structures 87 . The system-spanning flow structures formed at Re≥1.4 produce low k upturns in the energy spectra (see Fig. 6a and b) and the structure factors (see Fig. 6c and d).
Next, we show in Fig. 7a and b how the energy spectrum E(k) develops with time for off-symmetric and symmetric cases at Re=10.8. We can see E(k)~k −5/3 in the early stage before phase separation occurs, but the shape of E(k) starts to deviate from the power law in the late stage, reflecting the formation of coherent flow fields and macroscopic phase-separated domains. Such a crossover behaviour from the turbulent regime characterised by E(k)~k −5/3 to the condensate state is typical of 2D turbulence 87,88 . The presence of the inverse energy cascade is further evidenced by the negative energy flux Π E (k) (see Methods for the definition) 58,87 , as shown in Fig. 7c and d, corresponding to Fig. 7a and b respectively. We can clearly see the negative energy flux below the energy injection wavenumber (i.e., the particle size wavenumber), indicating that the energy is transported from a small to a large scale by the nonlinear convective term, i.e., the inverse energy cascade. This kinetic energy transfer occurs with the time scale of momentum diffusion over the system size L, i.e., L 2 /ν (see Supplementary Note 6 and Supplementary Fig. 7). Thus, the kinetic energy spectrum in the form of E(k)~k −5/3 is already established in this time scale. At this initial stage, local rotational flow fields created by individual rotors are correlated spatially but have yet to be self-organised. Later, with time, the flow field becomes self-organised on a larger and larger scale. During this process, the total kinetic energy keeps increasing since the spatial coherence of flow fields generated by rotors gradually increases   Fig. 6 Re-dependence on the final energy spectrum E(k) and structure factor S(k). a, b The k-dependence of the kinetic energy E(k) for off-symmetric (a) and symmetric (b) cases as a function of Re. For low Re, where phase separation is not induced, we can see approximately the −5/3 scaling law (see the dashed line) below the wavenumber corresponding to the particle size, which is the characteristic size of our system's energy input. The scaling law suggests the inverse cascade characteristic of two-dimensional (2D) turbulence (see below). We stress that this scaling is observed only before the formation of the condensate.
When the system-spanning condensate is formed by phase separation, we can no longer see such a scaling due to the formation of large-scale coherent flow structures 87 : The low wavenumber upturn reflects the macroscopic flow organisation accompanied by phase separation. We can confirm it by looking at the corresponding low-wavenumber upturns of the structure factor S(k) for asymmetric (c) and symmetric (d) cases. c, d The k-dependence of the structure factor S(k) for asymmetric (c) and symmetric (d) cases as a function of Re. Here we distinguish the rotation direction, as in Fig. 4 (see Materials and Methods).
(see Fig. 8). Consequently, the largest-scale coherent vortex or banded flow structure is eventually formed without a gradual domain coarsening process. The macroscopic spatial flow organisation also leads to the low-k upturn in the kinetic energy spectrum of the final state. The kinetic energy reaches a plateau when the final steady state is established (see Fig. 8).
We can also see the indication of the hexatic ordering accompanying the condensation for OSM (see Supplementary Note 2C and Supplementary Fig. 3d) from a sharp increase in the total kinetic energy upon ordering (see Fig. 8a). For SM, on the other hand, the kinetic energy gradually increases since the band structure is formed by self-organisation of hydrodynamic flow 'without condensation' (see Fig. 8a).
It is interesting to note that even for particle-scale Re < 1 (e.g., at Re=0.05), E(k) has the k −5/3~2 feature. As shown below, such turbulent behaviour is absent if we ignore the nonlinear convective acceleration term v ⋅ ∇ v.
Demixing induced by nonlinear, many-body long-range hydrodynamic interactions. Here, we consider how the inverse cascade characteristic of 2D turbulence is induced in oppositely rotating disk mixtures. The phase separation is intrinsically dynamic since the energy is constantly injected into the system via the torques applied to the rotors. First, we briefly review Onsager's argument on the interaction among 2D point vortices to gain an intuitive picture of this complex nonlinear problem. He considered a system composed of a finite number of 2D 'point' vortices (composed of an equal number of CW and CCW vortices) in an 'inviscid' fluid (i.e., no dissipation). He derived a potential energy function describing the gas of point vortices interacting with the other vortices. Then, he showed that under a finite phase volume constraint, the system's entropy, S, initially increases but decreases eventually with increasing the input energy E m kin . In other words, above some energy input E m kin , the quantity dS/dE kin , which is formally equal to the inverse of the temperature T, may become negative. In such a situation, vortices with the same sign cluster together, and vortices with opposite signs repel.
The original Onsager's point vortex model was introduced for an unbounded, infinite system. However, it was shown that the model also works for various geometries, including doubly periodic boundaries 89 . The final state consisting of two isolated CCW and CW vortices clusters, a dipolar situation, was also recently directly observed experimentally for a driven BEC fluid 70,90 . Our system is under a doubly periodic condition, which imposes zero-mean vorticity. A dipolar situation has also been observed in our system when the total area fraction of disks is small (see Supplementary Fig. 6).
However, the situation is entirely different for the higher total area fraction, as shown in Fig. 1b. There is no change in the local density for SM, and the two types of particles phase-separate to form an alternatively flowing band structure. Only the minority particles form a droplet-like condensate when the composition symmetry is broken. The majority of particles are dispersed in a spatiotemporally chaotic state. The crucial difference between the point-vortex and our rotating disk systems primarily originates from the absence and presence of the Magnus force, respectively, as discussed below. The above observation of a dipolar situation for SM at a low total area fraction may imply that the low total area fraction can be approximated by Onsager's point-vortex model at least qualitatively, which is the limit of the total area fraction → 0. However, due to the Magnus force and viscous dissipation, such dipolar condensation does not occur for a high area fraction. Now, we consider the turbulent phase demixing in our rotating disk mixtures intuitively. First, we should note that the macroscopic-scale Re increases with time, even for small particle-scale Re, due to the self-organisation of the flow field. Thus, the final macroscopic-scale Re determined by the spacespanning phase-separated structure is much larger than the particle-scale Re and large enough to induce turbulence. Then, it is natural to consider that the self-organisation of flow to form a  large-scale vortex for OSM or flow band structure for SM is a consequence of the inverse energy cascade. Many-body nonlinear hydrodynamic interactions make the motion of rotating disks chaotic in the early stage. However, they create large-scale coherent flow in the late stage through the inverse energy cascade mechanism that transports the kinetic energy from the energy injection to system-size wavenumbers. Thus, phase separation proceeds while enhancing flow coherency (see Fig. 3a and c), i.e., the local uniformity of flow besides enstrophy fluctuations. We can confirm this by the temporal increase of the total kinetic energy during phase separation under the constant energy input, as shown in Fig. 8a. Thus, we can view phase separation as the dynamic process of enhancing the spatial coherence of fluid flows produced by the torques imposed on individual rotors.
We also note that the characteristic time of the kinetic-energy increase, τ 1/2 , which is the time when the kinetic energy reaches half of its final value, steeply increases with decreasing Re towards a non-zero threshold value separating mixing and demixing states (see Fig. 8b). This again shows the importance of the inverse energy cascade in the phase separation of our system. Next, we consider the roles of the Magnus force and steric repulsive force arising from the finite size of disks, both of which are absent for a rotating point object (a = 0). For example, a system composed of only disks rotating in the same direction selforganises into a densely packed droplet at a low volume fraction 21 . However, a further increase in Re destabilises the condensed droplet state and transforms it into a space-spanning hexatic-order state. It is because the repulsive Magnus force between like-sign rotors of nonlinear origin eventually wins over the hydrodynamic attractions above a specific Re (see below). Thus, the Magnus effect induces repulsive interactions among rotors with the same chirality when Re and the total area fraction are large enough (see Fig. 7 of Ref. 21 ). This situation is markedly different from Onsager's point vortex model, in which interactions among point vortices with the same chirality are attractive. Now, we consider an interesting fundamental question of under which conditions like-sign rotors do not form a condensate (as seen in the majority phase of OSM and both phases of SM). For OSM, phase separation accompanies the condensation of the minority phase but is eventually stopped when the minority droplet can no longer decrease its volume due to steric repulsion. However, this is not necessarily the case for a low particle fraction, as mentioned above. For example, when the total area fraction is small enough, a mixture forms a dipolar configuration composed of two oppositely rotating condensed droplets (see Supplementary Note 5 and Supplementary Fig. 6).
On the other hand, for SM at high Φ, phase separation does not involve condensation due to the repulsive Magnus force between like-sign rotors under the nearly symmetric composition. When the area fraction of disks is high enough, the Magnus force induces repulsive interactions between rotors 21 . However, unlike a system composed of only rotors of the same sign, the turbulent dynamic fluctuations prevent the formation of stable hexatic order, resulting in a dynamically disordered structure.
Finally, we qualitatively consider what determines the SM-OSM boundary. When the composition (chirality) asymmetry increases, the difference in the expansion pressure produced by the repulsive Magnus force between the two phases becomes more significant, eventually leading to the compaction of the minority phase. It indicates that a specific asymmetry level of Φ and the resulting imbalance in the interparticle hydrodynamic interactions are necessary for destabilising band structures. As shown in Fig. 1a, the band structure can be formed for the more asymmetric composition at higher Re, which might be due to the stronger repulsion among the same type of rotors. We need a further detailed study to understand the Magnus force's roles in pattern formation.
A critical role of the nonlinear convective term in pattern evolution. To confirm the turbulent nonlinear origin of the phase separation without domain coarsening, we also study the behaviour with the term ρdv/dt but without the nonlinear convective acceleration term v ⋅ ∇ v, although artificial and not physically relevant. We show typical phase-separation behaviour for Φ = 0.5 in Fig. 9a. We can see a distinct difference in the particle distribution and flow time evolution from the corresponding SM case with the nonlinear convective term (see Fig. 2b). The system phase separates into droplets made of CCW and those made of CW, and their sizes keep growing until the finite-size effect comes into play. The hexatic order emerges around t = 1000 for both CW and CCW droplets. The streamline patterns and the spatial distribution of the velocity magnitude are also shown in Fig. 9b, t=1000 t=5000 t=10000 t=15000 t=200 where we can see two types of counter-rotating vortices in the final steady state. Unlike the corresponding patterns for simulations with the nonlinear convective term (see Fig. 2b), the streamlines are very smooth from the initial to the final stage, and there are few fluctuations in the velocity distribution. For offsymmetric compositions, the final sizes of CCW and CW droplets are different, reflecting the ratio between the total amounts of the two types of particles. Besides this difference, the basic behaviours are the same as in the symmetric case.
We can see the gradual domain coarsening during phase separation in the temporal changes in g(r) and S(k) in Fig. 10, unlike the case with the nonlinear convective term (see Fig 4). From the peak wavenumber of S(k) in Fig. 10c and d, we also measured the temporal change in the average droplet size R as a function of time t and found that R~t 1/2 . We confirm that the elementary growth process is due to collision and coalescence between the same-chirality droplets (see Supplementary Movie 6). However, since these simulations at large Re without the nonlinear convective term may not have physical relevance, we do not discuss the phase-separation behaviour in detail.
For large Re in the absence of the nonlinear convective term, we have confirmed that a pair of rotors with the same chirality attract, whereas those with the opposite chirality repel each other (see Supplementary Note 7 and Supplementary Figs. 8 and 9). This pair-level interaction can straightforwardly explain the phase-separation behaviour in Fig. 9: Particles with the same chirality gather and form rotating droplets while preserving the chirality, whereas droplets with different chirality repel each other. Thus, the droplet coarsening proceeds sequentially by collision and coalescence of droplets with the same chirality. The nonlinear convection term and the resulting inverse cascade mechanism entirely change the character of phase separation from with to without domain coarsening.
Since the nonlinear convective acceleration is critical for the inverse cascade in 2D turbulence, the above results, together with the negative energy flux Π E (k) (see Fig. 7c and d), strongly support that the phase separation without gradual domain coarsening is a consequence of the inverse energy cascade characteristic of 2D turbulence, even though particle-level Re is moderate (Re = O(1~10)).
Finally, we note that the pattern evolution depends on the system size determining the lowest wavenumber since the inverse cascade transfers energy towards the lowest wavenumber limited by the system size. (see Supplementary Note 8 and Supplementary Fig. 10). This finite-size effect is an interesting issue for future investigation.

Conclusion
A dilute system of binary mixtures of oppositely rotating disks in 2D fluid demixes due to interparticle nonlinear hydrodynamic interactions,. For an off-symmetric mixture, a rotating cluster droplet of the minority phase is first formed as a transient state and then transformed into a rotating crystalline droplet accompanied by its condensation. In contrast, a phase-separated structure with band-like morphology is formed for a nearly symmetric mixture. Both cases have no gradual domain coarsening process, and system-size phase-separated structures are directly formed from the mixed initial homogeneous state through a chaotic state. This unique phase-separation behaviour, i.e., the absence of domain coarsening, comes from the peculiar long-range nonlinear nature of hydrodynamic interactions in a 2D fluid, i.e., the inverse cascade in 2D turbulence. It makes the phase-separation behaviour of our system intrinsically different from the thermodynamic and corresponding dry disk mixtures.
The phase separation behaviour is partly controlled by the effective attraction between disks rotating in the same direction and repulsion between those rotating in the opposite directions, similar to the interaction of point vortices, induced by nonlinear hydrodynamic interactions. However, a critical difference from the point vortices' collective behaviour arises from the finite size of disks leading to the Magnus force and steric repulsion.
In the thermodynamic phase separation of passive systems, the ordering process proceeds in such a way that the free energy is reduced. Consequently, a late-stage coarsening occurs in order to lower the interface energy. In this case, the dissipation is closely related to the interface reduction. This is also true for dry active rotor mixtures 19 and the corresponding model system 91 , where the dissipation is related to interface reduction due to the shortrange nature of the interaction between rotors. On the other hand, in our system, the formation of phase-separated structures is a consequence of the kinetic energy transport from a small to a large scale via the inverse cascade, dissipating energy through enstrophy dissipation. Thus, our study reveals that dissipation plays an essential role in these three types of orderings, but its role is essentially different from each other.
Our study would shed fresh light on complex out-ofequilibrium pattern formation in the driven/active matter under the influence of long-range many-body hydrodynamic interactions with nonlinearity. It also provides deep insight into the 2D turbulence of finite-size vortices. In our study, we use an externally driven system. Elucidating the difference and similarity between driven and active systems in the nonlinear regime is an interesting topic for future research.

Materials and Methods
Simulation method. The FPD method treats solid colloidal particles as undeformable viscous liquid particles whose inside viscosity is much higher than solvent viscosity. A system of N particles whose centres of mass are located at r i (i = 1~N) immersed in the viscous liquid is represented by a smooth viscosity profile as where η ℓ is the solvent liquid viscosity and η c is the viscosity of the particle. Here the smooth particle field is given by where a is the particle radius and ξ is the interface thickness. Then, the equation of motion to be solved is where ρ is the mass density, v(r) is the velocity field, and Σ ¼ pI À η ∇v T þ ∇v À Á is the internal stress of the fluid. The pressure p is determined to satisfy the incompressibility condition ∇ ⋅ v = 0, and I is the unit tensor. Here f U (r) is the force density due to the interparticle interaction determined as f U ðrÞ ¼ À ∑ N i ðψ i ðrÞ=AÞ ∑ N j≠i ∂Uðjr ij jÞ=∂r ij , where U(r) is the interparticle potential (see below), r ij = r i − r j , and A = A i = ∫ drψ i (r) is the area of each particle. f T (r) is the sum of the force density due to the torque applied on particle i: where α i is the rotational force density applied to particle i (α i > 0 and < 0 imply counter-clockwise and clockwise particle rotation, respectively), and e i θ is the unit angular vector perpendicular to r − r i orienting to the counter-clockwise direction. The rotational force density amplitude α is the same for all particles: |α i | = α. The torque acting on particle i, T i , is obtained as Thus, the torque strength T = |T i | is given by T = α∫ dr|r − r i | 2 ψ i (r).
The activity of particle-based active/driven matter can be characterised by the particle-scale Reynolds number Re=a 2 ρΩ/η ℓ . Since Ω is a reactive quantity to the torque T and depends on each particle's environment, we define Re by the rotational force density α as follows. Since the torque strength should be equal to the rotational frictional force for a rotating disk 92 , we obtain T ¼ 4πη ' a 2 eff Ωλ R , where λ R is the correction factor due to the smooth interface and given by λ R = (∫ dr|r − r i | 2 ψ i (r))/(∫ dr|r − r i | 2 ψ i (r) 2 ) = 1.182 for our setting of ψ i (r). We confirm the above T − Ω relation by simulations and obtain the effective disk radius a eff = 1.08a. Thus, we have the following relation between Re and α for our system: Re = 0.0064/α. We use Re defined in this way as a control parameter.
In our FPD method, the particle rigidity is approximately expressed by introducing the smooth viscosity profile, η(r). Thus, the approximation is better for a larger viscosity ratio η c /η ℓ and a smaller ξ/a. By multiplying both sides of Eq. (3) by ψ i (r) and then performing its spatial integration, we can straightforwardly obtain an approximate equation of motion of particle i: M i dV i /dt = F i + K i , where M i = ρA = M and V i = ∫ drvψ i /A are the mass and the average velocity of particle i, respectively. This equation of motion is approximate since V i is the integrated velocity of the already updated grid velocity of particle i at t = t + Δt (not at time t). On the right hand side, F i = ∫ drψ i (f U + f T ) is the force arising from the interparticle interaction and the torque, and K i ¼ À R drψ i ∇ Á Σ ffi À R dS ini Á Σ is the force exerted by the fluid. Here we use the following approximate relation R dr∇ψ i Á Q ffi À R S in i dS i Á Q for an arbitrary tensor QðrÞ, where S i is the surface of particle i andn i is the unit outward normal vector to S i . In practical numerical calculations, the on-lattice velocity field, v(r, t + Δt), is evaluated from the physical quantities at time t by Eq. (3). Then we move particle i off-lattice as a rigid body by r i (t + Δt) = r i (t) + ΔtV i (t + Δt), where Δt is the time increment of the numerical integration. We note that our method preserves both translational and angular momentum (see below). The periodic boundary condition also imposes a zero mean vorticity.
In our simulation, the units of length ℓ and time τ are related as τ = ℓ 2 /(η ℓ /ρ), which sets both the scaled density and viscosity of the fluid region to unity. This τ is the time required for the fluid momentum to diffuse over a lattice size ℓ. The units of stress and energy are σ ¼ ρð'=τÞ 2 and ϵ ¼ σ' 3 , respectively. Furthermore, we set η c /η ℓ = 50, Δt = 0.0025, ℓ = 0.5, ξ = 0.5, and a = 6.4. This time step ensures the numerical stability for the viscosity ratio η c /η ℓ employed. The simulation box used typically was L 2 = 512 2 . We will use the same characters for the scaled variables to avoid cumbersome expressions. We solve the equation of motion [Eq.
(3)] by the Marker-and-Cell (MAC) method with a staggered lattice under the periodic boundary condition. The total number of particles is fixed to N = 320 besides the system-size dependence and the area-fraction dependence studies.
Validity of the FPD method. Here, we briefly mention the validity of our FPD method. For example, we confirmed that the FPD method reproduces the Stokes drag 79 . Although we study an athermal zero-temperature system in this study, we mention that the FPD method combined with fluctuating hydrodynamics has thermodynamic consistency 80,81 . Thus, we have confirmed the translational and rotational diffusion behaviours under thermal noise 80 . We also showed through the long-time tail analysis that our method satisfies the translational and rotational momentum conservations 80 . We also confirmed the validity of our method through the direct comparison of colloidal phase separation between confocal microscopy experiments and FPD simulations without arbitrary parameters 81 . Fujitani also mathematically proved that the FPD method, in the limit of the infinite colloid viscosity and the infinitely thin interface, reproduces the known results of the hydrodynamic motion of a sphere 93 . In our method, the colloidsolvent viscosity ratio η c /η ℓ and the finite interfacial thickness ξ for a colloid with radius a are introduced as approximations (see Methods). However, η c /η ℓ = 50 and a/ξ ≥ 6 provide a reasonably good approximation for our purpose 21,79 .
Moreover, in our simulations of rotating disks, we apply the torque directly to each particle; thus, the non-ideality coming from the finite colloid viscosity is negligible. We have also confirmed that a single rotor has the characteristics of the Rankin vortex (the increase of the velocity field in proportion to r inside and its 1/r decay outside) 21 and that the motion of rotating disk pairs reproduces the known results 20,65 (see Supplementary Fig. 1). We also confirmed the linear relation between the applied torque T and the rotation angular frequency Ω of an FPD disk with radius a: T ¼ 4πη ' a 2 eff Ω 92 with the effective radius a eff = 1.08a 21 .
Analyses. We measured sixfold hexatic bond orientational order parameter to characterise local structural order: ψ j 6 ¼ j 1 n j ∑ n j m¼1 e i6θ j m j, where n j is the number of nearest neighbours of particle j, i ¼ ffiffiffiffiffiffi ffi À1 p , and θ j m is the angle between (r m − r j ) and the x-axis, where particle m is a neighbour of particle j. Note that ψ j 6 ¼ 1 means the perfect hexagonal arrangement of six nearest-neighbour particles around particle j and ψ j 6 ¼ 0 means a random arrangement. The structure factor S(k, t) at time t is calculated by taking the circular average of S(k): where C k (t) is the Fourier transform of Cðr; tÞ ¼ ∑ i Θ 4 πσ 2 ftanh½a À r À r i ðtÞ =ξ þ 1=2g with r i being the centre of mass of particle i, Δk = 2π/L, and Θ = 1 for CCW and Θ = − 1 for CW. The parameter Θ is introduced to differentiate the two domains with different chiralities. This is similar to calculating concentration fluctuations around the average (distinguishing + and − domains) in ordinary phase separation 72 , a quantity typically used to study phase separation. We note that the structure factor of the particle density fluctuation without distinguishing chirality cannot detect phase separation for SM. The kinetic energy spectrum is obtained by taking the angular average of the velocity power spectrum 87 : where v(k) is the Fourier representation of the velocity field. The nonlinear energy flux Π E (k) across a sphere of radius k in Fourier space is calculated as 87 Here v <k is the velocity field v(r) low-pass filtered so that all wave numbers outside the sphere of radius k are set to zero:

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The algorithms and simulation codes used in this study are available from the corresponding author upon reasonable request.