Surface anchoring controls orientation of a microswimmer in nematic liquid crystal

Microscopic swimmers, both living and synthetic, often dwell in anisotropic viscoelastic environments. The most representative realization of such an environment is water-soluble liquid crystals. Here, we study how the local orientation order of liquid crystal affects the motion of a prototypical elliptical microswimmer. In the framework of well-validated Beris-Edwards model, we show that the microswimmer’s shape and its surface anchoring strength affect the swimming direction and can lead to reorientation transition. Furthermore, there exists a critical surface anchoring strength for non-spherical bacteria-like microswimmers, such that swimming occurs perpendicular in a sub-critical case and parallel in super-critical case. Finally, we demonstrate that for large propulsion speeds active microswimmers generate topological defects in the bulk of the liquid crystal. We show that the location of these defects elucidates how a microswimmer chooses its swimming direction. Our results can guide experimental works on control of bacteria transport in complex anisotropic environments. Interaction of active matter with geometrical constraints is an emerging topic of research due to its potential for design and control of flow patterns. By exploring various modes of swimming and strength of cell-liquid crystal interaction, the authors present a computational model of a microswimmer propulsion in a liquid crystal probing a range of non-trivial swimming dynamics.

U nderstanding self-propulsion in complex and anisotropic fluids is a growing area of research [1][2][3][4][5][6][7][8] . A realization of the experimental system combining motile bacteria with a nematic liquid crystal (LC) was introduced in refs. [2][3][4] . It was shown that bacteria in a LC demonstrate a wealth of intriguing phenomena, among which is the emergence of large-scale collective behavior. The injection of energy by microswimmers was also shown to result in the creation of moving topological defects (local singularities of molecular orientation) in the LC director field [9][10][11] . Manipulation of bacterial trajectories by the prescribed director patterning suggests novel approaches to control the transport of cargo 8 . Furthermore, many biological fluids, such as mucus, DNA solutions, or suspensions of long viruses exhibit a certain degree of liquid crystallinity [12][13][14][15][16] . Liquid crystallinity and mechanical shear also can control bacteria penetration in mucus 17 . Besides the theoretical importance, studies of bacteria motility in complex fluids provide clues in the propagation of bacteria-born diseases and colonization of tissues by pathogens.
The theoretical underpinning of these phenomena requires a better understanding of how an individual microswimmer navigates in the background flow. Although a theory for microswimmers in an isotropic Newtonian fluid is well-developed [18][19][20][21][22] , swimming in complex fluids lacks overall clarity. A distinguishing feature of swimming in an LC is that the LC exerts an elastic torque on microswimmers, forcing them to swim parallel to the LC director 2,23 . However, it was shown that under certain conditions a microswimmer with planar anchoring on its surface can overcome this torque and apparently swim across the preimposed director field 24 . Surface anchoring here means that LC molecules are oriented in a certain way to the surface. Two of the most common types of anchoring are parallel (planar anchoring) and perpendicular (homeotropic anchoring). In the recent study 5 , it was shown that the microswimmer may move parallel or perpendicular to LC director depending on its type, i.e., whether it is a pusher or a puller (a swimmer is called a pusher if the main source of propulsion is located behind, and a puller if otherwise 25 ). Both planar and homeotropic anchoring were tested. Results from ref. 5 suggest that spherical pushers swim parallel, whereas spherical pullers swim perpendicular to the LC director.
Here we investigate the role of the surface anchoring in the determination of stable swimming direction for various microswimmer's properties (aspect ratios of the microswimmer's shape) and their types (pusher or puller). To this end, we use a mathematical model for a microswimmer in LC based on the Beris-Edwards theory. Our main result is that there exists a critical value of the surface anchoring parameter for non-spherical microswimmers. Namely, the swimming of a puller occurs perpendicular in a sub-critical case and parallel in a super-critical case. We also investigate how the magnitude of the critical surface anchoring depends on the physical and geometrical parameters of the problem. Microswimmer's dynamics in LC is highly nonlinear and exhibits a wealth of interesting phenomena. For example, to elucidate the difference between pusher and puller microswimmers' behavior in LC, we consider large propulsion speeds at which swimming in LC is accompanied by formation of in-bulk topological defects. The location of these defects reveals why a microswimmer, whether a pusher or a puller, chooses to swim parallel or perpendicular to LC director. Our results may guide experimental works on control of transport of both biological and biomimetic microswimmers in anisotropic viscoelastic environments.

