1. Introduction
Microscopic biological organisms have adapted to a viscous world in which their inertia is inconsequential to their locomotion. For a typical micro-scale organism the Reynolds number 
$Re = \rho U L/\eta$ is small. For example, Escherichia coli has a characteristic length 
$L \sim 10\ \mathrm {\mu }\textrm {m}$ and a characteristic swimming speed 
$U\sim 10\ \mathrm {\mu }\textrm {m}\ \textrm {s}^{-1}$ in water (density 
$\rho \approx 10^{3}\ \textrm {{kg}}\ \textrm {m}^{-3}$ and dynamic viscosity 
$\eta \approx 10^{-3}\ \textrm {Pa}\ \textrm {s}$), leading to a negligibly small Reynolds number 
$Re = \rho U L/\eta \sim 10^{-5}\text {--}10^{-4}$ (Purcell Reference Purcell1977). At this scale, a swimmer reacts instantaneously to any forces, oblivious to any history of prior dynamics (Purcell Reference Purcell1977; Childress Reference Childress1981; Lauga & Powers Reference Lauga and Powers2009), and so the primary method of locomotion for macroscopic swimmers such as fish and humans, which relies on trading momentum with the fluid to generate a propulsive force, is ineffective here. Rather, the net translations of a microorganism are determined by the sequence of configurations it adopts to swim, independent of its deformation rate. Microswimmers must continually paddle or deform their bodies in a swimming pattern with non-reciprocal forward and reverse strokes to manipulate the drag forces for propulsion (Purcell Reference Purcell1977; Childress Reference Childress1981; Lauga & Powers Reference Lauga and Powers2009).
It is common for microorganisms to swim in geometrically confined media: in channels, near surfaces and interfaces and in films. A significant amount of theoretical and experimental work has been devoted to studying the effects of confinement on the motion of microscopic swimmers near solid walls (Pedley & Kessler Reference Pedley and Kessler1987; Lauga et al. Reference Lauga, DiLuzio, Whitesides and Stone2006; Berke et al. Reference Berke, Turner, Berg and Lauga2008; Drescher et al. Reference Drescher, Leptos, Tuval, Ishikawa, Pedley and Goldstein2009; Li & Tang Reference Li and Tang2009; Or & Murray Reference Or and Murray2009; Crowdy & Or Reference Crowdy and Or2010; Li et al. Reference Li, Bensson, Nisimova, Munger, Mahautmr, Tang, Maxey and Brun2011; Or, Zhang & Murray Reference Or, Zhang and Murray2011; Spagnolie & Lauga Reference Spagnolie and Lauga2012; Molaei et al. Reference Molaei, Barry, Stocker and Sheng2014; Ishimoto, Cosson & Gaffney Reference Ishimoto, Cosson and Gaffney2016), near fluid–fluid interfaces (Guasto, Johnson & Gollub Reference Guasto, Johnson and Gollub2010; Di Leonardo et al. Reference Di Leonardo, Dell'Arciprete, Angelani and Iebba2011; Wang & Ardekani Reference Wang and Ardekani2013; Lopez & Lauga Reference Lopez and Lauga2014; Masoud & Stone Reference Masoud and Stone2014; Stone & Masoud Reference Stone and Masoud2015) and in thin fluid layers atop a solid substrate (Lambert et al. Reference Lambert, Picano, Breugem and Brandt2013; Mathijssen et al. Reference Mathijssen, Doostmohammadi, Yeomans and Shendruk2016a,Reference Mathijssen, Doostmohammadi, Yeomans and Shendrukb; Ota et al. Reference Ota, Hosaka, Yasuda and Komura2018).
Motivated by recent experimental and theoretical studies of bacteria swimming in biofilms, in freely suspended thin films (Aranson et al. Reference Aranson, Sokolov, Kessler and Goldstein2007; Sokolov et al. Reference Sokolov, Aranson, Kessler and Goldstein2007) and on active proteins mimicking biological swimmers in lipid membranes (Huang, Chen & Mikhailov Reference Huang, Chen and Mikhailov2012), we study here the hydrodynamics of swimming microorganisms in a thin membrane. We treat the membrane as a continuous incompressible viscous fluid film of very small thickness. Flow fields in such a membrane are uniform throughout the thickness of the membrane. In contrast to a thin-film model that involves integration of three-dimensional (3-D) hydrodynamic equations across the thickness of the film, our membrane model is intrinsically two-dimensional (2-D), in the sense that the motion of molecules within the membrane in the direction normal to the plane of the membrane is forbidden. The membrane is embedded in a 3-D fluid of different viscosity. The motion of a swimmer in the membrane generates flows both in the membrane and in the surrounding fluid (see figure 1). While there have been some recent investigations of the hydrodynamics of swimmers in a thin layer of fluid sandwiched between fluids of a different viscosity (Leoni & Liverpool Reference Leoni and Liverpool2010; Rower, Padidar & Atzberger Reference Rower, Padidar and Atzberger2019), this problem is still largely unexplored.
Figure 1. Illustration of a flagellum confined to a plane of a thin membrane of 2-D viscosity 
$\eta _m$ sandwiched between two semi-infinite slabs of bulk fluid of 3-D viscosity 
$\eta$. The flagellum propagates transverse planar waves travelling with speed 
$c$ with respect to the flagellum.
As was demonstrated by Saffman & Delbrück (Reference Saffman and Delbrück1975) and Saffman (Reference Saffman1976), the amount of momentum imparted by the membrane to the bulk fluid is controlled by a hydrodynamic length scale, the so-called Saffman length 
$\ell _S$, given by the ratio of the 2-D membrane viscosity 
$\eta _m$ to the 3-D viscosity of the bulk fluid 
$\eta$, 
$\ell _S = \eta _m/(2\eta )$. If the membrane is perturbed by a localized force (applied in the plane of the membrane) at a point 
$\boldsymbol {x}$, for distances (measured from 
$\boldsymbol {x}$) much smaller than the Saffman length, the effect of the flows in the bulk on the membrane hydrodynamics is negligible. In this region, the fluid velocity in the membrane decays slowly (logarithmically) with distance, as in purely 2-D fluids. On the other hand, for distances much larger than the Saffman length, the contribution of the bulk fluid to the membrane dynamics is significant, and the membrane flow field decays inversely with the distance, a behaviour consistent with a 3-D dynamics.
Levine & MacKintosh (Reference Levine and MacKintosh2002) (LM) derived a Green function for a more general case of viscoelastic membranes. In the case of a purely viscous membrane that we consider here, there is no elastic response of the membrane, and a disturbance caused by a force results in the velocity field alone. The velocity of the membrane at position 
$\boldsymbol {x}^{\prime }$ due to an in-plane, localized force 
$\boldsymbol {f}(\boldsymbol {x}) = \boldsymbol {f}\delta (\boldsymbol {x})$ is determined by the LM response tensor 
$\boldsymbol {\alpha }(\boldsymbol {x} - \boldsymbol {x}')$,
\begin{equation} \boldsymbol{v}(\boldsymbol{x}) = \frac{1}{4{\rm \pi}\eta_m} \boldsymbol{\alpha}(\boldsymbol{x} - \boldsymbol{x}') \boldsymbol{\cdot} \boldsymbol{f}(\boldsymbol{x}'). \end{equation}
Here, 
$\boldsymbol {x}$ and 
$\boldsymbol {x}'$ are in-plane vectors with components 
$(x, y)$ and 
$(x^{\prime }, y^{\prime })$, respectively (refer also to figure 1 for our choice of the coordinate system). The response function 
$\boldsymbol {\alpha }(\boldsymbol {x} - \boldsymbol {x}')$ in (1.1) plays the role of the Oseen tensor in 
$3$-D hydrodynamics. The coupling between the membrane and the embedding bulk fluid is implicitly incorporated in 
$\boldsymbol {\alpha }(\boldsymbol {x} - \boldsymbol {x}')$.
As was shown by LM, the response function may be split into ‘parallel’ and ‘transverse’ contributions. In the component form we have
\begin{equation} \alpha_{\alpha\beta}(\boldsymbol{x}) = \alpha_{\|}(|\boldsymbol{x}|) \hat{x}_\alpha \hat{x}_\beta + \alpha_\perp(|\boldsymbol{x}|) [\delta_{\alpha\beta} - \hat{x}_\alpha\hat{x}_\beta], \end{equation}
where 
$\alpha ,\beta =x,y$, and 
$\hat {x}_{\alpha }$ is the 
$\alpha$ component of the unit vector 
$\hat {\boldsymbol {x}} =\boldsymbol {x}/|\boldsymbol {x}|$. In our notation, 
$\alpha _{\alpha \beta }$ corresponds to 
$-\textrm {i}\omega \alpha _{\alpha \beta }$ in the LM theory. The scalar functions 
$\alpha _{\|}$ and 
$\alpha _\perp$ are given by
\begin{equation} \left.\begin{gathered} \alpha_{\|}(\kappa)= \frac{\rm \pi}{\kappa} \boldsymbol{H}_1(\kappa) - \frac{2}{\kappa^{2}} - \frac{\rm \pi}{2}[Y_0(\kappa) + Y_2(\kappa)], \\ \alpha_\perp(\kappa) = {\rm \pi}\boldsymbol{H}_0(\kappa) - \frac{\rm \pi}{\kappa}\boldsymbol{H}_1(\kappa) + \frac{2}{\kappa^{2}} - \frac{\rm \pi}{2}[Y_0(\kappa)- Y_2(\kappa)], \end{gathered}\right\} \end{equation}
where 
$\boldsymbol {H}_\nu$ are Struve functions and 
$Y_\nu$ are Bessel functions of the second kind (Abramowitz & Stegun Reference Abramowitz and Stegun1965); 
$\kappa =|\boldsymbol {x}|/\ell _S$ is the non-dimensionalized distance between the point of application of the force and the point where the membrane velocity response is measured. Both 
$\alpha _\|(\kappa )$ and 
$\alpha _\perp (\kappa )$ diverge logarithmically as 
$\kappa \to 0$, while for large 
$\kappa$ we have 
$\alpha _\|(\kappa )\sim 1/\kappa$ and 
$\alpha _\perp (\kappa )\sim 1/\kappa ^{2}.$
In the small Reynolds number regime the inertia term in the Navier–Stokes equation can be neglected. We assume that the membrane has thickness 
$h$ and choose a coordinate system with the origin at the top surface of the membrane with 
$z=0$ (therefore, the bottom side of the membrane is at 
$z=-h$). The dynamics of the membrane embedded in a bulk fluid is governed by a modified Stokes equation and the incompressibility condition (Saffman Reference Saffman1976),
\begin{equation} {-}\boldsymbol{\nabla} p + \frac{\eta_m}{h} \nabla^{2} \boldsymbol{v}+ \frac{2\boldsymbol{f}}{h} = 0, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}=0, \end{equation}
where 
$p$ and 
$\boldsymbol {v}$ are the pressure and velocity fields of the membrane. The flows in the membrane set the bulk fluid into motion. The resulting flows in the embedding fluid, in their turn, exert traction on the membrane. Since the membrane is only a few molecular layers thick, the traction due to the bulk fluid produces a flow that is uniform throughout the thickness of the membrane, i.e. the fluid velocity does not depend on 
$z$ for 
$-h < z < 0$. In (1.4a,b) the coupling between the membrane and the bulk fluid is described by the force per unit volume 
$2\boldsymbol {f}/h$, with
\begin{equation} \boldsymbol{f} = \left.\eta \frac{\partial \boldsymbol{v}^{(3D)}}{\partial z}\right|_{z=0}, \end{equation}
where 
$\boldsymbol {v}^{(3D)}$ is the bulk fluid velocity. The factor of two in (1.4a,b) is due to an equal force 
$\boldsymbol {f}$ acting on the bottom side of the membrane. Equation (1.4a,b) can be written in a more compact form that we will use later,
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} = - \frac{2 \boldsymbol{f}}{h}, \end{equation}
where 
$\boldsymbol {\sigma }$ is the stress field of the membrane.
Being inertialess and swimming in the absence of external influences, a swimmer must maintain a zero net force 
$\boldsymbol {F}(t)$ and a zero torque 
$\boldsymbol {L}(t)$ on its body at every time instant,
\begin{gather} \boldsymbol{F}(t) = \int_{S} \boldsymbol{\sigma}\boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S, \end{gather}
\begin{gather}\boldsymbol{L}(t) = \int_{S} \boldsymbol{x} \times (\boldsymbol{\sigma}\boldsymbol{\cdot} \boldsymbol{n})\, \mathrm{d} S, \end{gather}
where the integration is over the surface of the swimmer and 
$\boldsymbol {n}$ is a unit vector normal to the surface and pointing away from the swimmer.
Many microorganisms such as spermatozoa, E. coli, and Caulobacter crescentus swim by moving thin extensions (flagella) on their bodies. Some creatures, like Paramecium, are covered in thousands of short hair-like appendages called cilia and propel themselves through a coordinated beating of these cilia. Since our primary goal is to study how the confinement to the plane of a membrane affects the swimming dynamics of a microorganism (rather than a detailed study of a particular microorganism), we consider here only minimal theoretical models of flagellated and ciliated microorganisms.
In § 2 we consider a headless, infinitely long ‘flagellum’ of infinitesimally small thickness propagating planar sinusoidal waves along its body. This is a 1-D analogue of the Taylor swimmer (Taylor Reference Taylor1951), an infinite plane in viscous fluid passing transverse sinusoidal waves. We recover Taylor's result for the swimming velocity in the limiting case of a pure 2-D hydrodynamics (the membrane in vacuum). We find that the membrane incompressibility condition imposes a constraint on the fluid dynamics that allows the flagellum to achieve much higher swimming speed than in pure 2-D and 3-D fluids for large ratios of the wavelength to the Saffman length.
In § 3 we study the propulsion of a flagellum of finite length and find its swimming speed and efficiency. In §§ 2 and 3 we apply the boundary-element method (BEM) that two of us have recently developed in work on hydrodynamic interaction of inclusions in freely suspended smectic films (Qi et al. Reference Qi, Nguyen, Park, Glaser, Maclennan, Clark, Kuriabova and Powers2014; Kuriabova et al. Reference Kuriabova, Powers, Qi, Goldfain, Park, Glaser, Maclennan and Clark2016; Qi et al. Reference Qi, Ferguson, Sechrest, Munsat, Park, Glaser, Maclennan, Clark, Kuriabova and Powers2017).
In § 4 we formulate the Lorentz reciprocal theorem for a quasi-2-D fluid and derive an equation for the swimming speed. We discuss the advantage of the method based on the Lorentz reciprocal theorem over the BEM for studying microscopic swimmers with swimming patterns that do not change the overall shape of the swimmers’ bodies (for example, swimmers propagating longitudinal compressive waves along their bodies). As an example of such a creature, we consider a simple model of a 2-D ‘squirmer’, a disk with a prescribed tangential velocity along its circumference. We find that, unlike its flagellated counterpart, a squirmer does not benefit from the presence of the bulk fluid: its swimming speed is lower than that in a purely 2-D fluid.
In § 5 we discuss our results and suggest further directions of investigation.
2. Infinitely long flagellum in a quasi-2-D membrane
The Taylor swimmer (Taylor Reference Taylor1951) is an infinite swimming 2-D sheet in a 3-D viscous fluid that propagates transverse waves of amplitude 
$b$ and wave speed 
$c = \omega /q$. In a frame moving with the swimming sheet (co-moving frame) the shape of the sheet is described by
\begin{equation} y = b \sin(q x - \omega t) = b\sin\xi, \end{equation}
where 
$\xi = q x - \omega t$ denotes the wave phase.
Taylor showed that the sheet with such a 1-D modulation travels, relative to the fluid at infinity, with a speed 
$U/c = (bq)^{2}/2 + O((bq)^{4})$ in the direction opposite to the wave velocity. We consider here a 1-D analogue of the Taylor swimmer: an infinitely long, infinitesimally thin flagellum confined to the plane of a viscous membrane embedded in bulk fluid (see figure 1) with prescribed motion given by (2.1), with 
$x$ and 
$y$ parametrizing the shape of the flagellum.
The swimming velocity of an infinitely long flagellum is time independent. Indeed, two snapshots of the waving flagellum taken at the same point on the 
$x$-axis differ only by a shift 
$\Delta x$ along the 
$x$-axis. Thus, a temporal shift at a fixed point 
$x$ is equivalent to a spatial displacement along the flagellum. A swimmer moving as a whole has the same translational velocity along its entire length. The swimming velocity, being invariant with respect to translations along the 
$x$-axis, must be invariant with respect to translations in time as well. Therefore, we can calculate the swimming speed for a single time instant and set 
$t=0$ in (2.1).
As in Taylor (Reference Taylor1951), we consider the case of an inextensible flagellum. In order to calculate the portion of the swimmer's velocity due to its distortion, we calculate the position of a material point of the flagellum as a function of time. In a frame moving with the wave (with speed 
$c$ relative to the co-moving frame) the shape of the flagellum does not change. The Cartesian coordinates for this frame are 
$(x',y)$, where 
$x'=x-ct$. In this reference frame a material particle of the flagellum travels a distance 
$\varLambda$ equal to the arclength of the flagellum spanned by one (linear) wavelength 
$\lambda$ during one period of oscillation 
$T = 2{\rm \pi} /\omega$,
\begin{equation} \varLambda = \frac{1}{q}\int_0^{2{\rm \pi}} \sqrt{1+ (bq)^{2} \cos^{2} \xi} \,\mathrm{d} \xi. \end{equation}
We will call 
$\varLambda$ the arcwise wavelength. The material particle's speed is, therefore,
\begin{align} C &= \frac{c}{2{\rm \pi}} \int_0^{2{\rm \pi}} \sqrt{1+ (bq)^{2} \cos^{2} \xi} \,\mathrm{d} \xi \end{align}
\begin{align} &\approx c\left(1+\frac{1}{4}b^{2}q^{2}-\frac{3}{64}b^{4}q^{4}\right). \end{align}
To determine the position of a material particle of the flagellum, we define the material coordinate 
$S$ to be the arclength coordinate 
$s$ of a material point at 
$t=0$. In the frame in which the nodes of the wave are fixed, arclength is related to the Cartesian coordinate 
$x'$ by
\begin{align} s&=\int_0^{x'}\sqrt{1+(bq)^{2}\cos^{2}(qx')}\,\mathrm{d} x' \end{align}
\begin{align} &\approx x'+\frac{b^{2}q}{8}[2qx'+\sin(2qx')]-\frac{b^{4}q^{3}}{256} [12qx'+8\sin(2qx')+\sin(4qx')]. \end{align}Reverting the series leads to
\begin{align} x'\approx s-\frac{b^{2}q}{8}[2qs+\sin(2qs)]+\frac{b^{4}q^{3}}{256}[28qs+16qs\cos(2qs)+16\sin(2qs)+5\sin(4qs)]. \end{align}
Using 
$y=b\sin qx'$ and 
$s=S-Ct$ leads to the position of the material point labelled by 
$S$ as a function of time 
$t$
\begin{align} x&\approx S-\frac{b^{2}q}{8}[2q S+\sin 2(qS-\omega t)]\nonumber\\ &\quad +\frac{b^{4}q^{3}}{256}[28qS+16qS\cos 2(qS-\omega t)+16\sin 2(qS-\omega t)+5\sin 4(qS-\omega t)],\end{align}
\begin{align} y&\approx b\sin(q S-\omega t)\nonumber\\ &\quad -\frac{b^{3}}{16}[4q^{3}S\cos(qS-\omega t)+q^{2}\sin(q S-\omega t)+q^{2}\sin 3(qS-\omega t)]. \end{align}
In the co-moving frame the components of a material particle's velocity 
$\boldsymbol {u}_S$ are given by
\begin{equation} {\boldsymbol{u}}_S=\left(\left.\frac{\partial{x}}{\partial t}\right|_S, \left.\frac{\partial{y}}{\partial t}\right|_S\right).\end{equation}
For an arbitrary value of 
$b$, the material particle's velocity can be calculated numerically using
\begin{gather} u_{S, x} =- C\cos\theta_S + c, \end{gather}
\begin{gather}u_{S, y} =-C\sin\theta_S, \end{gather}
with 
$\tan \theta _S = {\partial y}/{\partial x}|_S$. The total velocity of a material particle relative to the fluid at infinity (for 
$y\to \pm \infty$) is the sum of the surface disturbance and swimming velocities, 
${\boldsymbol {u}}_S + {\boldsymbol {U}}$.
The linearity of Stokes equations allows us to model the fluid velocity field in the membrane as a superposition of fluid velocities due to a (yet unknown) force density 
$\boldsymbol {f}(s)$ along the flagellum
\begin{equation} \boldsymbol{v}(\boldsymbol{x}) = \frac{1}{4{\rm \pi}\eta_m} \int_{\varGamma} \boldsymbol{\alpha}(\boldsymbol{x} - \boldsymbol{x}^{\prime}(s')) \boldsymbol{\cdot} \boldsymbol{f}(s')\,\mathrm{d} s^{\prime}, \end{equation}
where 
$\boldsymbol {v}(\boldsymbol {x})$ is the fluid velocity at an arbitrary point 
$\boldsymbol {x}$ on the membrane and the integration is along the (infinite) contour of the flagellum.
The spatial periodicity of the flagellum modulation implies the invariance of the flow field and the force density 
$\boldsymbol {f}(s)$ in (2.13) under translations along the 
$x$-axis by an integer multiple of the wavelength 
$\lambda =2{\rm \pi} /q$. Defining 
$x^{\prime }_m= x^{\prime } + m \lambda$, for all integers 
$m$, the integration on the right-hand side of (2.13) can be reduced to integration over a single wavelength,
\begin{equation} \boldsymbol{v}(\boldsymbol{x}) = \frac{1}{4{\rm \pi}\eta_m} \int_{\varGamma_0} \sum_{m=-\infty}^{\infty} \boldsymbol{\alpha}(\boldsymbol{x} - \boldsymbol{x}_m^{\prime}(s'))\boldsymbol{\cdot} \boldsymbol{f}(s^{\prime}) \,\mathrm{d} s^{\prime}, \end{equation}
where 
$\boldsymbol {x}_m^{\prime } = ( x^{\prime } + m \lambda , y^{\prime } )$, and 
$\boldsymbol {x}^{\prime }$ indicate the points on the flagellum that belong to a one-wavelength ‘window’ 
$\varGamma _0$.
We impose a no-slip boundary condition on the surface of the flagellum by setting the fluid velocity equal to the velocity of the material point on the flagellum, 
$\boldsymbol {v}(\boldsymbol {x}) =\boldsymbol {u}_S(\boldsymbol {x}) + \boldsymbol {U}$, with the surface disturbance velocity 
$\boldsymbol {u}_S(\boldsymbol {x})$ given by (2.10). Taking into account (2.14), we have
\begin{equation} \boldsymbol{u}_S(\boldsymbol {x}(s)) + \boldsymbol{U} = \frac{1}{4{\rm \pi}\eta_m} \int_{\varGamma_0} \sum_{m=-\infty}^{\infty} \boldsymbol{\alpha}(\boldsymbol{x}(s) - \boldsymbol{x}_m^{\prime}(s'))\boldsymbol{\cdot} \boldsymbol{f}(s^{\prime}) \,\mathrm{d} s^{\prime}, \end{equation}
where the points 
$\boldsymbol {x}(s)$ and 
$\boldsymbol {x}^{\prime }(s')$ belong to the flagellum contour 
$\varGamma _0$.
To close the system of equations for the force density 
$\boldsymbol {f}(s)$ and the swimming velocity 
$\boldsymbol {U}$, we also require the net force on the flagellum be equal to zero. The net force on the flagellum is the sum of all the forces applied along the flagellum contour and implicitly includes the traction due to the bulk fluid embedding the membrane,
\begin{equation} \int_{\varGamma_0} \boldsymbol{f}(s) \,\mathrm{d}s = 0. \end{equation} We solved (2.15) and (2.16) numerically in Matlab by splitting the integration path into 
$N$ straight-line segments of equal length 
$\Delta s$ and replacing the line integrals in (2.15) and (2.16) by summation over the segments,
\begin{gather} \boldsymbol{u}_S(\boldsymbol{x}_i)+ \boldsymbol{U} = \frac{1}{4{\rm \pi}\eta_m} \sum_{j = 1}^{N} \sum_{m=-M}^{M} \boldsymbol{\alpha}(\boldsymbol{x}_i - \boldsymbol{x}_{mj})\boldsymbol{\cdot} \boldsymbol{f}(\boldsymbol{x}_{j}) \Delta s, \end{gather}
\begin{gather}\sum_{j = 1}^{N} \boldsymbol{f}(\boldsymbol{x}_{j}) = 0. \end{gather}
In (2.17) and (2.18) 
$\boldsymbol {x}_{i(j)}$ are the coordinates of the segments’ midpoints. In (2.17) we introduced the truncation parameter 
$M$ for the infinite sum over the wavelengths.
For a starting value of parameter 
$M$ (usually 
$M =10$), we ran computations for five different values of parameter 
$N$ in the range 
$300$–
$1000$ and then extrapolated our results for the swimming velocity to 
$\Delta s \to 0$ (
$N\to \infty$). We then gradually increased the value of 
$M$ and repeated the computations until the solution for 
$\boldsymbol {U}$ converged, showing changes smaller than 
$0.5\,\%$ with further increase of the number of terms in the sum over 
$m$. The computations required increasingly more terms in the sum over 
$m$ for large amplitudes (
$bq > 1$) and large Saffman lengths (
$\lambda /\ell _S \ll 1$) due to strong long-range hydrodynamic interactions between segments of the flagellum in this (nearly 2-D) regime, and correspondingly large contribution to the flow field by the forces 
$\boldsymbol {f}(\boldsymbol {x}_{j})$ separated by multiple wavelengths along the flagellum.
We paid special attention to the diagonal term with 
$m=0$ and 
$i=j$ in (2.17). This term gives the fluid velocity in the close proximity of a localized force 
$\boldsymbol {f}(\boldsymbol {x})$. The response tensor 
$\boldsymbol {\alpha }(\boldsymbol {x})$ diverges due to logarithmic singularities in the functions 
$\alpha _{\|}(\boldsymbol {x})$ and 
$\alpha _\perp (\boldsymbol {x})$ in the limit 
$\boldsymbol {x}\to 0$ (see (1.3) and (1.4a,b)). In the close proximity of a localized force the fluid velocity is parallel to the force and is, therefore, determined by the parallel component of the response function 
$\alpha _\|(\boldsymbol {x})$. We expanded 
$\alpha _\|(\boldsymbol {x})$ about 
$\boldsymbol {x}=0$ and performed integration analytically over 
$\Delta s$ in the vicinity of 
$\boldsymbol {x}=0$. Therefore, for the diagonal term on the right-hand side of (2.17), which we denote as 
$\boldsymbol {{\mathcal {A}}}^{m=0}_{i=j}$, we have
\begin{align} \boldsymbol{{\mathcal{A}}}^{m=0}_{i=j} (\boldsymbol{x}_i) &= \frac{1}{4{\rm \pi}\eta_m} \boldsymbol{f}(\boldsymbol{x}_{i}) \int_{-\Delta s/2}^{\Delta s/2} \alpha_\|( z ) \,\mathrm{d}z \nonumber\\ &= \frac{1}{4{\rm \pi}\eta_m} \boldsymbol{f}(\boldsymbol{x}_{i}) \, 2 \lim_{\varepsilon \to 0} \int_{\varepsilon}^{\Delta s /2} \left [ \frac{1}{2} - \gamma + \frac{2z}{3\ell_S} + \log \frac{2\ell_S}{z} \right ] \mathrm{d}z \nonumber\\ &= \frac{1}{4{\rm \pi}\eta_m} \boldsymbol{f}(\boldsymbol{x}_{i}) \, \Delta s \left [ \frac{3}{2} + \frac{\Delta s}{6\ell_S} - \gamma + \log\left(\frac{4\ell_S}{\Delta s} \right) \right ], \end{align}
where 
$\gamma = 0.5772$ is the Euler constant.
Our computations confirm that an infinitely long flagellum has a non-vanishing component of the swimming velocity only along the 
$x$-axis, as expected by symmetry. In figure 2 we plot the swimming speed as a function of the dimensionless amplitude 
$bq$ for a range of wavelengths scaled by the Saffman length. In the limit of a pure 2-D hydrodynamics, which corresponds to large Saffman lengths (and small scaled wavelengths, 
$\lambda /\ell _S \ll 1$), the energy dissipation occurs primarily in the membrane, and the membrane's viscous drag on the flagellum makes the main contribution to the flagellum's propulsion. In this limit, our computations are in good agreement with the 2-D problem of the Taylor swimming sheet, as expected because a Taylor ‘string’ in a thin very viscous membrane is equivalent to the Taylor sheet in 3-D bulk fluid. For small amplitudes (
$bq\ll 1$), we recover Taylor's leading-order perturbative solution 
$U/c = (1/2) (bq)^{2}$. For larger amplitudes 
$bq$ (and 
$\lambda /\ell _S \ll 1$), our calculated swimming speed is in agreement with recent analytic and computational results of Sauzade and coworkers (Sauzade, Elfring & Lauga Reference Sauzade, Elfring and Lauga2011).
Figure 2. Calculated swimming speed vs. 
$bq$ for various ratios 
$\lambda /\ell _S=2{\rm \pi} /(q\ell _S)$ of an infinitely long inextensible flagellum (a) scaled by the wave speed 
$c$, and (b) scaled by 
$c (bq)^{2}$. The coloured (grey) curves are the results of our BEM computations described in the text. The black dashed curve corresponds to the swimming speed in a purely 2-D fluid, and the solid black curve is the local drag theory result of Gray and Hancock (Gray & Hancock Reference Gray and Hancock1955) for a flagellum in a 3-D unbounded fluid oscillating with moderately large amplitude. The solid curve in the inset in panel (b) represents an estimate for the ratio of the drag constants as a function of the scaled wavelength, as discussed in the text. The wave amplitude in the inset was set to 
$bq = 10^{-3}$, and the dotted line corresponds to 
$y = \ln x + \mbox {const.}$ for reference.
In figure 2 the black dashed curve corresponds to a flagellum swimming in a pure 2-D membrane (no bulk fluid surrounding the membrane). We briefly outline our computations for this limiting case in appendix A. The computations are similar to the boundary integral approach demonstrated in (Sauzade et al. Reference Sauzade, Elfring and Lauga2011) with the only difference that we neglected the double layer potential contribution. The solid black curve in figure 2 corresponds to the local drag theory for an infinitely long flagellum passing waves of moderately large amplitudes in a 3-D unbounded fluid (Gray & Hancock Reference Gray and Hancock1955).
As can be seen in figure 2, our BEM computations predict that, for wavelengths larger than the Saffman length (
$\lambda /\ell _S > 1$), the swimming speed in a quasi-2-D membrane exceeds that in purely 2-D and 3-D fluids. For qualitative explanation of this result we compare our BEM computations with the local drag model of Gray & Hancock (Reference Gray and Hancock1955). In the local drag approximation, one assumes that the viscous drag force on a small segment of the flagellum is proportional to the segment's velocity, and the total drag on the swimmer is a sum over these local drag forces. Thus, the local drag approximation does not explicitly take into account the long-range hydrodynamic interactions between distant segments of the flagellum. The local drag forces for the motion of a rod-like segment parallel and perpendicular to its geometrical axis are given by 
$F_\| = \zeta _\| v_\|$ and 
$F_\perp = \zeta _\perp v_\perp$, respectively, with the drag coefficients 
$\zeta _\|$ and 
$\zeta _\perp$.
We expect our BEM results to be in qualitative agreement with the local drag approximation in the limit of 
$bq\ll 1$ and 
$\lambda /\ell _S \gg 1$. For small amplitudes 
$bq$ the segments of an inextensible flagellum separated by large contour distances (
${>}\lambda$) do not come too close to each other, and for the wavelengths larger than the Saffman length the spatial decay of the flow field is faster (
${\sim }1/r$), in comparison with slower (logarithmic) decay rate for 
$\lambda /\ell _S \ll 1$. Thus, in the regime of 
$bq\ll 1$ and 
$\lambda /\ell _S \gg 1$, the cooperativity effect between distant segments of the flagellum is expected to be small. Gray & Hancock (Reference Gray and Hancock1955) obtained the swimming velocity of an infinitely long and thin flagellum to the leading order of amplitude 
$bq$,
\begin{equation} \frac{U}{c} = \frac{(bq)^{2}}{2} \left(\frac{\zeta_\perp}{\zeta_\|} - 1\right). \end{equation}
In three dimensions, the ratio of the drag coefficients for an infinitely thin rod is 
$\zeta _\perp /\zeta _\| = 2$. For inclusions in quasi-2-D membranes, the ratio 
$\zeta _\perp /\zeta _\|$ depends on the Saffman length. In the inset of figure 2 we plot our BEM results for 
$(\zeta _\perp /\zeta _\|)_{estim} \equiv 1 + 2U/(cb^{2}q^{2})$ as a function of 
$\lambda /\ell _S$. According to (2.20), 
$(\zeta _\perp /\zeta _\|)_{estim}$ should give us an estimate for the local drag anisotropy. For the plot in the inset we chose a small amplitude 
$bq = 10^{-3}$, when the comparison with the local drag calculation of Gray and Hancock is justified. As can be seen in the inset of figure 2, the effective ratio 
$(\zeta _\perp /\zeta _\|)_{estim}$ grows logarithmically with 
$\lambda /\ell _S$ for 
$\lambda /\ell _S \gg 1$.
This result is in qualitative agreement with the work of Levine and collaborators (Levine, Liverpool & MacKintosh Reference Levine, Liverpool and MacKintosh2004), some of which we summarize here. Levine et al. studied the drag coefficients for a rod-like inclusion of length 
$L$ moving in a quasi-2-D membrane, and showed that for rod-like inclusions of lengths smaller than the Saffman length (
$L/\ell _S \ll 1$), where the viscous dissipation occurs primarily in the membrane, the dependence of the drag coefficients on the size and orientation of the rod is weak: 
$\zeta _\perp /\zeta _\| \to 1$. For longer rods with 
$L/\ell _S \gg 1$, the dissipation is governed by the 3-D fluid surrounding the membrane, and the drag coefficients show a stronger dependence on the size of the rods. Levine et al. found that the drag coefficient 
$\zeta _\|$ for a thin rod in a quasi-2-D membrane is qualitatively similar to that in three dimensions and is given by
\begin{equation} \zeta_\| = \frac{2{\rm \pi} \eta L}{\ln\left(0.43 L/\ell_S\right)}. \end{equation}
However, the dependence of 
$\zeta _\perp$ on 
$L$ in a quasi-2-D membrane is very different from that in three dimensions
\begin{equation} \zeta_\perp=2{\rm \pi}\eta L; \end{equation}
$\zeta _\perp$ depends on 
$L$ linearly, without the logarithmic factor in the denominator. The linear dependence of 
$\zeta _\perp$ on 
$L$ indicates the local character of the drag and the effective absence of hydrodynamic interactions between different sections of the rod.
As emphasized by Levine et al. (Reference Levine, Liverpool and MacKintosh2004), this behaviour of 
$\zeta _\perp$ arises from the incompressibility of the membrane, 
$\boldsymbol {\nabla }_\perp \boldsymbol {\cdot }\boldsymbol {v}_\perp = 0$, where the symbol 
$\perp$ denotes differentiation ‘in plane’ and the components of the velocity field in the plane. In the case of a rod moving perpendicular to its long axis in a 3-D fluid, the fluid can flow past the rod by moving over and under it. In this flow pattern, the in-plane part of the incompressibility condition does not vanish: 
$\partial _x v_x + \partial _y v_y\neq 0$. A 2-D version of such a flow (its projection on the 
$xy$-plane) is impossible due to the membrane incompressibility. In a quasi-2-D membrane, the fluid moves the long way around the rod, and the flow extends over distances comparable to the largest rod dimension 
$L$. The membrane incompressibility constraint only affects the perpendicular drag on a filament. A segment of filament being dragged parallel to its long axis does not produce divergent flows in a simple fluid, and thus the flow character is unchanged by the presence of a membrane.
Therefore, for long wavelengths (
$\lambda /\ell _S \gg 1$), the membrane incompressibility is expected to lead to a logarithmic growth of the flagellum's effective drag anisotropy, 
$\zeta _\perp /\zeta _\| \propto \log (\lambda /\ell _S)$. An organism that relies on the drag anisotropy for propulsion would achieve greater swimming speeds in a quasi-2-D membrane than in pure 2-D or 3-D fluids, as is confirmed by our BEM computations.
For larger values of the amplitude, 
$bq >1$, the distant parts of the flagellum come closer to each other, and the long-range hydrodynamic interactions between the segments of the flagellum cannot be ignored. Thus, the drag force on a segment of the flagellum depends not only on the velocity of that segment but also on the motion of other parts of the flagellum. Our assumption is that locally the drag on a small segment of the flagellum is still anisotropic (the drag force is not strictly antiparallel to the segment velocity), and the anisotropy is more pronounced for larger values of the scaled wavelength, 
$\lambda /\ell _S \gg 1$. However, the drag force is now a complicated function of the flagellum's surface disturbance velocity field.
For a qualitative explanation of the swimming velocity results for larger values of 
$bq$, we can think about two competing contributions to the flagellum motility. On the one hand, larger values of 
$bq$ correspond to steeper angles between the segments of the flagellum and the direction of wave propagation and, respectively, a larger drag, which tends to increase the swimming speed. On the other hand, the long-range hydrodynamic interactions tend to reduce the viscous drag and the swimming speed.
In a quasi-2-D membrane in the regime of 
$\lambda /\ell _S \gg 1$, the membrane incompressibility condition leads to an enhanced drag anisotropy (in comparison with the motion in a 3-D fluid). At the same time, the hydrodynamic interactions tending to slow down the swimmer are less pronounced for 
$\lambda /\ell _S \gg 1$ due to a fast (
$1/r$) spatial decay rate of the fluid flow (similar to that in a 3-D fluid) and a reduced cooperativity between the distant parts of the flagellum. This qualitatively explains why the swimming speed in a membrane is greater than that in a pure 2-D or 3-D fluid for large values of 
$\lambda /\ell _S \gg 1$. For smaller values of 
$\lambda /\ell _S$, the drag anisotropy (helping the creature to swim faster) is less pronounced and the cooperativity effects (reducing the swimming speed) are stronger due to slower (logarithmic) spatial decay rate of the flow field. This explains why the swimming speed approaches that of a flagellum in a pure 2-D fluid for 
$\lambda /\ell _S \ll 1$ and arbitrary values of 
$bq$.
3. Finite-length flagellum in a quasi-2-D membrane
We also applied the BEM to the case of an inextensible, headless, infinitely thin flagellum of finite length. As in the case of an infinitely long flagellum, the motion of the swimmer is prescribed by a sinusoidal modulation, 
$y(s, t) = b\sin (q x(s) - \omega t + \phi _0)$. Here, 
$s$ is the arclength along the flagellum measured from the flagellum's hypothetical ‘head’ and 
$\phi _0$ is the initial phase. At every time instant, the shape of the flagellum is described by the curve 
$\boldsymbol {X}(s, t)$, where 
$\boldsymbol {X}(s, t) = (X(s), Y(s)) = (x(s), y(s, t) - y(0, t))$ (see figure 3). The unit tangent to the curve is 
$\boldsymbol {T}(s) = (\mathrm {d} X/\mathrm {d} s, \mathrm {d} Y/\mathrm {d} s)$.
Figure 3. In the body frame the ‘head’ of the flagellum is motionless and is placed at the origin of the coordinate system. The position of a material particle is determined by the arclength 
$s$ measured from the left end of the flagellum. The tangent vectors 
$\boldsymbol {T}(s)$ describe the instantaneous shape of the flagellum. The flagellum propagates planar sinusoidal waves to the right.
In the frame of the flagellum, a material point at position 
$s$ moves with velocity 
$\boldsymbol {u}_S(s, t) = \partial \boldsymbol {X}(s, t)/\partial t$. As was demonstrated by Higdon (Reference Higdon1979), in the case of transverse waves propagating along the flagellum, the material particle's velocity can be calculated in a different manner. In a reference frame moving with the wave, the shape of the flagellum is given by 
$\boldsymbol {X}^{w}(s - C t)$, where 
$C$ is the arcwise speed that we introduced in (2.3). Following Higdon, we note that
\begin{equation} X^{w}(s + \varLambda) = X^{w}(s) + \lambda, \quad Y^{w}(s + \varLambda) = Y^{w}(s), \end{equation}
where 
$\varLambda$ is the arcwise wavelength (see (2.2)) and 
$\lambda$ is the linear wavelength. The tangential vectors are identical in the body and wave frames since the wave frame simply translates with respect to the body frame and does not undergo rotation. Thus, 
$\boldsymbol {T}^{w}(s - C t) = \boldsymbol {T}(s- C t)$. The velocity of a material particle at 
$s$ in the wave frame is calculated as 
$\boldsymbol {u}^{w} (s, t) = \partial \boldsymbol {X}^{w}(s - C t)/\partial t = - C \partial \boldsymbol {X}^{w}/\partial s = - C \boldsymbol {T}(s- C t)$. The velocity of the ‘head’ in the wave frame is 
$- C \boldsymbol {T}(s - C t)|_{s=0} = - C \boldsymbol {T}( - C t)$. Therefore, the wave frame translates with respect to the body frame with velocity 
$C \boldsymbol {T}( - C t)$. Therefore, we can calculate the velocity of the material point 
$s$ in the body frame as
\begin{equation} \boldsymbol{u}_S(s, t) = C \boldsymbol{T}( - C t) - C \boldsymbol{T}(s - C t). \end{equation}The material velocity with respect to the fluid at infinity is then
\begin{equation} \boldsymbol{u}_S(s, t) + \boldsymbol{U}(t) + \boldsymbol{\varOmega}(t) \times \boldsymbol{X}(s, t), \end{equation}
where 
$\boldsymbol {U}(t)$ and 
$\boldsymbol {\varOmega }(t)$ are the translational and angular velocities of the flagellum.
Similarly to our treatment of an infinitely long flagellum in § 2, we model the fluid velocity due to the flagellum's motion as a linear superposition of velocities due to a force density (see (2.13)). Now the path of integration 
$\varGamma$ stands for the curve 
$\boldsymbol {X}(s, t)$ describing the instantaneous shape of the flagellum. The instantaneous swimming and angular velocities and the force density are calculated from the coupled integral equations for a no-slip boundary condition on the surface of the flagellum and the requirements of zero net force and torque on the flagellum,
\begin{gather} \boldsymbol{u}_S(s, t) + \boldsymbol{U}(t) + \boldsymbol{\varOmega}(t) \times \boldsymbol{X} = \frac{1}{4{\rm \pi}\eta_m} \int_\varGamma \boldsymbol{\alpha}(\boldsymbol{X} -\boldsymbol{X}^{\prime})\boldsymbol{\cdot} \boldsymbol{f}(s^{\prime}, t) \,\mathrm{d}s^{\prime}, \end{gather}
\begin{gather}\int_\varGamma \boldsymbol{f}(s, t) \,\mathrm{d}s = 0, \end{gather}
\begin{gather}\int_\varGamma \boldsymbol{X} \times \boldsymbol{f}(s, t) \,\mathrm{d}s= 0, \end{gather}
where 
$\boldsymbol {X} \equiv \boldsymbol {X}(s, t)$ and 
$\boldsymbol {X}^{\prime } \equiv \boldsymbol {X}(s^{\prime }, t)$.
Similar to the approach discussed in § 2, we solved the discretized version of these equations for the instantaneous velocities 
$\boldsymbol {\varOmega }(t)$, 
$\boldsymbol {U}(t)$ and the force density 
$\boldsymbol {f}(s, t)$. We normally calculated the angular and swimming velocities for approximately 60–80 snapshots per one period of oscillation and averaged them over one cycle of motion, 
$\langle \varOmega \rangle = ({1}/{T}) \int _0^{T} \varOmega (t) \,\mathrm {d}t \approx ({1}/{N_T}) \sum _{i=1}^{N_T} \varOmega (t_i)$ and 
$\langle \boldsymbol {U} \rangle = ({1}/{T} )\int _0^{T} \boldsymbol {R}(t)\boldsymbol {\cdot } \boldsymbol {U}(t) \,\mathrm {d}t \approx ({1}/{N_T}) \sum _{i=1}^{N_T}\boldsymbol {R}(t_i)\boldsymbol {\cdot } \boldsymbol {U}(t_i)$, where 
$\boldsymbol {R}(t)$ is the rotation operator that transforms the swimming velocity vector to the initial coordinate system 
$\boldsymbol {X}(s, 0)$, which is motionless with respect to the fluid at infinity, and 
$N_T$ is the number of snapshots per one period.
Our computations predict that 
$\langle \varOmega \rangle = 0$, and the flagellum swims in a straight line. In addition to an overall translation along the net swimming direction, 
$\langle \boldsymbol {U} \rangle$, the finite flagellum's motion involves pitching (rotation of the swimmer's centreline with respect to 
$\langle \boldsymbol {U} \rangle$ – note that the flagellum centreline is parallel to the wave centreline (about which actuation occurs) and passes through the centre of mass of the flagellum) and bobbing (motion of the flagellum's centre of mass perpendicular to 
$\langle \boldsymbol {U} \rangle$), see figure 4. During each cycle of motion, the flagellum goes through the same sequence of configurations. Two kinds of flagellum configuration, even and odd, are of particular interest. In the even configuration (illustrated in figure 5a) the flagellum has reflection symmetry with respect to the vertical line that passes through the flagellum's centre of mass. In the odd configuration the shape of the flagellum has point symmetry about the centre of the flagellum (marked by cross-hairs in figure 5b).
Figure 4. Propulsion of a flagellum during one period of oscillation. The time instants are separated by one eighth of the period. Flagella drawn with dashed lines correspond to the time moments 
$t=0$ and 
$t=T$. The thin dotted line is the trajectory of the flagellum's ‘head’ over one cycle of motion. The flagellum length is (a) 
$N_\varLambda =0.5$, (b) 
$N_\varLambda =1$. In both (a,b) the amplitude and the wavelength were set to 
$bq =1$ and 
$\lambda /\ell _S = 1$, respectively.
Figure 5. The swimming trajectory of the head (black) for a flagellum (a) starting its motion in an even configuration, and (b) starting its motion in an odd configuration. The arrow shows the direction of the average swimming velocity, and the dotted line is the centreline of the flagellum in the odd configuration shown. (c) The maximum pitching angle (in degrees) as a function of the flagellum's contour length 
$N_\varLambda$. We set 
$bq =1$.
Koehler, Spoor & Tilley (Reference Koehler, Spoor and Tilley2012), in work on the swimming of finite-length flagella in a Newtonian 3-D fluid, used symmetry arguments to prove that a flagellum in an even conformation has an instantaneous swimming velocity parallel to the flagellum's centreline. Koehler et al. (Reference Koehler, Spoor and Tilley2012) pointed out that the time reversal of the material particles’ velocities should lead to the reversal of the instantaneous swimming velocity, 
$\boldsymbol {U}_{-t} = -\boldsymbol {U}_t$, due to kinematic reversibility of Stokes equations. On the other hand, in an even configuration, the mirror reflection of the material particles’ velocities about the vertical line is equivalent to the time reversal. Therefore, the instantaneous swimming velocity in an even configuration must be identical to the mirror image of the time-reversed velocity. This condition requires the component of the swimming velocity normal to the centreline be equal to zero. Koehler et al. (Reference Koehler, Spoor and Tilley2012) also proved that a flagellum in an odd configuration has zero angular velocity. Thus, if a flagellum starts its motion from an even configuration, as shown in figure 5(a), its initial pitching angle is equal to zero. As the flagellum continues its motion, the pitching angle increases and reaches its maximum value at 
$t=T/4$, when the flagellum reaches an odd configuration. At this moment 
$\varOmega (t)$ = 0, and the pitching angle goes through a ‘turning point’. The grey curve in figure 5(a) shows the flagellum in one of its odd configurations. The magenta arrow shows the net direction of swimming, 
$\langle \boldsymbol {U} \rangle$. The dotted line is the flagellum's centreline in the odd configuration, and 
$\gamma _{max}$ denotes the maximum pitching angle.
If one aligns the 
$x$-axis of the laboratory frame with the flagellum's centreline when the flagellum is not in an even configuration, as in figure 5(b), the flagellum will appear to swim diagonally in such a coordinate system. A similar swimming pattern was reported by Peng, Pak & Elfring (Reference Peng, Pak and Elfring2016) in the work on flagella locomotion in granular media.
The maximum pitching angle depends on the length of the flagellum. In figure 5(c) we show the maximum pitching angles as a function of the flagellum's contour length, 
$N_\varLambda = L/\varLambda$ scaled by the arcwise wavelength 
$\varLambda$ for 
$bq = 1$ and various values of 
$\lambda /\ell _S$.
In figure 6(a) we plot the swimming speed averaged over one period of oscillations, 
$U = |\langle \boldsymbol {U} \rangle |$, as a function of a dimensionless parameter 
$bq$ for the flagellum contour length equal to one arcwise wavelength, 
$N_\varLambda =1$. Two competing mechanisms influence the swimming speed of the flagellum. On the one hand, larger values of 
$bq$ correspond to steeper angles between the flagellum and the direction of wave propagation and, therefore, a stronger propulsion force. On the other hand, for larger 
$bq$ values the segments of the flagellum come closer to each other. The hydrodynamic interactions between the segments tend to slow down the swimmer. The hydrodynamic interactions are stronger in the limiting case of a 2-D hydrodynamics (small 
$\lambda /\ell _S$ ratios) due to a slow, logarithmic spatial decay rate of the flow field. In the opposite limit of large 
$\lambda /\ell _S$ ratios, our calculations do not reproduce Higdon's results for the swimming speed in a purely 3-D fluid (Higdon Reference Higdon1979). Being qualitatively similar to Higdon's prediction, our calculations show much larger swimming speeds for 
$\lambda /\ell _S \gg 1$. As we discussed at the end of § 2, the incompressibility of the membrane sets a constraint on the fluid dynamics that leads to an effective drag anisotropy that grows logarithmically as a function of 
$\lambda /\ell _S$ for 
$\lambda /\ell _S \gg 1$. The enhanced drag anisotropy in a quasi-2-D membrane is responsible for larger swimming speeds in quasi-2-D membranes (in comparison with pure 2-D or 3-D fluids).
Figure 6. Calculated swimming speed scaled by the wave speed of an inextensible, headless finite-length flagellum (a) as a function of parameter 
$bq$ for a flagellum length equal to one arcwise wavelength, 
$N_\varLambda = 1$, (b) as a function of the scaled flagellum length 
$N_\varLambda$ for 
$bq = 1$. In (b) the circles represent calculations based on the Lorentz reciprocal theorem (see (4.9)). The horizontal dotted black lines are asymptotes for the swimming velocities in the limit 
$N_\varLambda \to \infty$ calculated using the method described in § 2.
In figure 6(b) we plot the swimming speed as a function of the scaled flagellum length 
$N_\varLambda$ for 
$bq=1$. For 
$N_\varLambda < 1$ the flagellum performs large yawing motion that is inefficient for swimming (see figure 4a). For larger values of 
$N_\varLambda$ the long-range hydrodynamic interactions taper off the growth of the swimming speed, and the speed approaches the values found for an infinitely long flagellum (shown as horizontal dotted lines in figure 6(b) for 
$\lambda /\ell _S = 100$ and 
$\lambda /\ell _S = 10$). The per cent difference between the swimming speeds of a finite flagellum of length 
$N_\varLambda =8$ and an infinitely long flagellum is 1.4 % for 
$\lambda /\ell _S = 100$ and 2.8 % for 
$\lambda /\ell _S = 10$. For smaller values of 
$\lambda /\ell _S$ the convergence of the swimming speed to that of an infinitely long flagellum is slower, due to greater influence of the long-range hydrodynamic interactions on the flagellum's dynamics. The ‘bumps’ in the curves reflect smaller yawing of the flagellum for the values of 
$N_\varLambda$ that are close (but not exactly equal) to the pitching angle minima in figure 5(c).
To find the flagellum motion that is optimal in terms of the power consumption, we calculated the swimming efficiency. As discussed in Koehler et al. (Reference Koehler, Spoor and Tilley2012), there are multiple efficiency metrics. Here, we calculated the efficiency as the ratio of the power required to pull the flagellum through the fluid at its average swimming speed, 
$F_{ pull} U = (1/ \mu _\| ) U^{2}$, to the average power 
$\langle P \rangle$ consumed by the swimmer over one period of motion,
\begin{equation} \eta = \frac{(1/ \mu_\| ) U^{2}}{\langle P \rangle}. \end{equation}
Here, 
$\mu _\|$ is the translational mobility of an inactive (straightened) flagellum for the motion along the flagellum's axis. The average power consumption is
\begin{equation} {\langle P \rangle} = \frac{1}{T} \int_0^{T} \mathrm{d} t \int_\varGamma \mathrm{d}s \boldsymbol{f}(s, t)\boldsymbol{\cdot} \boldsymbol{u}_S(\boldsymbol{X}, t). \end{equation}We calculated the mobility and the power consumption numerically using the BEM.
Our calculated flagellum efficiency is in qualitative agreement with Higdon (Reference Higdon1979). In figure 7(a) we plot the efficiency as a function of the amplitude 
$bq$ for a few values of the scaled wavelength 
$\lambda /\ell _S$ for a flagellum of length 
$N_\varLambda = 1$. For smaller values of 
$bq$ the segments of the flagellum have small angles with respect to the direction of wave propagation and produce a weak thrust. For larger values of 
$bq$ the flagellum ‘shrinks’ along the 
$x$-axis, and the stronger interference between the segments of the flagellum leads to a decrease in efficiency. The efficiency increases with 
$\lambda /\ell _S$ due to the reduced role of the long-range hydrodynamics on length scales exceeding the Saffman length 
$\ell _S$. The maximum efficiency falls at 
$bq \approx 1.1\text {--}1.2$.
Figure 7. Calculated flagellum efficiency of an inextensible finite-length flagellum (a) as a function of the flagellum amplitude 
$bq$ for 
$N_\varLambda = 1$, (b) as a function of 
$N_\varLambda$ for 
$bq = 1$.
In figure 7(b) we plot the efficiency as a function of the flagellum length 
$N_\varLambda$ for a for a few values of 
$\lambda /\ell _S$ for 
$bq=1$. For small values of 
$N_\varLambda$ the swimming of the flagellum is inefficient due to an excessive yawing motion and a weak overall thrust (see figure 4a). The efficiency reaches a maximum at 
$N_\varLambda \approx 1.3\text {--}1.4$, close to the pitching angle minima in figure 5(c). We note that this result is in accord with calculations of the efficiency of undulatory headless flagella in three dimensions (Higdon Reference Higdon1979; Dresdner, Katz & Berger Reference Dresdner, Katz and Berger1980; Koehler et al. Reference Koehler, Spoor and Tilley2012). The efficiency decreases with further growth of 
$N_\varLambda$ due to interference between the crests of the flagellum. The interference is stronger for larger amplitudes 
$bq$ since the crests are closer to each other, and the efficiency drops off more abruptly from its optimal value for larger values of 
$bq$. The secondary maxima correspond to a smaller yawing and larger propulsion for various values of 
$N_\varLambda$. Figure 4 demonstrates that the flagellum travels a considerable distance along the 
$y$-axis while making moderate overall progress along the 
$x$-axis.
In the following section we discuss an alternative computational approach to finding the swimming velocities using the Lorentz reciprocal theorem.
4. Lorentz reciprocal theorem for a quasi-2-D membrane
Finding the analytical solution for the swimming velocity can be a daunting task. One of the major difficulties is that one needs to solve Stokes equations with time dependent no-slip boundary conditions on the surface of the swimmer. Stone & Samuel (Reference Stone and Samuel1996) offered an elegant way to find the swimming velocity using the Lorentz reciprocal theorem (LRT) (Happel & Brenner Reference Happel and Brenner1965). For a fluid in three dimensions the LRT states that, if there are two solutions to Stokes equations and the incompressibility condition with the velocity fields and the stress tensors 
$(\boldsymbol {v}, \boldsymbol {\sigma })$ and 
$(\boldsymbol {v}^{\prime }, \boldsymbol {\sigma }^{\prime })$, respectively, that satisfy the same boundary conditions at infinity, then for a volume of fluid 
$V$ bounded by surface 
$S$, we have
\begin{equation} \oint_S \boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\sigma}^{\prime} \boldsymbol{\cdot} \boldsymbol{n}\, \mathrm{d} S= \oint_S \boldsymbol{v}^{\prime} \boldsymbol{\cdot} \boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S, \end{equation}
where 
$\boldsymbol {n}$ is the outward normal to the surface 
$S$.
Here we formulate the LRT approach for a finite swimmer confined to a quasi-2-D membrane. Let 
$\boldsymbol {v}$ and 
$\boldsymbol {\sigma }$ be the membrane velocity and the stress fields for the swimming problem. These velocity and stress fields are solutions of the Stokes equations and the incompressibility condition, (1.4a,b). They also satisfy the conditions of zero net force and torque on the swimming body, (1.7) and (1.8). For the reciprocal solution of (1.4a,b) we choose the membrane velocity 
${\boldsymbol {v}}^{\prime }$ and stress 
${\boldsymbol {\sigma }}^{\prime }$ fields due to an inactive object of the same shape as the swimmer and being dragged as a solid body with constant translational velocity 
${\boldsymbol {U}}^{\prime }$.
When the condition 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\sigma } = 0$ is relaxed (see (1.6)), a more general form of the LRT is (Kim & Karrila Reference Kim and Karrila1991)
\begin{equation} \oint_S {\boldsymbol{v}}^{\prime} \boldsymbol{\cdot} (\boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{n}) \,\mathrm{d} S - \int_V {\boldsymbol{v}}^{\prime} \boldsymbol{\cdot}(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma}) \,\mathrm{d} V = \oint_S \boldsymbol{v} \boldsymbol{\cdot} ({\boldsymbol{\sigma}}^{\prime} \boldsymbol{\cdot} \boldsymbol{n}) \,\mathrm{d} S - \int_V \boldsymbol{v} \boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\sigma}}^{\prime}) \,\mathrm{d} V, \end{equation}
where 
$V$ is the swimmer's volume bounded by the surface 
$S$. Here, we treat the volume occupied by the swimmer as being equivalent to the fluid domain of the same shape as the swimmer's and having the same velocity distribution as that of the material particles of the swimmer.
Let us consider the first term on the left-hand side of (4.2)
\begin{equation} \oint_S {\boldsymbol{v}}^{\prime} \boldsymbol{\cdot} (\boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{n}) \,\mathrm{d} S = {\boldsymbol{U}}^{\prime} \boldsymbol{\cdot} \int_{S_{w}} \mathrm{d} \boldsymbol{F}, \end{equation}
where we took into account that 
${\boldsymbol {v}}^{\prime }={\boldsymbol {U}}^{\prime }$ is a constant vector at the surface of the domain (uniform translation). Also, since 
$\boldsymbol {\sigma }$ does not have 
$z$-components, only 
$\boldsymbol {\sigma }\boldsymbol {\cdot } \boldsymbol {n}\, \mathrm {d} S = \mathrm {d} \boldsymbol {F}$ on the curvy wall of the domain, 
$S_w$ (see figure 8), will make a non-zero contribution.
Figure 8. A region in a quasi-2-D membrane of the same geometry as a swimmer. The integration surface 
$S$ in (4.5) is comprised of the curvy wall 
$S_{w}$ embedded in the membrane, the top (
$S_t$) and the bottom (
$S_b$) surfaces of the membrane that are in contact with the bulk fluid.
Taking into account (1.6), the second term on the left-hand side of (4.2) can be rearranged as
\begin{equation} -\int_V {\boldsymbol{v}}^{\prime} \boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} )\, \mathrm{d} V = {\boldsymbol{U}}^{\prime} \boldsymbol{\cdot} \int_V \frac{2\boldsymbol{f}}{h} \,\mathrm{d} V = {\boldsymbol{U}}^{\prime} \boldsymbol{\cdot} \int_{S_{t, b}} 2\boldsymbol{f} \, \mathrm{d} S, \end{equation}
where we take into account that the traction forces 
$2\boldsymbol {f}$ due to the fluid flows in the surrounding fluid act on the flat sides of the domain 
$S_{t,b}$ (see (1.6) and figure 8), and 
$\textrm {d}V = h \,\textrm {d}S$. Therefore, the left-hand side of (4.2) becomes
\begin{equation} {\boldsymbol{U}}^{\prime} \boldsymbol{\cdot} \left(\int_{S_{w}} \mathrm{d} \boldsymbol{F} + \int_{S_{t,b}} 2\boldsymbol{f} \, \mathrm{d} S \right)= {\boldsymbol{U}}^{\prime} \boldsymbol{\cdot} {\boldsymbol F} = 0, \end{equation}
where 
${\boldsymbol {F}}$ is the net force on the swimmer and is equal to zero.
Similarly, for the terms on the right-hand side of (4.2) we have
\begin{align} \oint_S \boldsymbol{v} \boldsymbol{\cdot} ({\boldsymbol{\sigma}}^{\prime} \boldsymbol{\cdot} \boldsymbol{n}) \,\mathrm{d} S - \int_V \boldsymbol{v} \boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\sigma}}^{\prime} ) \,\mathrm{d} V &= \int_{S_{w}}\boldsymbol{v}\boldsymbol{\cdot} ({\boldsymbol{\sigma}}^{\prime} \boldsymbol{\cdot} \boldsymbol{n}) \,\mathrm{d} S + \int_{S_{t, b}} \boldsymbol{v} \boldsymbol{\cdot} (2 {\boldsymbol{f}}^{\prime}) \, \mathrm{d} S \nonumber\\ &= \oint_S \boldsymbol{v} \boldsymbol{\cdot}\, \mathrm{d} {\boldsymbol F}^{\prime}, \end{align}
where in the last line of (4.6) we merged two terms into one integral over the total surface of the swimmer, and 
$\mathrm {d} {\boldsymbol {F}}^{\prime }$ denotes an elementary traction force on the inactive ‘swimmer’ being dragged with constant velocity. Thus, the Lorentz reciprocal relation, (4.2), assumes a compact form,
\begin{equation} 0 = \oint_S \boldsymbol{v} \boldsymbol{\cdot}\, \mathrm{d} {\boldsymbol{F}}^{\prime}. \end{equation}
Decomposing the surface velocity of the swimmer into the translational 
$\boldsymbol {U}(t)$ and the surface disturbance 
$\boldsymbol {u}_S(t)$ velocities, 
$\boldsymbol {v}(t) = \boldsymbol {U}(t)+ \boldsymbol {u}_S(t)$, we rewrite the Lorentz reciprocal relation in the form
\begin{equation} {\boldsymbol{F}}^{\prime}(t) \boldsymbol{\cdot} \boldsymbol{U}(t) = - \oint_{S(t)} \boldsymbol{u}_S \boldsymbol{\cdot}\, \mathrm{d} {\boldsymbol F}^{\prime}, \end{equation}similar to the equation derived by Stone & Samuel (Reference Stone and Samuel1996) for a swimmer in a 3-D fluid. In (4.8) the integration is performed over the instantaneous surface area of the swimmer. The generalization of (4.8) for the motion that involves rotation is
\begin{equation} {\boldsymbol{F}}^{\prime}(t) \boldsymbol{\cdot} \boldsymbol{U}(t) + {\boldsymbol{L}}^{\prime} (t) \boldsymbol{\cdot} \boldsymbol{\varOmega}(t)= - \oint_{S(t)} \boldsymbol{u}_S \boldsymbol{\cdot}\, \mathrm{d} {\boldsymbol{F}}^{\prime}, \end{equation}
where 
${\boldsymbol {L}}^{\prime }(t)$ is the torque applied to the inactive inclusion and 
$\boldsymbol {\varOmega }(t)$ is the swimmer's angular velocity.
The Lorentz reciprocal relation, (4.9), is particularly useful for computation of the swimmer's translational and rotational velocities, 
$\boldsymbol {U}(t)$ and 
$\boldsymbol {\varOmega }(t)$, when the stress tensor of the reciprocal problem (motion of an inactive body) is known. Unfortunately, it is also difficult to solve the reciprocal problem analytically for an inclusion of an arbitrary shape in a quasi-2-D membrane, since the coupling with the bulk fluid makes the problem essentially three-dimensional. When the reciprocal solution is not available, (4.9) can serve as an alternative computational path for finding the swimming velocity.
As a test of (4.9), we found the swimming velocities of a finite inextensible flagellum as described in § 3. We solved the reciprocal problem numerically by finding the force densities 
$\boldsymbol {f}(\boldsymbol {x}, t)$ for the uniform rotation of the flagellum about the 
$z$-axis and for the translational motion of the flagellum along the 
$x$- and 
$y$- axes for multiple flagellum conformations corresponding to various time instants of the swimming cycle. We then solved the resulting system of three equations (4.9) for the instantaneous angular velocity and 
$x$- and 
$y$-components of the translational velocity and found the average swimmer's speed over one period of oscillations. In figure 6(b) circles superimposed on the curves show the swimming velocities obtained using the Lorentz equation (4.9).
4.1. The 2-D squirmer example
The LRT can significantly simplify computations in the case of tangential deformations of the swimmer's body, when the overall shape of the creature remains unchanged. In this case the reciprocal problem can be solved for a single time instant, and the instantaneous swimming velocity can then be found from (4.9) by plugging in the time-dependent surface disturbance velocity 
$\boldsymbol {u}_S(t)$.
As an example, we consider a 2-D version of a squirmer, a creature that propels itself by beating its multiple hair-like appendages (cilia) in a periodic fashion. The periodic motion of the cilia carpet can be modelled by prescribing a velocity field on the surface of the squirmer. In a minimal model of a 2-D squirmer, a disk-like body of radius 
$a$ is propelled due to a tangential disturbance velocity 
$u_S(\phi ) = B_1\sin \phi + B_2 \sin (2 \phi )$ on the disk circumference (see figure 9), where 
$B_1$ and 
$B_2$ are constants. Here, 
$B_1$ corresponds to the strength of a source dipole flow field, and 
$B_2$ determines the strengths of the force dipole and source quadrupole fields (Blake Reference Blake1971; Papavassiliou & Alexander Reference Papavassiliou and Alexander2015). The material points in the interior of the disk are motionless in the frame of the squirmer, and they represent a passive body such as in a disk-shaped creature that only has cilia on its circular edge. When the constants 
$B_1$ and 
$B_2$ have the same sign, they describe a contractile swimmer. Otherwise, the model corresponds to an extensile swimmer.
Figure 9. The flow field due to a squirmer is modelled as a superposition of the flow fields due to point-like forces (blobs). The blobs on the circumference (red) move with tangential velocity 
$u_S(\phi ) = B_1\sin \phi + B_2 \sin (2 \phi )$ in the body frame of the squirmer. The blobs in the interior of the squirmer (blue) are motionless in the body frame.
In the limit of 
$a/\ell _S \ll 1$ the reciprocal problem for a disk was solved by Saffman (Reference Saffman1976) in the studies related to the Stokes paradox and a particle mobility in a quasi-2-D fluid. In appendix B we outline calculations for the swimming velocity of the squirmer in the limit 
$a/\ell _S \ll 1$ using Saffman's solution for the reciprocal stress tensor 
$\boldsymbol {\sigma }^{\prime }$. The LRT reproduces the known swimming speed 
$U = B_1/2$, for a 2-D squirmer in the limiting case of a pure 2-D hydrodynamics.
Since the analytical solution for the reciprocal problem for a disk of arbitrary radius 
$a/\ell _S$ is not readily available, we found the reciprocal stress 
$\boldsymbol {\sigma }^{\prime }$ numerically by adopting the method of regularized Stokeslets (RS) for a quasi-2-D membrane developed by Camley & Brown (Reference Camley and Brown2013) in their work on mobility of inclusions in a quasi-2-D membrane. In Camley & Brown (Reference Camley and Brown2013) the flow field due to a moving inclusion is modelled as a superposition of the flow fields due to point-like forces (blobs) tiling the disk area (see figure 9),
\begin{equation} u_\alpha (\boldsymbol{x}) = \sum_{i=1}^{N} \alpha_{\alpha\beta}(\boldsymbol{x} - \boldsymbol{x}_i) f^{\prime}_\beta(\boldsymbol{x}_i), \end{equation}
where 
$\alpha , \beta = x, y$; 
$N$ is the total number of blobs; 
$\boldsymbol {x}_i$ is the in-plane coordinate of the 
$i$th blob; 
$\alpha _{\alpha \beta }$ is the Levine–MacKintosh response function (see (1.2)); and 
$f^{\prime }_\beta (\boldsymbol {x}_i)$ is the unknown force distribution.
In the reciprocal problem the disk is being pulled through the membrane as a solid object with some given velocity 
$\boldsymbol {U}^{\prime }$. The force distribution 
$f^{\prime }_\beta (\boldsymbol {x}_i)$ over the blobs is found by imposing a no-slip boundary condition on each blob,
\begin{equation} U^{\prime}_\alpha= \sum_{j=1}^{N} \alpha_{\alpha\beta}(\boldsymbol{x}_i - \boldsymbol{x}_j) f^{\prime}_\beta(\boldsymbol{x}_j). \end{equation}
Due to the squirmer's reflection symmetry about the 
$x$-axis, the rotational motion of the squirmer in an unbounded domain is ruled out, and the swimming velocity 
$\boldsymbol {U}$ can be found from the discretized version of (4.8),
\begin{equation} \left( \sum_{i=1}^{N} \boldsymbol{f}^{\prime} (\boldsymbol{x}_i) \right) \boldsymbol{\cdot} \boldsymbol{U} = - \sum_{\textrm{{blobs}}\ \textrm{{on}}\ \textrm{{rim}}} \boldsymbol{u}_S(\boldsymbol{x}_j) \boldsymbol{\cdot} \boldsymbol{f}^{\prime} (\boldsymbol{x}_j). \end{equation}
The summation on the right-hand side of (4.12) is carried out only over the blobs on the squirmer's circumference since the inner blobs are motionless in the creature's frame. In the RS method the logarithmic singularity of the membrane response functions for 
$\kappa \to 0$ is eliminated by the regularization (smoothing) process that involves integration of the response function over the blob envelope function centred at 
$\kappa =0$. Camley & Brown (Reference Camley and Brown2013) selected a Gaussian function for the regularization. The width of the Gaussian is controlled by an auxiliary parameter 
$\varepsilon$. Camley and Brown set 
$\varepsilon = \delta /2$, where 
$\delta$ is the distance between the centres of adjacent blobs, calculated the inclusion mobilities for several values of 
$\delta$ in the range (0.03–0.07)a and extrapolated the results to the limit 
$\delta \to 0$. While the numerical calculations for the inclusion mobilities are only weakly dependent on the choice of 
$\varepsilon$, in our case of a swimming squirmer, the swimming velocities are more sensitive to the choice of the regularization parameter 
$\varepsilon$, since it effectively determines the thickness of the squirmer's deforming outer ring and its permeability, and therefore becomes a physical parameter.
In figure 10(a) we plot the LRT results for the scaled swimming speed of a squirmer, 
$U/(B_1/2)$, as a function of the squirmer radius for several values of parameter 
$\varepsilon$. As in Camley & Brown (Reference Camley and Brown2013), for a selected dependence of 
$\varepsilon$ on 
$\delta$ (e.g. 
$\varepsilon = \delta /6$), we calculated the swimming speed for a range of 
$\delta$ values and extrapolated it to 
$\delta \to 0$. As can be seen in figure 10(a), the regularization parameter 
$\varepsilon =\delta /12$ gives the swimming velocity that is close to the known value of 
$B_1/2$ in the limiting case of a pure 2-D membrane (membrane in vacuum). Our calculations also show that the scaled swimming velocity 
$U/(B_1/2)$ is independent of the ratio 
$B_2/B_1$.
Figure 10. LRT results for the squirmer swimming speed scaled by 
$B_1/2$ as a function of scaled squirmer radius 
$a/\ell _S$ (a) for various values of the regularization parameter 
$\varepsilon$, where 
$\delta$ is the distance between the centres of neighbouring blobs, (b) for 
$\varepsilon = \delta /12$ (circles). The black curve corresponds to the BEM calculations for 
$\varepsilon = \delta /12$. The swimming speed is independent of the ratio 
$B_2/B_1$.
In figure 10(b) we compare the results of calculations for the swimming speed obtained within the LRT and the BEM for 
$\varepsilon =\delta /12$. For the direct BEM calculation of the swimming speed and the force distribution we solved simultaneously the equations that impose no-slip boundary conditions and a zero net force on the swimmer,
\begin{gather} U_ \alpha = \sum_{\textrm{interior}\ \textrm{blobs}} \alpha_{\alpha\beta}(\boldsymbol{x}_j - \boldsymbol{x}_i) f_\beta(\boldsymbol{x}_i), \end{gather}
\begin{gather}U_\alpha + u_{S\alpha}(\boldsymbol{x}_j) = \sum_{\textrm{blobs}\ \textrm{{on}}\ \textrm{rim}} \alpha_{\alpha\beta}(\boldsymbol{x}_j - \boldsymbol{x}_i) f_\beta(\boldsymbol{x}_i), \end{gather}
\begin{gather}0=\sum_{i=1}^{N} \boldsymbol{f}(\boldsymbol{x}_i). \end{gather}
As can be observed in figure 10(b), the squirmer swimming velocity decreases with an increase of 
$a/\ell _S$ ratio. Larger values of 
$a/\ell _S$ correspond to a larger viscosity of the fluid embedding the membrane, which leads to a stronger traction on the ‘back’ and ‘belly’ of the creature.
In figure 11 we show our results for the squirmer efficiency as a function of the scaled radius 
$R/\ell _S$ for various 
$B_2/B_1$ ratios. The efficiency is defined as the ratio of the power required to drag an inactive disk through the membrane with the speed equal to the swimming speed of the squirmer, 
$(1/\mu ) U^{2}$, to the power expended by the squirmer 
$P$,
\begin{equation} \eta = \frac{(1/\mu) U^{2}}{P}. \end{equation}
Here, 
$\mu$ is the translational mobility of a disk. The analytical expression of the disk mobility for arbitrary disk radii 
$R/\ell _S$ was found by Hughes, Pailthorpe & White (Reference Hughes, Pailthorpe and White1981) (HPW). In our calculations, we used an accurate approximation of the complicated HPW mobility expression developed by Petrov & Schwille (Reference Petrov and Schwille2008). The power consumed by the squirmer is determined by
\begin{equation} P = \sum_{\textrm{{blobs}}\ \textrm{{on}}\ \textrm{{rim}}} \boldsymbol{u}_S(\boldsymbol{x}_i)\boldsymbol{\cdot} \boldsymbol{f}(\boldsymbol{x}_i). \end{equation}
Figure 11. Calculated squirmer efficiency as a function of squirmer radius.
The ‘neutral’ squirmer, with 
$B_2 = 0$, achieves the maximum efficiency. Similar to its 3-D and 2-D counterparts, a quasi-2-D squirmer uses only the 
$B_1$ term in its surface disturbance velocity for propulsion, while the 
$B_2$ term determines whether the swimmer is a ‘pusher’ (
$B_2<0$) or ‘puller’ (
$B_2>0$), and leads to additional energy costs. The swimming efficiency is independent of the sign of the 
$B_2/B_1$ ratio, which is consistent with Blake's results for a 2-D squirmer (Blake Reference Blake1971).
The peaks in the efficiency curves indicate two competing contributions to the squirmer's energy budget. Our model squirmer employs only the viscous drag generated in the membrane for propulsion. The traction due to the bulk fluid is detrimental to the squirmer's locomotion. In the limit 
$R/\ell _S \ll 1$, the viscous drag due to the bulk fluid is negligible. However, the long-range character of the flow field generated by the squirmer leads to a smaller thrust force, in comparison with the 3-D squirmer (Lighthill Reference Lighthill1952; Blake Reference Blake1971), leading to less efficiency than the case of 
$R/\ell _S\approx 1$. In the opposite limit of 
$R/\ell _S \gg 1$, the faster spatial decay rate of the flow field leads to a stronger thrust. However, the traction due to the bulk fluid is significant in this regime. From the dimensional analysis, we expect the thrust force generated by the squirmer in the membrane be proportional to 
$\eta _m U$, and the drag force due to the bulk fluid proportional to 
$\eta U R$. Thus, the ratio of the ‘harmful’ drag to the ‘useful’ drag is given by 
$R \eta /\eta _m \sim R/\ell _S$, again leading to a lower efficiency compared to the case of 
$R/\ell _S\approx 1$.
5. Conclusion
We have studied analytically and computationally the locomotion of microscopic organisms confined to a plane of a thin fluid membrane embedded in a bulk fluid of different viscosity. In our model the membrane is sufficiently thin, with material particles moving only in the plane of the membrane (motion in the perpendicular direction is forbidden).
The presence of the bulk fluid allows the introduction of a hydrodynamic length scale, the Saffman length, that controls the energy exchange between the membrane and the surrounding fluid. By varying the Saffman length, we make our model continuously vary between a pure 2-D system (large Saffman length) and a quasi-2-D system (small Saffman length). The hydrodynamic flows in the quasi-2-D membrane have features of both 3-D and 2-D hydrodynamics. We show that a flagellated swimmer in a viscous film (Saffman length smaller than swimmer characteristic length scale) swims faster than the same swimmer in a 3-D fluid. The speed up comes from the effectively larger perpendicular drag coefficient, which arises from the incompressibility of the membrane. It would be interesting to test this general prediction by measuring the swimming speed in both thin films and unconfined fluids of undulatory swimmers such as mammalian sperm or Caenorhabditis elegans. On the other hand, a circular squirmer, whose propulsion mechanism does not employ the local drag anisotropy, slows down for smaller Saffman lengths (in comparison with the squirmer's radius). Paramecium cells have recently been shown to exhibit sinusoidal trajectories in thick fluid films and frequent turns in thin fluid films (Jana et al. Reference Jana, Eddins, Spoon and Jung2015). In these films there is one solid bounding surface; our results could be used to help study how swimming behaviour depends on film thickness in films with two liquid interfaces.
The coupling of the membrane with the bulk fluid makes the problem three-dimensional and quite difficult for analytical treatment. We developed numerical schemes based on the boundary element method and the LRT. We show how the Lorentz reciprocal theorem can be used to simplify the computation of swimming speed, especially for swimmers such as the squirmer that do not change shape during a stroke. While we considered the minimal models of a flagellated swimmer and of a squirmer, our approach can be generalized to other swimmers’ geometries and swimming strokes, including catalytic particles such as Janus particles (Paxton et al. Reference Paxton, Kistler, Olmeda, Sen, St. Angelo, Cao, Mallouk, Lammert and Crespi2004; Michelin & Lauga Reference Michelin and Lauga2017).
Funding
This work has been supported by a Cottrell College Science Award (T.K.) and National Science Foundation grant no. 1437195 (T.R.P.). C.A. gratefully acknowledges support from a Frost Undergraduate Student Research Award.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Infinitely long flagellum in 2-D fluid
In the 2-D limit (membrane in vacuum) the system of (2.17) becomes
\begin{equation} \left. \begin{gathered} \boldsymbol{u}_S(\boldsymbol{x}_i)+ \boldsymbol{U} = \frac{1}{4{\rm \pi}\eta_m}\sum_{j = 1}^{N} \left( \int_{S_j} \boldsymbol{G}^{p}(\boldsymbol{x}_i - \boldsymbol{x}^{\prime}) \,\mathrm{d} \boldsymbol{x}^{\prime} \right) \boldsymbol{\cdot} \boldsymbol{f}(\boldsymbol{x}_{j}), \\ \sum_{j = 1}^{N} \boldsymbol{f}(\boldsymbol{x}_{j}) = 0, \end{gathered}\right\} \end{equation}
where 
$\boldsymbol {G}^{p}(\boldsymbol {x})$ is a 2-D periodic Stokeslet,
\begin{equation} \boldsymbol{G}^{p}(\boldsymbol{x} -\boldsymbol{x}^{\prime}) = \sum_{m = -\infty}^{\infty} -{ \boldsymbol{I}} \ln(q r_m) + \frac{\bar{\boldsymbol{x}}_{m} \bar{\boldsymbol{x}}_{m} }{r_m^{2}}, \end{equation}
with 
$\bar {\boldsymbol {x}} \equiv \{ \bar {x}, \bar {y}\} = \boldsymbol {x} - \boldsymbol {x}^{\prime }$, 
$\bar {\boldsymbol {x}}_{m} = \{\bar {x} + m ({2{\rm \pi} }/{q}), \bar {y} \}$ and 
$r_m=|\bar {\boldsymbol {x}}_{m}|$. In (A1) and (A2) all variables are dimensional. The periodic Stokeslets can be expressed in closed form (Pozrikidis Reference Pozrikidis1987; Sauzade et al. Reference Sauzade, Elfring and Lauga2011) using the analytic formula for the summation,
\begin{equation} A = \sum_{m=-\infty}^{\infty} \ln(|q r_m|) = \tfrac{1}{2} \ln [ 2 \cosh(q \bar{y}) - 2 \cos(q \bar x)]. \end{equation}The components of the periodic Stokeslet can be found as
\begin{gather} G_{xx}^{p} = - A - A_y + 1, \end{gather}
\begin{gather}G_{xy}^{p} = (q \bar{y}) A_x, \end{gather}
\begin{gather}G_{yy}^{p} = - A + (q \bar{y}) A_y , \end{gather}
where 
$A_x$, 
$A_y$ indicate the derivatives of 
$A$ with respect to 
$q\bar {x}$ and 
$q\bar {y}$ respectively. Similar to our treatment of the diagonal terms in (2.19) we eliminate the logarithmic singularity by analytic integration,
\begin{equation} \int_{S_i} \boldsymbol{G}^{p}(\boldsymbol{x}_i - \boldsymbol{x}^{\prime}) \,\mathrm{d} \boldsymbol{x}^{\prime} = { \boldsymbol{I}} \, 2 \lim_{\varepsilon \to 0} \int_{\varepsilon}^{\Delta s /2} (1 -\log(q z)) \,\mathrm{d} z = {\boldsymbol{I}} \Delta s ( 1 - \log(q\Delta s/2). \end{equation}Appendix B. Swimming velocity of a squirmer in the 2-D limit
We consider a tangential squirmer with a prescribed surface velocity of the form
\begin{equation} \boldsymbol{u}_S(\phi) = u_S(\phi) \hat{\boldsymbol{\phi}}= (B_1\sin\phi + B_2 \sin (2 \phi))\hat{\boldsymbol{\phi}}, \end{equation}
with free parameters 
$B_1$ and 
$B_2$. Since the disturbance velocity has only a 
$\hat {\boldsymbol {\phi }}$-component, the right-hand side of (4.9) becomes
\begin{align} \oint_{S(t)} {\boldsymbol{u}}_S \boldsymbol{\cdot}\, \mathrm{d} {\boldsymbol{F}}^{\prime} &= \int_{S_w} u_S(\phi) \hat{\boldsymbol{\phi}} \boldsymbol{\cdot} ({\sigma^{\prime}_{r\phi}} \hat{\boldsymbol{\phi}})\,\mathrm{d} S_w \nonumber\\ &= 2 h a \int_0^{\rm \pi} u_S(\phi) \sigma^{\prime}_{r\phi} \,\mathrm{d} \phi, \end{align}
where we took into account 
$\mathrm {d} S_w = h a \,\mathrm {d} \phi$, where 
$h$ is the thickness of the membrane. Therefore, (4.9) becomes
\begin{equation} {\boldsymbol{F}}^{\prime}(t) \boldsymbol{\cdot} \boldsymbol{U}(t) = - 2 h a \int_0^{\rm \pi} u_S(\phi) \sigma^{\prime}_{r\phi} \,\mathrm{d} \phi. \end{equation}
The membrane stress tensor element 
$\sigma _{r\phi }^{\prime }(r,\phi )$ is determined as
\begin{equation} \sigma_{r\phi}^{\prime}(r,\phi) = - \frac{\eta_m}{h}\left[ \frac{1}{r} \frac{\partial {u}^{\prime}_r}{\partial \phi} + \frac{\partial {u}^{\prime}_\phi }{\partial r} - \frac{{u}^{\prime}_\phi}{r} \right]. \end{equation}
For a special case of 
$a\ll \ell _S$ Saffman (Reference Saffman1976) found
\begin{gather} {\boldsymbol{F}}^{\prime} = \frac{4 {\rm \pi}\eta_m {\boldsymbol U}^{\prime}}{ \log(2 \ell_S/a) - \gamma }, \end{gather}
\begin{gather}\sigma_{r\phi}^{\prime}(r,\phi) = \frac{\eta_m}{h} \frac{4\alpha \sin\phi}{r^{3}}, \end{gather}with
\begin{equation} \alpha = \frac{a^{2} {U}^{\prime} }{2 ( \gamma - \log(2 \ell_S/a))}. \end{equation}
Here, 
$\gamma = 0.577$ is the Euler constant.
Plugging (B5), (B6) and (B7) in (4.9), after some simplifications we arrive at the squirmer swimming velocity in the limit of 
$a/\ell _S \ll 1$
\begin{equation} U = \int_0^{\rm \pi} u_S(\phi) \sin\phi \,\mathrm{d} \phi= \frac{B_1}{2}. \end{equation}
