Optimal slip velocities of micro-swimmers with arbitrary axisymmetric shapes

This article presents a computational approach for determining the optimal slip velocities on any given shape of an axisymmetric micro-swimmer suspended in a viscous fluid. The objective is to minimize the power loss to maintain a target swimming speed, or equivalently to maximize the efficiency of the micro-swimmer. Owing to the linearity of the Stokes equations governing the fluid motion, we show that this PDE-constrained optimization problem can be reduced to a simpler quadratic optimization problem, which we solve using a high-order accurate boundary integral method. We consider various families of shapes parameterized by the reduced volume and compute their swimming efficiency. We found that for a given reduced volume, prolate ellipsoids are the most efficient micro-swimmer shapes and that, irrespective of the shape, the optimal slip always corresponds to a neutral swimmer (as opposed to a pusher or a puller).


Introduction
The squirmer model (Lighthill 1952;Blake 1971) is widely adopted by mathematicians and physicists over the past decades to model ciliated micro-swimmers such as Opalina, Volvox and Paramecium (Lauga & Powers 2009). On a high level, this continuum model, sometimes referred to as the envelope model, effectively tracks the motion of the envelope formed by the tips of the densely-packed cilia, located on the swimmer body, while neglecting the motion below the tips. Individual and collective ciliary motions could be mapped to traveling waves of the envelope on the surface. Assuming no radial displacements of the surface and time-independent tangential velocity led to the simpler steady squirmer model (see Pedley 2016), wherein, a prescribed slip velocity on the boundary propels the squirmer. While the model was originally designed for spherical shapes, it has since been adapted to more general shapes and has recently been shown to capture realistic collective behavior of suspensions (Kyoya et al. 2015).
Shape is also a key parameter in the design of artificial micro-swimmers in applications such as targeted drug delivery. In particular, the squirmer model is often employed to study the propulsion of phoretic particles, which are micro-to nano-meter sized particles that propel themselves by exploiting the asymmetry of chemical reactions on their surfaces (Anderson 1989;Golestanian et al. 2007). A classical example is the Janus sphere (Howse et al. 2007), which consists of inert and catalytic hemispheres. When submerged in a suitable chemical solution, the asymmetry between the chemical reactions on the two hemispheres creates a concentration gradient. The gradient creates an effective steady slip velocity on the surface via osmosis that naturally suits the squirmer model. Besides the classical Janus spheres and bi-metallic nanorods (Paxton et al. 2004), more sophisticated shapes have also been proposed recently, such as two-spheres (Valadares et al. 2010;Palacci et al. 2015), spherocylinder (Uspal et al. 2018), matchsticks (Morgan et al. 2014) and microstars (Simmchen et al. 2017). Interestingly, Uspal et al. (2018) showed that special shapes of phoretic particles exhibit novel properties such as 'edgefollowing' when put close to chemically patterned surfaces.
Studying the efficiency of biological micro-swimmers is pivotal to understanding natural systems and designing artificial ones for accomplishing various physical tasks. The mechanical efficiency (Lighthill 1975) of the spherical squirmer can be directly computed, as its rate of viscous energy dissipation, or power loss, can be written in terms of the modes of the squirming motion. Michelin & Lauga (2010) found the optimal swimming strokes of unsteady spherical squirmers by employing a pseudo-spectral method for solving the Stokes equations that govern the ambient fluid and a numerical optimization procedure. Their approach, however, does not readily generalize to arbitrary shapes. On the other hand, Leshansky et al. (2007) analytically investigated the efficiency of microswimmers of prolate ellipsoidal shapes with a time-independent 'treadmilling' slip velocity and found that the efficiency increases unboundedly with the aspect ratio. Vilfan (2012) optimized the steady slip velocity and the shape at the same time, with constraints on its volume and maximum curvature. However, this work was restricted to the case where the power loss at each point is solely based on the local slip velocity, uncoupled from its neighboring points (e.g., when cilia are far apart).
In this note, we address the following broader question: Given an axisymmetric shape of a steady squirmer, what is the slip velocity that maximizes its swimming efficiency? The optimization problem, being quadratic, is reduced to a linear system of equations solved by a direct method, while forward exterior flow problems are solved using a boundary integral method. Those combined features produce a simple and efficient solution procedure. We introduce the optimization problem and our numerical solver in Section 2, present the optimal solution for various shape families in Section 3, followed by conclusions and a discussion on future research directions in Section 4.