Results
Here we examine how the stable swimming direction depends on physical and geometrical parameters. Swimmers typically orient themselves in LC either parallel or orthogonal to the far-field director alignment n ∞ . In our numerical studies, we consider n ∞ = i (horizontal, parallel to x-axis) and α(t) is the counterclockwise angle between direction of the principal axis of the microswimmer and n ∞ , as depicted in Fig. 1.
The microswimmer is represented by a rigid body which moves with velocity V(t) and rotates with angular velocity ω(t). Thus, a point of the microswimmer, occupying location x at time t, has instant velocity VðτÞ dτ is the location of the center of mass of the microswimmer at time t. Velocities V(t) and ω(t) are to be determined from the solution of the system coupling the microswimmer and the LC dynamics, described below.
To capture effect of microswimmer's shape anisotropy, the microswimmer is represented by an ellipse whose principal axis is given by unit vector p ¼ ðcosðαÞ; sinðαÞÞ and orientation angle α (Fig. 1). To model self-propulsion of the microswimmer, we impose the velocity slip boundary condition on the microswimmer's surface: Here, θ is the counterclockwise angle between x − x c (t) and vector p, unit vector τ is tangential pointing in the direction of increasing θ, parameter v 0 describes the velocity slip that determines the self-propulsion strength. Boundary condition (1) means that the activity (self-propulsion) enters the model through a prescribed slip velocity. This condition effectively describes a perturbation of the surrounding flow by, for example, rotating flagella used by a bacterium to swim or chemical reactions on the surface of a Janus particle leading to its motion. For a passive swimmer (no self-propulsion), v 0 = 0, Eq. (1) is the standard no-slip condition on solid surface. The expression for the slip velocity (the third term in the right-hand side of (1)) was first derived for a spherical swimmer in ref. 19 and was often used since then to model microswimmers in Stokes flow (e.g., see refs. 5,26 ). Swimming parameter β determines the microswimmer's type, whether it is a pusher or a puller. If β < 0, then the microswimmer is a pusher (e.g., Escherichia coli) such that the propulsion source is at the back. If β > 0, then the microswimmer is a puller (e.g., Chlamydomonas reinhardtii) such that the propulsion source is at the front. For β = 0 the microswimmer is a neutral one (e.g., Paramecium).
Conventionally, an LC is described by unit director field n(x) and scalar order parameter q(x), as in the classical Ericksen-Leslie The domain occupied by the microswimmer is modeled as an ellipse (light blue), oriented parallel to the unit vector p. The liquid crystal far-field director orientation is shown by n ∞ . The angle between the horizontal axis and the vector p is given by α. Small black segments represent the liquid crystal director field; near the microswimmer, these segments are tangential to the microswimmer's surface, thus illustrating planar anchoring. model [27][28][29] . Director n(x) gives the average direction of the LC molecules around location x, whereas order parameter q(x) can be thought of as the variation of LC molecular alignment from the average direction. As "head" and "tail" of an LC molecule are not distinguishable, directions n and −n are equivalent.
We employ the Beris-Edwards model 30 to describe LC and consider a symmetric traceless tensor Q(x) instead of fields n(x) and q(x): Here, I is the d × d identity matrix and d is the space dimension.
The flow of LC is described by velocity vector field v(x) and pressure p(x). The system for Q, v, and p constitutes the Beris-Edwards model, described in "Methods" and Supplementary Note 1.
To prescribe the director anchoring, we impose a Robin-type boundary condition on the tensor Q at the microswimmer's surface: Here, W is the strength of the surface anchoring, ν is the inward (relative to the microswimmer) normal vector, and tensor Q anchor describes, in the sense of relation (2), alignment in the direction n anchor ; q anchor is the corresponding scalar order parameter: Boundary condition (3) describes surface anchoring, i.e., how the microswimmer's surface aligns the nearby LC molecules (or, equivalently, director field n) in a certain direction. Unless otherwise specified, we impose planar surface anchoring, i.e., n anchor = τ, which means that LC director n orients itself tangential to the microswimmer's surface. One of the parameters describing LC is the averaged elastic constant K. From anchoring boundary condition (3), it follows that coupling between microswimmer and LC is characterized by effective anchoring strength W/K. Finally, translational and angular velocities, V(t) and ω(t), are determined via force and torque balances exerted on the microswimmer (see "Methods").
We also performed computational studies for pullers with various values of anchoring strength, see Fig. 2b. For various β and W/K, we computed the relaxation time as orientation angle α (t) converges to its limiting value α 1 ¼ lim t!1 αðtÞ, see Fig. 2c, d. In particular, our computations show that the convergence rate of α(t) for pullers to π/2 as t → ∞ decreases as anchoring strength parameter W increases with fixed K. This result indicates that the larger the anchoring between LC and the spherical puller, the slower the puller turns to swim with orientation angle α = π/2.
Our main finding is that for elongated (elliptical) microswimmers, the asymptotic behavior of the orientation angle α(t) is dramatically different from the spherical case. As an elongated passive rod with planar anchoring aligns along the nematic direction, the tendency to align parallel to LC director is reinforced for an elongated pusher with planar anchoring. Thus, no critical behavior is expected. In contrast, for a puller with planar anchoring, competition between propulsion (favoring perpendicular alignment) and the elongated shape (favoring parallel alignment) leads to a critical transition. Fig. 2 Orientation dynamics of spherical microswimmer. a, b Dynamics of the orientation angle α(t) for spherical microswimmer with planar anchoring for various swimming parameter β (a) and effective anchoring strength W/K for spherical puller with β = 5.0 (b). c, d Relaxation times τ rt , corresponding to a and b, respectively, of converging to the limiting angle, either α → 0 (parallel to far-field director n ∞ ) or α → π/2 (perpendicular to n ∞ ), for various swimming parameter β (c) and effective anchoring strength W/K (d). Markers, red crosses, and blue circles are from numerical integration, dashed lines sketch interpolation curves to visualize how the relaxation time depends on a parameter. As β → 0, the relaxation time diverges.
First, we considered dynamics of α(t) for an elliptic microswimmer with no anchoring, W = 0. We found that there is no symmetry with respect to the swimming parameter β observed for a spherical microswimmer in Fig. 2a, c. Namely, the relaxation time dependence is no longer close to an even function, as it can be seen in Fig. 3a, c. Most importantly, we found that for a puller with positive anchoring strength W > 0, the asymptotic orientation angle α ∞ depends on the value of anchoring strength W, as shown in Fig. 3b, d. We also found that the relaxation time apparently diverges as parameters β or W approach their critical values. For β = 0, one sees from Fig. 3a that α(t) increases and will eventually converge to π/2. On the other hand, for W/K = 0.2 no angular dynamics is clearly observable, as depicted in Fig. 3b. As it is unlikely that there is a stable steady state for π/4, we expect that in this case symmetry will eventually be broken as time evolves and α(t) will asymptotically converge to either 0 or π/2.
To further elucidate dependence of swimming orientation on anchoring strength, we considered a puller (β > 0) for various values of elastic constant K, microswimmer's aspect ratio ϵ, and swimming parameter β. The results of computational modeling are summarized in Fig. 4. We conclude that there exists a critical value of the surface anchoring W crit > 0, which depends, in particular, on K, ϵ, and β. In the sub-critical case, W < W crit , the microswimmer orients perpendicular to LC orientation n ∞ , α ∞ = π/2, whereas in super-critical case, W > W crit , the microswimmer eventually orients parallel to n ∞ , α ∞ = 0. We note that for some values of parameters K, ϵ, and β, finding the critical surface anchoring W crit is not straightforward. Specifically, in these cases, α(t) does not clearly converge neither to 0 nor to π/2 over the entire duration of computational study due to divergence of the relaxation time near the transition parameter value. On the other hand, ranges of W for which the large-time behavior of α(t) is not evidently sub-or super-critical are small intervals that can be effectively represented by a single value W crit .
It is important to note that results on critical anchoring obtained from our computational modeling are relevant for bacteria. Indeed, if the major axis of the domain occupied by a bacterium isL ¼ 10 μm (L is the original length corresponding to dimensionless length L = 1) andK ¼ 10 pN 24 , then the dimensionless parameterŴL=K ¼ WL=K ¼ 0:2 -close to critical value from Fig. 4 (note L = 1, so WL/K = W/K)-is obtained for W ¼ 0:2 Á 10 À6 J m À2 , which is a reasonable value for the anchoring strength parameter 24 . However, the anchoring strength can vary in a wide range depending on the surface treatment, e.g., see ref. 31 .
A crucial finding here is the existence of critical value of anchoring strength W crit . Namely, the elongated puller reorients parallel to LC director for small W as opposed to a spherical puller which reorients perpendicular to LC director for all W ≥ 0. It is caused by the aligning torque exerted by LC on an elongated object via planar surface anchoring condition (3). Indeed, for an elongated object with planar anchoring, the total LC energy is minimized at parallel-to-LC orientation α = 0 and thus this orientation is stable, e.g., see ref. 24 . It also implies that an elongated pusher with planar surface anchoring always reorients parallel to LC and does not possess a critical anchoring strength.
A somewhat opposite situation realizes for microswimmers with the homeotropic anchoring instead of planar one, i.e., nearby LC molecules orient perpendicular to the microswimmer's surface (n anchor = ν in (4)). In this case, the perpendicular orientation α = π/2 becomes more preferable than the parallel one with α = 0. Thus, an elongated puller reorients perpendicular to LC director, whereas an elongated pusher is expected to possess a critical anchoring: for small W, it reorients parallel to LC director and for large W it reorients perpendicular. Our numerical modeling for elongated pusher and homeotropic boundary conditions confirms this expectation, see Fig. 5a, b and Supplementary Note 4. In ref. 24 , apparent transition from parallel to perpendicular, to the director swimming was observed for pusher-like swimmers (Bacillus subtilis bacteria). However, this transition was a combined effect of homeotropic anchoring on the bounding plates and geometric confinement. Non-motile or slow bacteria aligned parallel to the director, i.e., perpendicular to the bounding plate. Sufficiently strong bacteria were able to distort homeotropic anchoring and make it local planar and thus swim parallel to the bounding plates. In a related study, ref. 31 , with Proteus mirabilis bacteria, no such transition was observed.

