Microswimmers in an axisymmetric vortex flow

Microswimmers are encountered in a wide variety of biophysical settings. When interacting with flow fields, they show interesting dynamical features such as hydrodynamical trapping, clustering, and preferential orientation. One important step towards the understanding of such features is to clarify the interplay of hydrodynamic flow with microswimmer motility and shape. Here, we study the dynamics of ellipsoidal microswimmers in a two-dimensional axisymmetric vortex flow. Despite this simple setting, we find surprisingly rich dynamics, which can be comprehensively characterized in the framework of dynamical systems theory. By classifying the fixed point structure of the underlying phase space as a function of motility and swimmer shape, we uncover the topology of the phase space and determine the conditions under which swimmers are trapped in the vortex. For spherical swimmers, we identify Hamiltonian dynamics, which are broken for swimmers of a different shape. We find that prolate ellipsoidal shaped microswimmers tend to align parallel to the velocity field, while oblate microswimmers tend to remain perpendicular to it. Additionally, we find that rotational noise allows microswimmers to escape the vortex with an enhanced escape rate close to the system's saddle point. Our results clarify the role of shape and motility on the occurrence of preferential concentration and clustering and provide a starting point to understand the dynamics in more complex flows.

Among the various mechanisms, shape and motility are the key parameters in quantifying microswimmer interaction with hydrodynamic flows [26][27][28][29]. For instance, motility is a crucial ingredient for the emergence of clustering of neutrally buoyant particles [13,16,27]. Shape, on the other hand, determines the dynamic reaction to hydrodynamic cues. Rod-like agents, for example, tend to align with the direction of local vorticity, while also spinning due to it [30][31][32][33][34]. Furthermore, it has recently been shown that self-propelled rod-shaped particles tend to align with the velocity field [29].
Given the complexity of turbulent and spatiotemporally varying flows, several investigations have focused on simple flows, which allow making contact with dynamical systems theory. For example, the (quasi-)periodic [35,36] and chaotic [37] motion of microswimmers in a Poiseuille flow can be understood through the underlying Hamiltonian dynamics. As another example, the ubiquity of vortices in natural environments and their dynamical impact on biological agents [38] continues to make the study of single vortex structures a good starting point for understanding microswimming in more complex, biologically relevant flows [39]. Isolating the impact of individual vortices on microswimming in a two-dimensional cellular flow has revealed barriers to particle transport as a function of shape and swimming speed [40,41]. Shape deformation, as another example, has been found to have a strong impact on the scattering dynamics of individual microswimmers in a single vortex structure [42]. This illustrates how the investigation of simple flow settings sheds light on the interplay between shape and swimming speed on microswimmer dynamics.
Here, we comprehensively characterize microswimmer dynamics in a single vortex structure by relating the observed physical phenomena to properties of the underlying dynamical system. In particular, we consider non-interacting swimmers in a twodimensional axisymmetric vortex flow and address their trapping properties as well as the occurrence of clustering (the spatially heterogeneous distribution of particles) and preferential orientation with respect to the velocity field. We idealize microswimmers as advected ellipsoidal particles, which additionally have a swimming direction and a constant self-propulsion speed. Furthermore, vorticity and shear induce particle tumbling, which alters the swimming direction.
This simple setting reveals surprising insights: we identify an effective swimming velocity, which takes into account shape, as a control parameter for this system. The effective anding of such features is to clarify the interplay of the hydrodynamic flow with microswimmer motility and shape. Here, we study the dynamics of ellipsoidal microswimmers in a two-dimensional axisymmetric vortex flow. Despite this simple setting, we find surprisingly rich dynamics, which can be comprehensively characterized in the framework of dynamical systems theory. By classifying the fixed-point structure of the underlying phase space as a function of motility and swimmer shape,swimming velocity allows to distinguish microswimmers that escape the vortex from the ones which remain trapped. Moreover, we find that spherical swimmers obey Hamiltonian dynamics, whereas other shapes break the symplectic structure of phase space. Hence phase-space contraction, and ultimately preferential concentration, can be set into the context of breaking of Hamiltonian dynamics by departure from a spherical shape. Finally, to quantify the robustness of our results, we investigate the impact of rotational noise. We find that a saddle point present in the relevant phase space plays an important role for the escape of microswimmers from the vortex core.

