Hydrodynamic Choreographies of Microswimmers

We unveil orbital topologies of two nearby swimming microorganisms using an artificial microswimmer, called Quadroar. Depending on the initial conditions of the microswimmers, we find diverse families of attractors including dynamical equilibria, bound orbits, braids, and pursuit–evasion games. We also observe a hydrodynamic slingshot effect: a system of two hydrodynamically interacting swimmers moving along braids can advance in space faster than non-interacting swimmers that have the same actuation parameters and initial conditions as the interacting ones. Our findings suggest the existence of complex collective behaviors of microswimmers, from equilibrium to rapidly streaming states.

submillimeter-scale Quadroar also a suitable candidate for various biological applications 35 such as drug delivery or autonomous surgery 34 . In macro scales, the Quadroar can be hired as a robotic swimmer 36,37 for inspection missions in highly viscous fluid reservoirs 38 .
In this study, we use the Quadroar swimmer and show that two model micro-swimmers in the Stokes regimethat generate flow fields resembling that of C. reinhardtii-have a rich two-body dynamics. Unlike other existing theoretical models that try to simulate the swimming mechanism of specific microorganisms 39 , the Quadroar is designed to induce an oscillatory flow field with anterior, side and posterior vortices in its surrounding 30 . Therefore, complex interactions that we find in the phase-space of two swimmers are generic characteristics of microorganisms generating anterior, side and posterior vortices.