Discussion
Our computations lead to a hypothesis that orientation dynamics can be effectively described by a one-dimensional ordinary differential equation _ α ¼ f ðαÞ. The form of scalar function f(α) depends on parameters β, ϵ, K, v 0 , etc. This hypothesis implies that steady swimming directions are determined by zeros of function f(α). For 0 ≤ α ≤ π/2, two such steady directions are given by symmetry of the problem: α = 0 and α = π/2.
We show that there exists the critical value of anchoring strength for fixed initial angle α| t=0 = π/4. We performed numerical simulations for various initial angles α| t=0 and we observed that effective dynamics is of the form where γ > 0 is a proportionality coefficient depending on problem parameters and γ can be found via γ = 1/(2τ rt |W crit − W|) where τ rt is the relaxation time, plotted in Figs. 2c, d, 3c, d, and 5b as a function of W or β with all other problem parameters fixed. One can easily verify that solutions α(t) of Eq. (5) converge to π/2 for W < W crit and to 0 for W > W crit regardless of initial conditions, and there is no intermediate steady state. To illustrate these observations, we depicted results of numerical integration of the model described in "Methods" for a sub-and super-critical values of W and various initial angles in Fig. 6. Thus, a microswimmer, in the framework of Beris-Edwards model described in "Methods", displays two steady directions with respect to LC director: parallel and perpendicular. Although the choice of the stable direction is the result of the competition between elastic and active hydrodynamic torques (as we explain below), it can be controlled by the anchoring strength W. Equation (5) and the existence of critical anchoring strength can also be explained as follows. In ref. 5 , it was justified that LC exerts the effective torque on spherical planar puller T eff / v 0 β sinð2αÞ, wherein T eff is independent of W. Thus, a puller, for any value β > 0 and W, eventually orients itself perpendicular to LC director, i.e., α| t≫1 ≈ π/2 for any initial angle 0 < α(0) ≤ π/2. In this work, we extend consideration to elongated microswimmers. If the surface anchoring is planar, it is known 24,32 that LC exerts stabilizing torque T stab on an elongated microswimmer, which rotates the microswimmer so that its major axis becomes parallel to the LC director. Stabilizing torque can also be explained by that the effective shear viscosity for velocities perpendicular to the director is larger than the viscosity parallel to the director, as discussed in ref. 33 . For a spherical particle, full hydrodynamic torque (including active torque in our notations), arising from the viscosity anisotropy leads to the direction selection 5 . In Supplementary Note 3 by combining arguments from ref. 5   non-spherical microswimmer, we show that T stab / ÀW sinð2αÞ. Therefore, the dynamics of an elongated puller microswimmer with planar anchoring in LC is the result of competition between active torque T eff (let us call it "active," as it is proportional to both self-propulsion v 0 and swimming parameter β) and stabilizing torque T stab . The active torque rotates the microswimmer perpendicular to LC director, whereas stabilizing one rotates the microswimmer parallel to LC director and the critical case occurs when the magnitudes of the torques coincide. Furthermore, if we add written above quantities to which both torques, T eff and T stab , are proportional, then we obtain the right-hand side of Eq. (5) for some γ and W crit .

