Emergence of lanes and turbulent-like motion in active spinner fluid

Assemblies of self-rotating particles are gaining interest as a novel realization of active matter with unique collective behaviors such as edge currents and non-trivial dynamic states. Here, we develop a continuum model for a system of fluid-embedded spinners by coarse-graining the equations of motion of the discrete particles. We apply the model to explore mixtures of clockwise and counterclockwise rotating spinners. We find that the dynamics is sensitive to fluid inertia; in the inertialess system, after transient turbulent-like motion the spinners segregate and form steady traffic lanes. At small but finite Reynolds number instead, the turbulent-like motion persists and the system exhibits a chirality breaking transition leading to a single rotation sense state. Our results shed light on the dynamic behavior of non-equilibrium materials exemplified by active spinners. Emergent collective behaviour has recently been addressed in systems of self-rotating particles, where motion, in particular, is an emergent phenomenon rather than a basic ingredient. Here, the authors derive a continuum model for mixtures of clockwise and counterclockwise Quincke spinners, demonstrating the emergence of same-spin phase separation, traffic lanes, sustained turbulent-like motion, and a chirality breaking transition depending on the fluid inertia of the system.

Simulations based on agent-based models, without hydrodynamics, for gear-shaped spinners 36,63 and rotating dimers 64 showed that a mixture of clockwise (CW) and counterclockwise (CCW) rotors segregate into same-spin phases. This behavior was actualized in an experiment 65 using a large number of small robots spinning on a table. The collective dynamics of fluidsuspended rotors, however, can be quite different from that observed in the dry systems due to hydrodynamic interactions. For example, unlike the frozen state of gear-shaped spinners at low density 36 , fluid-embedded spinners exhibit a gas-like phase, where the particles move randomly in the stirred fluid 38 . Furthermore, while both the dry and the fluid systems form similar large-scale dynamic patterns, such as lanes and vortices, their quantitative characteristics are different. Much insight has been gained by studying these systems using discrete-spinner models, however, this approach becomes computationally expensive when considering large systems with a high spinner density. In such cases, it becomes favorable to utilize continuum models that replace the discrete individual particles with a continuous distribution. Such continuum approach allows for the computationally feasible study of the collective dynamics of large systems 27,[66][67][68][69] .
The existing continuum models of self-rotating systems are phenomenological and typically based on the simplified hydrodynamic theory that includes either chiral components in the stress tensor or additional chiral forcing terms in the Navier-Stokes equation 45,46,70 . These models capture phenomena such as the surface waves in active spinner systems 46 , active coarsening, and the emergence of self-propelled vortex doublets 70 , etc. However, the phenomenological implementation of active rotation makes it challenging to assess if the observed dynamics are physical or dependent on the postulated continuum formulation. It also becomes difficult to directly compare experiment and theory.
Here we derive a continuum model for a monolayer of spherical spinners with CW-and CCW spins suspended in liquid in a three-dimensional domain. The model is motivated by experimental realizations such as spinners trapped at a fluid interface 39,71 or on a solid substrate 46 . We present numerical simulations assuming rotation due to an external DC (direct current) electric field, known as Quincke rotation. We explore the mixtures of CW and CCW spinners and study same-spin phase separation. We find that the dynamics is sensitive to fluid inertia: while in the inertialess system after a long transient turbulent-like motion the spinners settle and form steady traffic-like lanes, at small but finite Reynolds number the turbulent-like motion persists. For even larger Reynolds numbers, a population initially composed of an equal amount of CW/CCW spinners exhibits a chirality-breaking transition and collapses into a single rotation sense state. While our results are obtained for Quincke spinners, a similar continuum approach can be applied to model a monolayer of magnetic spinners energized by an applied AC magnetic field 39,71 . Our results advance a fundamental understanding of active systems, where large-scale chiral states emerge as a result of spontaneous symmetry breaking, as for example spinners, rotating bacteria, microtubules assays, etc.

