Optimal Ciliary Locomotion of Axisymmetric Microswimmers

Many biological microswimmers locomote by periodically beating the densely-packed cilia on their cell surface in a wave-like fashion. While the swimming mechanisms of ciliated microswimmers have been extensively studied both from the analytical and the numerical point of view, the optimization of the ciliary motion of microswimmers has received limited attention, especially for non-spherical shapes. In this paper, using an envelope model for the microswimmer, we numerically optimize the ciliary motion of a ciliate with an arbitrary axisymmetric shape. The forward solutions are found using a fast boundary integral method, and the efficiency sensitivities are derived using an adjoint-based method. Our results show that a prolate microswimmer with a 2:1 aspect ratio shares similar optimal ciliary motion as the spherical microswimmer, yet the swimming efficiency can increase two-fold. More interestingly, the optimal ciliary motion of a concave microswimmer can be qualitatively different from that of the spherical microswimmer, and adding a constraint to the ciliary length is found to improve, on average, the efficiency for such swimmers.

In this paper, we study the (hydrodynamic) swimming efficiency of ciliated microswimmers of an arbitrary axisymmetric shape. Specifically, the swimming efficiency is understood as the ratio between the 'useful power' against the total power. The useful power could be computed as the power needed to drag a rigid body of the same shape as the swimmer with the swim speed while the total power is the rate of energy dissipation through viscous stresses in the flow to produce this motion [Lighthill, 1952]. The goal of this paper is to find the optimal ciliary motion that maximizes the swimming efficiency for an arbitrary axisymmetric microswimmer.
In particular, in a seminal work, Michelin and Lauga [2010] studied the optimal beating stroke for a spherical swimmer using the envelope model. Specifically, the material points on the envelope were assumed to move tangentially on the surface in a time-periodic fashion, hence the swimmer retains the spherical shape. The flow field, power loss, swimming efficiency as well as their sensitivities, thereby, were computed explicitly using spherical harmonics. Their optimization found that the envelope surface deforms in a wavelike fashion, which significantly breaks the time-symmetry at the organism level similar to the metachronal waves observed in biological microswimmers.
Since most biological microswimmers do not have spherical shapes, there is a need for extending the previous work to more general geometries. Such an extension, however, is hard to carry out using semianalytical methods. Therefore, in this paper, we develop a computational framework for optimizing the ciliary motion of a microswimmer with arbitrary axisymmetric shape. We employ the envelope model, wherein, the envelope is restricted to move tangential to the surface so the shape of the microswimmer is unchanged during the beating period. We use a boundary integral method to solve the forward problem and derive an adjoint-based formulation for solving the optimization problem.
The paper is organized as follows. In Section 2, we introduce the optimization problem, derive the sensitivity formulas and discuss our numerical solution procedure. In Section 3, we present the optimal unconstrained and constrained solutions for microswimmers of various shape families. Finally, in Section 4, we discuss our conclusions and future directions.