Model
Consider an axisymmetric micro-swimmer whose boundary Γ can be obtained by rotating a curve γ about e 3 axis as shown in Fig. 1(a). Using the arc-length s ∈ [0, ] to parameterize the generating curve, its coordinate functions can be written as γ(s) = (x 1 (s), 0, x 3 (s)). Here, we restrict our attention to shapes of spherical topology, therefore, all shapes considered satisfy the conditions x 1 (0) = x 1 ( ) = 0 and x 1 (s) > 0, ∀ s ∈ (0, ). We assume that the micro-swimmer is suspended in an unbounded viscous fluid domain. The governing equations for the ambient fluid in the vanishing Reynolds number limit are given by the Stokes equations: where µ is the fluid viscosity, p and u are the pressure and flow field respectively. In the absence of external forces and imposed flow fields, the far-field boundary condition simply is A tangential slip u S defined on γ propels the micro-swimmer forward with a translational velocity U in the e 3 direction. Its angular velocity as well as the translational velocities image of a phototactic swimmer, which consists of a haematite particle extruded from a colloidal bead. (Aubret & Palacci 2018) in the e 1 and e 2 directions are zero by symmetry. Consequently, the boundary condition on γ is given by where τ is the unit tangent vector on γ. Note that, in order to avoid singularities, the slip must vanish at the end points: Due to the axisymmetry of Γ , the required no-net-torque condition on the freelysuspended microswimmer is automatically satisfied while the no-net-force condition reduces to one scalar equation where f is the active force density on the micro-swimmer surface (negative to fluid traction) and f 3 is its e 3 component. We quantify the performance of the micro-swimmer with slip velocity u S by its power loss while maintaining a target swimming speed U . The power loss is defined by (2.6) Note that P can be made arbitrarily small by lowering the swimming speed U . It is therefore necessary to compare the power loss of different swimmers that have the same swimming speed U . We note that a lower P with a fixed shape and swimming speed U corresponds to a higher efficiency, η = C D U 2 /P , as defined by Lighthill (1952), where C D is the drag coefficient of the given swimmer.

Boundary integral method for the forward problem
Before stating the optimization problem, we summarize our numerical solution procedure for (2.1) -(2.3). Again, we fix the swimming speed U , referred to from here onwards as the"target swimming speed", and assume that the tangential slip u S is given. In general, an arbitrary pair of u S and U does not satisfy the no-net-force condition (2.5). This condition will be treated as a constraint in our optimization problem. Therefore, the goal is to find the active force density f given the velocity on the boundary γ as in (2.3). We use the single-layer potential ansatz, which expresses the velocity as a convolution of an unknown density function µ with the Green's function for the Stokes equations G, from which the force density can be determined by convolution with the traction kernel T (Pozrikidis (1992)): where n is the unit normal vector pointing into the fluid. We can solve for µ by taking the limit of x → Γ in the above ansatz and substituting in (2.3). The boundary integrals in (2.7) become weakly singular on Γ , requiring specialized quadrature rules. Here, we use the approach of Veerapaneni et al. (2009) which performs an analytic integration in the θ−direction reducing the integrals to convolutions on the generating curve and applies a high-order quadrature rule designed to handle the log−singularity of the resulting kernels.

Optimization problem and its reformulation
The goal is to find a slip profile u S * (s) that minimizes the power loss P while maintaining the target swimming speed U of a given axisymmetrical micro-swimmer. Let J be the objective function, here equated to P defined in (2.6), and F be the net force functional: (2.8) They are slip velocity functionals as their values are completely determined by u S . The optimization problem can now be stated as follows: with U being the space of the all possible slip velocities satisfying (2.4). Notice that the no-net-force condition (2.5) is added as a constraint here. By (2.3) and linearity of the Stokes equation (2.1), the forward solution u and the net force F are affine in u S (u is linear in u S if F = 0). Consequently, J(u S ) is a quadratic functional and (2.9) is inherently a quadratic optimization problem. To make it more explicit, consider a discretized version of the slip optimization problem where u S is sought in the form for some set of m basis functions u S k satisfying (2.4). We adopt a B-spline formulation for these basis functions.Let (u 0 , p 0 , f 0 ) and (u k , p k , f k ) (with 1 k m) denote the solutions of the forward problem (2.1) with u = e 3 and u = u S k τ being their boundary conditions on γ, respectively.
The net force F (u S ) is then given by F (u S ) = 2πU F(ξ), where The elements of the m × m matrix A are given by A kj = γ f k · u S j τ x 1 ds. We have used the fact that γ f 0 · u S k τ x 1 ds = γ f k · e 3 x 1 ds for the linear term by the reciprocal theorem (Happel & Brenner 1973). We note that A is symmetric, also by the reciprocal theorem. Physically speaking, ξ T Aξ represents the scaled power loss of the swimmer being held still with its slip velocity parametrized by ξ, implying that A is positivedefinite. Now, the discretized optimization problem becomes min ξ∈R m J (ξ) subject to F(ξ) = 0. (2.13) Introducing the Lagrangian L(ξ, λ) := J (ξ) − 2λF(ξ), the slip optimization problem is reduced to solving the first-order stationarity equations for L given by (2.14) Note that forming the matrix requires (m + 1) solves of the forward problem (2.1) with appropriate boundary conditions. Since the micro-swimmer is assumed to be rigid, the single layer potential operator S as well as the traction operator T , required for forming A and F , are both fixed for a given shape. Therefore, we only need to form them once.