Results and discussion
Summary of the continuum model. In the "Methods" section, we derive a continuum model for a two-dimensional system of N active spinners of mass m and radius a suspended in a fluid with viscosity μ and density ρ confined to a plane of area A, see Fig. 1a. The spinner concentration, c = N/A, is assumed to be uniform. The system is governed by the following equations for the spinner angular velocity ω(r, t), fluid velocity u(r, t), and pressure p(r, t): where I p ¼ 2 5 ma 2 and D ¼ k B T 6πμa . The thermal energy of the system is given by k B T. Details of the model derivation are provided in the "Methods" section. Fig. 1 Schematics of the discrete system. a Illustration of spinners confined to a plane. Blue and yellow color corresponds to clockwise (CW) and counterclockwise (CCW) rotation. The axis of rotation for each spinner is shown as gray vertical line oriented in the z-direction. b The Quincke effect is an electrohydrodynamic instability that gives rise to a torque on a dielectric particle in a uniform DC (direct current) electric field 94 if the freecharge polarization is antiparallel to the applied field. Above a critical field strength E > E Q rotation about an axis in the plane perpendicular to the electric field (Ω ⋅ E = 0, where Ω is the individual spinner angular velocity) is generated by the misaligned induced dipole of the particle, P. The direction of rotation can be either CW or CCW, depending on the initial perturbation in the dipole direction.
Equation (1) is the spin angular momentum balance of the spinners; the three right-hand side terms account for the thermal noise, the torque on the spinners due to fluid drag, and the externally applied torque τ. Equation (2) is the linear momentum conservation for the fluid (the Navier-Stokes equations) with the addition of the last term on the right-hand side, which drives the fluid flow due to the rotation of the spinners. Equation (3) is the incompressibility condition since we assume constant fluid density ρ.
It is worth noting that Eqs. (1)-(3) are generic for any spherical spinners. The specifics of the particular physical system are introduced in the form of the externally applied torque τ, which is what causes the particles to spin. In this paper we consider the case of Quincke spinners, studied experimentally in ref. 72 , τ = c −1 (P×E) z , where P is the polarization vector of the individual spinner and E is the externally applied constant (DC) electric field, see Fig. 1b. The governing equations for the polarization vector P are provided in the "Methods" section, see Eqs. (28)- (29).
The complete model is provided in non-dimensional form in the Methods section, see Eqs. (31)- (36). All equations are nondimensionalized using the spinner radius a as the characteristic length scale and the Maxwell-Wagner polarization time t Q , which sets the magnitude of the rotation rate, as the characteristic time scale. The scaling of the physical variables and the definitions of the model parameters are summarized in Table 1. The scaling analysis shows that the most important parameter governing the behavior of the system is the Reynolds number for the fluid, Re ¼ ρa 2 =μt Q . For Re ( 1, the effect of inertial terms in Eq. (2) of is small; fluid inertia becomes important for Re ) 1. In the rest of the paper, we explore the system behavior as  show that the spinner fluid undergoes a phase separation into clusters of CCW (yellow) and CW spinners (blue), in agreement with the simulations using a discrete spinner model 38 . These clusters themselves are rotating, growing, separating, and reconnecting as also observed in the active coarsening continuum model in ref. 70 .
The flow is localized the interface between CW and CCW clusters, see Fig. 2b and Supplementary Movie 2. Physically, the fluid flow is canceled out in the region between the same spin spinners. Hence, even though the spinners might rotate very fast, the flow might be very weak, especially in the interior of the same spin regions; accordingly, while inertia might be very important in the angular conservation equation, it might be less so in the linear momentum conservation (Navier-Stokes equation for the fluid flow). At the boundary between clusters, opposite spin spinners pump fluid in the same direction causing bands of flow or edge currents, see an illustration of this in Fig. 2d.
Emergence of lanes depends on Reynolds number. An important question concerns the long-term behavior of a system composed of spinners that are initially randomly 50 CW : 50 CCW. Previous theoretical studies 38,41 of discrete micro-spinners have shown that the system exhibits separation into same-spin clusters that eventually evolve into emergent patterns, such as lanes or vortices. This work however was restricted to Stokes flow (Re ¼ 0), while recent colloidal experiments 39,44 suggest a nonnegligible Reynolds number.
Using the continuum model, if we consider Re ¼ 0 and an initial 50 CW : 50 CCW random configuration, lanes always eventually form, see Fig. 3a and Supplementary Movie 3, regardless of the initial configuration. Therefore, the continuum model is in agreement with the discrete spinner model. The formation of lanes can be qualitatively interpreted as a result of the fact that the Stokes flow corresponds to least dissipation. Segregation into the same spin regions is thus favored since the flow vanishes in the interior of the same spin region and the only remaining flow is confined to the interface between opposite spin regions. The flow (and thus the dissipation) is reduced if the interface between these same spin domains is minimized, which in the case of the periodic boundary conditions considered in our simulations, is a straight line. Conservation of mass further requires two interfaces with opposite flow directions resulting in a lane configuration.
A small, non-zero Reynolds number, such as Re ¼ 0:01, leads to distorted lanes, see Fig. 3b and Supplementary Movie 4. These quasi-stable lanes have turbulent-like flow near the interfaces between the two opposite-spin clusters. The lanes can also spontaneously dissolve into to a transient turbulent-like state and then reappear.
For small but finite Reynolds number, e.g., Re ¼ 0:1 and greater, the system only exhibits a turbulent-like motion with no emergent collective large-scale structure, as shown in Fig. 3c and Supplementary Movie 1. If the initial state is lanes, instead of a random configuration, the lanes become unstable and dissolve into turbulent-like motion, see Supplementary Movie 5.
For an even larger Reynolds number, such as Re ¼ 1, the ratio of CW to CCW spinners is no longer conserved and the system exhibits a chirality-breaking transition. As the inertial terms of the Navier-Stokes equation become more significant, the fluid vorticity is not solely determined by the spinner angular velocity distribution. Instead, the inertial-induced component of the fluid vorticity can cause the spinners to flip spin orientation. For this case, the fluid flow initially exhibits turbulent-like motion, but eventually begins to favor a single spin orientation, see Supplementary Movie 6. Clearly this behavior is not present in dry system of spinners.
We can quantitatively distinguish between the lane and the turbulent-like cases by studying the probability distribution function (PDF) of the fluid velocity, u. As shown in Fig. 4a, the PDF for lanes exhibits a clear bi-modal distribution due to the two sides of the lane. For the turbulent-like case, Fig. 4b, the PDF is Gaussian unlike the non-Gaussian behavior observed in experiments with colloidal spinners 71,73 . Overly populated tails in our simulation were found to be due to underresolved simulations. Velocity and vorticity correlation functions for turbulent-like regimes are displayed in Fig. 4c, d.
Finally, it is worth mentioning that the discrete spinner model 38,41 also observed the formation of vortices in the case of larger spinner concentrations, which is not captured by our continuum model. This deviation is likely because steric repulsion Table 1 The dimensionless parameters and their definitions in terms of the properties of the spinners and the environment.