and the expression for torque for
The above results were obtained for a relatively small selfpropulsion parameter v 0 . As v 0 increases, the steady-state swimming may become unstable and the microswimmer may generate topological defects in its wake. It is well-known that due to the anchoring surface condition, topological defects form near the microswimmer, e.g., so-called boojums appear for planar anchoring 34 . In our numerical studies, we have shown that if one increases self-propulsion parameter v 0 , then topological defects appear away from the microswimmer, see Figs. 7 and 8.
Specifically, we tested both puller and pusher microswimmers with stable horizontal (i.e., parallel to LC director) direction of motion for various values of self-propulsion parameter v 0 and fixed planar anchoring strength such that W/K = 0.35. When selfpropulsion is weak, v 0 = 0.1, we observed for both a pusher and a puller that there are two near-surface boojum defects, at the front and the rear. For a larger v 0 , the behavior of defects for pushers and pullers differs. Namely, the scalar order parameter for a puller displays two parabolic-shaped prolate valleys emanating from two side points of the microswimmer (blue zones in Fig. 7a-f). These two points are where self-propulsion velocity given by the additional slip component in flow velocity boundary condition (1) vanishes, i.e., where 1 þ β cosðθÞ ¼ 0. At these points, the no-slip boundary condition holds, so the relative fluid velocity with respect to moving microswimmer vanishes, but the velocity gradient is large due to large v 0 . Thus, the advection term S dominates in the equation for Q (see "Methods") at the point where fluid velocity vanishes, leading to a "defective" zone for LC, see Fig. 7g.
On the contrary, in computational modeling for a pusher, valleys with small scalar order parameter grow from the rear and the front of the microswimmer, not from sides, see Fig. 8a-f. All defects are located along the horizontal line containing the major axis of the microswimmer, see Fig. 8g.  Fig. 7 In-bulk defects around puller. a-f Scalar order parameter q(x) around horizontally translating puller microswimmer with planar anchoring and various values of self-propulsion parameter v 0 . Topological defects form as v 0 increases. For v 0 = 5.0, locations of defects are denoted by red and green circles, corresponding to defects of topological indexes "− 1/2" and "+1/2," respectively. g Director field for v 0 = 5.0; location of topological defects, denoted by red and green circles, are numerically determined as local minima of scalar order parameter which are less than q 0 = 0.05. The difference in geometry of topological defect locations can be explained by the velocity field generated by the microswimmer, depicted in Fig. 9a, d and sketched in Fig. 9b, e. Let us assume that defects are passive objects dragged along streamlines. Figure 9b, e shows that defects will likely end up moving along the red dashed curves that come out from side points for the flow generated by a puller, and from the front and rear for a pusher. This is consistent with our numerical modeling presented in Figs. 7 and 8. Moreover, the velocity field generated by the microswimmer explains why pullers, as opposed to pushers, with planar anchoring tend to swim perpendicular to the far-field LC orientation. We consider the side points on the microswimmer's surface, where selfpropulsion velocity vanishes (i.e., where 1 þ β cosðθÞ ¼ 0), and examine the velocity field at the vicinity of these points. From sketches in Fig. 9c, f, we conclude that at these points LC molecules tend to orient perpendicularly to a puller (in the sub-critical case, with not strong planar anchoring) and parallel to a pusher (this orientation configuration is consistent with our numerical modeling, see Figs. 7g and 8g). Therefore, if our consideration boils down to these side points, then we reproduce our main observation for the case of weak planar anchoring: pullers are perpendicular to LC, whereas pushers are parallel.
Overall, the phenomenon resembles the onset of the von Kármán vortex street in the flow past obstacle. In Fig. 7f, g, one can see a chain of defects with non-integer topological degrees inside these zones. Defects significantly affect other microswimmers' trajectories 9,10 and thus affect interactions between microswimmers in LC.
In conclusion, in this work we demonstrate that surface anchoring and shape anisotropy control the dynamics of microswimmers in a nontrivial way. For example, a puller may experience a transition from parallel to perpendicular orientation as the surface anchoring strength is decreased. Moreover, the increase in the self-propulsion speed may lead to unsteady motion and defect generation. Overall, these findings provide new insights into the control of biological and synthetic microswimmers in complex anisotropic fluids. With the rapid development of 3D microprinting technology, microswimmers with desired complex shape can be designed 35 . Patterning and functionalization of the microswimmer's surface may control the slip velocity distribution and thus the puller/pusher type (parameter β in our description). In addition, treating the microswimmers with surfactants or other surface-active molecules can control the anchoring strength and anchoring type (planar vs. homeotropic) 24,36 . Thus, by engineering shape and anchoring on the surface of the microswimmer, one enables control of swimming behavior.

