Birotor hydrodynamic microswimmers: From single to collective behaviour

A microswimmer composed of two oppositely rotating strongly coupled colloids in solution is here termed as birotor and investigated by means of hydrodynamic simulations. The related flow fields, swimmer velocities, and rotational diffusion are controlled by the properties of the fluid, the swimmer geometry, rotation frequency, and also by the substrate friction. Resulting from mutual hydrodynamic and steric interactions, birotor pairs might follow one another, or more frequently rotate around each other. For larger number of interacting swimmers the continuous formation and dissolution of small and rotating aggregates dominates the collective dynamics. The birotors motion is hydrodynamically enhanced at short distances, such that the average velocity of the swimmers shows to increase with density for the investigated range of densities. This is compensated by a decrease of rotational diffusive time, making that the overall effective diffusion decreases with density. These results constitute the first systematic analysis of the birotor microswimmer, which could be also further modified as an easy to manipulate active particle for various potential applications.

or cancer treatment [7]. For the design of artificial motile agents, several realisations have been suggested in the literature, ranging from self-thermophoresis [8][9][10] and selfdiffusiophoresis [11] to actuation via electric or magnetic fields [12,13]. The actuation via electro-magnetic radiation is of special interest, since the external energy necessary for locomotion can be continuously provided and typically does not require high configurational precision to the application of the radiation. A particular important challenge for the design of motile active agents is to include the possibility of controlling the direction and speed of motion in order to reorient, slow down, or speed up the swimmer.
A subclass of active matter refers to chiral activity which is composed of elements rotating in a fluid as a response to externally applied, or internally generated torques [14][15][16]. Resulting from such torques, the surrounding fluid flow is perturbed, which importantly affects neighbouring particles leading to specific interparticle interactions and consequent collective behaviour [17,18]. In these systems, symmetry breaking is an important ingredient to turn the rotational motion of the constituents, i.e., rotors, into translational motion. In this (e) Flow velocity in the propulsion direction uy calculated in the co-moving frame, in the two axes indicated in (b). Only half of the symmetric axes is shown. Axis 1, with lines in grey, is the propulsion direction, axis 2, with lines in orange, is the direction perpendicular to propulsion. Darker solid lines correspond to simulations with γ = 6.9Ω0, and lighter ones to γ = 2.1Ω0. Dotted lines correspond to the friction-free system, and dash-dotted lines correspond to eq. (1).
letter, we investigate a microswimmer composed of two coupled counter-rotating colloids, termed here as birotor, where rotors mutually break the hydrodynamically induced flow field symmetry, leading to the overall propulsion of the microswimmer. We describe the swimming behaviour of isolated swimmers first, and then in a large ensemble of collectively swimming birotors. The working mechanism could be related to ciliated microswimmers, such as Paramecium, which use cilia to generate a flow on the swimmer surface in order to propel [2,19,20]. We employ a minimal two-dimensional model by means of mesoscopic hydrodynamics simulations, and report how to overcome the well-known Jeffery paradox [21]. We show that the mean-square displacement of birotors can be mapped on that on the active Brownian particle model, while steric and hydrodynamic interactions show to significantly enhance the motion in the collective regime.
Model. -The considered birotor is a microswimmer composed by two linked equal sized colloids, as indicated in fig. 1(a), both rotating with angular velocities of equal magnitude Ω and opposite directions, while being suspended in a hydrodynamic solvent. Resulting from the no-slip boundary condition on the surface of the colloids, each colloid causes a rotational flow field in the surrounding solvent, such that the pair acts as a pump constantly pushing solvent and propelling into the opposite direction as a consequence.
The solvent is simulated by a method known as multiparticle collision dynamics (MPC) [22,23] which considers N point particles of mass m whose dynamics are updated in alternating ballistic streaming and collision steps. For the collision, the fluid particles are sorted into cells of length a and all particles in the same cell at the same time exchange momentum according to the following scheme. The relative velocity of each fluid particle with respect to the cell centre of mass velocity is rotated either clockwise or anti-clockwise with equal probability by an angle of π/2. A correction term is added after each collision in order to enforce angular momentum conservation in the fluid at cell level [24][25][26]. In order to ensure constant temperature T in the fluid, a cell-level thermostat [27] is applied. The streaming and collision schemes fulfil the hydrodynamic conservation laws and hence correctly reproduce the Navier-Stokes equations [28]. Here, the employed units are m, a, and k B T , with k B the Boltzmann constant, then the employed average fluid density is ρ = 10 m/a 2 , and the time between two collisions h = 0.02 ma 2 /(k B T ), which results in a fluid viscosity of η = 18 mk B T /a 2 , and a Schmidt number of Sc = 150.
The rotors are modelled as two-dimensional impenetrable circular no-slip boundaries. The interactions between solvent particles and colloids are taken into account via a bounce-back rule [29], i.e., fluid particles that penetrated the colloid in the ballistic streaming step are moved back onto the surface of the colloid and exchange momentum with the colloid before moving back into the opposite direction. Consequently, the fluid velocity on the colloid's surface is on average approximately equal to the velocity of the surface. In order to reduce the remaining finite slip, phantom solvent particles in the inner of the colloids with velocities drawn from a Maxwell-Boltzmann distribution centred around the colloids surface velocity take part in the collisions [30]. The considered colloid diameter is σ = 6a and the default rotational frequency Ω 0 = 0.01 k B T /(ma 2 ), resulting in a Reynolds number of Re 10 −1 . The two rotors are coupled with a harmonic potential U(r) = k(r − l 0 ) 2 /2 with a strong harmonic constant k = 1800k B T /σ 2 , which makes that the centre-to-centre distance between rotors r is on average l 0 = 1.5σ and maintained almost constant. Simulations are performed in a squared box of size 400a for the single swimmer case and 1000a for the collective.
This is a flow decaying with the inverse of the distance, which suggests that two coupled rotors are advected by the mutually generated flow, such that a first estimation for the birotor propulsion velocity can be obtained from the azimuthal flow created by a single rotor at distance r = l 0 from the rotor's centre, namely v p = σ 2 Ω/(4l 0 ). However in 1922, Jeffery studied the steady Stokes flow caused by two equally sized disks rotating with equal angular velocities and opposite directions in a resting and infinitely large two-dimensional fluid [21]. He showed that no solution to the Stokes equation can simultaneously fulfil the no-slip boundary condition on the rotors surface and at infinity. This scenario is nowadays known as the Jeffery paradox and implies that the flow velocity is not completely damped, such that given two close rotors in a two-dimensional setup, a finite flow velocity will remain at infinity [31]. In realistic setups, the induced two-dimensional flows can partially escape into the third dimension enforcing the outer boundary condition, i.e., that the fluid rests at infinity. On the other hand, solvent and rotors can be considered to have a non-negligible friction with a substrate, which could have for example a non-negligible rugosity. In this case, dissipation of energy has a similar effect, and the fluid also rests at infinity, for both the two-and the three-dimensional systems.
Here we consider a two-dimensional model of the rotors with an effective substrate friction γ, and in this case, the Stokes equation relating the variation of fluid velocity field u and pressure p, takes the form of a Brinkmann equation known from low-Re flow in porous media, with an additional linear friction term −γu [32] The friction adds then an additional cut-off factor to eq. (1), with a range decreasing for increasing friction.
In order to consider the effect of a substrate friction in the simulations in an effective manner, we consider additional phantom substrate particles in the collision step with velocities drawn from a Maxwell-Boltzmann distribution with zero mean velocity, as described in detail in [33]. These phantom substrate particles are randomly distributed over the whole simulation domain in each time step and take part in the collision step, such that momentum is transferred to the substrate and hence the outer boundary condition is fulfilled far away from the birotor. Unless otherwise stated, we consider a phantom substrate particle density of ρ subs = 5×10 −3 m/a 2 , which has shown to result in γ = 6.9 Ω 0 [33].
Birotor flow fields. -The averaged solvent velocity flow fields induced by a birotor can be calculated from the simulation data, and from these the solvent streamlines, as shown in fig. 1(b), (c), (d). The propulsion component of the average solvent velocity, u y , is obtained in two perpendicular axes centred at the birotor centre of mass (as indicated in fig. 1(c)), and displayed in fig. 1(e) as a function of s, the distance to the birotor centre of mass. Note that axis 1 does not cross the colloids, and that, for axis 2, the outer colloidal surfaces are placed at s = ±(l 0 + σ)/2. Simulations with no substrate friction show in fig. 1(b) that the birotor acts as a fluidic pump moving all the surrounding flow in a uniform direction. The flow at the swimmers surface and its direct vicinity is dictated by the no-slip boundary, i.e., the flow given by eq. (1) at r = σ/2, the rotors' surface, or u s = σΩ/2. At larger distances, the flow reaches a constant plateau value of u y 0.25u s in both directions, as can be seen in the dotted lines of fig. 1(e). The lack of flow decay far away from the birotor agrees qualitatively to the Jeffery paradox. Given the finite size domains employed in simulations, a non-vanishing flow in the steady state reaches an artificial plateau irrespective of the system size.
In contrast, simulations considering the substrate friction shown in fig. 1(c) clearly show that the created flow decreases for long enough distances. Close to the rotors surface the fluid is pumped against the propulsion direction, and due to the flow continuity, two circular flow patterns are formed, which is reminiscent of the field lines of a source dipole. This constitutes a proof of concept for the birotor behaviour as a microswimmer. The decay of the flow velocity in fig. 1(e) shows to occur in both directions, by at least two orders of magnitude at a distance from the birotor of approximately 10σ, which also ensures that the flow decays to zero sufficiently far away from the birotor, and it also shows to be faster the larger the friction. For the direction perpendicular to the propulsion, the flow change of sign is due to the presence of the stagnation points which appears earlier for larger friction values. The decay in eq. (1) shows to be a reasonable estimation, which even agrees quantitatively with the decay in the direction perpendicular to the propulsion for the case with γ = 2.1Ω 0 . This value of γ ensures the decay of the flows far away, but has only little impact on the near field, which is dominated by the flow created by the nearest rotor, agreeing therefore with eq. (1). The velocity in the propulsion direction includes a small interval between the colloids where the flow goes in the propulsion direction, but due to continuity this quickly reverses, and after a short distance, the flow decays again.
In order to understand the effect of hydrodynamic interactions when two or more particles are close to each other it is important to describe also the flow fields in the lab frame. In fig. 1(d), we approximate this flow by subtracting the particle average velocity v p from the flow field of the fixed birotor in fig. 1(c). The resulting flow is reminiscent of that created by a 3D source dipole, also known as a sliplet, which is a minimal model for a forcefree microswimmer [8]. From the results in fig. 1(d) we can infer that the birotor has a front-repulsive and rearattractive interactions. Sideways close interacting neighbours are dragged against each other's motions, which hinders clustering, and simultaneously enhances the particle relative motion.
Single swimmer dynamics. -Besides the induced advective movement on the birotors, the solvent thermal fluctuations lead to a diffusive behaviour of the birotors centre of mass, and of the propulsion direction, this is to a diffusive spreading of the direction of propulsion. This combination gives rise to the typical motion of active Brownian particles, as can be seen in the trajectory of fig. 2(a). The birotor centre of mass time-normalised mean-square displacements (MSD) Δr 2 /t for systems with different values of the substrate friction γ are shown in fig. 2(b). Three regimes can be distinguished. At short times, the system is superdiffusive instead of diffusive, which is a well-known effect for diffusion in twodimensional solvents [34]. At intermediate times, the system becomes clearly ballistic since the birotor advection dominates the motion. At longer times, the timenormalized MSD reaches a plateau indicating the diffusive regime. The last points at very long times, indicate the end of the statistical accuracy. The MSD are larger at all times for lower values of the substrate friction, which is a first indication that the effective birotor velocity decreases with increasing friction. The MSD curves can furthermore be fitted with the prediction from the active Brownian particle model [11], as shown in fig. 2(b), which nicely agrees with the birotor behaviour. Moreover these fits provide an estimation of both the swimmer velocity v p and the rotation diffusive time τ R . Fits are performed for tΩ ∈ (0.7, 3000) which consider the ballistic and diffusive regimes, and with the value of the self-diffusion of a single rotating colloid D = 3.73 × 10 −4 σ 2 k B T /(ma 2 ) as estimated for the single rotating disks [14]. The resulting velocity v p is depicted in fig. 2(c) as a function of the input rotational frequency Ω, and the increase shows to be mainly linear, as expected from eq. (1). Furthermore, v p is expected to be inversely proportional to l 0 , since changing the distance between the coupled rotors changes the mutually experienced advection. Increasing friction diminishes the birotor velocity v p , as shown fig. 2(d), and slightly increase τ R as shown in the inset of fig. 2(d). The propulsion velocity v p estimated from eq. (1) considers the flow induced by just one colloid, and would only be valid in the limit of γ → 0. The measured values of v p > v p result from the additive contributions of each of the two rotors conforming the birotor. Also the substrate friction diminishes the created flows, which explains that with increasing γ, the actual propulsion speed v p decreases as shown in fig. 2(d). For the obtained values v p and τ R , the actual statistical inaccuracy of the data appears larger than the one provided by the error bars obtained from the fits, such that a more precise functional dependence is not possible. The data for the related Péclet number, Pe = v p τ R /σ, is shown in fig. 2(e) and shows that Pe increases with increasing with Ω, and although the dependence with γ is not clear, we assume that the trend is consistent with a constant or slight increase. Note that, for the default parameters used in this letter, the propulsion velocity is v p = 0.66v p , and the rotational diffusion of the birotor orientation is τ R Ω/(2π) = 18.65. This means that the influence of the substrate friction with γ = 6.9Ω on the created flows reduces the propulsion velocity to 66% of the predicted velocity, and the swimmers orientation is on average reoriented after approximately 20 rotor revolutions of the intrinsic rotation. Consequently, the birotor's rotational Péclet number is Pe rot 13, which says that the birotor travels roughly a length of 13σ before the velocity direction is reoriented.
Collective dynamics. -When two or more birotors interact, it is not only the direct steric encounters but also the short-and long-range hydrodynamic interaction which play a crucial role in determining the system dynamics [10]. Two representative snapshots are shown in figs. 3(a), (b) at low and medium densities, respectively, which can also be seen in the two corresponding Supplementary Movies movie1.mp4 and movie2.mp4 (SM). In both cases, multiple encounters are frequent, which makes the configurations clearly distinguishable from those of passive Brownian systems. There is also no evident stable clustering which makes the structures significantly different from those of phase-separated Brownian active systems [35,36]. For both densities, a number of single swimming birotors can be seen, and also a number of metastable birotor pair configurations which can be understood in terms of the flow field shown in fig. 1(d). When two birotors are close and perfectly aligned, as sketched in fig. 3(c), then the rear-attraction of the front birotor is compensated by the front-repulsion of the second birotor and they effectively form a channel and move together. When the alignment is not perfect though, the two birotors quickly rotate relative to each other by the torque induced by the flow, see fig. 3(d), such that the configuration is highly unstable. For birotor pairs with opposite propulsion directions, the steric interactions lead to a configuration where the four respective rotor positions lay on the corners of a parallelogram, as in fig. 3(e), where the thrust centres do not lay on the same line as the centre of mass of the two birotors, such that the pair rotates. These rotating pairs are more stable than the aligned ones, but typically thermal fluctuations of the surrounding solvent, or collisions with other particles, destabilize the arrangement after one to two rotations. Encounters of three to six birotors in the case of lower densities, or up to twenty for higher densities, occur frequently as can be seen in figs. 3(a), (b) (see the SM). These encounters are very dynamic, where propulsion, relative, and collective rotations, as described above for the pairs, can be observed for all particles interacting in smaller groups within the encounter, which are also frequently recombining. Some other particle collective arrangements like the flow-trail structure in fig. 3(f) might eventually also appear. This makes that these multi-particle encounters are short-lived, producing structures that are continuously grouping, dissolving, and quickly reforming again with other neighbouring birotors. Note that even at higher densities, we expect that stable aggregation will be largely avoided, although flow-trail ordered structures will become more frequent and stable eventually providing the system with a more turbulent-like dynamics.
The described dynamic behaviour can be partially quantified by the time-normalised mean-square displacement of the birotors ensembles Δr 2 /t, as shown in fig. 4(a). At low density, the isolated birotors dominate the collective dynamics, such that the mean-square displacement approaches that of the isolated single birotor, and progressively deviates for increasing density values. Two interesting effects occur as the density increases. At shorter times, the average quadratic displacement increases with density. This means that birotors become faster due to the more frequent interactions with other birotors. This is a purely hydrodynamic effect, since the isolated rotors motion get enhanced by the influence of the flow induced by all neighbouring birotors. This differs from many systems of self-propelled particles where motion is mostly hampered by the interactions, and is related to similar observations for non-coupled rotors where hydrodynamic interactions of rotors lead to translational actuation [14]. The enhancement of the velocity is expected to change trend for very large densities where steric interactions would essentially impede the birotors motion. A second effect is that steric and hydrodynamic interactions change constantly the propulsion direction as described above, leading to a decrease of τ R . This can be already observed by the shift to earlier times of crossover to diffusive behaviour of the MSD with increasing density. Finally, the overall effective diffusion behaviour, shown as the plateau at longer times for the time-normalized MSD in fig. 4(a), decreases with density. Although the overall dynamic behaviour of the system cannot be explained with the active Brownian particle model, the MSD curves can still be fitted by prediction from the active Brownian particle model [11], as shown by the lines in fig. 4(a). The fits are performed with the same considerations as for the single swimmer case, and also provide an estimation for v p and τ R as an average over time and over all particles in the ensemble. The results in fig. 4(b) confirm the increase of average particle velocity v p which almost doubles in the investigated density range. Simultaneously the rotational diffusion time τ R in fig. 4(c) shows to decrease for more than a factor ten. Although these results could in principle result counter-intuitive, the combination of both effects result in the expected decrease of both, the effective Péclet number in fig. 4(d), and the enhanced effective diffusion in fig. 4(e).