Parameter
Symbol Definition Spinner diffusion D * Dt Q /a 2 Inverse particle Reynolds number κ 8πμa 3 t Q /I p Normalized electric field γ E/E c Reynolds number Re ρa 2 /μt Q Area filling fraction α Nπa 2 /A Here a is the particle radius, t Q is the Maxwell-Wagner relaxation time t Q , D is the diffusion constant, E is the electric field, E c is the critical electric field, N is the number of spinners, A is the area covered by spinners, μ is the fluid dynamic viscosity, I p is the particle moment of inertia, ϵ is given by Eq. (30).
is neglected in the continuum model but becomes non-negligible at higher spinner concentrations. Steric repulsion can be included in the continuum model by introducing a particle stress tensor 27 .
Turbulent-like flow has a power-law energy spectrum. In many quasi-two-dimensional active systems exhibiting turbulent-like behavior, such as magnetic rollers 39 , or swimming bacteria 74 , the energy is injected at the microscopic scale. For example, rotation of bacterial flagella results in the onset of large-scale collective behavior and so-called "active turbulence" 6 , characterized by the inverse cascade 75 . It is tempting to assume that in the Quincke system, the inverse cascade should be present as well due to the microscopic rotation of individual spinners. However, the situation is less clear in our model because the microscopic energy injection scale is coarse-grained in the continuum approximation and replaced by the effective driving torque τ in Eq. (1). This torque may slowly vary in space and does not produce the inverse cascade.
As the simulations show turbulent-like behavior, it is of interest to compute the energy spectrum, E(k), and the energy flux in kspace, Π E . The energy spectrum E(k) is obtained by taking the Fourier transform of the spatial auto-correlation of u: such that r 0 ¼ x 0x þ y 0ŷ and 〈. . . 〉 |k|=k is the average over all wavenumbers k such that |k| = k. Correspondingly, the energy flux is obtained as follows Here u <k is the low-pass filtered velocity field with the wave numbers outside |k| < k set to zero: andũðk; tÞ is the Fourier transform of u(r, t) 76 . A hallmark of the inverse cascade in 2D turbulence is the energy scaling E(k)~k −5/3 and negative value of the energy flux Π E (k) for small wavenumbers k. As shown in Fig. 5, the energy cascade scales as k −5/3 . However, the energy flux is negative only in a very narrow range, Π E (k) < 0 for k < k 0 ≈ 0.016, which is not consistent with the energy cascade observed in 2D classical turbulence [76][77][78] and colloidal experiments 44,71 . We suspect that this behavior is due to the relatively small system size in the computations and coarse-gaining of the microscopic energy infection scale. In the cases where lanes form at long times, the fluid flow initially exhibits transient turbulence-like motion with the same k −5/3 scaling, prior to settling into lanes.
The turbulence inertial range is characterized by two lengthscales, the Taylor microscale 76 , and the integral scale 79 . The Taylor microscale falls in between the large-scale eddies and the small-scale eddies. The integral scale L I can be extracted from the velocity autocorrelation function in Fig. 4c. It gives L I ≈ 100 for particle Reynolds number Re = 0.1. The corresponding r.m.s. (root mean square) velocity u rms can be determined from Fig. 4b and is about 5, resulting in the integral scale Reynolds number of the order of 50. The Taylor microscale λ is several times smaller, and is determined from the fluctuations of velocity and velocity gradients, ð∂ r huiÞ 2 ¼ u 2 rms =λ 2 . It roughly can be estimated from the velocity gradients (vorticity) autocorrelation function, Fig. 4d, which gives λ ≈ 15.
Initially, the spinners begin to phase separate into larger and larger same spin clusters. However, the energy spectrum exhibits a peak that occurs at the length scale of the domain,~L. This demonstrates the cluster size becomes limited by the size of the domain. The spatial correlation functions for the fluid flow 73 and fluid vorticity are another useful metric for studying turbulent flows. As shown in Fig. 4, the turbulent flow case exhibits no characteristic length scale as seen in meso-scale turbulence. This behavior is consistent with colloidal spinner experiments 71,73 .
Activity enhances transport. The capability of the spinner system for active transport can be quantified by tracking the position of point particles s i that are being passively advected by the fluid flow, ∂s i /∂t = u(s i , t) and computing the mean-square-displacement (MSD), js i ðt þ t 0 Þ À s i ðt 0 Þj 2 . In addition, the mean square velocity (MSV) for the fluid flow, juj 2 , is another useful metric for the studying the active transport of the system.
As shown in Fig. 6a, the MSD initially exhibits a transient ballistic regime for all shown Re values. Counterintuitively, the MSD and MSV decrease as Re increases, Fig. 6b. For increasing Re, the spinner-driven flow begins to compete with the effects from the inertial terms. Also, the force term in Eq. (2), accounting for the spinners disturbing the fluid, is proportional to viscosity. Therefore, the disturbance in the fluid generated by the spinners does not get damped down with higher fluid viscosity.
The curves for different Re values begin to qualitatively diverge in the long term. The lane formation cases, Re ¼ 0 and Re ¼ 0:01, show sustained ballistic motion in MSD (~t 2 ) and a relatively constant MSV, as all the motion is confined to the interface of the lanes. Yet for sustained turbulent-like behavior, Re ¼ 0:1 and Re ¼ 1, the MSD shifts away from a ballistic regime closer towards a diffusive regime (~t). However, a reliable evaluation of the diffusion coefficient would require a much longer integration time. For Re = 1, the MSV also exhibits an eventual decrease. This is due to the chirality breaking of the spinner orientation.