Model
Consider an axisymmetric microswimmer whose boundary Γ is obtained by rotating a generating curve γ of length about e 3 axis, as shown in Figure 1(a). We adopt the classic envelope model [Lighthill, 1952] and assume that the ciliary tips undergo time-periodic tangential movements along the generating curve. Let s = α(s 0 , t) be the ciliary tip's arclength coordinate on the generating curve γ at time t for a cilium rooted at s 0 . The tangential slip velocity of this material point in its body-frame is thus u S (s, t) = u S (α(s 0 , t), t) = ∂ t α(s 0 , t). (1) In addition to the time-periodic condition, the ciliary motion α needs to satisfy two more conditions to avoid singularity [Michelin and Lauga, 2010]. First, the slip velocities should vanish at the poles α(0, t) = 0 and α( , t) = , ∀ t ∈ R + ,  Figure 1: (a) Schematic of the microswimmer geometry. The shape is assumed to be axisymmetric, obtained by rotating the generating curve γ about the e 3 axis. The tip of the cilium rooted at s 0 at time t is given by s = α(s 0 , t). (b) Illustration of the algorithm for computing the slip velocity at the quadrature points u S (s q , t). We first compute the "tip" position and the corresponding tip velocities (open blue circles) of cilia rooted at the N q quadrature points s q (closed blue circles). We then obtain the slip velocities at sample points uniformly distributed along the generating curve (open red squares) by a cubic interpolation. The slip velocity at any arclength (black curve) are then obtained by a high-order B-spline interpolation from the sample points. We have reduced the number of quadrature and sample points in this figure (compared to values used in the numerical experiments) to avoid visual clutter. and second, α should be a monotonic function, that is, The last condition ensures the slip velocity is unique at any arclength s; in other words, crossing of cilia is forbidden. While in reality, cilia do cross, this condition is enforced to ensure validity of the continuum model. In the viscous-dominated regime, the flow dynamics is described by the incompressible Stokes equations at every instance of time − µ∇ 2 u + ∇p = 0, ∇ · u = 0, where µ is the fluid viscosity, p and u are the fluid pressure and velocity fields respectively. In the absence of external forces and imposed flow field, the far-field boundary condition is simply The free-swimming microswimmer also needs to satisfy the no-net-force and no-net-torque conditions. Owing to the axisymmetric assumption, the no-net-torque condition is satisfied by construction, and the no-net-force condition is reduced to one scalar equation where x 1 is the e 1 component of x, f is the active force density the swimmer applied to the fluid (negative to fluid traction) and f 3 is its e 3 component. Given any ciliary motion α(s 0 , t) that satisfies (2) & (3), there is a unique tangential slip velocity u S (s, t) defined by (1). Such a slip velocity propels the microswimmer at a translational velocity U (t) in the e 3 direction, determined by (6). Its angular velocity as well as the translational velocities 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 γ. Thereby, the instantaneous power loss P (t) can be written as The second term on the right-hand-side is zero provided that the no-net-force condition (6) is satisfied. Following Lighthill [1952], we quantify the performance of the microswimmer by its swimming efficiency , defined as where P = P (t) and U = U (t) are the instantaneous power loss and swim speed, · denotes the time-average over one period, and C D is the drag coefficient defined as the total drag force of towing a rigid body of the same shape at a unit speed along e 3 direction. The coefficient C D depends on the given shape γ only; for example, C D = 6πµa in the case of a spherical microswimmer with radius a.
In our simulations, we normalize the radius of the microswimmer to unity, and the period of the ciliary motion to 2π. It is worth noting that the swimming efficiency (9) is size and period independent, thanks to its dimensionless nature. The Reynolds number of a ciliated microswimmer of radius 100µm and frequency 30Hz submerged in water can be estimated as Re ∼ 10 −4 , confirming the applicability of Stokes equations.

Numerical algorithm for solving the forward problem
Before stating the optimization problem, we summarize our numerical solution procedure for the governing equations (4) -(7). By the quasi-static nature of the Stokes equation (4), the flow field u(x, t) can be solved independently at any given time, and the time-averages can be found using standard numerical integration techniques (e.g., trapezoidal rule). Here we adopt a boundary integral method (BIM) at every time step. A similar BIM implementation was detailed in our recent work Guo et al. [2021] which studied the optimization of time-independent slip profiles. The main procedures are summarized below.
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: The force density can then be evaluated as a convolution of µ with the (negative of) traction kernel: We convert these weakly singular boundary integrals into convolutions on the generating curve γ by performing an analytic integration in the orthoradial direction, and apply a high-order quadrature rule designed to handle the log−singularity of the resulting kernels [Veerapaneni et al., 2009]. The Stokes flow problem defined at any time t by equations (4) -(7) is then recast as the BIM system for the unknowns µ and U (t) obtained by substituting (10) in (7) and (11) in (6). The numerical solution method consists in discretizing γ into N p non-overlapping panels, each panel supporting the nodes of a 10-point Gaussian quadrature rule. The single-layer operator is approximated in Nyström fashion, by collocation at the N q = 10N p quadrature nodes, while the values of µ are sought at the same quadrature nodes. The resulting BIM system is where the vectors µ = µ(s q , t) and u S = u S (s q , t) are the unknown density and the given slip velocity at all quadrature nodes s q , S is the axisymmetric single-layer potential operator (which is fixed for a given shape γ), B is the column vector reproducing e 3 at each quadrature node, C is the row vector such that C[µ] = Γ f (x) · e 3 dΓ is the total traction force in the e 3 direction. The algorithm to obtain the slip velocity at the quadrature nodes at a given time u S (s q , t) is summarized in Figure 1(b). Specifically, we start by computing the corresponding ciliary tip position s = α(s q , t) and the slip velocity u S (s, t) from (1). These tip positions s can be highly nonuniform, depending on the form of α, which could be difficult for the forward solver. To circumvent this difficulty and to find a smooth representation of the slip velocities on the quadrature points, we first find the slip velocities at N s sample points uniformly distributed along the generating curve by interpolating u S (s, t) (we use the routine PCHIP in MATLAB); the slip velocities at the quadrature nodes u S (s q , t) are then in turn interpolated from the N s sample points using high-order B-spline bases. An alternative approach could be to follow the position and the slip velocity of each material point. In other words, one can use u S (s, t) directly on the right-hand-side of (12), which will bypass the interpolation steps mentioned above. However, it requires re-assembly of the matrix S at every time step, significantly increasing the computational cost.