Summary and conclusions. -
The dynamics of birotors, coupled pairs of counter-rotating colloids in a twodimensional low-Re flow induce an unphysical flow known as the Jeffery paradox, which we show can be overcome by the consideration of substrate friction. The intensity of the substrate friction quantitatively changes the created hydrodynamic flow field and decreases the propulsion velocity; however, the qualitative behaviour is independent of the value of the substrate friction. The flows mutually induced on the rotating surfaces advect the rotor pair forward, with a propulsion velocity linearly proportional to the rotors angular velocity and inversely proportional to the centre-to-centre distance between the rotors. The flows induced by the birotors are reminiscent of a sliplet flow and the mean-squared displacements can be effectively mapped to the active Brownian particle model. For large birotors ensembles, hydrodynamic and steric interactions among pairs or triplets of birotors lead to aligned swarming, or to anti-aligned rotating groups, together with the frequent formation and dissolution of aggregates. For increasing densities, birotors show to become clearly faster at short times, which is in contrast to the behaviour of many other microswimmers where the motion is mostly prevented by the interactions with other swimmers. The velocity enhancement of the birotor system is due to the frequent hydrodynamic encounters with other birotors, which is also related to the translation velocity induced in systems of synchronously rotating colloids [14]. However, systems of birotors and of rotors are significantly different, mainly because birotors also behave as swimmers in very dilute conditions and do not form large rotating clusters in concentrated situations. Note that hydrodynamic birotors are intrinsically different from other existing microswimmers such at self-phoretic colloids, flagellated bacteria, or squimmer models.
The experimental construction of these birotors might be relatively involved, but possible in practice for example by the construction of mechanically rotating components, by the use of phoretic multicompounds [37], or by building the birotors out of one particle carrying a magnetic dipole and another one an electric dipole, and applying two rotating fields one magnetic and one electric. In this case, not only could the speed of the birotor be tuned by the rotation frequency, but also the direction of propulsion could be varied by modifying the rotational frequency of the two fields independently. The birotor could then alternate from an intrinsic rotation motion into a straight propelling one. The manoeuvrability of such structure would certainly require specific detail investigation varying for example the colloid sizes, or also dynamically modifying each of the angular velocities. This promising working mechanism might trigger the development of micromachines with interactively controlled trajectories, which would provide them with promising potential applications [6,[38][39][40]. Machine learning algorithms might help then for the selfregularisation of such systems, in order for example to bring one of such birotors directly to a target [41][42][43]. * * * The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC) and the Helmholtz Data Federation (HDF) for funding this work by providing services and computing time on the HDF Cloud cluster at the Jülich Supercomputing Centre (JSC).

Data availability statement:
The data that support the findings of this study are available upon reasonable request from the authors.