Conclusions
To summarize, we develop a continuum coarse-grained model for a suspension of discrete microspinners. Assuming Quincke rotation, we then present numerical simulations that demonstrate the emergence of lanes of same-spin fluid and turbulent-like behavior depending on the Reynolds number. Our work sheds light on active spinner materials and makes testable predictions for experiments on the effects of fluid inertia and chirality symmetry breaking.
For the case of the Stokes flow, the formation of lanes of spinners is favored as it minimizes the interfaces between counter-spin clusters, thus minimizing the energy of the system. This behavior is in agreement with the existing discrete spinner models 36,38 . Yet, even small inertia causes turbulent-like behavior, which becomes more pronounced with increasing Re: turbulent lanes for Re $ 0:01 and sustained turbulent-like flow for Re ≥ 0:1. We also quantitatively characterize lane formation and the turbulent-like flow through the mean square displacement, MSV, and the energy spectrum statistics. Furthermore, for non-zero Reynolds numbers we observe that initially 50:50 CW/CCW population of spinners exhibit a chirality symmetry breaking transition and consequent condensation to a single sense of rotation state.
While this work used periodic boundaries, future work could apply this model with both fixed, impenetrable boundaries and soft, deformable boundaries, as well as arrays of artificial obstacles. Another possible avenue is the case of non-uniform spinner densities. As the spinner densities, in experiments, are often nonconstant in space and time, it would be of interest to see how this would affect the collective and turbulent behavior seen in this study. The model can be modified to analyze other experimental systems, e.g., magnetic spinners at water-air interface 25,39 and could also be easily generalized for a three-dimensional suspension of neutrally buoyant micro-rotors.