Model equations
We model microswimmers as inertialess particles advected by a velocity field u. The particles are additionally capable of self-propulsion with swimming velocity v s in directionp [43]. The swimmer position x obeys the equation of motioṅ where the flow field is evaluated at the Lagrangian position of the swimmer u(t, x(t)). A simple but effective way to introduce shape is to consider ellipsoidal microswimmers [44]. Particle orientation can be described the the particle's symmetry axisp. Additionally, in many relevant settings, microswimmers are smaller than the smallest hydrodynamic scales (such as the Kolmogorov length scale in turbulence). In this limit, the spinning and tumbling of the particles' orientation can be described by Jeffery's equation [45], which takes the forṁ Here, S ij = (∂ i u j + ∂ j u i )/2 is the strain tensor and ω = ∇ × u the vorticity. The ratio between the ellipsoid's major and minor axes λ defines the shape parameter as α = (λ 2 − 1)/(λ 2 + 1). The parameter α interpolates shapes between an oblate ellipsoid (0 < α < 1), a sphere (α = 1), or a prolate ellipsoid (−1 < α < 0).

Two-dimensional vortex flow
In two spatial dimensions, (1) and (2) have three degrees of freedom: two coordinates for position and one swimming angle It is also worth noting that in two dimensions ellipsoids can only tumble. In the following, we restrict ourselves to an axisymmetric vortex flow of the form u(t, x) = u(r)ê φ , where r = x 2 + y 2 is the radial coordinate andê φ a unit vector along the angular coordinate. By inserting (3) into (1) and (2) we obtain the equations of motioṅ Inspection of these equations in polar coordinates reveals rotational invariance. Without loss of generality, the dynamics of the system can therefore be reduced by one degree of freedom by transforming into a co-rotating frame. This is achieved by introducing new coordinates In the co-rotating coordinate system {X, P }, the swimmers are rotated so that they swim along the positive X axis. That is, microswimmers only differ by a rotation from the laboratory frame {x, y}. Hence, the radial coordinate is identical in both coordinate systems, Introducing an angle variable defined through tan ψ := P/X helps in the description of the dynamics. The equations of motion for X and P can be obtained from (4)-(6) aṡ where h(ψ) := cos(2ψ) and f (r) := r 2 d dr The function h(ψ) is a geometric term stemming from the relative orientation of an ellipsoid with respect to the flow. The function f (r) contains the dependence on the velocity field u(r). The angle θ, characterizing the orientation, is a slave variable to X and P and evolves according to the equatioṅ Hence, this co-rotating representation reduces the three-dimensional dynamics of the swimmers in the laboratory frame {x, y, θ} to an effective two-dimensional dynamics in the co-rotating frame {X, P }. As a concrete example, we consider a stationary Lamb-Oseen vortex, a prototypical vortex structure which is representative for a large class of hydrodynamic vortices [46]. The essential feature of this field is that the velocity profile interpolates between linear growth near the core (corresponding to a solid body rotation) and a Gaussian decay far from the core, leading to differential rotation. Its velocity field is given by where u 0 and r 0 are the maximum azimuthal velocity and its corresponding radial coordinate. The vortex is stationary, and hence its width r 0 is kept fixed. The constant σ is the nontrivial solution to the equation exp(σ) = (1 + 2σ) [47], which is obtained by fixing u 0 to r 0 . Using r 0 and r 0 /u 0 as length and time scales, respectively, (4)-(6) and (8)-(11) can be non-dimensionalized: As a result, we have a non-dimensional swimming velocityṽ s = v s /u 0 . In the following, we drop the tildes and work with the non-dimensionalized swimming speed. For the numerical results, we integrate the ordinary differential equations (4)-(6) using a fourth-order Runge-Kutta method with a time step dt = 0.05. For each of the simulations shown in figure 1 and figure 3 we initialized 5 × 10 5 trajectories. For the radial distribution function in figure 4 we initialized 8 × 10 6 microswimmers. Typical dynamics of microswimmers following (4)-(6) and the corresponding representation in the co-rotating frame {X, P } are shown in figure 1. Here, the quasiperiodic trajectories observed in the laboratory frame can be explained by the coupling of the typical angular velocity of the microswimmers in the co-rotating frame with their rotation frequency given by (11).
In the next section, we explore the features of (8) and (9) to uncover the fixed-point structure of the underlying dynamics and precisely characterize the microswimmer observed phenomena.