Results
We tested the convergence of our numerical solvers rigorously; the boundary discretization for all the numerical examples presented here is chosen so that at least 6-digit solution accuracy is attained (determined via self-convergence tests). Here we focus on analysis of the optimal solutions for various micro-swimmer shape families. Let V be the volume enclosed by the swimmer. We normalize lengths by the radius of a sphere of equivalent volume i.e., by R = (3V /4π) 1/3 , and velocities by the swimming speed U . A simple calculation shows that, for a micro-swimmer submerged in water of size R = 5 µm and the speed of one body-length per second, the Reynolds number (Re) ≈ 5 × 10 −5 ; thereby, confirming the validity of the Stokes equation (2.1). We will use the dimensionless reduced volume, defined by ν = 6 √ πV /A 3/2 where A is the surface area of the given shape, to characterize each shape family. The largest possible value of ν, attained by spheres, is ν = 1, while for example ν decreases monotonically for ellipsoids as the aspect-ratio is increased.
We first consider six different micro-swimmer shapes and plot their optimal slip profiles obtained by solving (2.14) in Figure 2. In each case, we also show the flow fields in both the body and lab frames. The optimal slip velocities plotted against the arclength, measured from north pole to south pole, are shown in the insets. In the case of a sphere (Fig. 2(a)), we recover the standard result that the optimal profile is a sine curve (Michelin & Lauga 2010). The optimal slip velocity of the prolate swimmer, shown in Fig. 2(b), 'flattens' the sine curve in the middle while that of the oblate swimmer, shown in Fig. 2(c), 'pinches' the sine curve. Additionally, the peak value of the optimal slip velocity is low for the prolate swimmer, and high for the oblate swimmer, compared to the spherical swimmer.
Next, we consider three shapes corresponding to different shape families. In Fig. 2(d), we consider the 'wavy' configuration obtained by adding high-order axisymmetric modes to the spherical shape. The optimal slip velocity follows the general trend for that of (a), while lower slip velocities are observed at the troughs, qualitatively consistent to those obtained in Vilfan (2012). The spherocylinder (Fig. 2(e)) resembles closely the prolate ellipsoid of Fig. 2(b) with the same aspect ratio, its optimal slip velocity being nearly the same (albeit with a slightly narrower plateau and higher peak slip velocity). Finally, we investigate the optimal slip velocity of the stomatocyte shape ( Fig. 2(f)), which is the only non-convex shape among those considered here. Similar to that of the oblate swimmer, the general slip velocity is like a pinched sine wave. However, one distinguishing feature is that slip velocity is nearly zero over part of its surface, namely the cup-like region in its posterior.
We note that judging by the flow field, all optimal swimmers studied here appear to be neutral swimmers (as opposed to pushers or pullers), even though some of them (e.g., the stomatocyte) have front-back asymmetric shapes. Additionally, the optimal slip velocity is proportional to the target swimming speed U due to linearity of the Stokes equations. As a consequence, while the results only showcase micro-swimmers propelling themselves in the positive e 3 direction, the optimal solution u S * for swimming in the opposite direction is merely a change of sign.
Next, we study the optimal active force density f corresponding to the same shapes. Its normal and tangential components are plotted in Fig. 3. We note that by the no-net-force condition (2.5), the power loss reduces to P = 2π γ f · (u S τ )x 1 ds, implying that only the tangential component contributes to the power loss. The change in tangential forces as a function of arclength loosely resembles that of the optimal slip velocity, mediated by the local curvature along the generating curve. Qualitatively, a low local curvature suppresses the traction and a high local curvature amplifies it. Slip velocities scaled by their local curvatures are shown in black dotted curves for a reference. In Fig. 4, we plot the minimal power loss as a function of the reduced volume for various shape families, scaled by that of the spherical swimmer with the same volume. The minimal power loss for prolate ellipsoids monotonically decreases as the shape gets more slender; in contrast, it is well-known that the shape with the minimal fluid drag is one with approximately 2:1 aspect ratio (Pironneau 1973). By slender body theory, the power loss of a prolate ellipsoid scales as ∼ µα 2/3 U 2 , where α is the aspect ratio (see Leshansky et al. (2007)). On the other hand, the minimal power loss for oblate ellipsoids grows rapidly as the reduced volume is increased. Shapes of the spherocylinder family behave similarly to the prolate ellipsoids, and converge to the spherical case when the length of the cylinder reduces to 0, as expected. It is however worth pointing out that spherocylinder costs more power loss than prolate ellipsoids with the same reduced volume; this relates to the fact that the peak slip velocity for spherocylinder is higher than that of the prolate ellipsoid ( Fig. 2 (b)&(e)). The stomatocyte family is constructed by 'pulling' the rim of the shape, effectively making the shape 'taller' and curls deeper and deeper inside. We find that 'taller' shapes requires lower power loss for this shape family, which is qualitatively consistent with the ellipsoid family. Finally, we note that the power loss of the snowman family (two spheres attaching with each other) is quite robust to the relative sizes of the two spheres. The power loss is only about 25% higher than that of a single sphere in the limit case where the two spheres are of the same size.
A few other examples that take more generic shapes are also shown in Figure 4. The optimal slip velocities are colored on their surfaces while their power loss is shown in the form of scatter points. The generating curve of these shapes are formed by spherical harmonics. We note that the optimal performance of shapes that appear similar can be very different. For example, the difference in power loss between examples 6 and 8 is about 150% of the spherical swimmer, or 60% of example 6. This result is a strong indicator that the slip velocity of the artificial swimmer, as well as its shape, must be carefully designed to achieve good performance. We note that the minimal power loss for all the shape familites considered here are bounded from below by the curve for prolate ellipsoids, agreeing with the analysis made by Leshansky et al. (2007).