Kinematics and Numerical Framework
The Quadroar consists of an I-shaped frame including an active chassis and two axles of length 2b (see Fig. 1). Each axle at its two ends is connected to two disks of radii a. The length of the chassis is variable and is equal to 2l + 2s(t) where s(t) is the contribution from the expansion/contraction of a linear actuator installed in the middle of the chassis. The angular position of each disk D n (n = 1, …, 4) with respect to the leg of its axle is denoted by − π ≤ ϑ n ≤ π. We define a body-fixed Cartesian coordinate system (x 1 , x 2 , x 3 ) with its origin at the geometrical center of the frame. The (x 1 , x 2 )-plane lies in the plane of the swimmer and the x 1 -axis is along the chassis. We also define a global Cartesian coordinate system (X 1 , X 2 , X 3 ) as is shown in Fig. 1(b). The body-fixed coordinates x i are related to the global coordinates X i (i = 1, 2, 3) through a transformation matrix R that depends on three orientation (Euler's) angles of the swimmer. We assume that the influence of each disk on its surrounding environment can be modeled as a point force and a point torque 40 .
For each of the two swimmers, j = A, B, and for each of their four disks n = 1, …, 4, the point forces (f jn ) and torques (τ jn ), expressed in the global coordinate frame are given by 30 jn j j n j n ,body where μ denotes dynamic viscosity of the surrounding fluid; v jn and ω jn are the linear and angular velocities of each disk with respect to the swimmer's hydrodynamic center and body-fixed coordinate frame, respectively; v j,c is the absolute velocity of the hydrodynamic center and ω j,body is the angular velocity of the jth swimmer expressed in terms of Euler's angles; u jn and 2Ω jn are the velocity and vorticity fields of the fluid at the center of each disk; G is the isotropic rotation tensor and K jn is the translation tensor corresponding to disk n of swimmer j, defined as 33 : Since a self-propelled swimmer in the Stokes regime is force-free and torque-free, we must have ∑ = = f 0 n jn 1 4 and also τ ∑ T  jn  jn  jn  1  4 for j = A, B. These four sets of vectorial equations (two vectorial equations for each swimmer) require the values of velocities u jn and spins Ω jn to be complete and solvable for v c and ω body . The linear nature of the Stokes equation allows us to invoke superposition and obtain 30,40 : where . The scalars z kn,j and z kn,ij are the magnitudes of the vectors X kn,j and X kn,ij , respectively, and r jn denotes the position vector of the nth disk in the swimmer's local coordinate frame. In all expressions we have i,j = A, B with the condition i ≠ j in each expression.
We assume that disks D n (n = 1, …, 4), of each swimmer j = A, B are spinning with angular velocities ω ϑ = ϑ =˙c where δω is a detuning parameter, and the length of the linear actuator at the middle of the chassis varies according to s(t) = s 0 [1 − cos(ω s t)]/2. Throughout our simulations we set a = 1, b/a = 4, l/a = 4, s 0 /l = 1/2, and ω s = 1. The characteristic time scale of the two-body system is T s = 2π/ω s . The parameter c 0 affects both the swimmer's dynamics and flow field around it. For c 0 ≈ 0.5, it has been shown 30 that the flow field induced by the Quadroar closely resembles that of C. reinhardtii alga 29 . Our numerical experiments show that the resemblance holds for almost any c 0 ≥ 1. This is similar to the recent experimental observations 41 that swimming speed or beat frequency do not have a considerable effect on bioconvection behavior of C. reinhardtii cells. The similarity in behaviors for this broad range of c 0 , which even includes the single-frequency case (i.e., c 0 = 1), adequately addresses the concern about whether any of emerging dynamical regimes is affected by the presence of two different frequencies. This further highlights the significance of having an oscillatory flow field. To speed up numerical simulations, we set c 0 = 50. For individual swimmers, nonzero values of δω significantly increase the number of orbital families, and in some cases lead to densely interwoven quasi-periodic rosette-shaped trajectories capable of inducing chaotic mixing in the surrounding environment 30 . Here, to focus on the basics of mutual interactions, we consider δω = 0. The results of this study are still valid for small δω, but start to deviate and become more involved as δω increases. Note that a = 1 μm leads to a Quadroar of ∼8-12 μm, which is similar to the size of a C. reinhardtii cell. Moreover, setting c 0 = 50 results in the frequency of 50 Hz for the disks, reminiscent of the flagella beat frequency for a C. reinhardtii cell 29 . For a more detailed description of the flow field around an isolated Quadroar the reader is referred to earlier publications on the Quadroar dynamics 30 .

Results
Two of our swimmers, depending on their relative initial locations (dX 1 , dX 3 ) portray a range of various trajectories as a result of their mutual hydrodynamic interactions. These trajectories range from converging, diverging, and oscillatory motions (which are also seen in other artificial microswimmers 23 ), to forming braids [ Fig. 2(b)], and even dynamical equilibria [Figs 2(e)and 3(b)] which, to the best of our knowledge, have never been observed in low-Reynolds-number swimming. We also report capture into bound orbits [ Fig. 3(c)] for two interacting microorganisms swimming in an infinite unbounded fluid. Interestingly, a similar behavior is observed in the lab for two Volvox colonies attracted by the chamber ceiling 42 .
In order to systematically study different possibilities of two-swimmer wiggling, induced by hydrodynamic interactions, we simultaneously consider the effects of swimming direction and relative initial locations, which also cover phase shift effects. Since there is no explicit time dependency in the Stokes equation, swimmers with an arbitrary phase shift between them (as a result of being launched at different times) can be considered as two swimmers with initial locations described at the moment that the second swimmer is turned on. For simplicity, we focus on the planar phase space. Nevertheless, our findings can be inherently generalized to 3D space. Our study has also been conventionally arranged in two general categories: (i) the two swimmers are released in the same direction such that their initial x 3 -axes are parallel and both aligned with the positive X 3 -axis (cf. Fig. 1), and (ii) the two swimmers are initially facing opposite directions such that at t = 0 the following conditions hold: , where the hat sign denotes unit vector. The resulting parameter space for each of these general cases is still valid for small perturbations. For larger perturbations, however, the parameter space starts to deviate from the presented plot and gradually tends to that of the other general case. For example, by changing the relative angle between the swimmers' initial x 3 -axes from zero to π, the corresponding parameter space diagram will gradually transform from Fig. 2a (swimming in the same direction) to Fig. 3a (swimming in opposite direction).
The parameter space for the trajectories of two swimmers released parallel and in the same direction [case (i)] is displayed in Fig. 2(a) with sample trajectories demonstrated in Fig. 2(b-e). In these figures, the swimmers would follow dashed lines in the absence of hydrodynamic interactions. If the two swimmers are released close to each other, and depending on their relative locations, they form a variety of braids with different shapes [ Fig. 2(b)]. Interestingly, we find that forward translational motion along a braid is faster, sometimes by a factor of two, than the motion of individual swimmers in the absence of hydrodynamic interactions. This phenomenon, which we refer to as hydrodynamic slingshot effect, can be easily deduced from Fig. 2(b): two hydrodynamically interacting swimmers advance along braids, and therefore in space (motion along colored lines), faster than non-interacting swimmers (moving along dashed lines) whose actuation parameters and initial conditions exactly match those of the interacting ones. The slingshot effect is caused by a synergistic process: each swimmer induces an advection field that sums with the relative velocity of its companion swimmer with respect to the  Figure 4 demonstrates a snapshot of the flow field and corresponding streamlines induced by the system of two swimmers advancing along a braid-like trajectory. It shows how each swimmer induces an advection field at the geometric center of the other swimmer, boosting its absolute velocity. The net flow field could also be described as a constructive interference of the two swimmers' flow fields (see Fig. 4). The resulting net flow field is similar to the one induced by a single Quadroar swimmer 33 but with a more powerful propellers (higher values of c 0 ).
Other families of trajectories that we observe for interacting swimmers belong to a general family of non-orbiting paths including diverging and converging trajectories [ Fig. 2(c)]. Non-orbiting paths may occur as pursuit-evasion games when one of the swimmers chases the other one [ Fig. 2(d)]. The most interesting non-orbiting path that we have found happens when the swimmers get in a reverse motion [ Fig. 2(e), colored dark blue in Fig. 2(a)], eventually reaching to a dynamical equilibrium. In dynamical equilibrium states, the swimmers' propellers are working continuously and their flow fields form a saddle structure (Fig. 5). The net flow of the saddle structure is zero, so follows the equilibrium state. In the space between the swimmers, fluid is pumped out in a direction almost parallel to the chassis of both swimmers, and is sucked back normal to the chassis. Four prominent vortices are formed around the propellers of the swimmers. These vortices are enclosed by a large-scale hyperbolic structure. Our long-term simulations show that dynamical equilibria are stable to small perturbations. This is a counter-intuitive property because the existence of hyperbolic structures usually implies local instability. The existence of dynamical equilibria for N > 2 swimmers is an unsolved problem, whose solution can sharpen our understanding of bacterial clustering and motile cell accumulations 43,44 .
If the initial distance of the swimmers is large enough, their hydrodynamic interaction will be very small, drifting the swimmers slightly off their straight trajectories. We have also observed a switch between different trajectories as the two-body system evolves: a converging motion may end up in an equilibrium state, or pursuit-evasion game may bifurcate to either of converging or diverging paths. In the parameter space, we have marked these cases with two or more colors [ Fig. 2(a)].
Two swimmers starting their motions in opposite directions [case (ii)] exhibit different orbital topologies from what we observed for co-directional ones. Their two-body dynamics depends on the impact parameter dX 1 [ Fig. 3(a)]. When the impact parameter is relatively small (  dX a / 5 1 ) and the swimmers initially move towards   Fig. 3(c)]. It is noted that capture into a quasi-periodic orbit is a transitional state between dynamical equilibria and deflecting trajectories. Such transitional states fill a complex fractal-shaped region of the parameter space, showing high degree of sensitivity to initial conditions [see the zoomed-in box in Fig. 3(a)] with the dominant length-scale of a disk radius. This result suggests the existence of highly chaotic N-body systems of swimming microorganisms. Although details of trajectories in an orbiting motion can be complex, the bounded nature of the overall two-body motion in an infinite fluid domain is a unique physical process, for which many applications can be sought. Examples include mixing by microswimmers and trapping microorganisms by artificial microswimmers.