Fixed-point analysis
Visual inspection of the microswimmer dynamics in the co-rotating frame in figure 1 reveals an intricate behavior. These observations can be made precise by a fixed-point analysis. Because the dynamical system (8)- (9) is two-dimensional, the Poincaré-Bendixson theorem [48] applies, and the topology of the dynamics is completely determined by the fixed points. Already in figure 1 (d)-(f) we can visually identify two fixed points in phase space. These fixed points are marked as blue and green dots in figure 2 (a) and (b), and are studied in the following. Calculating the fixed points {X FP , P FP } of (8)-(9) trivially leads to X FP = 0 for any axisymmetric flow profile. As a consequence, r = |P FP | and ψ = sign(P FP ) π/2. The subsequent analysis depends on the specific flow. For the Lamb-Oseen vortex (12) we find that the inequality f (r) ≤ 0 is valid for all radii. This means that fixed points exist only for P FP > 0 and ψ = π/2. The fixed points of (8)- (9) are then given by the solutions to The type of fixed points of a two-dimensional system is determined by the determinant and the trace of the Jacobian matrix J [48]. At the fixed-point coordinates (13) we obtain Tr(J) = 0, These quantities reveal the nature of the fixed points and their dependence on swimmer shape. Because det(J) ∈ R and Tr(J) = 0, this type of system can only have either center or saddle points. As a consequence, trajectories remaining bound to the vortex core correspond to areas in phase space enclosed by the homoclinic or heteroclinic orbits of the saddle points. All other areas in phase space lead to and come from infinity. Moreover, (13) motivates the definition of an effective swimming velocity as As the trace and determinant are independent of the swimming speed v s , and the term (1 − α) 2 is non-negative, swimmers with the same effective swimming velocity have identical types of fixed points located at the same coordinates. This implies that the effective swimming velocity can be used to classify swimmers according to their shape and swimming speed. Therefore, v eff plays the role of a control parameter for the topology of the phase space. To obtain the bifurcation diagram of this system the equation can be solved graphically. The solution to this equation for the Lamb-Oseen vortex (12) is shown in figure 2 (c). In the case of the Lamb-Oseen vortex (figure 2) the condition 0 < v eff < v max ensures the existence of a solution pair to (13), which corresponds to a center-saddle pair. As long as this solution exists, so does the homoclinic orbit, enclosing a region of trapped microswimmers. Figure 2 (a) shows the fixed-point pair and the homoclinic orbit. By choosing a larger value for v eff , as in figure 2 (b), the fixed points merge and undergo a saddle-node bifurcation at v max . Additionally, the condition v eff < v max defines a region in the {v s , α} parameter space for which swimmers are trapped. For the Lamb-Oseen vortex we numerically obtain v max ≈ 0.726 as the maximum of (16). Hence, by taking into account the role of shape, microswimmers with much lower swimming velocity v s than the maximum fluid velocity u 0 can escape the vortex. More concretely, for constant v s thin elongated microswimmers have a divergent v eff in the limit α → 1, and hence always escape the vortex in this limit.