Optimization problem
The goal of this work is to find the optimal ciliary motion for a given arbitrary axisymmetric shape, that is, the ciliary motion α (s 0 , t) that maximizes the swimming efficiency : where A is the space of all possible time-periodic ciliary motion satisfying (2) & (3). It is, however, not easy to define and manipulate finite-dimensional parametrizations of α that remain in that space. To circumvent this difficulty, we follow the ideas in Michelin and Lauga [2010] and represent α in terms of a time-periodic function ψ(x, t), such that where is the total length of the generating curve γ. Note that α is also (implicitly) a function of time t, through ψ = ψ(x, t). It is easy to verify that α given by (14) satisfies the boundary conditions (2) and the monotonicity requirement (3) for any choice of ψ. Conversely, for any α satisfying (2) and (3), there is at least one ψ that provides α. As a result, the optimization problem is recast as finding where ψ(·, t) is only required to be square-integrable over [0, ] for any t.
We use a quasi-Newton BFGS method [Nocedal and Wright, 2006] to optimize the ciliary motion via ψ, which requires repeated evaluations of efficiency sensitivities with respect to perturbations of ψ. The sensitivities of power loss and swim speed are derived using an adjoint-based method, while the efficiency sensitivity is found using the quotient rule thereafter. The adjoint-based method exhibits a great advantage against the traditional finite difference method when finding the sensitivities, as regardless of the dimension of the parameter space, the objective derivatives with respect to all design parameters can here be evaluated on the basis of one solve of the forward problem for each given ciliary motion α. The derivations are detailed below.