Methods
Computational model. The system for Q, v, and p constitutes the Beris-Edwards model: Equation (6) describes the dynamics of tensor Q, and S = S(∇v, Q) measures how the flow gradient ∇v rotates and stretches the order parameter. The expression for S(∇v, Q) can be found in Supplementary Note 1. The parameter Γ is the director relaxation rate and H is the variation of the Landau-de Gennes free energy F : where K is the averaged elastic constant and parameters a, b, and c are the Landaude Gennes coefficients. The external force F exter is introduced for two-dimensional case (d = 2) only. It approximately describes the effect of upper and lower bounding plates confining the LC film. Due to planar anchoring on the plates, the aligning torque F exter is imposed on Q. We define the term F exter by 9 : where ζ an is the alignment strength, R π/2 is the matrix of rotation by angle π/2, Q f ¼ ff À 1 2 I is the tensor corresponding to the constant aligning direction f = n ∞ , |f| = 1. The term F exter does not change the magnitude of Q but forces Q to be collinear with Q f . Director f = n ∞ represents the far-field alignment, i.e., n ≈ n ∞ far from the microswimmer, as shown in Fig. 1. Equality n(x) ≡ n ∞ corresponds to the case when there is no microswimmer and LC is in equilibrium. Equation (7) describes dynamics of velocity field v and is derived from the linear momentum balance equation. The term σ H = η(∇v + ∇v T ) − pI is the hydrodynamic stress tensor with viscosity η. The term σ LC takes into account how anisotropy of LC affects flow v, for its definition see Supplementary Note 1. Equation (8) is the incompressibility condition.
System ((6), (7), (8)) is supplemented with boundary conditions: selfpropulsion slip condition (1) and surface anchoring condition (3). It is noteworthy that if no external force F exter , such as the one defined by (10), is exerted, and flow v is negligible, then the right-hand side of equation (6) and boundary condition (3) on the microswimmer's surface S can be recovered by variation of the LC energy functional E Q , where The second term in the definition of E Q accounts for the surface anchoring energy and its form is similar to the classical Rapini-Papoular condition 37 . We write the force and torque balance for a microswimmer in an over-damped limit (i.e., microswimmer's inertia is neglected): Here, F ν (x) = (σ H + σ LC )ν is the force exerted by the LC on the microswimmer at point x on the microswimmer's surface S. The second term in the left-hand side of torque balance (13) originates from an asymmetric part of the stress tensor. Thus, the total torque exerted by the media (LC) on the rigid microswimmer is the sum of torques F ν and additional force couples ℓ. Using the energy principle, we show (see Supplementary Note 2) that ℓ can be computed by Here, α is the orientation angle of the microswimmer and e z = (0, 0, 1) T . Provided that boundary condition (3) is taken into account, Eq. (14) is equivalent to the expression for the couple stress tensor in ref. 38 and how the torque is computed in the lattice Boltzmann computational model in ref. [39][40][41][42] . Thus, to simulate microswimmer dynamics in LC, we solve Eqs. (6), (7), and (8) subject to boundary conditions (1) and (3), as well as force and torque balances (12) and (13).
Numerical implementation. The computational model is simulated by applying finite element method (FEM) and using the open-source finite element platform FreeFEM++ 43 . Equations are solved in the microswimmer's frame, i.e., in the frame moving with velocity V(t), so the center of the microswimmer is always in origin. It is noteworthy that the frame does not rotate with the microswimmer, i.e., the microswimmer does not translate but rotates in this frame. The long axis of the microswimmer L is chosen to be 1. We simulate the model in the box (−R, R) × (−R, R), R = 20, imposing periodic boundary conditions on Q, v, and p for x = ±R and y = ±R. We assume that the Reynolds number is small and thus we neglect the left-hand side of (7) by equating it to 0.
For time-dependent problem, given initial orientation angle α(0), we compute Q(0) by solving stationary problem ΓH + F exter = 0 with boundary conditions (3). To this end, we consider auxiliary problem ∂ tQ ¼ ΓH þ F exter withQð0Þ ¼ q 0 ðn 1 n 1 À I=2Þ such that q 0 ≠ 0 and H = 0 for Q ¼Qð0Þ. Then initial condition Q(0) is obtained as the large-time limit ofQ. At each time step, we first update Q from (6) using semi-implicit scheme: ðQ nþ1 À Q n j rÀv n dt Þ=dt À ΓKΔQ nþ1 ¼ S n þ ΓðH n À KΔQ n Þ þ F n exter (here n is the time step number). Here, we used FEM with piecewise quadratic elements (P 2 elements) for Q. Then, v is updated by inverting Stokes operator in (7) where the right-hand side is with Q = Q n+1 . This is done by FEM with piecewise quadratic elements for flow v (P 2 elements) and piecewise linear elements for pressure p (P 1 elements). Finally, V and ω are computed using force and torque balances with the updated v and Q.
We use an adaptive triangular mesh for FEM so that the mesh is finer close to the microswimmer. Specifically, the mesh is such that there are 200 grid points along the microswimmer's surface and 400 on the periodic box boundary. Near the microswimmer, the ratio between average diameter of mesh (triangle) and the diameter of microswimmer (its long axis) is~1 : 200. Convergence in spatial discretization was checked using finer mesh with ratio from 1 : 200 to 1 : 1600 for some selected cases. At each time step, the mesh is re-generated due to the rotation of microswimmer. We also checked that the results of numerical simulations are stable with respect to the size of the periodic cell R ≫ 1 by performing simulations for R = 20 and R = 40. In the study of topological defects, steady swimming direction was obtained by simulating the time-dependent problem for sufficiently large number of time steps.

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

Code availability
The code to carry out the simulations is available from the corresponding author on reasonable request.