Hamiltonian dynamics and phase-space contraction
Next, we explore the effect of shape on the microswimmer dynamics and its relation to phase-space contraction. Initializing microswimmers homogeneously inside the homoclinic orbit leads to shape-dependent stationary distributions, as shown in figure 3.
For equal effective swimming velocity v eff , changing shape leads to a variety of density  Figure 3.
Microswimmer distribution (gray density plot) for different shapes and effective swimming velocities, namely v eff = 0.25 in (a)-(c) and v eff = 0.5 in (d)-(f). For spherical shape (b,e) we clearly observe the phase-space-conserving Hamiltonian dynamics in {X, P } space. A homogeneously filled homoclinic orbit (red line) at initial time will remain so throughout time evolution. Other shapes break the Hamiltonian structure, and phase-space contraction Γ (23) leads to preferential alignment of microswimmers with respect to the velocity field. Oblate microswimmers (a,d) tend to concentrate along the X-axis while prolate microswimmers (c,f) tend to concentrate along the P -axis. This corresponds to preferential swimming perpendicular or parallel to the flow, respectively, as illustrated by the denser regions in the {X, P } plots for non-spherical microswimmers.The black lines show sample trajectories.
An even stronger statement can be made in our case. By using the {X, P } coordinates of the co-rotating frame, we can reveal the Hamiltonian structure of the equations of motion, i.e.Ẋ = ∂H(X, P ) ∂P , where the Hamilton function is given by for arbitrary axisymmetric velocity fields. As a consequence of the Hamiltonian dynamics, phase-space volume is conserved and does not contract or expand. Therefore, starting from a homogeneous distribution, spherical swimmers maintain a homogeneous distribution as time evolves. While spherical swimmers obey Hamiltonian dynamics, other swimmer shapes break the Hamiltonian structure. This can be seen by recasting (8) and (9) aṡ As a result, phase space can contract or expand for non-spherical swimmers. The phasespace contraction rate is given by Non-spherical microswimmers develop denser accumulations in some regions of phase space and depart from an initially homogeneous distribution. The fact that for the Lamb-Oseen vortex (12) the inequality f (r) ≤ 0 holds, together with (23) implies a constant sign of Γ inside each quadrant in {X, P } space. As trapped microswimmers traverse the different quadrants, they periodically switch between expanding and contracting quadrants. Microswimmers exiting an expansion quadrant will show a minimum in density, whilst those exiting a contraction quadrant will show maximum density. Hence, by considering the sign of α and f (r), we conclude that denser regions are formed along the X-axis for oblate ellipsoids and along the P -axis for prolate ellipsoids.
In figure 3 the different density regions for oblate and prolate microswimmers can be identified, as well as the homogeneous distribution for spheres. Recall that in the corotating coordinate frame {X, P }, the swimmers are rotated so that they always point in the positive X direction. The velocity field, on the other hand, is rotationally invariant. Therefore contraction along the different axes reveals that oblate ellipsoids (α < 0) swim predominantly perpendicular to the velocity field (denser regions along the Xaxis) while prolate ellipsoids (α > 0) mostly remain parallel to it (denser regions along the P -axis). That means that phase space contracts in such a way that, starting from random initial conditions, trapped microswimmers show shape-dependent preferential orientation parallel or perpendicular to the flow. Similar effects have been observed in chaotic, mildly turbulent flows [29]. Interestingly, the dynamical features of the microswimmers in the co-rotating frame lead to clustering in the laboratory frame. To characterize the spatial distribution of microswimmers, we consider the radial distribution function. As the laboratory frame {x, y} and the co-rotating frame {X, P } differ only by a rotation, the radial distribution of microswimmer ensembles is identical in both cases. That means that integrating the distribution function in the co-rotating frame along the angle variable exactly corresponds to the radial distribution function in the laboratory frame for any ensemble configuration. In the case of trapped spherical microswimmers (α = 0), starting from homogeneous initial conditions which fill out the homoclinic orbit (as in figure 1), the radial distribution function is unaltered as time evolves. In this case, integrating a constant density inside the homoclinic orbit over the angle variable yields the radial distribution function. Here, it is not necessary to integrate the equations of motion for an ensemble of microswimmers to obtain the radial distribution function; it can be obtained from the shape of the homoclinic orbit alone. For other microswimmer shapes this approach is not feasible, as phase-space contraction sets in under time evolution and a non-trivial stationary distribution develops. Microswimmers with equal v eff have the same type of fixed points at the same coordinates. Nevertheless by changing shape we observe a variety of quasi-periodic orbits and density distributions. These differences are rooted in the shape of the homoclinic orbit as well as the phase-space contraction rate Γ, which is shown in the background of figure 3, normalized by the maximum contraction rate Γ max := 2 max(|f (r)|). Therefore, the radial distribution function of bound microswimmers is a function of both phase-space contraction and the shape of the homoclinic orbit. Both of these differ for swimmers of different shape and swimming speed, even if their effective swimming speed is identical (see figure 4). However, the maximum swimming radius r 0 of trapped microswimmers, beyond which the radial distribution function is zero, is a common property of swimmers with identical effective swimming speeds. This can be explained by the fact that the saddle point is the point on the homoclinic orbit with the largest radius. Hence the saddle point determines the maximum swimming radius of trapped microswimmers, which is a constant for microswimmers of equal effective swimming speed.