Sensitivity analysis
We start by finding the sensitivities in terms of the slip profile u S . The sensitivities in terms of the auxiliary unknown ψ will be found subsequently by a change of variable. As the concept of adjoint solution in general rests on duality considerations, we recast the forward flow problem in weak form for the purpose of finding the sought sensitivities of power loss and swim speed, even though the numerical forward solution method used in this work does not directly exploit that weak form. Specifically, we recast the forward problem (4) - (7) in mixed weak form (see, e.g., Brezzi and Fortin [1991, Chap. 6 where the bilinear forms a and b are defined by and D[u] := (∇u + ∇ T u)/2 is the strain rate tensor. ·, · Γ is a short-hand for the inner product on Γ. For Similarly, with a slight abuse of notation, the power loss functional could be written as The Dirichlet boundary condition (7) is (weakly) enforced explicitly through (16 b), rather than being embedded in the velocity solution space V, as this will facilitate the derivation of slip derivative identities; this is in fact our motivation for using the mixed weak form (16). Condition (16 c) is the no-net-force condition (6).
First-order sensitivities of functionals at u S are defined as directional derivatives, by considering perturbations of u S of the form for some ν in the slip velocity space and η ∈ R. Then, the directional (or Gâteaux) derivative of a functional J(u S ) in the direction ν, denoted by J (u S ; ν), is defined as For the power loss functional, we obtain (since the derivative of u S in the above sense is ν) where f and U are the derivatives of the active force f and swim speed U solving problem (16), considered as functionals on the slip velocity u S : Differentiating the weak formulation (16) of the forward problem with respect to u S leads to the weak formulation of the governing problem for the derivatives (u , f , p , U ) of the solution (u, f , p, U ) Here we have assumed without loss of generality that the test functions in (16) verify v = 0, g = 0, and q = 0, which is made possible by the absence of boundary constraints in V. At first glance, evaluating P (u S ; ν) in a given perturbation ν appears to rely on solving the derivative problem (22). However, a more effective approach allows to bypass the actual evaluation of f . Let the adjoint problem be defined by i.e. (û,p) are the flow variables induced by prescribing a unit velocity e 3 on Γ. For later convenience, we let F 0 denote the (nonzero) net force exerted on Γ by the adjoint flow: Problem (23) in strong form is defined by equations (4) - (7) with U = 1, u S = 0. In fact, F 0 takes the same value as the drag coefficient C D in (9). Then, combining the derivative problem (22) with the forward problem (16) or the adjoint problem (23) with appropriate choices of test functions allows to derive expressions of P (u S ; ν) and U (u S ; ν) which do not involve the forward solution derivatives.
The derivation and the explicit expression of each term in (33) are given in the Appendix. Finally, the efficiency sensitivity in terms of ψ readily follows by the quotient rule

Constraints on surface displacement
The unconstrained optimization problem (15) introduced above has the tendency to converge to unphysical/unrealistic strokes, where each cilium effectively 'covers' the entire generating curve. For a more realistic model, we should add a constraint on the length of the cilium. To this end, and again following Michelin and Lauga [2010], we replace the initial unconstrained optimization problem (15) with the penalized optimization problem where the (non-negative) penalty term C(ψ), defined as serves to incorporate the kinematical constraint A(ψ) ≤ c in the optimization problem. The functional A(ψ) in (36) is a measure of the amplitude of the displacement of individual material points for the stroke (through α), and c is a threshold parameter to bound A(ψ) (a smaller c corresponding to a stricter constraint). H is a smooth non-negative penalty function defined by which for large enough Λ 2 approximates u → 2Λ 1 u 2 Y (u) (Y being the Heaviside unit step function). The multiplicative parameter Λ 1 then serves to tune the severity of the penalty incurred by violations of the constraint A(ψ) ≤ c. We use Λ 1 = 10 4 and Λ 2 = 10 4 in our numerical simulations unless otherwise mentioned. The optimization results are not sensitive to the choice of Λ 1 and Λ 2 . A small caveat of the penalty function (37) is that it has a (small) bump at Λ 2 u ≈ −1.109. This bump would occasionally trap the optimizations into local extrema that have significantly lower efficiencies, depending on the initial guesses. Perturbing Λ 2 for such cases helps to alleviate the problem. The physically most relevant definition of A would be the actual displacement amplitude of an individual point, i.e., ∆s = [α max (s 0 )−α min (s 0 )]/2. The strong nonlinearity of this measure, however, is not appropriate for the computation of the gradient. Following Michelin and Lauga [2010], we measure the displacement by its variance in time: The maximum displacement ∆s max = max s0 (∆s) will be found post-optimization for the optimal ciliary motion α to better illustrate our results in Section 3. Like the initial problem (15), the penalized problem (35) is solvable using unconstrained optimization methods, and we again adopt a quasi-Newton BFGS algorithm to optimize the ciliary motion. Applying the chain rule to the penalty functional C(ψ), we obtain the derivative of the penalty term in the direction ofψ as The derivative of the penalized objective functional E(ψ) is therefore where and C are given by equations (34) and (39), respectively.
3 Results and discussion

Parametrization
We parametrize ψ(s 0 , t) such that where B k are the 5th order B-spline basis functions and their coordinates ξ k (t) are expanded as trigonometric polynomials ξ k (t) = a 0k /2 + n j=1 [a jk cos jt + b jk sin jt] to ensure time-periodicity. Taken together, we have ψ(s 0 , t) = m k=1 a 0k 2 + n j=1 (a jk cos jt + b jk sin jt) B k (s 0 ) so that the finite-dimensional optimization problem seeks optimal values for the m(2n + 1) coefficients a 0k , a jk and b jk . The initial guesses are chosen to be low frequency waves with small wave amplitudes. To obtain such initial waves, the coefficients of the zeroth Fourier mode a 0k /2 are randomly chosen from a uniform distribution within [0, 1], the first Fourier modes a 1k and b 1k are randomly chosen from a uniform distribution within [0, 0.01], and the coefficients for higher order Fourier modes j > 1 are set to 0. To evaluate the gradient of E(ψ) with respect to the design parameters a 0k , a jk and b jk , we use (40) withψ taken as the basis functions of the adopted parametrization (42), i.e.ψ(s 0 , t) = B k (s 0 )/2,ψ(s 0 , t) = B k (s 0 ) cos jt andψ(s 0 , t) = B k (s 0 ) sin jt, respectively. In terms of parametrization, local minima are multiple in the parameter space, since multiplying optimal parameters by a constant factor yields the same optimum for α.

Spheroidal swimmers
By way of validation, we start with optimizing the ciliary motion of a spherical microswimmer. The efficiency as a function of iteration number for the unconstrained optimization (15) is shown in Figure 2(a) in blue. The maximum efficiency is about 35%. The ciliary motion of the optimal spherical microswimmer is shown in Figure 2(b). Each curve follows the arclength coordinate of a cilium tip over one period. We observe, similar to the results of Michelin and Lauga [2010], clearly distinguished strokes within the beating period. In particular, cilia travel downward 'spread out' during the effective stroke (corresponding to a stretching of the surface), but travel upward 'bundled' together during the recovery stroke in a shock-like structure (corresponding to a compression of the surface). This type of waveform is known as an antiplectic metachronal wave [Knight-Jones, 1954, Blake, 1972. We note that this optimal ciliary motion produces an efficiency higher than the 23% efficiency obtained numerically by Michelin and Lauga [2010, Fig. 11]. This is due to a larger maximum displacement ∆s max ≈ 0.45 in our optimizations (translated to a maximum angle of 81 degrees vs 53 degrees). Our optimization result aligns well with their results using the analytical ansatz [Michelin and Lauga, 2010, Fig. 14]. Additionally, we found that increasing the number of Fourier mode n increases the maximum displacement as well as the efficiency; the optimal ciliary motion of higher n also exhibits a higher slope for the shock-like structures (results not shown here). This is again consistent with their analytical ansatz, which shows that the efficiency approaches 50% in the limit of the maximum displacement approaches 90 degrees, and the corresponding 'width' of the shock in this limit is infinitely small. The mean slip velocity of the Eulerian points within each period are almost identical to the optimal time-independent slip velocity scaled by the swim speed, as shown in Figure 2(d).
The optimal unconstrained prolate spheroidal microswimmer with a 2:1 aspect ratio has an efficiency ≈ 69%, about twice as high as the spherical microswimmer as shown in Figure 2(a). This roughly two-fold increase in efficiency is also observed in the optimal time-independent microswimmers [Guo et al., 2021]. The optimal ciliary motion is very close to that of the spherical swimmer ( Fig. 2(b)&(c)) , while the mean slip velocity of the Eulerian points are between the optimal time-independent slip velocity of the same shape and those of a spherical swimmer, as shown in Figure 2(e). As a sanity check, swapping the ciliary motions obtained from optimizing the spherical swimmer and the prolate swimmer leads in both cases to lower swimming efficiencies. Specifically, a spherical swimmer with the ciliary motion shown in Figure 2(c) has 34% swimming efficiency and a prolate swimmer with the ciliary motion shown in Figure 2(b) has 65% swimming efficiency (compared to 35% and 69% using the 'true' optimal ciliary motions, respectively).
We then turn to the case in which the cilia length is constrained by prescribing a bound on the displacement variance (38). We control the maximum variance by tuning c in (36), and the efficiencies are plotted against the maximum displacement ∆s max scaled by the total arclength in Figure 3. Three different random initial guesses are used for each c. The unconstrained optimization results for the spherical and prolate spheroidal swimmers are also shown in the figure for reference. Notably, for both the unconstrained swimmers, the length of the cilia is roughly half the total arclength of the generating curve (∆s max ≈ /2). In other words, a cilium rooted at the equator would be able to get very close to both poles during the beating cycle. In general, a smaller variance (tighter constraint) leads to a lower efficiency, as expected. The efficiency results of spherical microswimmers closely match those reported by Michelin and Lauga [2010]. The efficiencies of the prolate spheroidal microswimmer under constraints are also shown in Figure 3. Similar to the spherical microswimmer, the efficiency increases roughly linearly with the scaled cilia length ∆s max / , and converges to the kinematically unconstrained optimal microswimmer as the maximum variance c is increased. It is noteworthy that adding a constraint in the cilia length not only limits the wave amplitudes, but also breaks the single wave with larger amplitude into multiple waves with smaller amplitudes (Fig. 4(a)), which resemble the metachronal waves of typical ciliated microswimmers such as Paramecium. More interestingly, the mean slip velocity in the constrained case can be qualitatively different from the time-independent optimal slip velocity, as shown in Figure 4(b). In particular, the mean slip velocity around the equator is significantly higher than the time-independent slip velocity, while the mean slip velocity near the poles are closer to zero. This can be inferred from the ciliary motions, as the cilia only move slightly near the poles, whereas multiple waves with significant amplitudes travel around the equator within one period.

Non-spheroidal swimmers
We then investigate the effects of shapes on the optimal ciliary motions and the swimming efficiencies. In particular, we examine whether a single wave travelling between north and south poles always maximizes the swimming efficiency, and whether adding a constraint in the cilia length is always detrimental to the swimming efficiency.
We consider a family of shapes whose generating curves are given by: (x, z) = (R(θ) sin θ, R(θ) cos θ), where R(θ) = (1 + δ cos 2θ) is a function that makes the radius non-constant, and θ ∈ [0, π] is the parametric coordinate. For 0 < δ < 1, the radius is the smallest at θ = π/2, corresponding to a 'neck' around the equator. In the limit δ = 0, the generating curve reduces to a semicircle and the swimmer reduces to the spherical swimmer.
The optimization results are depicted in Figure 5 for 0 ≤ δ ≤ 0.8. Some corresponding shapes are  Figure 4: Ciliary motion (a) and mean slip velocity (b) for the optimal spherical swimmer with constraint (∆s max / ≈ 5.0%). The efficiency is ≈ 6.9%, and the swim speed is U ≈ 0.091. The swimmer forms multiple waves in the equatorial region, leading to a high slip velocity at s ≈ 0.5 . The motion close to the poles is nearly zero. The dashed curve in (b) is the time-independent optimal slip velocity of the spherical swimmer, scaled by the swim speed. The video of the optimal ciliary motion can be found in the online supplementary material (Movie 3). Blue dashed curves are the optimal time-independent slip velocities scaled by the swim speed. In these simulations, we increase the number of panels N p = 40 to resolve the sharp shape change. The videos of the optimal ciliary motions can be found in the online supplementary material (Movie 4 & 5) shown as insets. The median efficiencies of ten Monte Carlo simulations are plotted for each δ value, and compared against the time-independent efficiencies. For all three cases (constrained, unconstrained, and time-independent), the efficiencies increase as δ increases from 0 to 0.3. This is because increasing δ in this regime makes the shape more elongated. Increasing δ further reduces the efficiencies as the 'neck' at the equator becomes more and more pronounced. Additionally, the unconstrained microswimmers, on average, have better efficiencies than the microswimmers with kinematic-constraints for 0 ≤ δ ≤ 0.6.
Interestingly, unconstrained optimization may result in worse ciliary motions on average when the shape is highly curved, compared to its kinematically-constrained counterpart. Specifically, the constrained microswimmers have higher median efficiencies for δ ≥ 0.7. We note that the unconstrained optimizations are likely to be trapped in local optima where the ciliary motion forms a single wave ( Fig. 5(b)), whereas the constrained optimizations are 'forced' to find the ciliary motion with multiple waves split at the equator (Fig. 5(c)), because of the constrained cilia length. Additionally, our numerical results show that a single wave travelling between the north and south poles is not as efficient as two separate waves travelling within each hemisphere for this shape. Figures 5(d)&(e) show that the single wave generates a high mean slip velocity at the position where the generating curve bends inward (the equator), whereas the two separate waves generate a mean slip velocity similar to that obtained from the time-independent optimization. In a way, the constraint in cilia length is helping the optimizer to navigate the parameter space.
To better understand the effects of constraints on the highly curved shapes, we present the statistical results of the thin neck microswimmer (δ = 0.8) with various constraints in Figure 6. In general, the highest efficiency from the Monte Carlo simulations increases with the constraint for c ≤ 0.8, similar to the case of spheroidal swimmers (Figure 3). Keep increasing c has limited effect on the highest efficiencies, indicating that the constraint is no longer limiting the optimal ciliary motion. The median efficiencies (red horizontal lines), on the other hand, decreases with the constraint if c ≥ 0.8, consistent with the observation from Figure 5. It is worth noting that the constrained optimization is more likely to get stuck in very low efficiencies (e.g., the lowest outlier for c = 0.8), possibly due to the secondary bump of the penalty function C mentioned earlier.
All data points from the optimization are plotted in Figure 6(b) as function of the maximum displacement ∆s max . The efficiencies grow almost linearly until ∆s max ≈ 0.25 , as in the case of spheroidal swimmers, and decrease for larger ∆s max . This is another evidence that the optimal ciliary motion for this shape consists of two separate waves traveling within each hemisphere. We want to emphasize that unconstrained optimization can still reach the optimal ciliary motion, as shown in the box of c = ∞. However it is more likely to reach the sub-optimal ciliary motion compared to the constrained cases.

Conclusions and Discussions
In this paper, we extended the work of Michelin and Lauga [2010] and studied the optimal ciliary motion for a microswimmer with arbitrary axisymmetric shape. In particular, the forward problem is solved using a boundary integral method and the sensitivities are derived using an adjoint-based method. The auxiliary function ψ is parametrized using high-order B-spline basis functions in space and a trigonometric polynomial in time. We studied the constrained and unconstrained optimal ciliary motions of microswimmers with a variety of shapes, including spherical, prolate spheroidal, and concave shapes which are narrow around the equator. In all cases, the optimal swimmer displays (one or multiple) traveling waves, reminiscent of the typical metachronal waves observed in ciliated microswimmers. Specifically, for the spherical swimmer with limited cilia length ( Fig. 4(a)), the ratio between the metachronal wavelength close to the equator and the cilia length could be estimated as λ M W /∆s max ≈ 0.2 /0.05 = 4. This ratio lies in the higher end of the data collected in Velho Rodrigues et al. [2021, Table 9] for biological ciliates, which reports ratio ranging between 0.5 to 4. Our slightly high ratio estimate may not be surprising after all, as the envelope model prohibits the crossing between neighboring cilia. We showed that the optimal ciliary motions of prolate microswimmer with a 2:1 aspect ratio are very close to the ones of spherical microswimmer, while the swimming efficiency can increase two-fold. The mean slip velocity of unconstrained microswimmers also tend to follow the optimal time-independent slip velocity, which can be easily computed using our recent work [Guo et al., 2021].
Most interestingly, we found that constraining the cilia length for some shapes may lead to a better efficiency on average, compared to the unconstrained optimization. It is our conjecture that this counterintuitive result is because the constraint effectively reduces the size of the parameter space, hence lowering the probability of being trapped in local optima during the optimization. Although the concave shapes studied in Section 3.3 are somewhat non-standard, they allows us to gain insights into the effect of local curvature on optimal waveform. Incidentally, these shapes are also observed for ciliates in nature (e.g. during the cell division process).
It is worth pointing out that works on sublayer models (explicitly modeling individual cilia motions) have reported swimming or transport efficiencies in the orders of 0.1 ∼ 1% (see, e.g., Elgeti and Gompper [2013], Ito et al. [2019], Omori et al. [2020]), much lower than the optimal efficiency reported here and others using the envelope models. This large difference can possibly be attributed to the fact that the envelope model we adopted here considers only the energy dissipation outside the ciliary layer (into the ambient fluid), while sublayer models in general considers energy dissipation both inside and outside the ciliary layer. Research has shown that the energy dissipation inside the layer could be as high as 90 ∼ 95% of the total energy dissipation, due to the large shear rate inside the layer (see, e.g., Keller and Wu [1977], Ito et al. [2019]). We note that it is possible to incorporate energy dissipation inside the ciliary layer in the envelope model, as previously done in Vilfan [2012], albeit for a time-independent slip profile. Additionally, the difference could also be due to modeling assumptions on the cilia length and the number of cilia. In particular, the cilia length considered in sublayer models are usually below 1/10 of the body length. Omori et al. [2020] showed that the swimming efficiency increases with the cilia length as fast as powers of 3 in the short cilia limit, and the number of cilia also has a significant positive effect on the swimming efficiency (the envelope model assumes a ciliary continuum). Factoring all three factors (energy inside/outside, cilia length, number of cilia) could bridge the gap between the results obtained from these two types of models.
It is without a doubt that maximizing the hydrodynamic swimming efficiency is not the sole objective for biological microswimmers. Other functions such as generating feeding currents [Riisgård andLarsen, 2010, Pepper et al., 2013] and creating flow environment to accelerate mixing for chemical sensing [Supatto et al., 2008, Shields et al., 2010, Nawroth et al., 2017 are also important factors to consider as a microswimmer. The effect of such multi-tasking on the ciliary dynamics is not well understood. Nevertheless, our work provides an efficient framework to investigate the hydrodynamically optimal ciliary motions for microswimmers of any axisymmetric shape, and could provide insights into designing artificial microswimmers.
This completes our derivation of (33). For the ciliary motion (14) used here, introducing the shorthand notation I(f, g; s) := s 0 f (x)g(x)dx, we