Methods
Derivation of the continuum model. Let us consider N spherical spinners of radius a and mass m confined to a 2D plane of area A. The translational motion of a spinner is the result of being passively advected by the surrounding fluid flow, u, and Brownian diffusion. The translational velocity of an individual spinner located at position r is then: The thermal noise of the system is modeled by a random variable ξ(t) that is Wiener process governed by the standard normal distribution. In the low-Reynolds number limit, the diffusion constant is D = k B T/6πμa, where k B T is the thermal energy and μ is the fluid viscosity. We also assume low spinner concentration such that the effect of steric repulsion is negligible. We introduce a probability density, Ψ(r, Ω, t), of finding a spinner at position r and with angular velocity Ω. Since there are N total spinners, the normalization for Ψ is given as Using Eqs. (7)- (14) and the Fokker-Planck equation 80 , it can be shown that the governing equation for Ψ is given by the following Smoluchowski equation assuming an incompressible fluid (∇ ⋅ u).
From this point on the derivation of our continuum model follows ref. 67 . It is useful to define a function h:i Ω such that hXi Ω ðr; tÞ ¼ Z 1
Therefore Eq. (12) gives the evolution equation for the spinner density of the system. It is worth noting that, with appropriate boundary conditions, c is conserved and ∫ A cdA = N according to Eq. (8).
Spin angular momentum equation. If we multiply Eq. (9) by Ω and integrate with respect to Ω, it is straightforward to show that We can again ignore the boundary term as spinners with |Ω| → ∞ are non-physical and have zero probability. The conservation of spin angular momentum of a single spinner is I p _ Ωðr; tÞ ¼ ζ p 1 2 ð∇ uÞ z ðr; tÞ À Ωðr; tÞ þ τ e ðΩ; tÞ ð 14Þ where I p ¼ 2 5 ma 2 is moment of inertia and ζ p = 8πμa 3 is the friction coefficient for a sphere and the externally applied torque on the spinner is τ e (Ω, t). The effect of thermal noise is ignored in Eq. (14) because the spinners rotation rate is very high and thus the rotational Peclet number Ω/D r ≫ 1, where D r = k B T/(8πμa 3 ).
Using Eq. (14), Eq. (13) simplifies to Furthermore, if we assume a constant spinner density c and define ω ¼ ω=c,  Computations performed with L = 300a to increase resolution, L is the integration domain size, a is the particle radius.
The above derivation is based on the assumption that a spinner is not selfpropelled. In other words, uðr 0 ; tÞ is dependent only on spinners located at r≠r 0 .
Fluid flow. Assuming Stokes flow, the fluid flow generated by a single spinner can be approximated as a point rotlet 81 , with torque 8πμa 3 Ωẑ. The resulting fluid flow is given by the solution to the following forced Stokes equation: For distribution of N spinners with center of mass position r 0i and angular velocity Ω i , the fluid flow can be written as The pressure distribution, p(r, z, t) can be computed similarly. By multiplying Eq. (17) by Ψ(r 0 ) and integrating, the Stokes equation for the spinner distribution can be shown to be For flows with non-negligible inertia, Eq. (20) can be approximated for low Reynolds flow as Since the majority of experiments and previous theoretical work involve an effectively two-dimensional system 25,39,60,61 , we will therefore consider a twodimensional fluid flow in the plane of the spinners. For the case of two-dimensional flow, a frictional drag term, −βu, can be added to the right side of Eq. (21). For small β values, this term does not qualitatively change the observed behaviors. The effect of larger β is more extensively studied in ref. 70  Spinner's torque. A common spinner system consists of ferromagnetic colloids in either a rotating or oscillating magnetic field 42 . The external torque is then given by where μ 0 permeability of vacuum, magnetization is M, and the magnetic field is H. An oscillating magnetic field allows for spontaneous-symmetry breaking in the direction of rotation. Therefore, spinners can have either a CW or CCW spin 44,85 . However, a rotating magnetic field predetermines the spinner's direction of rotation, causing the spinners to be either all CW or all CCW 46 . We, instead, consider a system of Quincke spinners driven by a DC electric field. Quincke rotation has the advantage of being well studied [86][87][88][89] and it allows for simultaneous populations of CW and CCW spinners.
Quincke rotation. For a fluid with conductivity and permittivity (σ 1 ,ϵ 1 ) and a suspended spherical particle with conductivity and permittivity (σ 2 ,ϵ 2 ), Quincke rotation occurs when and the external electric field |E| is above a threshold given by where is the Maxwell-Wagner relaxation time. The characteristic angular velocity for a Quincke spinner is Neglecting the higher-order electromagnetic interactions 89 between spinners, we instead focus on the spinners' hydrodynamic interactions. If R is the distance between two rotors, the hydrodynamic forces scale as~R −2 while the electromagnetic forces 86 scale as~R −4 . The external torque from the electric field is then given as with electric polarization density P(r, t) (dipole moment per unit area) governed by where the equilibrium polarization density is P eq is and Eq. (28)-Eq. (29) are provided in ref. 83 with the inclusion an artificial diffusion term for numerical stability.
Continuum model of a Quincke spinner fluid. We will assume a constant, uniform electric field E parallel to x-y plane. Without loss of generality, we can define E ¼ Ex. We nondimensionalize the system of equations in order to reduce the number of parameters. Equations (1)-(3) and Eqs. (26)- (29) can be nondimensionalized as Eqs. (31)-(36) using the scalings and dimensionless parameters shown in Table 1. X * is the corresponding dimensionless value of a variable X.
Re The system involves seven dimensionless parameters: the dimensionless spinner diffusion coefficient (D * ), the particle Reynolds number (κ −1 ), the strength of the electric field relative to critical field (γ), the Reynolds number (Re) for the fluid, and the percentage of area covered by spinners (α), and a small artificial diffusion constant (D Ã P ) for the polarization evolution equation. Since Quincke rotation only occurs for E > E c , we will operate for γ > 1. For the dimensionless form of the equations see Supplementary Note 1.
The parameters γ and α control how fast the spinners rotate and the corresponding effect on the fluid flow. D * determines the characteristic length for the interface between two opposite-spinning clusters. Therefore, the only parameters that could significantly change the qualitative behavior of the system are Re and κ.
In the initialization of ω * , each grid point is randomly assigned a small positive value (CW) or a small negative value (CCW).
Scaling and dimensionless parameters. We scale the length by the particle radius a, time by the t Q is the Maxwell-Wagner relaxation time t Q . The electric field E is normalized by the critical electric field E c . Here N is the number of spinners, A is the area covered by spinners, and μ is the fluid dynamic viscosity. Correspondingly, the polarization P is scaled by Nϵ/A, where N is the number of spinners, A is the area covered by spinners, and ϵ is given by Eq. (30). Finally, the pressure p is normalized by μ/t Q , where μ is the fluid dynamic viscosity. Definitions of the model dimensionless parameters are given in Table 1.
Numerical methods. The reported results are numerical simulations of the continuum model using a square domain of length L, double periodic boundary conditions. The code is highly parallelized using CUDA to run on a NVIDIA graphics card. Spatial derivatives are computed using a pseudo-spectral method 90 . All Fourier and inverse transforms are computed using the fast Fourier transform (FFT). Equations (31)-(34) are then time-stepped using first-order exponential time differences 91 . For Stokes flow, the fluid velocity is solved using the streamfunction formulation 92 . For non-Stokes flow, the fluid velocity is time-stepped using the method of Chorin Projection 93 . Details about the numerical implementation can be found in Supplementary Note 2.

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