Impact of rotational noise
So far, we have considered only deterministic swimmers. In realistic biophysical settings, fluctuations play an important role. We explore this by considering the impact of rotational fluctuations on microswimmers. We introduce a Wiener process dW into (2) and obtain dp =ṗ dt +ê z ×p g dW, where the deterministic partṗ is given by (2). Here we are using the non-dimensionalized quantities, such that a Péclet number can be defined as Pe := 1/g 2 . This equation is understood in the Stratonovich sense so thatp remains normalized. With the parametrization (3) the equations of motion for microswimmers with rotational noise correspond to additive noise on the swimming angle dx =ẋ dt, dθ =θ dt + g dW, (27) where again the deterministic partsẋ,ẏ, andθ are given by (4)- (6). For the numerical results presented in this section, we solve (25)- (27) using the Euler-Maruyama method with a time step dt = 0.0005 and a total of 5 × 10 5 trajectories. Switching to the corotating {X, P }-frame, we obtain the stochastic equations (in the Stratonovich sense) dX =Ẋ dt + g P dW, withẊ andṖ given by (21) and (22). The effect of noise in the co-rotating frame corresponds to rotational diffusion. This can be seen from the Fokker-Planck equation for the density distribution function ρ where Γ is the phase-space contraction rate (23) and µ is a rotational drift operator Additionally D is a diffusion operator given by Setting g = 0 in (30) yields a Liouville equation for the deterministic part of the dynamics (8)- (9). For v max < v eff the drift term due to swimming dominates, and microswimmers always escape the vortex core. Without swimming (v eff = 0) only rotational drift is present, i.e. the microswimmers behave as passive tracers. In the regime 0 < v eff < v max the presence of the fixed-point pair and the homoclinic orbit leads to a stationary solution to the Liouville equation, corresponding to stationary distributions as shown in figure 3. The addition of rotational noise on the swimming direction leads to rotational diffusion (32) in the co-rotating frame. This induces microswimmer transfer across the homoclinic orbit. Therefore, starting from a stationary distribution of microswimmers in the homoclinic orbit, all microswimmers will eventually escape the vortex core. Figure 5 illustrates this phenomenon for various swimmer shapes. The saddle point plays an important role in this context: the maximum expansion direction at this fixed point leads to enhanced diffusion, allowing microswimmers to escape the vortex core faster than they would do with just rotational difussion. For a constant v eff we observe that oblate ellipsoids escape the homoclinic orbit faster than prolate ellipsoids. This is due to the fact that for constant v eff prolate microswimmers swim slower than oblate microswimmers.

Summary and conclusions
We have studied self-propelled ellipsoidal particles as idealized microswimmers in a twodimensional axisymmetric vortex flow. In particular, we have investigated under which conditions swimmers are trapped by the vortex, and whether they exhibit preferential orientation. This simple setting reveals surprising insights: due to the axisymmetry of the problem, the phase space is two-dimensional and can be parameterized by the swimmer's radial position and an orientation angle (relative to the position vector). Topologically, the phase space features a saddle point and a center. Microswimmers bound to the vortex core follow closed orbits inside a homoclinic orbit. Clustering in the laboratory frame occurs as a consequence of phase-space contraction. Shape plays a decisive role: for spherical swimmers, we have showe that the dynamics is Hamiltonian, excluding clustering as a result of phase space conservation. However, non-spherical particles break the Hamiltonian structure, hence enabling phase-space contraction and shape-induced clustering.
To determine whether microswimmers are trapped, we identified the effective swimming velocity as the central control parameter. The effective swimming velocity depends both on the swimming velocity and on the shape: at a constant swimming velocity, prolate ellipsoids have a larger effective swimming velocity than oblate ellipsoids. This allowed us to map different swimmer shapes to topologically equivalent phase spaces. Using a bifurcation analysis, we determined the maximum velocity for a given flow profile such that swimmers faster than this velocity are fast enough to escape the vortex core.
Notably, this maximum velocity is lower than the maximum azimuthal flow velocity, implying that microswimmers with a smaller swimming velocity than the advecting flow field can escape the vortex. Shape plays a role also here as prolate swimmers can more easily escape than oblate swimmers. As the shapes of many plankton and bacteria species can be approximated as thin rods [44], this effect may have implications for motile species in aquatic environments. In particular, a prolate shape may yield advantages such as avoiding hydrodynamic trapping while grazing or escaping predation, without having to dedicate additional energy to swim faster. Finally, we investigated the impact of rotational noise. We find that the inclusion of rotational noise allows for initially trapped microswimmers to escape the vortex core. The presence of the saddle point leads to enhanced escape rates as the maximum expansion direction quickly drives microswimmers away from the vortex core. While we focused on the simple case of an axisymmetric flow, our results might help to also better understand the impact of shape and motility on microswimmer dynamics in more complex flows.

Acknowledgments
This work was supported by the Max Planck Society. MW also gratefully acknowledges a Fulbright-Cottrell Award grant. The authors thank Ramin Golestanian for helpful discussions.