Conclusions and Discussion
Here we have shown that two microswimmers in Stokes regime can stop each other by forming a dynamical equilibrium in an infinite fluid domain. Furthermore, depending on where the two swimmers are released, they may also get trapped into bounded orbits and revolve about each other indefinitely. We have systematically studied the entire phase space of a hydrodynamically interacting two-swimmer system, and identified the basins of dynamical equilibria and periodic orbits in the parameter space. We also found other diverse sets of orbits including closely winding braids, and pursuit-evasion dynamics. Sensitivity to initial conditions, slingshot effect for motions along braids, dynamical equilibria, and capture into bound orbits, as demonstrated in this study, can have unexpected implications to motion of microorganisms. Nonlocal models of passive and active stresses due to hydrodynamic and steric interactions 9 will then need modifications as diffusion in the phase space cannot be modeled only as a function of macroscopic streaming velocity.
cos sin sin cos cos Since there is a coordinate-type singularity in T, when θ = ±π/2, all computations have been carried out in the space of unit quaternions q, and then outputs are mapped back onto the space of Euler's angles α 33 : Three-dimensional Beads Model Simulation. In order to validate our models of the point forces and torques of the disks, we develop a full three-dimensional beads realization of the disks 45,46 . We first briefly explain the concept of beads model using Fig. 6a, then compare the results of our numerical method with those of beads model simulation.
A single spherical bead, moving with velocity v 0 , in Stokes regime induces a well-known velocity field in the surrounding fluid. For an arbitrary point in cylindrical coordinate system, this velocity field in radial and tangential directions is given by: where R 0 is the radius of the bead, and re r is the position vector of the arbitrary point with respect to the bead's center. Figure 6a demonstrates two beads B 1 and B 2 , which are moving with absolute velocities v 1 and v 2 in the stationary frame. With respect to the background fluid, B 1 (B 2 ) has a hydrodynamic velocity v 1H (v 2H ), and thus induces a velocity field v 1H,2 (v 2H,1 ) at the position of B 2 (B 1 ). So, the hydrodynamic velocity of each bead is given by the following implicit formula: Generalization of this simple idea, in order to formulate hydrodynamics of a system composed of N beads, results in the following system of linear algebraic equations:  = , the general hydrodynamic relations between beads is defined here as:  where r mn = r n − r m = (x, y, z). Applying this representation to the general formulation of the system of N beads presented in equation (8), leads to the following implicit system of linear equations: This system of equations can then be solved using standard linear algebra methods. The inputs of the system are absolute velocities, v i , which are assigned to individual beads that assemble a rigid body, and outputs are hydrodynamic velocities.
We now validate our numerical method of modeling hydrodynamic interactions using a 3D beads model realization of two nearby rotating disks, where each disk is composed of a large number of beads (see Fig. 6b). The optimum number of beads required to model each disk is determined through the convergence of results. For the presented case in Fig. 6b, as an example, the optimum number of beads is 331, which corresponds to R 0 /a = 1/21 ≈ 0.05. It should be noted that the thickness (2R 0 ) of each disk can be neglected compared to its diameter (2a), as expected by the swimmer model. Then, the general set of equations (14) must be solved for the entire system of the beads. Drag element exerted on each bead is then determined by multiplying the translational drag coefficient, Figure 7. Magnitudes of the total force (a) and torque (b) exerted on two interacting rotating disks of radius a as a function of their distance. The schematics of our setup has been shown in Fig. 6b. We have set ω = |ω 1 | = |ω 2 | = 0.5ω s . Stars and circles respectively correspond to the results of our theoretical (point force and torque) and 3D beads models.
Scientific REPORTS | (2018) 8:3670 | DOI:10.1038/s41598-018-21832-w 6πμR 0 , to the bead's consequent hydrodynamic velocity. At the end, the final results of this three-dimensional simulation is compared to our numerical results of the point force and torque models. Figure 6b demonstrates the schematics of the problem setup, where two disks of radius a and located at a distance of d are rotating with angular velocities ω 1 and ω 2 . The interaction of these two disks is modeled using: (i) our point force and torque models, and (ii) the full three-dimensional beads simulation. Our results displayed in Fig. 7 show a good agreement between the two models, even for the smallest possible distance, d = 2a, between the disks. Figure 7 represents the magnitudes of the total force and torque exerted on disks as a function of their distance. The disks are counter-rotating with |ω 1 | = |ω 2 | = 0.5ω s , and the total exerted force and torque (on each) are computed instantaneously for different distances between them.