1. Introduction
Many swimming microorganisms propel themselves by periodically beating the active slender appendages on their cell surfaces. These slender appendages are known as cilia or flagella depending on their lengths and distribution density. Eukaryotic flagella, such as those found in mammalian sperm cells and algae cells, are often found in small numbers, whereas ciliated swimmers such as Paramecium and Opalina present many hundreds of cilia densely packed on the cell surfaces (Brennen & Winet Reference Brennen and Winet1977; Witman Reference Witman1990). Besides the locomotion function for microswimmers, cilia inside mammals serve various other functions, such as mucociliary clearance in the airway systems and transportation of egg cells in fallopian tubes (see Satir & Christensen (Reference Satir and Christensen2007), and references therein). Cilia are also found to be critical in transporting cerebrospinal fluid in the third ventricle of the mouse brain (Faubel et al. Reference Faubel, Westendorf, Bodenschatz and Eichele2016) and in creating active flow environments to recruit symbiotic bacteria in a squid-vibrio system (Nawroth et al. Reference Nawroth, Guo, Koch, Heath-Heckman, Hermanson, Ruby, Dabiri, Kanso and McFall-Ngai2017).
Owing to the small length scale of cilia, the typical Reynolds number is close to zero. In this regime, inertia is negligible and the dynamics is dominated by viscous effects. As a result, many effective swimming strategies familiar to our everyday life become futile. For example, waving a rigid tail back-and-forth will not generate any net motion over one period. This is known as time reversibility, or the ‘scallop theorem’, which states that a reciprocal motion cannot generate net motion (Purcell Reference Purcell1977). Microswimmers therefore need to go through non-time-reversible shape changes to overcome and exploit drag (Lauga & Powers Reference Lauga and Powers2009).
Ciliated microswimmers break the time reversibility on two levels. On the individual level, each cilium beats in an asymmetric pattern: during the effective stroke, the cilium pushes the fluid perpendicular to the cell surface like a straight rod, and then moves almost parallel to the cell surface in a curly shape during the recovery stroke, in preparation for the next effective stroke. On the collective level, neighbouring cilia beat with a small phase difference that produces travelling waves on the cell surface, namely the metachronal wave. Existing evidence suggests that the optimal ciliated swimmers exploit the asymmetry on the collective level more than that on the individual level (Michelin & Lauga Reference Michelin and Lauga2010; Guo et al. Reference Guo, Nawroth, Ding and Kanso2014).
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 ‘useful power’ against ‘total power’. Useful power can be computed as the power needed to drag a rigid body of the same shape as the swimmer at the swim speed, while the total power is the rate of energy dissipation through viscous stresses in the flow to produce this motion (Lighthill Reference Lighthill1952). The goal of this paper is to find the optimal ciliary motion that maximises the swimming efficiency for an arbitrary axisymmetric microswimmer.
Studies of ciliated microswimmers can be loosely classified into two types of models. One type is known as the sublayer model in which the dynamics of each cilium is explicitly modelled, either theoretically (Blake & Sleigh Reference Blake and Sleigh1974; Brennen & Winet Reference Brennen and Winet1977) or numerically (Gueron & Liron Reference Gueron and Liron1992, Reference Gueron and Liron1993; Guirao & Joanny Reference Guirao and Joanny2007; Osterman & Vilfan Reference Osterman and Vilfan2011; Eloy & Lauga Reference Eloy and Lauga2012; Elgeti & Gompper Reference Elgeti and Gompper2013; Guo et al. Reference Guo, Nawroth, Ding and Kanso2014; Ito, Omori & Ishikawa Reference Ito, Omori and Ishikawa2019; Omori, Ito & Ishikawa Reference Omori, Ito and Ishikawa2020). The other type is known as the envelope model (commonly known as the squirmer model if the slip profile is time independent), which takes advantage of the densely packed nature of cilia, and traces the continuous envelope formed by the cilia tips. The envelope model has been extensively applied to study the locomotion of both single and multiple swimmers (e.g. see Lighthill Reference Lighthill1952; Blake Reference Blake1971; Ishikawa, Simmonds & Pedley Reference Ishikawa, Simmonds and Pedley2006; Ishikawa & Pedley Reference Ishikawa and Pedley2008; Michelin & Lauga Reference Michelin and Lauga2010; Vilfan Reference Vilfan2012; Brumley et al. Reference Brumley, Polin, Pedley and Goldstein2015; Elgeti, Winkler & Gompper Reference Elgeti, Winkler and Gompper2015; Guo et al. Reference Guo, Zhu, Liu, Bonnet and Veerapaneni2021; Nasouri, Vilfan & Golestanian Reference Nasouri, Vilfan and Golestanian2021), as well as the nutrient uptake of microswimmers (e.g. Magar, Goto & Pedley Reference Magar, Goto and Pedley2003; Magar & Pedley Reference Magar and Pedley2005; Michelin & Lauga Reference Michelin and Lauga2011, Reference Michelin and Lauga2013). While originally developed for spherical swimmers, the envelope model has been generalised to include spheroidal swimmers (e.g. Ishimoto & Gaffney Reference Ishimoto and Gaffney2013; Theers et al. Reference Theers, Westphal, Gompper and Winkler2016).
In particular, in a seminal work, Michelin & Lauga (Reference Michelin and Lauga2010) studied the optimal beating stroke of 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 optimisation showed that the envelope surface deforms in a wave-like 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. However, such an extension is hard to carry out using semi-analytical methods. Therefore, in this paper, we develop a computational framework for optimising the ciliary motion of a microswimmer with an arbitrary axisymmetric shape. We employ the envelope model, wherein the envelope is restricted to move tangentially to the surface, such that the shape of the microswimmer is unchanged during the beating period. We use a boundary-integral method (BIM) to solve the forward problem and derive an adjoint-based formulation for solving the optimisation problem.
The paper is organised as follows. In § 2, we introduce the optimisation problem, derive the sensitivity formulas and discuss our numerical solution procedure. In § 3, we present the optimal unconstrained and constrained solutions for microswimmers of various shape families. Finally, in § 4, we discuss our conclusions and future directions.
2. Problem formulation
2.1. Model
Consider an axisymmetric microswimmer whose boundary 
$\varGamma$ is obtained by rotating a generating curve 
$\gamma$ of length 
$\ell$ about the 
$\boldsymbol {e}_3$ axis, as shown in figure 1(a). We adopt the classic envelope model (Lighthill Reference Lighthill1952) and assume that the ciliary tips undergo time-periodic tangential movements along the generating curve. Let 
$s=\alpha (s_0,t)$ be the ciliary tip arc length coordinate on the generating curve 
$\gamma$ at time 
$t$ for a cilium rooted at 
$s_0$. The tangential slip velocity of this material point in its body-frame is thus
\begin{equation} u^{S}(s, t)= u^{S}(\alpha(s_0, t), t)= \partial_t \alpha(s_0, t). \end{equation}
Figure 1. (a) Schematic of the microswimmer geometry. The shape is assumed to be axisymmetric, obtained by rotating the generating curve 
$\gamma$ about the 
$\boldsymbol {e}_3$ axis. The tip of the cilium rooted at 
$s_0$ at time 
$t$ is given by 
$s=\alpha (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 arc length (black curve) is 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.
In addition to the time-periodic condition, the ciliary motion 
$\alpha$ needs to satisfy two more conditions to avoid singularity (Michelin & Lauga Reference Michelin and Lauga2010). First, the slip velocities should vanish at the poles
\begin{equation} \alpha(0,t) = 0 \quad\text{and} \quad \alpha(\ell,t) = \ell, \quad \forall\, t\in\mathbb{R}^+, \end{equation}
and second, 
$\alpha$ should be a monotonic function, that is,
\begin{equation} \partial_{s_0} \alpha(s_0,t) > 0, \quad \forall\, (s_0,t)\in [0,\ell]\times\mathbb{R}^+. \end{equation}
The last condition ensures the slip velocity is unique at any arc length 
$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:
\begin{equation} {-}\mu\nabla^2\boldsymbol{u} + \boldsymbol{\nabla} p = \boldsymbol{0},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{equation}
where 
$\mu$ is the fluid viscosity, and 
$p$ and 
$\boldsymbol {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
\begin{equation} \lim_{\boldsymbol{x}\rightarrow\infty}\boldsymbol{u}(\boldsymbol{x},t) = \boldsymbol{0}. \end{equation}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
\begin{equation} \int_\varGamma \boldsymbol{f}(\boldsymbol{x},t)\boldsymbol{\cdot}\boldsymbol{e}_3\mathrm{d}\varGamma=2{\rm \pi}\int_\gamma {f}_3(\boldsymbol{x},t)\, x_1\,\mathrm{d}s= 0, \end{equation}
where 
$x_1$ is the 
$\boldsymbol {e}_1$ component of 
$\boldsymbol {x}$, 
$\boldsymbol {f}$ is the active force density the swimmer applied to the fluid (negative to fluid traction) and 
$f_3$ is its 
$\boldsymbol {e}_3$ component.
Given any ciliary motion 
$\alpha (s_0,t)$ that satisfies (2.2) and (2.3), there is a unique tangential slip velocity 
${u}^{{S}}(s,t)$ defined by (2.1). Such a slip velocity propels the microswimmer at a translational velocity 
$U(t)$ in the 
$\boldsymbol {e}_3$ direction, determined by (2.6). Its angular velocity as well as the translational velocities in the 
$\boldsymbol {e}_1$ and 
$\boldsymbol {e}_2$ directions are zero by symmetry. Consequently, the boundary condition on 
$\gamma$ is given by
\begin{equation} \boldsymbol{u}(\boldsymbol{x}(s),t) = {u}^{{S}}(s,t)\boldsymbol\tau(s)+{U}(t)\boldsymbol{e}_3, \end{equation}
where 
$\boldsymbol \tau$ is the unit tangent vector on 
$\gamma$. Thereby, the instantaneous power loss 
$P(t)$ can be written as
\begin{align} P(t) &= \int_\varGamma \boldsymbol{f}(\boldsymbol{x},t) \boldsymbol{\cdot} \boldsymbol{u}(\boldsymbol{x}, t) \, \mathrm{d}\varGamma \nonumber\\ &= 2{\rm \pi}\left[ \int_\gamma \boldsymbol{f}(s,t) \boldsymbol{\cdot} \boldsymbol{\tau}(s) u^{S}(s,t)\, x_1 \, \mathrm{d}s + U(t)\int_\gamma \boldsymbol{f}(s,t) \boldsymbol{\cdot} \boldsymbol{e}_3\, x_1 \, \mathrm{d}s \right]. \end{align}The second term on the right-hand side is zero provided that the no-net-force condition (2.6) is satisfied.
Following Lighthill (Reference Lighthill1952), we quantify the performance of the microswimmer by its swimming efficiency 
$\epsilon$, defined as
\begin{equation} \epsilon = \frac{C_D \langle{U}\rangle^2}{\langle{P}\rangle}, \end{equation}
where 
${P}=P(t)$ and 
${U}=U(t)$ are the instantaneous power loss and swim speed, 
$\langle \cdot \rangle$ 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 the 
$\boldsymbol {e}_3$ direction. The coefficient 
$C_D$ depends on the given shape 
$\gamma$ only; for example, 
$C_D = 6{\rm \pi} \mu a$ in the case of a spherical microswimmer with radius 
$a$.
In our simulations, we normalise the radius of the microswimmer to unity, and the period of the ciliary motion to 
$2{\rm \pi}$. It is worth noting that the swimming efficiency (2.9) is size and period independent, thanks to its dimensionless nature. The Reynolds number of a ciliated microswimmer of radius 
$100\ \mathrm {\mu } \text {m}$ and frequency 30 Hz submerged in water can be estimated as 
${Re} \sim 10^{-4}$, confirming the applicability of the Stokes equations.
2.2. Numerical algorithm for solving the forward problem
Before stating the optimisation problem, we summarise our numerical solution procedure for the governing equations (2.4)–(2.7). By the quasi-static nature of the Stokes equation (2.4), the flow field 
$\boldsymbol {u}(\boldsymbol {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). We use a BIM at every time step. A similar BIM implementation was detailed in our recent work Guo et al. (Reference Guo, Zhu, Liu, Bonnet and Veerapaneni2021) in which we studied the optimisation of time-independent slip profiles. The main procedures are summarised below, with key equations highlighted in boxes.
We use the single-layer potential ansatz, which expresses the velocity as a convolution of an unknown density function 
$\boldsymbol {\mu }$ with Green's function for the Stokes equations:
\begin{equation} \boldsymbol{u}(\boldsymbol{x}) = \frac{1}{8{\rm \pi}}\int_\varGamma \left(\frac{1}{|\boldsymbol{r}|} \boldsymbol{I} + \frac{\boldsymbol{r} \otimes \boldsymbol{r} }{|\boldsymbol{r}|^3}\right) \, \boldsymbol{\mu}(\boldsymbol{y}) \, \mathrm{d}\varGamma(\boldsymbol{y}), \quad\text{where}\ \boldsymbol{r} = \boldsymbol{x} - \boldsymbol{y}.\end{equation}
The force density can then be evaluated as a convolution of 
$\boldsymbol {\mu }$ with the (negative of) traction kernel:
\begin{equation} \boldsymbol{f}(\boldsymbol{x}) =\frac{1}{2}\boldsymbol{\mu}\left(\boldsymbol{x}\right) + \frac{3}{4{\rm \pi}}\int_{\varGamma} \left(\frac{\boldsymbol{r} \otimes \boldsymbol{r} }{|\boldsymbol{r}|^5}\right) (\boldsymbol{r}\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x}))\boldsymbol{\mu}\left(\boldsymbol{y}\right)\mathrm{d}\varGamma\left(\boldsymbol{y}\right). \end{equation}
We convert these weakly singular boundary integrals into convolutions on the generating curve 
$\gamma$ by performing an analytic integration in the orthoradial direction, and then apply a high-order quadrature rule designed to handle the log singularity of the resulting kernels (Veerapaneni et al. Reference Veerapaneni, Gueyffier, Biros and Zorin2009). The Stokes flow problem defined at any time 
$t$ by (2.4)–(2.7) is then recast as the BIM system for the unknowns 
$\boldsymbol {\mu }$ and 
$U(t)$ obtained by substituting (2.10) into (2.7) and (2.11) into (2.6). The numerical solution method consists in discretising 
$\gamma$ 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=10 N_p$ quadrature nodes, while the values of 
$\boldsymbol {\mu }$ are sought at the same quadrature nodes. The resulting BIM system is
\begin{equation} \begin{bmatrix}
\mathcal{S} & -\mathcal{B} \\ \mathcal{C} & 0
\end{bmatrix} \begin{bmatrix} \boldsymbol{\mu}\\
{U(t)} \end{bmatrix} = \begin{bmatrix}
{\boldsymbol{u}^{{S}}}\\ {0} \end{bmatrix},
\end{equation}
where the vectors 
$\boldsymbol {\mu }=\boldsymbol {\mu }(s_q,t)$ and 
$\boldsymbol {u}^{{S}}=\boldsymbol {u}^{{S}}(s_q,t)$ are the unknown density and the given slip velocity at all quadrature nodes 
$s_q$, 
$\mathcal {S}$ is the axisymmetric single-layer potential operator (which is fixed for a given shape 
$\gamma$), 
$\mathcal {B}$ is the column vector reproducing 
$\boldsymbol {e}_3$ at each quadrature node and 
$\mathcal {C}$ is the row vector, such that 
$\mathcal {C}[\boldsymbol {\mu }] = \int _\varGamma \boldsymbol {f}(\boldsymbol {x})\boldsymbol {\cdot }\boldsymbol {e}_3 \mathrm {d}\varGamma$ is the total traction force in the 
$\boldsymbol {e}_3$ direction.
The algorithm to obtain the slip velocity at the quadrature nodes at a given time 
$\boldsymbol {u}^{S}(s_q,t)$ is summarised in figure 1(b). Specifically, we start by computing the corresponding ciliary tip position 
$s=\alpha (s_q,t)$ and the slip velocity 
$u^{S}(s,t)$ from (2.1). These tip positions 
$s$ can be highly non-uniform, depending on the form of 
$\alpha$, which can 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 
$\boldsymbol {u}^{S}(s,t)$ directly on the right-hand side of (2.12), which will bypass the interpolation steps mentioned above. However, it requires reassembly of the matrix 
$\mathcal {S}$ at every time step, significantly increasing the computational cost.
2.3. Optimisation problem
The goal of this work is to find the optimal ciliary motion for a given arbitrary axisymmetric shape; that is, the ciliary motion 
$\alpha ^{\star }(s_0,t)$ that maximises the swimming efficiency 
$\epsilon$:
\begin{equation}
\alpha^{{\star}}=\mathop{\textrm{arg max}}\limits_{\alpha\in\mathcal{A}}\epsilon(\alpha),
\end{equation}
where 
$\mathcal {A}$ is the space of all possible time-periodic ciliary motion satisfying (2.2) and (2.3). It is, however, not easy to define and manipulate finite-dimensional parametrisations of 
$\alpha$ that remain in that space. To circumvent this difficulty, we follow the ideas in Michelin & Lauga (Reference Michelin and Lauga2010) and represent 
$\alpha$ in terms of a time-periodic function 
$\psi (x,t)$, such that
\begin{equation} \alpha(s_0,\psi) = \frac{\ell\int_0^{s_0}{[\psi(x,t)]^2\,\mathrm{d}\kern0.04em x}}{\int_0^{\ell}{[\psi(x,t)]^2\,\mathrm{d}\kern0.04em x}}, \end{equation}
where 
$\ell$ is the total length of the generating curve 
$\gamma$. Note that 
$\alpha$ is also (implicitly) a function of time 
$t$, through 
$\psi = \psi (x,t)$. It is easy to verify that 
$\alpha$ given by (2.14) satisfies the boundary conditions (2.2) and the monotonicity requirement (2.3) for any choice of 
$\psi$. Conversely, for any 
$\alpha$ satisfying (2.2) and (2.3), there is at least one 
$\psi$ that provides 
$\alpha$. As a result, the optimisation problem is recast as finding
\begin{equation} \psi^{{\star}} =
\mathop{\textrm{arg max}}\limits_{\psi} \epsilon(\psi),
\end{equation}
where 
$\psi (\cdot ,t)$ is only required to be square-integrable over 
$[0,\ell ]$ for any 
$t$.
We use a quasi-Newton Broyden–Fletcher–Goldfarb–Shanno (BFGS) method (Nocedal & Wright Reference Nocedal and Wright2006) to optimise the ciliary motion via 
$\psi$, which requires repeated evaluations of efficiency sensitivities with respect to perturbations of 
$\psi$. 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 over the traditional finite-difference method when finding the sensitivities, because 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 
$\alpha$. The derivations are detailed below.
2.4. Sensitivity analysis
We start by finding the sensitivities in terms of the slip profile 
$u^{S}$. The sensitivities in terms of the auxiliary unknown 
$\psi$ will be found subsequently by a change of variable. As the concept of an 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 (2.4)–(2.7) in mixed weak form (see e.g. Brezzi & Fortin Reference Brezzi and Fortin1991, chapter 6). That is, to find 
$(\boldsymbol {u}, p, \boldsymbol {f}, U) \in \boldsymbol {\mathcal {V}}\times \mathcal {P}\times \boldsymbol {\mathcal {F}}\times \mathbb {R},$ such that
$$\begin{gather} a(\boldsymbol{u}, \boldsymbol{v}) - b(\boldsymbol{v},p) - b(\boldsymbol{u},q) - \langle\, \boldsymbol{f}, \boldsymbol{v}\rangle_{\varGamma} = 0 \quad \forall\, (\boldsymbol{v},q) \in \boldsymbol{\mathcal{V}}\times\mathcal{P}, \end{gather}$$
$$\begin{gather}\langle \boldsymbol{g}, \boldsymbol{e}_3\rangle_\varGamma U + \langle \boldsymbol{g}, u^{S} \boldsymbol{\tau} \rangle_\varGamma - \langle \boldsymbol{g}, \boldsymbol{u}\rangle_\varGamma = 0 \quad \forall \boldsymbol{g} \in \boldsymbol{\mathcal{F}}, \end{gather}$$
$$\begin{gather}\langle\, \boldsymbol{f}, \boldsymbol{e}_3\rangle_\varGamma = 0, \end{gather}$$
where the bilinear forms 
$a$ and 
$b$ are defined by
\begin{equation} a(\boldsymbol{u}, \boldsymbol{v}) := \int_\varOmega 2\mu \boldsymbol{D}[\boldsymbol{u}]:\boldsymbol{D}[\boldsymbol{v}]\, \mathrm{d} V, \quad b(\boldsymbol{v},q) := \int_\varOmega q\, \text{div}\,\boldsymbol{v}\, \mathrm{d}V, \end{equation}
and 
$\boldsymbol {D}[\boldsymbol {u}] := (\boldsymbol {\nabla } \boldsymbol {u}+\boldsymbol {\nabla }^{\textrm {T}}\boldsymbol {u})/2$ is the strain-rate tensor; we use 
$\langle \cdot , \cdot \rangle _\varGamma$ as shorthand for the inner product on 
$\varGamma$. For example, 
$\langle \boldsymbol {f}, \boldsymbol {v}\rangle _\varGamma = \int _\varGamma \boldsymbol {f} \cdot \boldsymbol {v}\, \mathrm {d}\varGamma$. Similarly, with a slight abuse of notation, the power loss functional can be written as 
$P(u^{S}) := \langle \boldsymbol {f}, u^{S}\boldsymbol {\tau } + U\boldsymbol {e}_3 \rangle _\varGamma$, where 
$U := U(u^{S})$ is the swim speed functional.
The Dirichlet boundary condition (2.7) is (weakly) enforced explicitly through (2.16b), rather than being embedded in the velocity solution space 
$\boldsymbol {\mathcal {V}}$, as this will facilitate the derivation of slip-derivative identities; this is in fact our motivation for using the mixed weak form (2.16). Condition (2.16c) is the no-net-force condition (2.6).
First-order sensitivities of functionals at 
$u^{S}$ are defined as directional derivatives, by considering perturbations of 
$u^{S}$ of the form
\begin{equation} u^{S}_{\eta} = u^{S} + \eta \nu \end{equation}
for some 
$\nu$ in the slip velocity space and 
$\eta \in \mathbb {R}$. Then, the directional (or Gâteaux) derivative of a functional 
$J(u^{S})$ in the direction 
$\nu$, denoted by 
$J'(u^{S};\nu )$, is defined as
\begin{equation} J'(u^{S};\nu) = \lim_{\eta\to0} \frac{1}{\eta}\left( J[u_{\eta}^{S}]-J[u^{S}] \right). \end{equation}
For the power-loss functional, we obtain (because the derivative of 
${u}^{S}$ in the above sense is 
$\nu$)
\begin{equation} P'(u^{S};\nu) = \langle\, \boldsymbol{f}',{u}^{S}\boldsymbol{\tau} + U\boldsymbol{e}_3 \rangle_{\varGamma} + \langle\, \boldsymbol{f},\nu\boldsymbol{\tau} \rangle_{\varGamma} + \langle\, \boldsymbol{f},\boldsymbol{e}_3 \rangle_{\varGamma}U', \end{equation}
where 
$\boldsymbol {f}'$ and 
$U'$ are the derivatives of the active force 
$\boldsymbol {f}$ and swim speed 
$U$ solving problem (2.16), considered as functionals on the slip velocity 
$u^{S}$:
\begin{equation} \boldsymbol{f}' = \lim_{\eta\to0} \frac{1}{\eta}\left( \boldsymbol{f}[u_{\eta}^{S}]-\boldsymbol{f}[u^{S}] \right), \quad U' = \lim_{\eta\to0} \frac{1}{\eta}\left( U[u_{\eta}^{S}]-U[u^{S}] \right). \end{equation} Differentiating the weak formulation (2.16) of the forward problem with respect to 
$u^{S}$ leads to the weak formulation of the governing problem for the derivatives 
$(\boldsymbol {u}', \boldsymbol {f}', p', U')$ of the solution 
$(\boldsymbol {u}, \boldsymbol {f}, p, U)$:
$$\begin{gather} a(\boldsymbol{u}',\boldsymbol{v}) - b(\boldsymbol{u}',q) - b(\boldsymbol{v},p') - \langle\, \boldsymbol{f}',\boldsymbol{v} \rangle_{\varGamma} =0 \quad\, \forall(\boldsymbol{v},q)\in\boldsymbol{\mathcal{V}}\times\mathcal{P}, \end{gather}$$
$$\begin{gather}\langle \nu\boldsymbol{\tau},\boldsymbol{g}\rangle_{\varGamma} + U'\langle \boldsymbol{e}_3,\boldsymbol{g}\rangle_{\varGamma} - \langle \boldsymbol{u}',\boldsymbol{g} \rangle_{\varGamma} = 0 \quad \forall\, \boldsymbol{g}\in\boldsymbol{\mathcal{F}}, \end{gather}$$
$$\begin{gather}\langle\, \boldsymbol{f}',\boldsymbol{e}_3 \rangle_{\varGamma} = 0 . \end{gather}$$
Here we have assumed without loss of generality that the test functions in (2.16) verify 
$\boldsymbol {v}' = \boldsymbol {0}$, 
$\boldsymbol {g}' = \boldsymbol {0}$, and 
$q' = 0$, which is made possible by the absence of boundary constraints in 
$\boldsymbol {\mathcal {V}}$.
At first glance, evaluating 
$P'(u^{S};\nu )$ in a given perturbation 
$\nu$ appears to rely on solving the derivative problem (2.22). However, a more effective approach allows us to bypass the actual evaluation of 
${\boldsymbol {f}}'$. Let the adjoint problem be defined by
$$\begin{gather} a(\hat{\boldsymbol{u}},\boldsymbol{v}) - b(\hat{\boldsymbol{u}},q) - b(\boldsymbol{v},\hat{p}) - \langle\hat{\boldsymbol{f}},\boldsymbol{v}\rangle_{\varGamma} = 0 \quad \forall(\boldsymbol{v},q)\in\boldsymbol{\mathcal{V}}\times\mathcal{P}, \end{gather}$$
$$\begin{gather}\langle \boldsymbol{e}_3,\boldsymbol{g}\rangle_{\varGamma} - \langle \hat{\boldsymbol{u}},\boldsymbol{g}\rangle_{\varGamma} = 0 \quad \forall\boldsymbol{g}\in\boldsymbol{\mathcal{F}}, \end{gather}$$
i.e. 
$(\hat {\boldsymbol {u}},\hat {p})$ are the flow variables induced by prescribing a unit velocity 
$\boldsymbol {e}_3$ on 
$\varGamma$. For later convenience, we let 
$F_0$ denote the (non-zero) net force exerted on 
$\varGamma$ by the adjoint flow:
\begin{equation} F_0 := \langle\, \hat{\boldsymbol{f}},\boldsymbol{e}_3 \rangle_{\varGamma}. \end{equation}
Problem (2.23) in strong form is defined by (2.4)–(2.7) with 
$U=1,\,u^{S}=0$. In fact, 
$F_0$ takes the same value as the drag coefficient 
$C_D$ in (2.9).
Then, combining the derivative problem (2.22) with the forward problem (2.16) or the adjoint problem (2.23) with appropriate choices of test functions allows us to derive expressions of 
$P'(u^{S};\nu )$ and 
$U'(u^{S};\nu )$ which do not involve the forward solution derivatives.
Specifically, we set the test functions to 
$(\boldsymbol {v},q,\boldsymbol {g})=(\boldsymbol {u}',p',\boldsymbol {f}')$ in ((2.16a) and (2.16b)) of the forward problem and 
$(\boldsymbol {v},q,\boldsymbol {g})=({\boldsymbol {u}},{p},{\boldsymbol {f}})$ in ((2.22a) and (2.22b)) of the derivative problem. Then, the combination (2.22a) + (2.22b) 
$-$ (2.16a) 
$-$ (2.16b) is evaluated, to obtain
\begin{equation} \langle\, \boldsymbol{f}',{u}^{S}\boldsymbol{\tau} + U\boldsymbol{e}_3 \rangle_{\varGamma} = \langle\, \boldsymbol{f},\nu\boldsymbol{\tau} \rangle_{\varGamma} + \langle\, \boldsymbol{f},\boldsymbol{e}_3 \rangle_{\varGamma}U'. \end{equation}Substituting (2.25) into (2.20), and recalling the no-net-force condition (2.6), we have
\begin{equation} \boxed{ P'(u^{S};\nu) = 2\langle\, \boldsymbol{f},\nu\boldsymbol{\tau} \rangle_{\varGamma} ={4{\rm \pi}}\int_\gamma (\boldsymbol{f} \boldsymbol{\cdot} \boldsymbol{\tau}) \, \nu x_1 \, \mathrm{d}s.} \end{equation} Likewise, setting the test functions to 
$(\boldsymbol {v},q,\boldsymbol {g})=(\boldsymbol {u}',p',\boldsymbol {f}')$ in the adjoint problem (2.23) and 
$(\boldsymbol {v},q,\boldsymbol {g})=(\hat {\boldsymbol {u}},\hat {p},\hat {\boldsymbol {f}})$ in ((2.22a) and (2.22b)) of the derivative problem (2.22), then evaluating (2.22a) + (2.22b) 
$-$ (2.23a) 
$-$ (2.23b), yields
\begin{equation} 0 = \langle\, \hat{\boldsymbol{f}},\nu\boldsymbol{\tau} \rangle_{\varGamma} + \langle\, \hat{\boldsymbol{f}},U'\boldsymbol{e}_3 \rangle_{\varGamma} - \langle\, \boldsymbol{f}',\boldsymbol{e}_3 \rangle_{\varGamma} = \langle\, \hat{\boldsymbol{f}},\nu\boldsymbol{\tau} \rangle_{\varGamma} + F_0U'. \end{equation}
Note that 
$\langle \boldsymbol {f}',\boldsymbol {e}_3 \rangle _{\varGamma } = 0$ according to (2.22c). Rearranging terms in (2.27), we have
\begin{equation} \boxed{ U'(u^{S};\nu) ={-}\frac{1}{F_0}\langle\, \hat{\boldsymbol{f}},\nu\boldsymbol{\tau} \rangle_{\varGamma} ={-} \frac{2{\rm \pi}}{F_0} \int_\gamma (\hat{\boldsymbol{f}}\boldsymbol{\cdot} \boldsymbol{\tau}) \, \nu x_1\, \mathrm{d}s.} \end{equation} The sensitivity formulas (2.26) and (2.28) are not practically applicable in this form to the current optimisation problem, because the constraints (2.2) and (2.3) are not easy to enforce on parametrisations of the unknown slip profiles 
$u^{S}$. For this reason, we rewrite the quantities of interest as functionals of 
$\psi$, and the connection between 
$\psi$ and 
$\alpha$ is given by (2.14). Specifically, the slip profile is
\begin{equation} u^{S}(s,t) = \partial_{t}\alpha(s_0,\psi) = \partial_\psi\alpha(s_0,\psi;\dot{\psi}) = \partial_\psi\alpha\left( \beta(s,\psi),\psi;\dot{\psi} \right) = v^{S}(s,\psi), \end{equation}
where 
$\dot {\psi } := \partial _t \psi$ and 
$\beta (s,\psi )$ is the inverse function of 
$\alpha$, i.e. 
$s_0 = \beta (s,\psi )$. The average power-loss and swim-speed functionals are written as
\begin{equation} \langle\mathbb{P}\rangle(\psi) := \langle P \rangle(u^{S}), \quad \langle \mathbb{U}\rangle(\psi):= \langle U\rangle(u^{S}) \quad \text{with}\ u^{S}(s,t) = v^{S}(s,\psi). \end{equation} On applying the change of variables 
$s = \alpha (s_0, \psi )$ in the integrals (2.26) and (2.28) and averaged over one period, we obtain
$$\begin{gather} \boxed{\langle\mathbb{P}\rangle'(\psi;\hat{\psi}) = 2 \int_0^{2{\rm \pi}} \int_\gamma \boldsymbol{f}(\alpha)\boldsymbol{\cdot}\boldsymbol{\tau}(\alpha)\,x_1(\alpha)\, v^{S}{}'(s,\psi;\hat{\psi}) \, \partial_s\alpha \, \mathrm{d}s_0 \,\mathrm{d}t,} \end{gather}$$
$$\begin{gather}\boxed{ \langle \mathbb{U}\rangle'(\psi;\hat{\psi}) ={-}\frac{1}{F_0} \int_0^{2{\rm \pi}} \int_\gamma \hat{\boldsymbol{f}}(\alpha)\boldsymbol{\cdot}\boldsymbol{\tau}(\alpha)\,x_1(\alpha)\,v^{S}{}'(s,\psi;\hat{\psi}) \, \partial_s\alpha\, \mathrm{d}s_0 \, \mathrm{d}t,} \end{gather}$$
where 
$v^{S}{}'(s, \psi ; \hat {\psi })$ is the directional derivative of 
$u^{S}$ with respect to 
$\psi$ and in the direction 
$\hat {\psi }$. Specifically, we can show that
$$\begin{gather} v^{S}{}'(s,\psi;\hat{\psi}) \, \partial_s\alpha(s_0,\psi) \mathrm{d}s_0 = \left\{ \partial_s\alpha(s_0,\psi)\,\left[ \partial^2_\psi \alpha \left( s_0,\psi;\hat{\psi},\dot{\psi}\right) + \partial_\psi \alpha \left( s_0,\psi;\dot{\hat{\psi}} \right)\,\right] \right. \nonumber\\ \left. -\partial_{\psi s}\alpha\left( s_0,\psi;\dot{\psi} \right)\; \partial_{\psi}\alpha\left( s_0,\psi;\hat{\psi} \right) \right\} \,\mathrm{d}s_0. \end{gather}$$
The derivation and the explicit expression of each term in (2.33) are given in Appendix A. Finally, the efficiency sensitivity in terms of 
$\psi$ readily follows by the quotient rule
\begin{equation} \epsilon'(\psi;\hat\psi) = C_D\frac{2\langle\mathbb{U}\rangle \langle\mathbb{U}\rangle'\langle\mathbb{P}\rangle - \langle\mathbb{U}\rangle^2\langle\mathbb{P}\rangle'}{\langle\mathbb{P}\rangle^2}. \end{equation}2.5. Constraints on surface displacement
The unconstrained optimisation problem (2.15) introduced previously 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 & Lauga (Reference Michelin and Lauga2010), we replace the initial unconstrained optimisation problem (2.15) with the penalised optimisation problem
\begin{equation} \psi^{{\star}} = \textrm{arg max}_{\psi} E(\psi), \quad E(\psi) = \epsilon(\psi) - C(\psi), \end{equation}
where the (non-negative) penalty term 
$C(\psi )$, defined as
\begin{equation} C(\psi) = \int_0^{\ell} H({A}(\psi) - c) \,\mathrm{d}s_0, \end{equation}
serves to incorporate the kinematic constraint 
$A(\psi )\leq c$ in the optimisation problem. The functional 
${A}(\psi )$ in (2.36) is a measure of the amplitude of the displacement of individual material points for the stroke (through 
$\alpha$) and 
$c$ is a threshold parameter to bound 
$A(\psi )$ (a smaller 
$c$ corresponding to a stricter constraint). We use 
$H$ as a smooth non-negative penalty function defined by
\begin{equation} H(u) = \varLambda_1\left[1+\tanh \left(\varLambda_2{u}\right)\right]u^2, \end{equation}
which for a large enough 
$\varLambda _2$ approximates 
$u\mapsto 2\varLambda _1 u^2 Y(u)$ (
$Y$ being the Heaviside unit step function). The multiplicative parameter 
$\varLambda _1$ then serves to tune the severity of the penalty incurred by violations of the constraint 
$A(\psi )\leq c$. We use 
$\varLambda _1 = 10^4$ and 
$\varLambda _2 = 10^{4}$ in our numerical simulations unless otherwise mentioned. The optimisation results are not sensitive to the choice of 
$\varLambda _1$ and 
$\varLambda _2$. A small caveat of the penalty function (2.37) is that it has a (small) bump at 
$\varLambda _2 u \approx -1.109$. This bump can occasionally trap the optimisations into local extrema that have significantly lower efficiencies, depending on the initial guesses. Perturbing 
$\varLambda _2$ for such cases helps to alleviate the problem.
The most relevant physical definition of 
${A}$ would be the actual displacement amplitude of an individual point, i.e. 
$\Delta s= [\alpha _{max}(s_0) - \alpha _{min}(s_0)]/2$. However, the strong nonlinearity of this measure is not appropriate for the computation of the gradient. Following Michelin & Lauga (Reference Michelin and Lauga2010), we measure the displacement by its variance in time:
\begin{equation} {A}(\psi) = \langle (\alpha(s_0,\psi) - \langle\alpha\rangle(s_0))^2\rangle. \end{equation}
The maximum displacement 
$\Delta s_{max} = \max _{s_0}(\Delta s)$ will be found post-optimisation for the optimal ciliary motion 
$\alpha ^{\star }$ to better illustrate our results in § 3.
Like the initial problem (2.15), the penalised problem (2.35a,b) is solvable using unconstrained optimisation methods, and we again adopt a quasi-Newton BFGS algorithm to optimise the ciliary motion. Applying the chain rule to the penalty functional 
$C(\psi )$, we obtain the derivative of the penalty term in the direction of 
$\hat \psi$ as
\begin{equation} C'(\psi; \hat\psi) = \int_0^\ell H'({A}(\psi) - c) {A}'(\psi; \hat\psi) \,\mathrm{d}s_0. \end{equation}
The derivative of the penalised objective functional 
$E(\psi )$ is therefore
\begin{equation} E'(\psi;\hat{\psi}) = \epsilon'(\psi;\hat{\psi}) - C'(\psi;\hat{\psi}), \end{equation}
where 
$\epsilon '$ and 
$C'$ are given by (2.34) and (2.39), respectively.
3. Results
3.1. Parametrisation
We parametrise 
$\psi (s_0,t)$, such that
\begin{equation} \psi(s_0,t) = \sum_{k=1}^m \xi_k(t) B_k(s_0), \end{equation}
where 
$B_k$ are the fifth-order B-spline basis functions and their coordinates 
$\xi _k(t)$ are expanded as trigonometric polynomials 
$\xi _k(t) = {a_{0k}}/{2} + \sum _{j=1}^n [a_{jk} \cos jt + b_{jk} \sin jt]$ to ensure time periodicity. Taken together, we have
\begin{equation} \psi(s_0,t) = \sum_{k=1}^m \left[\frac{a_{0k}}{2}+ \sum_{j=1}^n (a_{jk} \cos jt + b_{jk} \sin jt)\right] B_k(s_0) \end{equation}
so that the finite-dimensional optimisation 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. A brief discussion about the initial guesses can be found in Appendix B. To evaluate the gradient of 
$E(\psi )$ with respect to the design parameters 
$a_{0k}$, 
$a_{jk}$ and 
$b_{jk}$, we use (2.40) with 
$\hat {\psi }$ taken as the basis functions of the adopted parametrisation (3.2), i.e. 
$\hat {\psi }(s_0,t)=B_k(s_0)/2$, 
$\hat {\psi }(s_0,t)=B_k(s_0)\cos jt$ and 
$\hat {\psi }(s_0,t)=B_k(s_0)\sin jt$, respectively. In terms of parametrisation, local minima are multiple in the parameter space, because multiplying optimal parameters by a constant factor yields the same optimum for 
$\alpha$.
3.2. Spheroidal swimmers
By way of validation, we start with optimising the ciliary motion of a spherical microswimmer. The efficiency 
$\epsilon$ as a function of iteration number for the unconstrained optimisation (2.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 arc length coordinate of a cilium tip over one period. We observe, similar to the results of Michelin & Lauga (Reference Michelin and Lauga2010), 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 Reference Knight-Jones1954; Blake Reference Blake1972). We note that this optimal ciliary motion produces an efficiency higher than the 
$23\,\%$ efficiency obtained numerically by Michelin & Lauga (Reference Michelin and Lauga2010, figure 11). This is owing to a larger maximum displacement 
$\Delta s_{max} \approx 0.45\ell$ in our optimisations (translated to a maximum angle of 
$81^{\circ }$ vs 
$53^{\circ }$). Our optimisation result aligns well with their results using the analytical ansatz (Michelin & Lauga Reference Michelin and Lauga2010, figure 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\,\%$ at the limit of the maximum displacement approaches 
$90^{\circ }$, 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).
Figure 2. Unconstrained optimisation history of a spherical swimmer and a prolate swimmer with a 
$2\,{:}\,1$ aspect ratio. The optimal spherical swimmer has an efficiency 
$\epsilon \approx 35\,\%$ and swim speed 
$\langle U\rangle \approx 1.2$. The optimal prolate swimmer has an efficiency 
$\epsilon \approx 69\,\%$ and swim speed 
$\langle U\rangle \approx 1.5$. The results for the spherical and prolate swimmers are depicted in blue and green colours, respectively. (a) The efficiency is shown as a function of iteration number. (b,c) The ciliary motions of the optimal swimmers. (d,e) The time-averaged slip velocities (at Eulerian points) are shown in solid curves. Dashed curves are the time-independent optimal slip velocities of the given shape scaled by the swim speed (Guo et al. Reference Guo, Zhu, Liu, Bonnet and Veerapaneni2021). Parameters used in the optimisation: 
$m = 25, n = 2$. The number of panels 
$N_p = 20$, number of sample points 
$N_s = 80$ and number of time steps per period 
$N_t = 50$. We use the same in the remainder of the paper unless otherwise stated. Note that the vertical axes of panels (b) and (c) are flipped so that the north pole (
$s=0$) appear on the top of the figure. The corresponding waveforms are known as antiplectic metachronal waves (tips are spread out during the effective stroke and close together during the recovery stroke). Videos of the optimal ciliary motions can be found in the online supplementary material (Movies 1 and 2) available at https://doi.org/10.1017/jfm.2021.744.
The optimal unconstrained prolate spheroidal microswimmer with a 
$2\,{:}\,1$ aspect ratio has an efficiency 
$\epsilon \approx 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. Reference Guo, Zhu, Liu, Bonnet and Veerapaneni2021). The optimal ciliary motion is very close to that of the spherical swimmer (figure 2b,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 check, swapping the ciliary motions obtained from optimising 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 (2.38). We control the maximum variance by tuning 
$c$ in (2.36), and the efficiencies are plotted against the maximum displacement 
$\Delta s_{max}$ scaled by the total arc length 
$\ell$ in figure 3. Three different random initial guesses are used for each 
$c$. The unconstrained optimisation 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 arc length of the generating curve (
$\Delta s_{max}\approx \ell /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 & Lauga (Reference Michelin and Lauga2010). 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 
$\Delta s_{max}/\ell$, and converges to the kinematically unconstrained optimal microswimmer as the maximum variance 
$c$ is increased.
Figure 3. Efficiency as a function of maximum displacement of ciliary tips. Blue and green symbols represent spherical and prolate spheroidal swimmers (
$2\,{:}\,1$ aspect ratio), respectively. Diamond symbols are the optimal unconstrained case. Open symbols are optimisation results of spherical swimmers taken from Michelin & Lauga (Reference Michelin and Lauga2010, figure 11).
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 (figure 4a), 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.
Figure 4. Ciliary motion (a) and mean slip velocity (b) for the optimal spherical swimmer with constraint (
$\Delta s_{max}/\ell \approx 5.0\,\%$). The efficiency is 
$\epsilon \approx 6.9\,\%$ and the swim speed is 
$\langle U\rangle \approx 0.091$. The swimmer forms multiple waves in the equatorial region, leading to a high slip velocity at 
$s\approx 0.5\ell$. 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).
3.3. 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 the north and south poles always maximises 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(\theta )\sin \theta , R(\theta )\cos \theta )$, where 
$R(\theta ) = (1+\delta \cos 2\theta )$ is a function that makes the radius non-constant, and 
$\theta \in [0,{\rm \pi} ]$ is the parametric coordinate. For 
$0<\delta <1$, the radius is the smallest at 
$\theta = {\rm \pi}/2$, corresponding to a ‘neck’ around the equator. In the limit 
$\delta =0$, the generating curve reduces to a semicircle and the swimmer reduces to a spherical swimmer.
The optimisation results are depicted in figure 5 for 
$0\le \delta \le 0.8$. Some corresponding shapes are shown as insets. The median efficiencies of 10 Monte Carlo simulations are plotted for each 
$\delta$ value, and compared against the time-independent efficiencies. For all three cases (constrained, unconstrained and time-independent), the efficiencies increase as 
$\delta$ increases from 
$0$ to 
$0.3$. This is because increasing 
$\delta$ in this regime makes the shape more elongated. Increasing 
$\delta$ 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\le \delta \le 0.6$.
Figure 5. Constrained optimisations could lead to more efficient ciliary motions for microswimmers with a thin ‘neck’ on average. (a) Efficiencies of the microswimmers with various neck widths. The median efficiencies of the time-dependent optimisations across 10 randomised initial conditions are shown for each shape in cross symbols ‘
$\times$’. Unconstrained and constrained optimisations (
$c=1$) are depicted in blue and green, respectively. Efficiencies of the microswimmers with time-independent slips are shown, using black circle symbols ‘
$\circ$’, as a reference. (b,c) Ciliary motions of microswimmers with 
$\delta =0.8$ from unconstrained and constrained optimisations from the same initial guess. The swimming efficiencies are 
$20\,\%$ and 
$29\,\%$, respectively. (d,e) Mean slip velocity corresponding to the ciliary motions in (b,c). 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 (Movies 4 and 5).
Interestingly, unconstrained optimisation 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 
$\delta \ge 0.7$. We note that the unconstrained optimisations are likely to be trapped in local optima where the ciliary motion forms a single wave (figure 5b), whereas the constrained optimisations are ‘forced’ to find the ciliary motion with multiple waves split at the equator (figure 5c), 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. Figure 5(d,e) shows 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 optimisation. In a way, the constraint in cilia length is helping the optimiser 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 (
$\delta = 0.8$) with various constraints in figure 6. In general, the highest efficiency from the Monte Carlo simulations increases with the constraint for 
$c\le 0.8$, similar to the case of spheroidal swimmers (figure 3). Further 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 \ge 0.8$, consistent with the observation from figure 5. It is worth noting that the constrained optimisation is more likely to get stuck in very low efficiencies (e.g. the lowest outlier for 
$c=0.8$), possibly owing to the secondary bump of the penalty function 
$C$ mentioned earlier.
Figure 6. Statistical results of thin neck microswimmer of 
$\delta = 0.8$ with various constraint 
$c$ for 10 Monte Carlo simulations. The unconstrained simulation is denoted by 
$c=\infty$. (a) Efficiencies grouped by the constraint 
$c$. For each box, the central mark indicates the median of the 10 random simulations, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The outliers are denoted by red 
$+$ symbols. (b) Efficiencies plotted against the maximum displacement 
$\Delta s_{max}/\ell$. The numerical parameter 
$\varLambda _2$ is set to be 
$10^4$ by default. Occasionally the optimisation might stop within merely a few iterations, making the ciliary motion stuck in a very inefficient local minimum. Setting 
$\varLambda _2$ to 
$10^3$ for these cases (most of the time) cures the problem.
All data points from the optimisation are plotted in figure 6(b) as function of the maximum displacement 
$\Delta s_{max}$. The efficiencies grow almost linearly until 
$\Delta s_{max} \approx 0.25 \ell$, as in the case of spheroidal swimmers, and decrease for larger 
$\Delta s_{max}$. This is further evidence that the optimal ciliary motion for this shape consists of two separate waves travelling within each hemisphere. We want to emphasise that unconstrained optimisation can still reach the optimal ciliary motion, as shown in the box of 
$c=\infty$. However, it is more likely to reach the suboptimal ciliary motion compared with the constrained cases.
4. Conclusions and discussions
In this paper, we extended the work of Michelin & Lauga (Reference Michelin and Lauga2010) and studied the optimal ciliary motion for a microswimmer with an arbitrary axisymmetric shape. In particular, the forward problem is solved using a BIM and the sensitivities are derived using an adjoint-based method. The auxiliary function 
$\psi$ is parametrised 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) travelling waves, reminiscent of the typical metachronal waves observed in ciliated microswimmers. Specifically, for the spherical swimmer with limited cilia length (figure 4a), the ratio between the metachronal wavelength close to the equator and the cilia length could be estimated as 
${\lambda _{MW}}/{\Delta s_{max}}\approx {0.2\ell }/{0.05\ell } = 4$. This ratio lies in the higher end of the data collected in Velho Rodrigues, Lisicki & Lauga (Reference Velho Rodrigues, Lisicki and Lauga2021, table 9) for biological ciliates, which reports a 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 neighbouring 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. Reference Guo, Zhu, Liu, Bonnet and Veerapaneni2021).
Most interestingly, we found that constraining the cilia length for some shapes may lead to a better efficiency on average, compared with the unconstrained optimisation. 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 optimisation. Although the concave shapes studied in § 3.3 are somewhat non-standard, they allow 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 modelling individual cilia motions) have reported swimming or transport efficiencies in the orders of 
$0.1\sim 1\,\%$ (see, e.g. Elgeti & Gompper Reference Elgeti and Gompper2013; Ito et al. Reference Ito, Omori and Ishikawa2019; Omori et al. Reference Omori, Ito and Ishikawa2020), 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 consider 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\sim 95\,\%$ of the total energy dissipation, owing to the large shear rate inside the layer (see e.g. Keller & Wu Reference Keller and Wu1977; Ito et al. Reference Ito, Omori and Ishikawa2019). We note that it is possible to incorporate energy dissipation inside the ciliary layer in the envelope model, as previously done in Vilfan (Reference Vilfan2012), albeit for a time-independent slip profile. Additionally, the difference could also be due to modelling assumptions on the cilia length and the number of cilia. In particular, the cilia length considered in sublayer models is usually below 
$1/10$ of the body length. Omori et al. (Reference Omori, Ito and Ishikawa2020) showed that swimming efficiency increases with the cilia length as fast as powers of three 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 variables (energy inside/outside, cilia length, number of cilia) could bridge the gap between the results obtained from these two types of models.
Clearly, maximising the hydrodynamic swimming efficiency is not the sole objective for biological microswimmers. Other functions such as generating feeding currents (Riisgård & Larsen Reference Riisgård and Larsen2010; Pepper et al. Reference Pepper, Roper, Ryu, Matsumoto, Nagai and Stone2013) and creating flow environment to accelerate mixing for chemical sensing (Supatto, Fraser & Vermot Reference Supatto, Fraser and Vermot2008; Shields et al. Reference Shields, Fiser, Evans, Falvo, Washburn and Superfine2010; Ding et al. Reference Ding, Nawroth, McFall-Ngai and Kanso2014; Nawroth et al. Reference Nawroth, Guo, Koch, Heath-Heckman, Hermanson, Ruby, Dabiri, Kanso and McFall-Ngai2017) are also important factors to consider as a microswimmer. The effect of such multitasking 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.
A straightforward extension of our work is to allow more general ciliary motions, e.g. including deformations normal to the surface. Such a swimmer will display time-periodic shape changes and the optimisation will require the derivation of shape sensitivities. Additionally, the computational cost would also increase significantly because the matrix in (2.12) needs to be updated at every time step. Our framework is also open to many generalisations and could, for example, help in accounting for the multiple factors mentioned previously, such as mixing for chemical sensing, in the study of optimal ciliary dynamics.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2021.744.
Funding
The authors gratefully acknowledge support from NSF under grants DMS-1719834, DMS-1454010 and DMS-2012424.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivations of sensitivities
Here, we include the detail derivations that lead to (2.33) and the explicit expressions of the terms therein.
Recall that the power loss and the swim speed can be written as functionals of 
$\psi$, as shown in (2.30). The sensitivities of 
$\langle \mathbb {P}\rangle$ and 
$\langle \mathbb {U}\rangle$ can thus be formulated by considering perturbed versions of 
$\psi$ as in
\begin{equation} \psi_{\eta}(x,t) = \psi(x,t) + \eta \hat{\psi}(x,t), \quad \text{i.e. } \psi_{\eta} = \psi + \eta \hat{\psi}, \end{equation}
so that the perturbed location 
$s_{\eta }$ at time 
$t$ of the material particle initially located at 
$s_0$ is given by
\begin{equation} s_{\eta} = \alpha(s_0,\psi_{\eta}), \end{equation}
with the functional 
$\alpha$ unchanged. Similar to (2.29), the perturbed slip velocity 
$u^{S}_{\eta }(s,t)$ satisfies
\begin{equation} u^{S}_{\eta}(s,t) = \partial_\psi\alpha\left( \beta(s,\psi_{\eta}),\psi_{\eta};\dot{\psi}_{\eta} \right) = \upsilon^{S}(s,\psi_{\eta}), \end{equation}
where 
$\beta$, the inverse function of 
$\alpha$, is also unchanged.
Notice that 
$u^{S}$ and 
$u^{S}_{\eta }$ given by (2.29) and (A3) are evaluated at the same time 
$t$ and current location 
$s$ (the latter being thus reached from different initial positions 
$\beta (s,\psi )$ and 
$\beta (s,\psi _{\eta })$). This allows us to define the directional derivative 
$\upsilon ^{S}{}'(s,\psi ;\hat {\psi })$ of 
$u^{S}$ with respect to 
$\psi$ in the direction 
$\hat {\psi }$, as a total derivative with respect to 
$\eta$:
\begin{equation} \upsilon^{S}{}'(s,\psi;\hat{\psi}) := \lim_{\eta\to0} \frac{1}{\eta} \left[ u^{S}_{\eta}(s,t) - u^{S}(s,t) \right] = \frac{\text{d}}{\text{d}\eta} \partial_\psi\alpha \left.\left( \beta(s,\psi_{\eta}),\psi_{\eta};\dot{\psi}_{\eta} \right) \right|_{\eta=0}. \end{equation}Carrying out the above differentiation in a straightforward way, we find
$$\begin{gather} \upsilon^{S}{}'(s,\psi;\hat{\psi}) = \partial_{\psi s}\alpha\left( \beta(s,\psi),\psi;\dot{\psi} \right)\; \partial_{\psi}\beta\left( s,\psi;\hat{\psi} \right) \nonumber\\ + \partial_{\psi\psi}\alpha\left( \beta(s,\psi),\psi;\dot{\psi}\,,\, \hat{\psi} \right) + \partial_\psi\alpha\left( \beta(s,\psi),\psi;\dot{\hat{\psi}} \right). \end{gather}$$
Moreover, for any 
$\psi$, the functions 
$\alpha$ and 
$\beta$ are linked through
\begin{equation} s = \alpha\left( \beta(s,\psi),\psi \right), \end{equation}
which, upon taking the directional derivative in the direction 
$\hat {\psi }$ and using the chain rule, yields
\begin{equation} 0 = \partial_s\alpha\left( \beta(s,\psi),\psi \right) \partial_{\psi}\beta\left( s,\psi;\hat{\psi} \right) + \partial_{\psi}\alpha\left( \beta(s,\psi),\psi;\hat{\psi} \right). \end{equation}
The above equality allows us to eliminate 
$\partial _{\psi }\beta$ from (A5), to obtain
$$\begin{gather} \upsilon^{S}{}'(s,\psi;\hat{\psi}) ={-}\partial_{\psi s}\alpha\left( \beta(s,\psi),\psi;\dot{\psi} \right)\; \frac{\partial_{\psi}\alpha\left( \beta(s,\psi),\psi;\hat{\psi} \right)} {\partial_s\alpha\left( \beta(s,\psi),\psi \right)} \nonumber\\ + \partial_{\psi\psi}\alpha\left( \beta(s,\psi),\psi\;;\;\dot{\psi}\,,\, \hat{\psi} \right) + \partial_\psi\alpha\left( \beta(s,\psi),\psi;\dot{\hat{\psi}} \right). \end{gather}$$In practice, the slip velocity derivative 
$\upsilon ^{S}{}'$ given by (A8) is more conveniently expressed in the initial arc length variable 
$s_0=\beta (s,\psi )$. Moreover, in the event that 
$\psi (s_0,t) = 0$ for some 
$s_0$ and 
$t$, 
$\upsilon ^{S}{}'$ given by (A8) blows up because 
$\partial _s \alpha (\beta (s,\psi ),\psi ) = 0$ in this case, whereas 
$\upsilon ^{S}{}'\,\mathrm {d}s$ remains finite if expressed in terms of 
$s_0$ (because 
$\mathrm {d}s = \partial _s\alpha (s_0,\psi ) \,\mathrm {d}s_0$). Upon effecting the change of variable 
$s=\alpha (s_0,\psi )$ in the integrals (2.26) and (2.28), we obtain
$$\begin{gather} \langle\mathbb{P}\rangle'(\psi;\hat{\psi}) = 4{\rm \pi} \left\langle\int_\gamma R(\alpha(s_0,\psi))\, \boldsymbol{f}(\alpha(s_0,\psi),t)\cdot\boldsymbol{\tau}(\alpha(s_0,\psi))\,\upsilon^{S}{}'(s,\psi;\hat{\psi}) \, \partial_s\alpha(s_0,\psi) \,\mathrm{d}s_0\right\rangle , \end{gather}$$
$$\begin{gather}\langle\mathbb{U}\rangle'(\psi;\hat{\psi}) = \frac{-2{\rm \pi}}{F_0}\left\langle \int_\gamma R(\alpha(s_0,\psi))\, \hat{\boldsymbol{f}}(\alpha(s_0,\psi))\cdot\boldsymbol{\tau}(\alpha(s_0,\psi))\,\upsilon^{S}{}'(s,\psi;\hat{\psi}) \, \partial_s\alpha(s_0,\psi) \,\mathrm{d}s_0\right\rangle , \end{gather}$$where, thanks to (A8), we have used
\begin{align}
\upsilon^{S}{}'(s,\psi;\hat{\psi}) \,\mathrm{d}s &=
\upsilon^{S}{}'(s,\psi;\hat{\psi}) \,
\partial_s\alpha(s_0,\psi) \,\mathrm{d}s_0,
\nonumber\\ &= \left\{ \partial_s\alpha(s_0,\psi)\,\left[
\partial^2_\psi \alpha \left(
s_0,\psi;\hat{\psi},\dot{\psi}\right) + \partial_\psi
\alpha \left( s_0,\psi;\dot{\hat{\psi}} \right)\,\right]
\right.\nonumber\\ &\left.\quad -\partial_{\psi
s}\alpha\left( s_0,\psi;\dot{\psi} \right)\;
\partial_{\psi}\alpha\left( s_0,\psi;\hat{\psi} \right)\right\} \mathrm{d}s_0. \end{align}
This completes our derivation of (2.33).
For the ciliary motion (2.14) used here, introducing the shorthand notation 
$I(f,g;s):=\int _0^s f(x)g(x) \,{\textrm {d}\kern0.04em x}$, we have
\begin{align} \alpha(s_0,\psi)&= \frac{\ell I(\psi,\psi;s_0)}{I(\psi,\psi;\ell)}, \end{align}
\begin{align} \partial_s\alpha(s_0,\psi)&= \frac{\ell\psi^2(s_0)}{I(\psi,\psi;\ell)}, \end{align}
\begin{align} \partial_{\psi}\alpha\left( s_0,\psi;\hat{\psi} \right)&= \frac{2\ell I(\psi,\hat{\psi};s_0)}{I(\psi,\psi;\ell)} - 2\alpha(s_0,\psi)\frac{I(\psi,\hat{\psi};\ell)}{I(\psi,\psi;\ell)}, \end{align}
\begin{align} \partial_{s\psi}\alpha\left( s_0,\psi;\dot{\psi} \right) &= \frac{2\ell\psi(s_0)\dot{\psi}(s_0)}{I(\psi,\psi;\ell)} - 2\ell\frac{I(\psi,\dot{\psi};\ell)\, \psi^2(s_0)}{\left( I(\psi,\psi;\ell) \right)^2}, \end{align}
\begin{align} \partial^2_\psi \alpha \left( s_0,\psi\;;\;\hat{\psi},\dot{\psi} \right) &= \frac{2\ell I(\hat{\psi},\dot{\psi};s_0)}{I(\psi,\psi;\ell)} -2\alpha(s_0, \psi) \frac{I(\hat{\psi},\dot{\psi};\ell)}{I(\psi,\psi;\ell)} \nonumber\\ &\quad -\frac{2I(\psi,\hat{\psi};\ell)}{I(\psi,\psi;\ell)} \partial_{\psi}\alpha\left( s_0,\psi;\dot{\psi} \right) -\frac{2I(\psi,\dot{\psi};\ell)}{I(\psi,\psi;\ell)} \partial_{\psi}\alpha\left( s_0,\psi;\hat{\psi} \right). \end{align}Appendix B. Initial coefficient sensitivity
In our optimisations, the initial guesses are chosen to be low-frequency waves with small wave amplitudes. These are obtained by choosing the coefficients of the first Fourier modes from a uniform distribution within 
$[0,0.01]$ (to restrict the initial wave amplitudes), and setting the coefficients of the higher modes to 0 (to discourage high-frequency waves).
Restricting our attention to low-frequency waves effectively sets a time scale in our problem. That is, it helps us to focus on the ciliary motion within sone beating cycle which is given by the base Fourier mode. Note that there is a danger of confusing the (spatial) Legendre modes used in Blake (Reference Blake1971) and the (temporal) Fourier modes studied here. While the swim speed is determined by the first Legendre mode, introducing higher order Fourier modes would affect the swim speed. Specifically, cilia beating twice as fast (beating two cycles in the same time span) could double the swim speed. However, the efficiency would remain unchanged because of the simultaneous increase of the power loss.
Owing to the high-dimensional nature of the problem (hundreds of degrees of freedom), many local optima exist. As shown in figure 7(a), a large initial range of the Fourier coefficient (e.g. 
$[0,1]$) increases the risk of the optimiser getting stuck close to an unsuitable local optimum. For example, an initial waveform as shown in figure 7(c) can only be optimised to a waveform shown in figure 7(e), which has a swimming efficiency as low as 
$2\,\%$. On the other hand, the initial wave with small amplitudes (as shown in figure 7b) could almost always be optimised to the waveform with swimming efficiency 
$\epsilon \approx 35\,\%$.
Figure 7. Sensitivity to the initial Fourier coefficient. (a) Optimised efficiencies for the unconstrained spherical swimmer with the initial first Fourier mode chosen from 
$[0,0.01]$ (blue), 
$[0,0.1]$ (black), 
$[0,1]$ (red), respectively. (b,d) The initial and final waveforms of the case where the range is 
$[0,0.01]$. (c,e) The initial and final waveforms of the case where the range is 
$[0,1]$. Results with different initial conditions are highlighted by their corresponding colours.
Open access