Conclusions
In this work, we provided a solution procedure for the PDE-constrained optimization problem of finding the optimal slip profile on an axisymmetric micro-swimmer that minimizes the power loss required to maintain a target swimming speed. While it can be extended to other objective functions, we exploited the quadratic nature of the power loss functional in the control parameters to simplify and streamline the solution procedure. In the general case, an adjoint formulation and iterative optimization algorithms can be employed. Regardless of the formulation, however, the use of boundary integral method to solve the Stokes equations greatly reduces the computational cost due to dimensionality reduction. Solving any of the examples presented in this work, for example, required only a few seconds on a standard laptop. Extending our procedure to fully three-dimensional (non-axisymmetric) shapes is straightforward; the key technical challenge is incorporating a high-order boundary integral solver, for which open-source codes are now available (e.g., see Gimbutas & Veerapaneni (2013)).
The optimal slip velocities of different families of bio-relevant shapes and common phoretic particle shapes were studied and compared against idealized shapes. We found that optimal slip velocity for all shapes, even with front-back asymmetric shapes, could be viewed as neutral swimmers. This mirrors the findings in Stone & Samuel (1996) that swimmers that create no vorticity in the surrounding fluid spend less energy than those which do. Additionally, the power loss of micro-swimmers seem to be bounded from below by prolate ellipsoids with the same reduced volume.
The optimization procedure developed in this work can directly be employed in the design pipeline of autophoretic particles. For example, in the case of diffusiophoresis, the computed optimal slip profile for a given shape can be used to formulate the chemical coating pattern of the phoretic particles. Another natural extension of this work is to relax the steady slip assumption and consider time-periodic squirming motion as done in Michelin & Lauga (2010). This would be particularly useful for studying the ciliary locomotion of micro-organisms with arbitrary shapes. Furthermore, building on the recent work of Bonnet et al. (2020), we are developing solvers for the shape optimization problem of finding the most efficient microswimmer shapes under specified area, volume or other physical constraints.
Authors gratefully acknowledge support from NSF under grants DMS-1719834 and DMS-1454010.