# The hydrodynamics of slender swimmers near deformable interfaces

Sankalp Nambiar  Nordita, KTH Royal Institute of Technology and Stockholm
University, Stockholm 10691, Sweden  J. S. Wettlaufer<EMAIL_ADDRESS>Nordita, KTH Royal Institute of Technology and Stockholm University, Stockholm
10691, Sweden  Yale University, New Haven, CT 06520-8109, USA

###### Abstract

We study the coupled hydrodynamics between a motile slender microswimmer and a
deformable interface that separates two Newtonian fluid regions. From the
disturbance field generated by the swimming motion, we quantitatively
characterize the interface deformation and the manner in which the coupling
modifies the microswimmer translation itself. We treat the role of the swimmer
type (pushers and pullers), size and model an interface that can deform due to
both surface tension and bending elasticity. Our analysis reveals a strong
dependence of the hydrodynamics on the swimmer orientation and position. Given
the viscosities of the two fluid media, the interface properties and the
swimmer type, a swimmer can either migrate towards or away from the interface
depending on its configurations. When the swimmer is oriented parallel to the
interface, a pusher-type swimmer is repelled from the interface at short times
if it is swimming in the more viscous fluid. At long times however, pushers
are always attracted to the interface, and pullers are always repelled from
it. On the other hand, swimmers oriented orthogonal to the interface exhibit a
migration pattern opposite to the parallel swimmers. In consequence, a host of
complex migration trajectories emerge for swimmers arbitrarily oriented to the
interface. We find that confining a swimmer between a rigid boundary and a
deformable interface results in regimes of attraction towards both surfaces
depending on the swimmer location in the channel, irrespective viscosity
ratio. The differing migration patterns are most prominent in a region of
order the swimmer size from the interface, where the slender swimmer model
yields a better approximation to the coupled hydrodynamics.

##  I Introduction

Swimming microorganisms are often found in fluid environments near interfaces
that can be either rigid Rothschild (1963); Frymier et al. (1995); DiLuzio et
al. (2005); Berke et al. (2008); Di Leonardo et al. (2011); Ferracci et al.
(2013); Molaei et al. (2014) or compliant Montecucco and Rappuoli (2001);
Moriarty et al. (2008); Lee et al. (2008); Sekirov et al. (2010). The
hydrodynamics in such systems involves a coupling between the intrinsic
swimming motion of the microswimmer and the boundary Trouilloud et al. (2008);
Spagnolie and Lauga (2012); Dias and Powers (2013); Ledesma-Aguilar and
Yeomans (2013); Bechinger et al. (2016); Shaik and Ardekani (2017); Ishikawa
(2019); Daddi-Moussa-Ider et al. (2019); Škultéty and Morozov (2020); Dalal et
al. (2020), due to which a host of rich and complex dynamical responses emerge
Sokolov et al. (2009); Kurtuldu et al. (2011); Woodhouse and Goldstein (2012);
Lushi et al. (2014); Costanzo et al. (2014); Wioland et al. (2016); Mathijssen
et al. (2016); Bechinger et al. (2016); Theillard et al. (2017). For instance,
observations near rigid interfaces include aspects of confinement induced
microswimmer migration such as upstream swimming and boundary accumulation
Berke et al. (2008); Li and Tang (2009); Rusconi et al. (2014); Ezhilan et al.
(2015); Bearon and Hazel (2015); Barry et al. (2015); Ezhilan and Saintillan
(2015); Mathijssen et al. (2016); Wioland et al. (2016); Manabe et al. (2020);
Vennamneni et al. (2020), changes in the confinement pressure and fluctuation
forces Yan and Brady (2015); Ezhilan et al. (2015); Lee et al. (2017);
Škultéty and Morozov (2020), as well as boundary-induced changes in the
swimmer trajectory Molaei et al. (2014); Hu et al. (2015); Wioland et al.
(2016); Manabe et al. (2020). When swimming near complaint interfaces however,
hydrodynamic effects drive interfacial deformation. The resulting coupled
dynamics depends on surface tension and/or elasticity lea (2007); Bickel
(2007); Pozrikidis (2007); Daddi-Moussa-Ider et al. (2016); Rallabandi et al.
(2018). Depending on the properties of the deformable interface, the fluid
medium and the swimmer model adopted, simulations and continuum theories have
shown enhanced pumping of the disturbance flow field Dias and Powers (2013),
as well as deformation-induced enhancements Ledesma-Aguilar and Yeomans
(2013); Daddi-Moussa-Ider et al. (2019) and retardation J. and Amy (1995);
Dalal et al. (2020) in the swimmer translation.

Analytical studies that have investigated the hydrodynamics of microswimmers
and passive particles near deformable interfaces, have treated the bodies as
singularities (point forces, force-dipoles, quadropoles, etc.) Bickel (2007);
Trouilloud et al. (2008); Daddi-Moussa-Ider et al. (2016); Rallabandi et al.
(2018); Daddi-Moussa-Ider et al. (2019), two-dimensional swimming sheets Dias
and Powers (2013), or spherical particles/squirmers Berdan and Leal (1982);
Shaik and Ardekani (2017). Each of these unique approaches can explore
observations in the regime of their applicability. For example, Dias and
Powers (2013) showed that a two-dimensional sheet swimming near a deformable
interface can generate a pumping velocity that can be either attractive or
repulsive to the interface depending on the viscosity ratio of the two fluid
regions. Moreover, they found that an increase in the bending stiffness
resulted in enhanced swimmer translation but reduced the fluid pumping. Daddi-
Moussa-Ider et al. (2019) considered higher order terms of a multipole
analysis for determining the effect of the coupled hydrodynamics on the
translational and rotational motion of a microswimmer near a deformable
interface. By considering an interface which exhibits resistance to both
bending and shear, they established that the swimmer translation was enhanced
due to interface bending and suppressed due to shear resistance. On the other
hand Shaik and Ardekani (2017) modeled the swimmer as a sphere with a
prescribed surface slip velocity (the spherical squirmer approach Blake
(1971)) and quantitatively characterized both the interface deformation and
the swimmer translation for different squirmer configurations. They found
that, depending on the initial orientation of the squirmer relative to the
interface, the squirmer either swam towards or away from the interface.

Analytical progress in such approaches has been contingent upon a small
parameter, such a small capillary number $Ca$ (the ratio of viscous fluid
stress to surface tension) Trouilloud et al. (2008); Shaik and Ardekani
(2017), or small sheet deformations in a two-dimensional domain Dias and
Powers (2013). Alternatively, the small parameter emerges from the separation
of scales of the swimmer size relative to its distance from the interface;
when it is sufficiently far away from the interface to evoke a multipole
analysis Trouilloud et al. (2008); Daddi-Moussa-Ider et al. (2019). There are,
however, systems where the microswimmers are at distances of order their size
$L$ from the interface, that aren’t necessarily constrained as in a quasi-two-
dimensional setting Montecucco and Rappuoli (2001); Moriarty et al. (2008).
Moreover, several bacteria, such as $H.pylori$, $B.subtilis$, have a strongly
orientable geometry. The orientability can be characterized in terms of an
aspect ratio $\kappa$, which is defined as the total swimmer length
(head+flagellar bundle combined) to its lateral extent. Typical values of the
aspect ratio are of the order $\kappa\sim 10$, implying that $\kappa\gg 1$
Montecucco and Rappuoli (2001); Berg (2008); Elgeti et al. (2015); Rusconi et
al. (2014). One would therefore like to relax some of the constraints and
develop a more general framework to model microscopic swimmers near deformable
interfaces. These could include aspects of the finite swimmer size and
geometry, characterizing the hydrodynamics for swimmer distances of $O(L)$
from the interface, and going beyond the limit where the ratio of the viscous
fluid stress to the deformation driving stress is small.

Here, we consider the motion of a slender microswimmer moving in the proximity
to a deformable interface that separates two Newtonian fluid regions of
different viscosities. By modulating the distribution of forces along the
swimmer length, we are able to model both pusher- and puller-type swimmers
Kasyap et al. (2014); Nambiar et al. (2021), hereafter referred to as pushers
and pullers, respectively. Importantly, using slender-body theory enables us
to resolve the hydrodynamics in regions of order the swimmer size $L$ owing to
the weak inverse logarithmic scaling of the disturbance flow field in $\kappa$
Batchelor (1970). In particular, this weak scaling also allows us to avoid the
smallness of the stress ratios discussed above and thereby access a wider
parameter regime of interfacial properties. We characterize the nature of
interface deformation and swimmer translation due to the coupled hydrodynamics
for different swimmer orientations. We treat an interface that can deform due
to both surface tension and bending elasticity, and our study extends the
framework developed to model such systems by Bickel (2007) and Daddi-Moussa-
Ider et al. (2016).

The manuscript is organized as follows. In Sec. II, we formulate the governing
equations and the boundary conditions for a slender swimmer disturbing the
fluid medium in the vicinity of a deformable interface. Specifically, in Sec.
II.1 we describe the Stokes equations of the two fluid region along with the
relevant velocity and stress boundary conditions at the interface, and in Sec.
II.2 we describe the non-dimensional system, approximations to the field
variables and the boundary conditions that exploit slender swimmer model.
These non-dimensional equations and boundary conditions are solved in Appendix
A, for swimmers oriented parallel to the interface and Appendix B for swimmers
oriented orthogonal to the interface. In Appendix C, we validate the
approximations to the boundary condition for swimmers oriented orthogonal to
the interface. We discuss the results for the swimmer translation and
interface deformation, first for swimmers oriented parallel to the interface
in Sec. III.1, and for swimmers oriented orthogonal to the interface in Sec.
III.2. In each case, we describe the results pertaining to the interface
deformation and of the effect on the swimming motion due to the hydrodynamic
coupling. In Sec. III.3 we generalize our analysis for arbitrary swimmer
orientations, focusing only on the swimmer translation and rotation. In Sec.
IV we analyze the change in the swimmer translation when it is confined
between a deformable interface and a rigid boundary, and oriented parallel to
both. Lastly, in Sec. V we summarize the results and present concluding
remarks on our analysis.

##  II Slender swimmers near a deformable interface: a coupled hydrodynamics
framework

In this section, we derive the coupled set of differential equations that
characterize the disturbance flow field and the interface deformation due to a
force- and torque-free microswimmer translating in the vicinity of a
deformable interface. As described in Sec. I we model a slender microswimmer
moving near a deformable interface at distances of order its own body length
away from it. The interface separates two density-matched Newtonian fluid
regions that have different viscosities, and is modeled such that it can
deform due to surface tension and bending elasticity as is the case for
typical lipid bilayer membranes or vesicles and cells Faucon et al. (1989);
Fygenson et al. (1997); Rosen et al. (2005); Freund (2014). A schematic is
given in figure 1.

Figure 1: A schematic representation of a slender swimmer translating near a
deformable interface with speed $V_{s}$ along its director vector $\bm{p}$.
The swimmer is characterized by an aspect ratio $\kappa$ defined as the ratio
of its total length $L$ to its lateral extent $W$. The circles along the axial
length of the fore-aft symmetric swimmer represent a line distribution of
stokeslets characterizing the head and tail for pushers. The disturbance flow
field generated by the swimming motion deforms the interface, and $u_{z}$ is
the interface deformation (solid gray line) relative to its initially flat
undeformed (dotted line) state $z_{0}$.

###  II.1 Equations for the fluid velocity and the boundary condition

The equations governing the two fluid regions are the Stokes equations and the
continuity equations which characterize the disturbance flow field $\bm{v}$
and the pressure field $P$, and are:

$\displaystyle-\nabla P_{\alpha}+\eta_{\alpha}\nabla^{2}\bm{v}_{\alpha}$
$\displaystyle=$
$\displaystyle\int_{-\frac{L}{2}}^{\frac{L}{2}}\bm{f}_{\alpha}\bm{\delta}(\bm{x}-\bm{x}_{s}-V_{s}\bm{p}t-s\bm{p})\mathrm{d}s,\;\;\text{and}$
(1a) $\displaystyle\bm{\nabla}\cdot\bm{v}_{\alpha}$ $\displaystyle=$
$\displaystyle 0.$ (1b)

Here, the subscript $\alpha\in[1,2]$ represents the two fluid regions, $\eta$
is the fluid viscosity, $\bm{\delta}(\bm{x})$ is the Dirac-delta function
Lighthill (1958), $L$ is the total swimmer length, $V_{s}$ is the swimming
speed, and $s$ refers to the distance along the slender swimmer axial
coordinate, whose center-of-mass is at position $\bm{x}_{s}$ and is oriented
along the direction $\bm{p}$. The orientation vector $\bm{p}$ is characterized
by a polar angle $\theta$ relative to the vertical and azimuthal angle $\phi$
in the plane of the undeformed interface. The kernel in the integral on the
right-hand side of Eq. (1a), $\bm{f}_{\alpha}$, represents the forcing due to
the slender swimmer and characterizes a line distribution of Stokeslets along
the axial coordinate of the swimmer. For a fore- and aft-symmetric swimmer,
which disturbs in the surrounding fluid medium as it moves, one may write
$\bm{f}_{\alpha}$ as Kasyap et al. (2014); Nambiar et al. (2021):

$\bm{f}_{\alpha}=\begin{cases}D\,\eta_{\alpha}V_{s}\bm{p}\,\text{sgn}(s)/(\ln\kappa)\\!\\!&;\alpha=1\\\
0\\!\\!&;\alpha=2,\\\ \end{cases}\\\ \\\ $ (2)

where $\kappa$ refers to the slender swimmer aspect ratio, and is defined as
the ratio of the total swimmer length to its lateral extent, as shown in Fig.
1. The above form of $\bm{f}_{\alpha}$ describes two swimming mechanisms,
namely, pushers (rear propelled) and pullers (fore propelled). In the far-
field both swimming mechanisms exhibit a disturbance flow field akin to a
point-force dipole with opposing dipole strength Lauga and Powers (2009);
Subramanian and Nott (2011). This opposing character is encoded in the
parameter $D=-1\,(+1)$ in Eq. (2) which describes the specific nature of the
force-dipole swimming mechanism, that is extensile for pushers (contractile
for pullers). Typically, pusher-type bacteria such as $E.coli$, $B.subtilis$
are fairly slender with $\kappa\sim 10$ when the swimmer length is measured
based on the cell body and the flagellar bundle length combined Brennen and
Winet (1977); Patteson et al. (2016); Kasyap et al. (2014); Rusconi et al.
(2014).

Equation (1) must be solved subject to appropriate boundary conditions at the
interface. We consider an impenetrable no-slip interface of infinitesimal
thickness within the continuum framework lea (2007). This translates to the
following velocity boundary conditions:
$\bm{v}_{1}\cdot\bm{n}=\bm{v}_{2}\cdot\bm{n}$ and
$\bm{v}_{1}\cdot(\bm{I}-\bm{nn})=\bm{v}_{2}\cdot(\bm{I}-\bm{nn})$, where
$\bm{n}=\bm{\nabla}F/|\bm{\nabla}F|$ is the unit normal to the interface
$F=z-z_{0}-u_{z}$ lea (2007), and $\bm{I}$ is the identity tensor. The
$z$-component of the interface deformation relative to its undeformed planar
state is $u_{z}$ (see Fig. 1). As described in Sec. I, the interface can
deform when there is a flow in the fluid medium surrounding it, since it has a
finite surface tension and can to bend elastically. Therefore, the normal
component of the stress undergoes a jump across the interface:
$\sigma_{1}^{zz}|_{u_{z}^{+}}-\sigma_{2}^{zz}|_{u_{z}^{-}}=\gamma\bm{\nabla}\cdot\bm{n}+\delta
F_{\text{bend}}$, where
$\bm{\sigma}_{\alpha}=-P\bm{I}+\eta_{\alpha}(\bm{\nabla}\cdot\bm{v}+\bm{\nabla}\cdot\bm{v}^{t})$
is the stress tensor, $\gamma$ the isotropic surface tension and $\delta
F_{\text{bend}}$ represents the stress jump due to bending. Here, we consider
linear elastic bending, and use the Helfrich model Helfrich  to relate the
interface deformation to the bending stress. This model has been used
extensively in the literature for analyzing swimming sheets, Stokeslets and
other higher-order singularities near deformable interfaces (see Dias and
Powers, 2013; Daddi-Moussa-Ider et al., 2016; Rallabandi et al., 2018; Daddi-
Moussa-Ider et al., 2018, 2019, and also references therein). For simplicity
we consider the tangential component of the stress jump across the interface
to be continuous, implying that
$\sigma_{1}^{xz}|_{u_{z}^{+}}=\sigma_{2}^{xz}|_{u_{z}^{-}}$ and
$\sigma_{1}^{yz}|_{u_{z}^{+}}=\sigma_{2}^{yz}|_{u_{z}^{-}}$. However, we note
that the present analytical framework allows incorporation of a shear
resistance on the interface Daddi-Moussa-Ider et al. (2019). Now, in addition
to these fluid equations and boundary conditions, we also have the kinematic
condition at the interface, which provides a direct relation between the
interface deformation and the disturbance flow field measured at the interface
Berdan and Leal (1982); lea (2007); Shaik and Ardekani (2017), and is:

$\frac{\partial F}{\partial t}+\bm{v}\cdot\bm{\nabla}F=0.$ (3)

With the interface deformation and the disturbance field characterized, we
express the swimmer translation velocity and rate of rotation due to the
hydrodynamic interaction with the interface as follows. From viscous slender-
body theory Batchelor (1970); Kim and Karrila (2005) the velocity $\bm{V}$ for
a slender swimmer is:

$\bm{V}(s\bm{p})-\bm{v}(s\bm{p})|_{\text{swim
location}}=\frac{\bm{f}_{1}}{4\pi\eta_{1}}\cdot(\bm{I}+\bm{pp})\ln\kappa,$ (4)

where the second term on the left-hand side of Eq. (4) is the disturbance flow
field at the location of the swimmer. For a rigid slender swimmer,
$\bm{V}=\bm{V}^{T}+\bm{\omega}\wedge\bm{p}s$, with $\bm{V}^{T}$ being the
translational velocity and $\bm{\omega}$ the angular velocity of the swimmer
emanating from the coupled hydrodynamics. The force-free
($\int_{-\frac{1}{2}}^{\frac{1}{2}}\bm{f}_{1}\mathrm{d}s=0$) and torque-free
($\int_{-\frac{1}{2}}^{\frac{1}{2}}s\bm{p}\wedge\bm{f}_{1}\mathrm{d}s=0$)
swimming constraints translate to the following expressions for the swimmer
translational velocity $\bm{V}^{T}$ and rotation rate
$\dot{\bm{p}}=\bm{\omega}\wedge\bm{p}$ Nambiar et al. (2021),

$\displaystyle\bm{V}^{T}$ $\displaystyle=$
$\displaystyle\int_{-\frac{1}{2}}^{\frac{1}{2}}\bm{v}(s\bm{p})\mathrm{d}s,\;\;\text{and}$
(5a) $\displaystyle\dot{\bm{p}}$ $\displaystyle=$ $\displaystyle
12\int_{-\frac{1}{2}}^{\frac{1}{2}}(\bm{I}-\bm{pp})\cdot\bm{v}(s\bm{p})s\mathrm{d}s.$
(5b)

Therefore, in addition to the swimmers intrinsic motion, a slender swimmer
translates along the length averaged disturbance field at its center and
rotates due to the first moment of the disturbance field along its axial
coordinate.

###  II.2 Non-dimensionalization

We non-dimensionalize the equations and the boundary conditions using $L$ as a
length scale, $V_{s}$ as a velocity scale and $\eta_{1}V_{s}/L$ as a scale for
the fluid stress. It is convenient to solve the problem in a reference frame
moving with the swimmer centre:
$\bm{r}=\bm{x}-\bm{x}_{s}-(V_{s}t_{c}/L)\bm{p}\,t$;
$\bm{\nabla}\equiv\bm{\nabla}_{\bm{r}}$ and $\partial(\cdot)/\partial
t=\partial(\cdot)/\partial t-(V_{s}t_{c}/L)\bm{p}\cdot\nabla_{\bm{r}}(\cdot)$,
where we have allowed the time scale $t_{c}$ to remain arbitrary for the
moment. This yields the following non-dimensional form of the Stokes equations
and the continuity equations, Eq. (1), in the two fluid regions:

$\displaystyle-\nabla P_{1}+\nabla^{2}_{r}\bm{v}_{1}$ $\displaystyle=$
$\displaystyle\frac{D\,\bm{p}}{\ln\kappa}\int_{-\frac{1}{2}}^{\frac{1}{2}}\text{sgn}(s)\bm{\delta}\left(\bm{r}-s\bm{p}\right)\mathrm{d}s,$
(6a) $\displaystyle\bm{\nabla}_{\bm{r}}\cdot\bm{v}_{1}$ $\displaystyle=$
$\displaystyle 0,$ (6b)
$\displaystyle-\nabla_{\bm{r}}P_{2}+\lambda\nabla^{2}_{r}\bm{v}_{2}$
$\displaystyle=$ $\displaystyle 0,\;\text{and}$ (6c)
$\displaystyle\bm{\nabla}_{\bm{r}}\cdot\bm{v}_{2}$ $\displaystyle=$
$\displaystyle 0,$ (6d)

where Eqs. (6a) and (6c) correspond to Eq. (1a), and Eqs. (6b) and (6d) to Eq.
(1b) for $\alpha\in[1,2]$. Note that, for brevity, we retain the original
notation of the field variables in the non-dimensional Eqs. (6a)-(6d). Now,
the leading order flow field generated by a slender-body is
$O(\ln\kappa)^{-1}$ Batchelor (1970); Kim and Karrila (2005); Lauga and Powers
(2009); Kasyap et al. (2014), where we emphasize that the aspect ratio
$\kappa\gg 1$. As the disturbance field $\bm{v}$ in the Stokes equations (6a)
is linear in the forcing, it will also scale as $O(\ln\kappa)^{-1}$, and is
therefore asymptotically small. Hence, we can consider small interface
deformations without having to invoke the smallness of the ratio of viscous
fluid stresses to the deformation driving stresses on the interface (surface
tension and bending). In other words, owing to the weak flow field of a
slender-body, the slender swimmer model enables simplifications of the
boundary conditions.

For small deformations, we can linearize the expressions for the unit normal
and the curvature at the interface as: $\bm{n}\approx\mathbf{1}_{z}$ and
$\bm{\nabla}_{\bm{r}}\cdot\bm{n}\approx-\Delta_{r_{\|}}u_{z}$, where
$\mathbf{1}_{z}$ is the $r_{z}$-component of the unit normal Bickel (2007);
Daddi-Moussa-Ider et al. (2016). Within the linear response framework the
height function, $u_{z}$, depends only upon the $r_{x}$ and $r_{y}$
coordinates, and an analytical expression for the stress jump due to bending
using the Helfrich elastic energy function dimensionally is: $\delta
F_{\text{bend }}=\kappa_{\beta}\Delta_{\|}^{2}u_{z}$, where $\bm{\Delta}_{\|}$
represents the Laplacian in the $r_{x}-r_{y}$ plane and $\kappa_{\beta}$ the
bending modulus of the interface Daddi-Moussa-Ider et al. (2016). Thus, the
non-dimensional form of the boundary conditions specified above Eq. (3) are:

$\displaystyle v_{1z}|_{r_{z_{0}}^{+}}$ $\displaystyle=$ $\displaystyle
v_{2z}|_{r_{z_{0}}^{-}},$ (7a)
$\displaystyle\bm{v}_{1}\cdot\left.(\bm{I}-\bm{1}_{z}\bm{1}_{z})\right|_{r_{z_{0}}^{+}}$
$\displaystyle=$
$\displaystyle\bm{v}_{2}\cdot\left.(\bm{I}-\bm{1}_{z}\bm{1}_{z})\right|_{r_{z_{0}}^{-}},$
(7b) $\displaystyle\left.\left(\frac{\partial v_{1z}}{\partial
r_{x}}+\frac{\partial v_{1x}}{\partial r_{z}}\right)\right|_{r_{z_{0}}^{+}}$
$\displaystyle=$ $\displaystyle\lambda\left.\left(\frac{\partial
v_{2z}}{\partial r_{x}}+\frac{\partial v_{2x}}{\partial
r_{z}}\right)\right|_{r_{z_{0}}^{-}},$ (7c)
$\displaystyle\left.\left(\frac{\partial v_{1z}}{\partial
r_{y}}+\frac{\partial v_{1y}}{\partial r_{z}}\right)\right|_{r_{z_{0}}^{+}}$
$\displaystyle=$ $\displaystyle\lambda\left.\left(\frac{\partial
v_{2z}}{\partial r_{y}}+\frac{\partial v_{2y}}{\partial
r_{z}}\right)\right|_{r_{z_{0}}^{-}},\;\text{and}$ (7d)
$\displaystyle-\left(P_{1}|_{r_{z_{0}}^{+}}-P_{2}|_{r_{z_{0}}^{-}}\right)+2\left(\left.\frac{\partial
v_{1z}}{\partial r_{z}}\right|_{r_{z_{0}}^{+}}-\lambda\left.\frac{\partial
v_{2z}}{\partial r_{z}}\right|_{r_{z_{0}}^{-}}\right)$ $\displaystyle=$
$\displaystyle-\frac{\gamma}{\eta_{1}V_{s}}\Delta_{r_{\|}}u_{z}+\frac{\kappa_{\beta}L^{2}}{\eta_{1}V_{s}}\Delta^{2}_{r_{\|}}u_{z},$
(7e)

where $\lambda=\eta_{2}/\eta_{1}$ is the viscosity ratio. We are studying an
interfacial deformation process and to non-dimensionalize the kinematic
condition Eq. (3), we can choose between the scale for the elastic bending and
surface tension, corresponding to $t_{c1}=\eta_{1}L^{3}/\kappa_{\beta}$ and
$t_{c2}=\eta_{1}L/\gamma$, respectively. For ultra-soft interfaces, the
surface tension can be $\gamma\sim O(10)^{-5}N/m$ or smaller Faucon et al.
(1989); de Hoog and Lekkerkerker (2001); Rosen et al. (2005), and for example
the bending rigidity of vesicles is $\kappa_{\beta}\sim O(10)^{-19}J$ Faucon
et al. (1989); Song and Waugh (1990, 1993); Fygenson et al. (1997). In a fluid
medium with a viscosity approximately that of water, $\eta_{1}\sim
O(10)^{-3}\mathrm{Ns}/\mathrm{m}^{2}$, and a swimmer of size $L\sim
O(10)^{-5}m$, this implies that $t_{c1}\sim O(1)-O(10)^{2}$ and $t_{c2}\sim
O(10)^{-3}-O(1)$. Independent of the time scale chosen, the relative
importance of the surface tension and bending stress on the right-hand side of
Eq. (7e) can then be characterized by the ratio of the two scales, namely
$\gamma L^{2}/\kappa_{\beta}$. Depending on the stiffness of the interface,
for typical values of the time scales, $\gamma L^{2}/\kappa_{\beta}$ can vary
from $O(1)$ to $\gamma L^{2}/\kappa_{\beta}\gg 1$. Therefore, we choose
$t_{c1}$ as a time scale noting that choosing $t_{c2}$ is equally plausible.
In the literature, $\gamma L^{2}/\kappa_{\beta}$ is also treated as an elasto-
capillary length Style et al. (2017); Andreotti and Snoeijer (2020), to
characterize the length scales over which surface tension or elasticity
dominate. The dimensionless kinematic boundary condition Eq. (3) becomes

$\frac{\partial u_{z}}{\partial
t}+\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\bm{v}\cdot\bm{\nabla}_{\bm{r}_{\|}}u_{z}=\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\Big{[}\left.v_{z}\right|_{u_{z}}+p_{z}\Big{]},$
(8)

where $p_{z}=\cos\theta$, and the right-hand side of Eq. (8) applies at the
location of the deformed interface $u_{z}$.

In principle, there is a third time scale in the problem if we account for the
intrinsic reorientation time scale of a microswimmer. For instance,
$t_{c3}=\tau$ or $t_{c3}=D_{r}^{-1}$, where $\tau$ is the mean run duration of
a run-and-tumble particle (RTP) Berg (2008); Koch and Subramanian (2011);
Elgeti et al. (2015); Ezhilan et al. (2015); Lauga (2016); Berg (2018) and
$D_{r}$ the rotary diffusivity coefficient of an active Brownian particle
(ABP) Koch and Subramanian (2011); Rusconi et al. (2014); Takatori et al.
(2014); Berg (2018); Vennamneni et al. (2020). For simplicity, we neglect this
third scale and treat only straight-swimmers, that is, the class of swimmers
that do not tend to intrinsically reorient. In Sec. V, we briefly discuss the
relevance of the swimmer reorientation time and how it compares to the
response scale of the interface deformation.

For small deformations, one can use the method of domain perturbation lea
(2007); Daddi-Moussa-Ider et al. (2018) to represent this term at the location
of the planar undeformed interface $z_{0}$ instead. Namely, $v_{z}|_{u_{z}}$,
can be perturbed about $r_{z_{0}}$ as: $v_{z}|_{u_{z}}\approx
v_{z}|_{r_{z_{0}}}+u_{z}(\partial v_{z}/\partial
r_{z})|_{r_{z_{0}}}+O(u_{z})^{2}+\ldots$ , and therefore considering terms up
to $O(u_{z})$ in the velocity, one may rewrite Eq. (8) as:

$\frac{\partial u_{z}}{\partial
t}+\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\bigg{[}\bm{v}\cdot\bm{\nabla}_{\bm{r}_{\|}}u_{z}-u_{z}\frac{\partial
v_{z}}{\partial
r_{z}}\bigg{]}=\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\Big{[}\left.v_{z}\right|_{r_{z_{0}}}+p_{z}\Big{]}.$
(9)

Note that the second and third nonlinear terms on the left-hand side of Eq.
(9) involve products of the field variable ($u_{z}$, $\bm{v}$) and their
associated gradients. These terms are usually neglected if the forcing is
sufficiently far away from the interface so that terms involving their
products are asymptotically small in the far-field Bickel (2007); Daddi-
Moussa-Ider et al. (2016). Given that we consider swimmers at $O(1)$ distances
from the interface, these terms can not be neglected a-priori. However, as
noted below Eq. 6, the disturbance field $\bm{v}$ is $O(\ln\kappa)^{-1}$.
Therefore, if we express the interface deformation as a perturbation series in
$(\ln\kappa)^{-1}$,
$u_{z}=u_{z1}\times(\ln\kappa)^{-1}+u_{z2}\times(\ln\kappa)^{-2}+\ldots$, the
second and third terms on the left-hand side of Eq. (9) will be
$O(\ln\kappa)^{-2}$, and hence asymptotically smaller than the disturbance
field $\bm{v}$. A similar argument, owing to the asymptotically weak flow
field of a slender swimmer, has recently been used to treat the tracer
diffusivity and the velocity variance in a suspension of interacting slender
swimmers in bulk Nambiar et al. (2021). Thus, to a leading order in
$(\ln\kappa)^{-1}$, one can simplify Eq. (9) as:

$\frac{\partial u_{z}}{\partial
t}\approx\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\Big{[}\left.v_{z}\right|_{r_{z_{0}}}+p_{z}\Big{]},$
(10)

for a slender-body near a deformable interface. We can solve the above system
of linear equations (6a)-(6d) and the boundary conditions (7a)-(7e) in Fourier
space Lighthill (1958); Bickel (2007); Daddi-Moussa-Ider et al. (2016), and we
define the two-dimensional Fourier transform of any variable $A(r_{x},r_{y})$
as:
$\hat{A}(\bm{k})=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}x\mathrm{d}yA(r_{x},r_{y})\exp(-2\pi{i\mkern
1.0mu}\bm{k}\cdot\bm{r}_{\|})$. Determining the Fourier transformed
disturbance field $\hat{v}_{z}$ is straightforward, as shown in Appendix A, B
for swimmers oriented parallel and perpendicular to the interface,
respectively. However, here we show the Fourier transformed kinematic boundary
condition Eq. (10):

$\frac{\text{d}\hat{u}_{z}}{\text{d}t}=\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\left[\hat{v}_{z}|_{r_{z_{0}}}+p_{z}\frac{\delta(k)}{\pi
k}\right],$ (11)

which is solved in the next section for different swimmer configurations
relative to the interface. For swimmers translating parallel to the interface
$p_{z}=0$, in which case the right-hand side of Eq. (11) is only the Fourier
transformed disturbance velocity at the undeformed interface. This
approximation holds as long as the ratio of the viscous to bending stress:
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}$ remains $O(1)$. However, if
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}\gg 1$, the nonlinear boundary condition
given by Eq. (9) would have to be used. In Appendix C, we validate the above
approximation for the interface deformation by comparing $u_{z}$ derived from
solving Eq. (9) with that from Eq. (10), for the specific case of swimmers
oriented orthogonal to the interface over a range of
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}$. In the following sections, we fix the
aspect ratio $\kappa=10$ and the ratio of the viscous to bending stresses
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.

##  III  Swimmers near a single interface

In general, simultaneously solving Eq. (5) and Eq. (11) results in an
inherently unsteady problem, with the coupling involving a dependence on both
the distance of the swimmer to the interface as well as its orientation.
Therefore, in what follows, we first discuss the results for a swimmer
oriented parallel to the interface, then for a swimmer orthogonal to the
interface, and finally for a swimmer with an arbitrary orientation. In Sec. IV
we consider a swimmer confined between a rigid boundary and a deformable
interface.

###  III.1  Microswimmers swimming parallel to the interface

To understand the coupled hydrodynamics, it is useful to refer to the equation
for the interface deformation and for the swimmer motion for swimmers oriented
parallel to the interface. For the former, we use the expression for
$\hat{v}_{1z}|_{r_{z_{0}}}$ from Eq. (39) in Appendix A.3 and obtain the
following equation for the time evolution of $\hat{u}_{z}$:

$\displaystyle\frac{\text{d}\hat{u}_{z}}{\text{d}t}$ $\displaystyle+$
$\displaystyle\frac{\pi
k}{(1+\lambda)}\left(4\pi^{2}k^{2}+\Gamma\right)\hat{u}_{z}=\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\left[\frac{D}{\pi
k(1+\lambda)\ln\kappa}\sin^{2}\left(\frac{\pi}{2}kp_{l}\right)r_{z_{0}}\exp(2\pi
kr_{z_{0}})\right],$ (12)

where $\Gamma\equiv\gamma L^{2}/\kappa_{\beta}$ characterizes the relative
importance of surface tension to bending stress. Note that the interface
deformation is linked to the swimmers instantaneous configuration via
$r_{z_{0}}$ and $\bm{p}$, both of which can evolve in time. Here,
$\dot{\bm{p}}=0$, as is also true for a generic force-dipole swimmer near a
planar or weakly deforming boundary Dunstan et al. (2012); Spagnolie and Lauga
(2012); Daddi-Moussa-Ider et al. (2019). However, the same is not true for
$\bm{V}^{T}$, since the $z$-component of the translation is responsible for
the coupling with the interface deformation in Eq. (12) via $r_{z_{0}}$. We
only focus on $V_{z}^{T}$, which is,

$\displaystyle V_{z}^{T}\equiv\frac{\mathrm{d}r_{z_{0}}}{\mathrm{d}t}$
$\displaystyle=$ $\displaystyle-\int\mathrm{d}\bm{k}\frac{\sin(\pi
kp_{l})}{\pi
kp_{l}}\left[\frac{2Dr_{z_{0}}^{2}}{\ln\kappa}\left(\frac{1-\lambda}{1+\lambda}\right)\sin^{2}\left(\frac{\pi}{2}kp_{l}\right)\exp(4\pi
kr_{z_{0}})\right.$ (13)
$\displaystyle\left.\qquad+\left(\frac{\kappa_{\beta}}{\eta_{1}V_{s}L^{2}}\right)\frac{\pi
k}{(1+\lambda)}\left(4\pi^{2}k^{2}+\Gamma\right)(1-2\pi
kr_{z_{0}})\hat{u}_{z}\exp(2\pi kr_{z_{0}})\right],$

where we have used the definition of the inverse Fourier transform:
$A(r_{x},r_{y})=\int\mathrm{d}\bm{k}\hat{A}$ $\exp(2\pi{i\mkern
1.0mu}\bm{k}\cdot\bm{r}_{\|})$; $\bm{k}=(k_{x},k_{y})$ and
$\bm{r}_{\|}=(r_{x},r_{y})$. Thus, characterizing the coupled hydrodynamics
due to a swimmer translating parallel to the interface reduces to
simultaneously solving Eqs. (12) and (13) to obtain $\hat{u}_{z}$ and
$r_{z_{0}}$, and then the disturbance flow field given by Eq. (36) and Eq.
(39) in Appendix A.3.

We now interpret the nature of the deformation influenced by the disturbance
field generated from the swimmer. In Fig. 2a we plot the three-dimensional
contour of the interface deformation due to a pusher in fluid region 1 for the
parameters $\lambda=0.5$ and $\Gamma=0.1$. As expected, for pushers (or
pullers) oriented parallel to the interface, the deformation is not radially
symmetric Shaik and Ardekani (2017). The qualitative nature of the deformation
remains the same for all values of $\lambda$ (not shown), with the magnitude
of the deformation decreasing as $\lambda$ increases, as expected from Eq.
(12), noting that for $\lambda\gg 1$, $u_{z}\propto 1/\lambda$.

(a) Three-dimensional contour of the interface deformation

(b) Two-dimensional contour of the interface deformation

Figure 2: The non-dimensional interface deformation plotted for a swimmer of
unit length oriented along the $r_{x}-$direction (a) viewed at an angle of 30∘
from the undeformed interface and (b) top view. The black line in (b) spanning
from $r_{x}\in[-1/2,1/2]$ signifies the slender swimmer. The colorbar given on
the right-hand side of (b) applies to both figures. The plot is for
$\Gamma=0.1$, $\lambda=0.5$ at $t\approx 2$, $\kappa=10$ and
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.

To characterize the interface response with $\Gamma$, the deformation is
plotted as a function of the radial distance parallel and orthogonal to the
swimmer axis in Fig. 3 for $\lambda=0.5$ and $\Gamma=$ $0.1,1,5$. It is clear
that increasing the importance of surface tension, decreases the magnitude of
the deformation. In the far-field, however, at any finite time, the
deformation is independent of $\Gamma$ and scales as $O(r_{\|})^{-3}$, as
shown in the insets of Fig. 3. In fact, for any small but finite time, one can
obtain an analytical expression for the interface deformation assuming
$r_{z_{0}}$ to be nearly constant; this constraint of a constant $r_{z_{0}}$
is reasonable up to an $O(1)$ change in time, as will be seen below while
interpreting the swimmer translation. The resulting expression for $u_{z}$ is:

$\left.u_{z}\right|_{r_{\|}\gg
1}\approx\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\frac{D}{(1+\lambda)\ln\kappa}\frac{r_{z_{0}}t}{16\pi}\;\frac{r_{\|}^{2}-3r_{x}^{2}}{r_{\|}^{5}},$
(14)

highlighting independence from $\Gamma$, and the $O(r_{\|})^{-3}$ scaling in
the far-field. The approximate form of $u_{z}$ given by Eq. (14) is shown in
the insets of Fig. 3, and it agrees well with the numerically determined
values.

(a) Along swimmer orientation

(b) Orthogonal to swimmer orientation

Figure 3: The interface deformation, $u_{z}$, plotted for a swimmer oriented
along the $r_{x}-$coordinate for three values of $\Gamma=$ 0.1, 1 and 5 (a)
along $r_{x}$ (b) along $r_{y}$. In both the figures $\lambda=0.5$,
$r_{z_{0}}(0)=1$, $t\approx 2$, $\kappa=10$ and
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$. Legends under $\circ$ refer to $u_{z}$
obtained from solving the pair of equations (12), (13), and those under
$\diamond$ from solving the approximate equations (15), (16). The (black)
dotted line in the insets is the far-field approximation from Eq. (14).

Now, the interface deformation is also a function of the distance of a swimmer
from the interface: At distances of order the swimmer size, one needs to
numerically solve Eqs. (12) and (13) simultaneously, whereas simplifications
can be made when the swimmer is farther away. We rescale the wavevector in the
interface deformation, Eq. (12), as $\bar{k}\equiv k\bar{r}_{z_{0}}$ and time
as $\bar{t}=t/\bar{r}_{z_{0}}$, which is valid in the regime
$\bar{r}_{z_{0}}\equiv\left|r_{z_{0}}(0)\right|\gg 1$. Seeking simplifications
appropriate for $\bar{k}$, $\bar{t}\sim O(1)$ Hinch (1991) yields the
following reduced equation for the interface deformation:

$\displaystyle\frac{\text{d}\hat{u}_{z}}{\text{d}\bar{t}}+\frac{\pi\bar{k}}{(1+\lambda)}\Gamma\hat{u}_{z}=\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\left[\frac{\pi
Dr_{z_{0}}}{4(1+\lambda)\ln\kappa}\bar{k}p_{l}^{2}\exp\left(2\pi\bar{k}\frac{r_{z_{0}}}{\bar{r}_{z_{0}}}\right)\right].$
(15)

The simplified expression for the vertical component of the swimmer
translation is

$V_{z}^{T}=-\frac{3D}{256\pi\ln\kappa}\left(\frac{1-\lambda}{1+\lambda}\right)\frac{1}{r_{z_{0}}^{2}}+\left(\frac{\kappa_{\beta}}{\eta_{1}V_{s}L^{2}}\right)\int\mathrm{d}\bm{k}\frac{\pi
k\Gamma}{(1+\lambda)}\left(1-2\pi\bar{k}\frac{r_{z_{0}}}{\bar{r}_{z_{0}}}\right)\hat{u}_{z}\exp\left(2\pi\bar{k}\frac{r_{z_{0}}}{\bar{r}_{z_{0}}}\right).$
(16)

Eq. (15) highlights two important qualitative changes in the interface
deformation for swimmers that are far from the interface. First, the interface
deformation is controlled by surface tension alone, as the term proportional
to bending is $O(\bar{r}_{z_{0}})^{-2}$ smaller than $\Gamma$ on the left-hand
side of Eq. (12), and hence is absent from the leading order approximation.
Second, the above rescalings of $k$ and $t$ recover the original dependence of
the rescaled variables on $r_{z_{0}}$, as in Eq. (12). This implies that, for
$r_{z_{0}}\gg 1$, the time evolution slows down in proportion to $r_{z_{0}}$.
In other words, swimmers at an $O(1)$ distance from the interface influence an
$O(1)$ growth of the interface deformation at $O(1)$ times. However, a swimmer
that is farther away ($r_{z_{0}}\gg 1$) will influence an $O(1)$ growth of the
interface deformation only when $t\sim O(r_{z_{0}})$. In the insets of Fig. 3
we plot the interface deformation determined from simultaneously solving Eqs.
(15) and (16). The agreement of both equations is excellent for $\Gamma=5$ for
$r_{\|}\sim O(1)$, but the deviation of the approximate solution is large for
$\Gamma=0.1$. This is expected as Eq. (15) relies on $\Gamma$ being larger
than the bending term, and hence, for $\bar{r}_{z_{0}}=1$ it relies solely on
the largeness of $\Gamma$. At large $r_{\|}$, the agreement is good for the
same reason underlying Eqs. (15) and (16), but now for $r_{\|}\gg 1$ instead
of $r_{z_{0}}\gg 1$. For large $\bar{r}_{z_{0}}$ (not shown), the
approximations in Eqs. (15) and (16) agree well even for $\Gamma<O(1)$. We
note that the far-field approximation corresponds to modeling the swimmer as a
point force-dipole, and the disagreement at small to moderate values of
$\Gamma$, particularly for $r_{z_{0}}\sim O(1)$ is indicative of the dipolar
swimmer overestimating the interface deformation.

A result arising from the neglect of the bending term in Eq. (15) is a subtle
qualitative change in the deformation behavior when bending dominates
$(\Gamma\ll 1)$ versus when surface tension dominates $(\Gamma\gg 1)$. In
particular, following a zero-crossing at an $O(1)$ distance away from the
swimmer, when $\Gamma\gg 1$, the deformation monotonically transitions to an
$O\left(r_{\|}\right)^{-3}$ far-field character, whereas for $\Gamma$ up to
$O(1)$, the deformation exhibits an intermediate scaling before transitioning
to the far-field decay. As shown in the insets of Fig. 3, there is a
characteristic maxima in $\left|u_{z}\right|$ for $\Gamma=0.1$ (see the red
curve corresponding to the solution of Eqs. (12) and (13)). These features are
absent in the approximate Eqs. (15) and (16), and are most evident in Fig. 3a
(see orange and dark green curves in the inset). This qualitative difference
emerges since the bending mode, which dominates in this regime, depends on
higher order derivatives in the curvature as specified by the stress boundary
condition, Eq. (7e).

In Fig. 4, we plot the vertical translation velocity $V_{z}^{T}$ of a pusher
oriented parallel to the interface for three specific viscosity ratios
$\lambda=$ 0.5, 1 and 1.5. We note two interesting features of the vertical
migration. First, as shown in Fig. 4a, when $\lambda<1$, $V_{z}^{T}$ is
positive for short times and negative for long times. That is, a pusher is
repelled from the interface at short times, while at long times it is
eventually attracted to the interface. Second, the crossover time from
repulsion to attraction is shorter for swimmers closer to the interface. For
$r_{z_{0}}\sim O(1)$, the swimmer undergoes a change in its vertical
translational velocity at a $t_{\text{transition}}<O(1)$, whereas, when
$r_{z_{0}}\gg 1$, $t_{\text{transition}}\sim O(r_{z_{0}})$; here,
$t_{\text{transition}}$ is the time taken for $V_{z}^{T}$ to vanish. For
instance, when $r_{z_{0}}\in[1,2]$, $t_{\text{transition}}<0.5$ and when
$r_{z_{0}}\geq 4$, $t_{\text{transition}}>1$. This is unlike the case of
$\lambda\geq 1$ as shown in Fig. 4b and c, where the pusher is monotonically
attracted to the interface. Pullers (not shown) behave in the exact opposite
manner, that is, for $\lambda<1$ they are attracted to the interface at short
times and repelled at long times, with the crossover time for attraction to
repulsion having the same character; for $\lambda\geq 1$ pullers are
monotonically repelled from the interface.

(a) $\lambda=0.5$

(b) $\lambda=1$

(c) $\lambda=1.5$

Figure 4: The vertical swimmer translation velocity $V_{z}^{T}$ of pushers
plotted as a function of the initial distance of the swimmer from the
interface $r_{z_{0}}$ at different time instants for the viscosity ratios (a)
$\lambda=0.5$, (b) $\lambda=1$ and (c) $\lambda=1.5$. The plot is to be
interpreted as follows: at each $r_{z_{0}}$, the vertical black dotted lines
trace the time evolution of $V_{z}^{T}$. In all the plots $\kappa=10$,
$\Gamma=1$ and $\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.

In Fig. 5, we plot the relative swimmer trajectories of pushers as a function
of time for swimmers starting from different initial locations $r_{z_{0}}(0)$,
for $\lambda=$ 0.5, 1 and 1.5 as well. Owing to the repulsion ($V_{z}^{T}>0$)
at short times when $\lambda<1$ as highlighted in Fig. 4a, the relative
swimmer trajectory: $r_{z_{0}}(t)-r_{z_{0}}(0)\geq 0$ in Fig. 5a. Moreover, as
expected from the discussion surrounding Fig. 4, the time it takes for the
swimmer to return to its initial location depends sensitively upon its initial
distance from the interface. For instance, when $r_{z_{0}}(0)=1$, the return
time is about $t\sim 0.25$, whereas, when $r_{z_{0}}(0)=5$, return time is
$t\sim 3.7$. Such a non-monotonic vertical swimmer translation for $\lambda<1$
stands in contrast to what has been observed for swimmers near rigid
boundaries Berke et al. (2008); Li and Tang (2009); Spagnolie and Lauga
(2012); Elgeti et al. (2015); Ezhilan et al. (2015), where the pushers are,
hydrodynamically, only attracted to the interface. Even for model swimmers
near deformable interfaces, the swimmer migration has been reported to behave
monotonically for a given $\lambda$ Dias and Powers (2013); Shaik and Ardekani
(2017). For instance, Dias and Powers (2013) noted that the average pumping
velocity between a swimming sheet and the interface was of negative for
$\lambda<1$ and positive for $\lambda>1$. Shaik and Ardekani (2017) calculated
the vertical translational velocity of a spherical squirmer over a range of
$\lambda$, and did not find a qualitative change in the swimmer migration
across $\lambda=1$ (see Figs. 2 and 8 therein).

Figure 5: The relative vertical swimmer trajectory $r_{z_{0}}(t)-r_{z_{0}}(0)$
of pushers plotted as a function of time, for the viscosity ratios
$\lambda=0.5$, 1 and 1.5. In all the plots $\kappa=10$, $\Gamma=1$ and
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.

The first of the above observations regarding the non-monotonic swimmer
translational for $\lambda<1$ can be explained from a closer inspection of
$V_{z}^{T}$ given by Eq. (13). The translational velocity has two principal
contributions: (1) the instantaneous contribution from Stokes flow that is
proportional to $(1-\lambda)/(1+\lambda)$, (2) the time dependent term
associated with the interface deformation $\hat{u}_{z}$. In fact, at $t=0$,
$V_{z}^{T}$ (red circles in Fig. 4) is precisely the instantaneous
contribution before the interface deforms. The time dependent term remains
negative at all times $t>0$, and its magnitude gradually increases as the
interface deformation grows. However, the instantaneous term changes sign
depending on whether $\lambda$ is greater than or less than unity, although
its magnitude remains nearly constant. For pushers, when $\lambda<1$ ($>1$),
the instantaneous term is positive (negative), and it vanishes when
$\lambda=1$. Its magnitude remains nearly constant because it changes only via
the hydrodynamic interaction-induced change in the swimmer vertical motion
$r_{z_{0}}(t)$, which does not change appreciably over $O(1)$ times implying
that $\Delta r_{z_{0}}(t)\equiv r_{z_{0}}(t)-r_{z_{0}}(0)\ll r_{z_{0}}(0)$.
Hence, among the two fluid regions, if a pusher is present in the more viscous
fluid ($\lambda<1$), it starts off being repelled at short times due to the
dominance of the repulsive instantaneous term. At $t=t_{\text{transition}}$,
the growing time dependent deformation term exactly matches the nearly
constant instantaneous Stokes term, whence $V_{z}^{T}$ vanishes, and this is
characterized by the peak in the relative swimmer trajectory in Fig. 5 (left
panel). For $t>t_{\text{transition}}$ the swimmer experiences an attraction to
the interface owing to the dominance of the time dependent term. This explains
the repulsion of pushers from the interface at short times and attraction at
long times only when it is in the more viscous fluid.

To explain the second observation of the time dependence of the swimmer
translation as a function of its distance from the interface, we note that the
time dependence in the problem originates from the kinematic boundary
condition of the interface deformation given by Eq. (12). The swimmer
translation velocity $V_{z}^{T}$(and thus, $r_{z_{0}}$) obeys the same time
dependence as $u_{z}$ due to the linear relationship between the two field
variables in Eq. (13), independent of $\lambda$. Therefore, such a delay in
the time to return back to its original position when $\lambda<1$ is expected
following the arguments explained in the context of the interface deformation.
That is, when $r_{z_{0}}\sim O(1)$, the interface deformation, and hence,
$V_{z}^{T}$ evolves on $O(1)$ time units. However, when $r_{z_{0}}\gg 1$, the
interface deformation evolves on a time scale of $O(r_{z_{0}})$. Therefore, as
shown in Fig. 5 for $\lambda=0.5$, the farther the pusher (puller) is from the
interface, the longer it will sense a repulsion (an attraction), albeit with
reduced intensity.

(a)

(b)

Figure 6: The relative swimmer trajectory $r_{z_{0}}(t)-r_{z_{0}}(0)$ of a
pusher plotted as a function of time. In (a) the initial swimmer distance to
the interface is fixed: $\left|r_{z_{0}}\right|=1$ and the ratio of surface
tension to bending stress is varied: $\Gamma=$ 0.1, 1, 10. In (b) both
$r_{z_{0}}(0)$ and $\Gamma$ are varied so as to yield similar swimmer
translation, with $r_{z_{0}}(0)=1.5,2$ and $\Gamma=1,4$, respectively. In both
figures $\lambda=$ 0.5, $\kappa=10$ and $\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.
In the inset the corresponding swimmer translation velocities $V_{z}^{T}$ are
plotted for the same parameters.

To show the relative importance of surface tension to bending stress on the
swimmer motion, in Fig. 6a we plot the relative swimmer translation as a
function of time for different values of $\Gamma(\equiv\gamma
L^{2}/\kappa_{\beta})$. Clearly, increasing the surface tension results in a
more rapid translation, with a dramatic reduction in the time spent in the
repulsive state. For instance, as $\Gamma$ is increased from 0.1 to 10, the
return time decreases from $t\approx 0.28$ to $t\approx 0.06$. Moreover, as
evident from the inset of Fig. 6a, the swimmer translation velocity approaches
a quasi-steady state when $\Gamma\gg 1$. Note the emphasis on ‘quasi-steady
state’, as the problem is inherently an unsteady one, even though $V_{z}^{T}$
appears to be almost steady for $\Gamma=10$ following the initial transient.
This is because as the swimmer is moving closer to the interface, it will
continue to accelerate, given that the interface undergoes a stronger
deformation from swimmers closer to it. In turn, the translation velocity will
continue to grow. Nevertheless, this acceleration remains small after
sustaining an initial growth in the deformation, as the driving force is still
small (see the scales along the ordinate of Fig. 6a).

This analysis shows that the swimmer translation must be analyzed in terms of
both distance to the interface and interface properties. That is, both
$r_{z_{0}}$ and $\Gamma$ have comparable effects on the swimmer translation
when these quantities are large. Note, however, that owing to the exponential
damping of $u_{z}$ and $V_{z}^{T}$ with respect to $r_{z_{0}}$ even an $O(1)$
change in the distance of the swimmer to the interface results in a large
reduction of the swimmer translational velocity. Therefore, a swimmer near an
interface with relatively low surface tension can exhibit similar
translational dynamics to that of a swimmer that is slightly farther away from
the interface with higher surface tension. An example of such a trajectory is
highlighted in Fig. 6b, where the initial locations of the swimmers are
displaced by 0.5 units of $r_{z_{0}}$, but $\Gamma$ is more than doubled. For
an $O(1)$ change in time, the trajectories in two scenario are similar.

###  III.2  Microswimmers swimming orthogonal to the interface

We now consider swimmers oriented orthogonal to the interface ($p_{z}=\pm 1$
and $p_{x}=p_{y}=0$), since the results are easier to interpret relative to
swimmers with arbitrary orientation. Here too, the rate of rotation
$\dot{\bm{p}}=0$, and therefore characterizing the swimmer motion requires
simultaneously solving for the interface deformation and the vertical
translation.

The equation governing the interface deformation is obtained using the
expression for $\hat{v}_{1z}|_{r_{z_{0}}}$ from Eq. (46) (Appendix B) in Eq.
(11), and is:

$\displaystyle\frac{\partial\hat{u}_{z}}{\partial t}+\frac{\pi
k}{1+\lambda}\left(\Gamma+4\pi^{2}k^{2}\right)\hat{u}_{z}$ $\displaystyle=$
$\displaystyle\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\frac{D}{(1+\lambda)\ln\kappa}\frac{\exp(2\pi
kr_{z_{0}})}{4\pi^{2}k^{2}}\big{[}2(1-\pi kr_{z_{0}})\left(\cosh(\pi
k)-1\right)-\pi k\sinh(\pi k)\big{]}$ (17)
$\displaystyle+\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)p_{z}\frac{\delta(k)}{\pi
k}.$

The expression for the vertical component of the swimmer translation velocity
$V_{z}^{T}$ is:

$\displaystyle V_{z}^{T}$ $\displaystyle=$
$\displaystyle\left(\frac{2\pi}{1+\lambda}\right)\left(\frac{\kappa_{\beta}}{\eta_{1}L^{2}V_{s}}\right)\int_{0}^{\infty}\mathrm{d}k\;k\left(\Gamma+4\pi^{2}k^{2}\right)\exp(2\pi
kr_{z_{0}})\left[\pi k\cosh(\pi k)-2(1-\pi kr_{z_{0}})\sinh(\pi
k)\right]\hat{u}_{z}$ (18)
$\displaystyle+\frac{D}{4\pi\ln\kappa}\left(\frac{1-\lambda}{1+\lambda}\right)\Bigg{[}\frac{1}{2-40r_{z_{0}}^{2}+128r_{z_{0}}^{4}}\bigg{(}12r_{z_{0}}^{2}-\left(5-4r_{z_{0}}-100r_{z_{0}}^{2}+80r_{z_{0}}^{3}+320r_{z_{0}}^{4}\right.$
$\displaystyle\left.\qquad-256r_{z_{0}}^{5}\right)\ln[-2+4r_{z_{0}}]+\left(5-8r_{z_{0}}-100r_{z_{0}}^{2}+160r_{z_{0}}^{3}+320r_{z_{0}}^{4}-512r_{z_{0}}^{5}\right)\ln[-1+4r_{z_{0}}]$
$\displaystyle\qquad+\left(5+8r_{z_{0}}-100r_{z_{0}}^{2}-160r_{z_{0}}^{3}+320r_{z_{0}}^{4}+512r_{z_{0}}^{5}\right)\ln[1+4r_{z_{0}}]$
$\displaystyle\qquad+\left(-5-4r_{z_{0}}+100r_{z_{0}}^{2}+80r_{z_{0}}^{3}-320r_{z_{0}}^{4}-256r_{z_{0}}^{5}\right)\ln[2+4r_{z_{0}}]$
$\displaystyle\qquad-24r_{z_{0}}^{2}+3(1-20r_{z_{0}}^{2}+64r_{z_{0}}^{4})\ln\left[\frac{4-16r_{z_{0}}^{2}}{1-16r_{z_{0}}^{2}}\right]\bigg{)}+\frac{1}{2(4r_{z_{0}}^{2}-1)}\Bigg{]};\left|r_{z_{0}}\right|>\frac{1}{2}.$

Here, the interface deformation and the swimmer translation exhibit a radial
symmetry in the plane of the interface, which enables a convenient analytical
formulation of the instantaneous part of $V_{z}^{T}$ as given in Eq. (18).
Such a radial symmetry is specific to swimmers orthogonal to the interface,
and is a consequence of the linearity of the problem. Note that we have not
incorporated the self-swimming term in $V_{z}^{T}$, which would involve
addition of a constant factor $p_{z}$ to Eq. 18. Therefore, the discussion
here pertains to shakers - swimmers that disturb the fluid with a flow field
$V_{s}/\ln\kappa$ but do not self-propel Stenhammar et al. (2017).

We now analyze the interface deformation created by the swimming motion. In
Fig. 7, we plot the interface deformation due to a pusher as a function of the
radial distance along the plane of the interface, $r_{\|}$, for $\Gamma=0.1$,
1 and 10. As defined previously, $r_{\|}=0$ is the swimmer location on the
interface. For large $r_{\|}$ the deformation $u_{z}$ is independent of
$\Gamma$, similar to that of swimmers parallel to the interface (see the
collapse in the inset of Fig. 7 for $r_{\|}\gg 1$). This far-field scaling is
$O\left(r_{\|}\right)^{-3}$, and can be identified by solving Eq. (17)
analytically, applying the approach used to obtain Eq. (14) in Sec. III.1,
yielding:

$\left.u_{z}\right|_{r_{\|}\gg
1}\approx\frac{1}{2\pi\Gamma}\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\frac{D}{\ln\kappa}\bigg{[}\frac{r_{z_{0}}^{2}}{4(r_{\|}^{2}+r_{z_{0}}^{2})^{\frac{3}{2}}}+\frac{r_{z_{0}}(1+\lambda)\Big{(}\Gamma
t-2r_{z_{0}}(1+\lambda)\Big{)}}{\Big{(}4r_{\|}^{2}(1+\lambda)^{2}+(\Gamma
t-2r_{z_{0}}(1+\lambda))^{2}\Big{)}^{\frac{3}{2}}}\bigg{]}.$ (19)

The above approximation also holds for $r_{z_{0}}\gg 1$ at $r_{\|}=O(1)$, and
is plotted in the inset of Fig. 7b.

However, the interface deformation for this swimmer configuration exhibits an
opposite trend to that of swimmers oriented parallel to the interface; namely,
for a pusher the deformation occurs away from the swimmer (into fluid region
2) at short separations, whereas, its positive for $r_{\|}\gg 1$, as also
noted previously in Lee et al. (2008) for point force-dipoles and Shaik and
Ardekani (2017) for spherical squirmers. This is expected as the disturbance
field generated by a pusher exits along the ends of its long-side and is
thereby pushing onto the interface when it is oriented orthogonal to it. Given
that the interface deformation also contributes to the evolution of the
vertical translation, this opposing trend will have important consequences for
the nature of the migration of swimmers oriented orthogonal to the interface,
as will be discussed below.

(a) $r_{z_{0}}(0)=1$

(b) $r_{z_{0}}(0)=5$

Figure 7: The interface deformation, $u_{z}$, due to a pusher plotted as a
function of the radial distance along the interface $r_{\|}$ for $\Gamma=$
0.1, 1 and 10 at time $t=2$ (a) for $r_{z_{0}}(0)=1$ and (b) $r_{z_{0}}(0)=5$.
The dotted lines in the inset represent the far-field $O(r_{\|})^{-3}$
scaling. Legends under $\circ$ refer to $u_{z}$ obtained from solving the pair
of equations (17), (18), and those under $\diamond$ from solving (20), (21).
The ‘Far-field’ in the inset of (b) refers to $u_{z}$ from Eq. (19) for
$\Gamma=1$.

The radial symmetry in the plane of the interface enables simplifications to
both $u_{z}$ and $V_{z}^{T}$ for swimmers positioned far away from the
interface ($r_{z_{0}}\gtrsim O(1)$). Following the discussion in Sec. III.1
for swimmers oriented parallel to the interface, we rescale the wavevector
$\bar{k}=k\bar{r}_{z_{0}}$ and time $\bar{t}=t/\bar{r}_{z_{0}}$ in Eq. 17 to
find

$\frac{\partial\hat{u}_{z}}{\partial\bar{t}}+\frac{\pi\Gamma\bar{k}}{1+\lambda}\hat{u}_{z}=\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\frac{\pi
Dr_{z_{0}}}{4(1+\lambda)\ln\kappa}\bar{k}\exp\left(2\pi\bar{k}\frac{r_{z_{0}}}{\bar{r}_{z_{0}}}\right),$
(20)

with the corresponding expression for vertical component of the translation
velocity becoming

$V_{z}^{T}=-\left(\frac{2\pi^{2}\Gamma}{(1+\lambda)\bar{r}_{z_{0}}^{3}}\right)\left(\frac{\kappa_{\beta}}{\eta_{1}L^{2}V_{s}}\right)\int_{0}^{\infty}\mathrm{d}\bar{k}\;\bar{k}^{2}\exp\left(2\pi\bar{k}\frac{r_{z_{0}}}{\bar{r}_{z_{0}}}\right)\hat{u}_{z}+\frac{3D}{128\pi\ln\kappa}\left(\frac{1-\lambda}{1+\lambda}\right)\left(\frac{1}{r_{z_{0}}^{2}}\right).$
(21)

In the insets of Fig. 7, we plot $|u_{z}|$ using both the full equations (17)
and (18), and the approximate equations (20) and (21). Similar to the
discussion in Sec. III.1, for $r_{z_{0}}(0)=1$ in Fig. 7a, the agreement of
both pair of equations is excellent for $\Gamma=10$, but the deviation of the
approximate solution is large for $\Gamma=1$ and 0.1. For $r_{z_{0}}=5$
however, this holds for $\Gamma<O(1)$, as seen in the inset of Fig. 7b.

(a) $\lambda=0.5$

(b) $\lambda=1$

(c) $\lambda=1.5$

Figure 8: The vertical swimmer translation velocity $V_{z}^{T}$ of pushers
plotted as a function of the initial distance of the swimmer from the
interface $r_{z_{0}}$ at different time instants for the viscosity ratios (a)
$\lambda=0.5$, (b) $\lambda=1$ and (c) $\lambda=1.5$. The plot is to be
interpreted as follows: at each $r_{z_{0}}$, the vertical black dotted lines
trace the time evolution of $V_{z}^{T}$ as indicated by the black arrow. In
all the plots $\kappa=10$, $\Gamma=1$ and
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.

In Fig. 8, we plot the time trace of the vertical component of the swimmer
translation velocity $V_{z}^{T}$ of a pusher at different initial swimmer
distances from the interface $r_{z_{0}}(0)$, for $\lambda=$ 0.5, 1, and 1.5.
Similar to swimmers parallel to the interface, for $\lambda<1$, both repulsive
and attractive regimes are exhibited by $V_{z}^{T}$ as a function of time.
Here, the deformation term is repulsive owing to the opposite character of the
interface deformation, as discussed below Eq. (18). The instantaneous term in
$V_{z}^{T}$ is attractive (repulsive) when $\lambda<1$ $(\lambda>1)$.
Therefore, when $\lambda<1$, for short times, $V_{z}^{T}$ is attractive
(negative), and is repulsive (positive) at long times. In turn, the relative
swimmer trajectory $r_{z_{0}}(t)-r_{z_{0}}(0)$ shown in Fig. 9 is the
qualitative mirror image about $r_{z_{0}}(t)-r_{z_{0}}(0)=0$, to that of a
swimmer parallel to the interface in Fig. 5.

The transition time from attraction to repulsion, $t_{\text{transition}}$,
exhibits a character similar to that of swimmers parallel to the interface,
and is expected from the similar scaling’s of Eq. (15) and Eq. (20). Namely,
for $r_{z_{0}}\sim O(1)$, $t_{\text{transition}}\sim O(1)$, whereas for large
$r_{z_{0}}$, $t_{\text{transition}}\sim O(r_{z_{0}})$. However, the magnitude
of $V_{z}^{T}$ is larger for swimmers oriented orthogonal to the interface
compared to swimmers oriented parallel to the interface. At $r_{z_{0}}(0)=1$,
for a swimmer parallel to the interface $V_{z}^{T}$ changes from $+0.5\times
10^{-3}$ to about $-10^{-3}$ in a time interval $t\in[0,4]$, as can be seen in
Fig. 4a (see the dotted trace of the first vertical line). In contrast, for
swimmers oriented orthogonal to the interface, $V_{z}^{T}$ transitions from
about $-10^{-3}$ to about $+2.5\times 10^{-3}$ in the same time interval, as
shown in Fig. 8a. Since the velocity at $t=0$, is that prescribed by the
instantaneous component, and that at large times is dominated by the
deformation, such character shows that both terms of $V_{z}^{T}$ are larger
for swimmers oriented orthogonal to the interface. Therefore, in a given time
interval, the relative swimmer vertical translation of an orthogonally
oriented swimmer as shown in Fig. 9 is larger than that of a swimmer parallel
to the interface (see Fig. 5). In other words, the coupled hydrodynamics more
strongly repels a pusher oriented orthogonal to the interface than it attracts
a pusher oriented parallel to the interface.

Figure 9: The relative vertical swimmer trajectory $r_{z_{0}}(t)-r_{z_{0}}(0)$
of pushers plotted as a function of time, for the viscosity ratios $\lambda=$
0.5, 1 and 1.5. In all the plots $\kappa=10$, $\Gamma=1$ and
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.

The relative importance of surface tension to bending stress on the swimmer
motion (not shown), is similar to that of swimmers parallel to the interface,
as discussed in Sec. III.1 and as is evident in Eqs. (17) and (18). Increasing
the surface tension results in a more rapid translation, and for $\lambda<1$
this implies a reduction in the time spent being attracted to the interface.
Thus, this monotonic response in $\Gamma$ implies that swimmers farther away
from an interface with larger $\Gamma$ can exhibit migration similar to those
closer to an interface with lower $\Gamma$. This is qualitatively similar in
character to Fig. 6b, but again of an opposite trend.

###  III.3  Microswimmers swimming arbitrarily oriented to the interface

We now generalize the coupled hydrodynamics to arbitrary swimmer orientations.
As explained in Appendix B, the derivation procedure for any non-parallel
swimmer orientation remains the same, and the changes to the fluid equations
due to a non-trivial in-plane orientation component are discussed in Appendix
D. Unlike the previous cases, here the vertical component of the rotation rate
is non-zero; $\dot{p}_{z}\neq 0$. For brevity we omit the large equations for
$\hat{u}_{z}$, $V_{z}^{T}$ and $\dot{p}_{z}$, the girth of latter two
associated with the contributions from all three orientation components. In
the following, we only focus on the role of the coupled hydrodynamics on the
swimmer vertical translation and rotation. We study how the change in the
swimmer orientation affects the nature of the swimmer migration, and the final
swimmer orientation given an arbitrary initial orientation. As in Sec. III.2,
here too we do not account for the self-swimming term in $V_{z}^{T}$, which is
instantaneously shifted by constant factor $p_{z}$.

As in the case of swimmer translation, the rotation rate has two parts; one
from the time dependent interface deformation and the other from the
instantaneous Stokes flow field. However, because the time dependent interface
deformation term in $\dot{p}_{z}$ is smaller than the instantaneous term, the
rotation rate is dominated by the latter. This quasi-steady behavior allows us
to construct a phase portrait in orientation space, as shown in Fig. 10a
Strogatz (2015). For pushers, only trajectories in the first and third
quadrants are admissible and hence any initial orientation $p_{z}\neq 0$
approaches $p_{z}=\pm 1$ depending on whether its directed away from
($p_{z}\to 1$ or $\theta\to 0$) or towards the interface ($p_{z}\to-1$ or
$\theta\to 180^{\circ}$). Now, from the analysis of swimmers parallel to the
interface in Sec. III.1, we know that $\dot{p}_{z}=0$, for $p_{z}=0$.
Therefore, $p_{z}=0$ is an unstable fixed point for pusher, and hence a stable
fixed point for pullers. Similar rotational preferences have been reported by
Shaik and Ardekani (2017) for spherical squirmers near a deformable interface.

(a)

(b)

Figure 10: (a) The rotation rate $\dot{p}_{z}$ plotted as a function of
$p_{z}(\equiv\cos\theta)$. (b) The vertical component of the translation
velocity $V_{z}^{T}(t)$ of pushers plotted as a function of time, for a range
of initial swimmer orientation relative to the $r_{z}$ axis $\theta(0)=$
$30^{\circ}$, $45^{\circ}$, $60^{\circ}$, $75^{\circ}$, $85^{\circ}$ and
$89.5^{\circ}$. In the inset, $V_{z}^{T}$ is shifted by the time it takes the
$\theta(0)=89.5^{\circ}$ curve to attain the corresponding $\theta(0)$. In all
the plots $\lambda=$ 0.5, $\kappa=10$, $\Gamma=1$ and
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.

In Fig. 10b the time dependence of the vertical component of the swimmer
translation velocity $V_{z}^{T}$ is shown for pushers starting from different
initial orientations $\theta(0)$. When a swimmer is nearly aligned with the
interface ($\theta(0)=89.5^{\circ}$), $V_{z}^{T}$ decreases until $t\sim
O(1)$, owing to the dominant contribution from the parallel configuration, as
anticipated from Sec. III.1. On the other hand, $V_{z}^{T}$ grows at large
times, owing to the dominant contribution of the orthogonal configuration as
$\theta(t)\to 0$. Therefore, the translation velocity $V_{z}^{T}$ of swimmers
nearly aligned with the interface goes through a minimum, which is negative
for nearly parallel swimmers (see curves of $\theta(0)=$ $85^{\circ}$ and
$89.5^{\circ}$ ) and positive for those that start at a smaller initial
orientation, say $\theta(0)\lesssim 75^{\circ}$ and $\gtrsim 45^{\circ}$. At
$\theta(0)\lesssim 45^{\circ}$, $V_{z}^{T}$ only increases in time with a
$\theta(0)$ dependent $t=0$ intercept. Therefore, the early time trajectories
are qualitatively and quantitatively dependent upon $\theta(0)$ but all of the
curves collapse at late times.

Owing to the sensitive dependence of $V_{z}^{T}(t)$ on $\theta(0)$ at short
times, the swimmer migration has a rich set of trajectories up to times of
$O(1)$. In Fig. 11a, we plot the relative vertical swimmer trajectories
$r_{z_{0}}(t)-r_{z_{0}}(0)$ for the same values of $\theta(0)$. The inset
shows short time oscillations for $\theta(0)=$ $85^{\circ}$ and
$89.5^{\circ}$, the latter of which briefly becomes weakly attractive, whereas
the curvature gradually changes as $\theta(0)$ decreases. At long times the
slopes asymptote for all $\theta(0)$ and the overall displacement increases as
$\theta(0)$ decreases.

(a)

(b)

Figure 11: (a) The relative vertical swimmer trajectory
$r_{z_{0}}(t)-r_{z_{0}}(0)$ of pushers plotted as a function of time for a
range of initial swimmer orientation relative to the $r_{z}$ axis $\theta(0)=$
$30^{\circ}$, $45^{\circ}$, $60^{\circ}$, $75^{\circ}$, $85^{\circ}$ and
$89.5^{\circ}$. The viscosity ratio is $\lambda=0.5$, and the inset is a
zoomed version of the abscissa for $t\in[0,1]$. (b) The relative vertical
swimmer trajectory $r_{z_{0}}(t)-r_{z_{0}}(0)$ of pushers plotted as a
function of time, for the viscosity ratios $\lambda=$ 0.5, 1 and 1.5. In the
inset the corresponding vertical swimmer translation velocity $V_{z}^{T}$ of
pushers plotted as a function of time. In all the plots $\kappa=10$,
$\Gamma=1$ and $\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.

In Fig. 11b, we plot the relative vertical swimmer trajectory
$r_{z_{0}}(t)-r_{z_{0}}(0)$, for three values of the viscosity ratio
$\lambda=$ 0.5, 1 and 1.5 for a pusher that starts off initially almost
parallel to the interface ($\theta=85^{\circ}$). Although the three
trajectories are qualitatively similar, here a swimmer exhibits regimes of
attraction and repulsion for $\lambda\geq 1$. The underlying principle for
such a response can be explained as follows. For swimmers that start off
nearly aligned with the interface, at short times $V_{z}^{T}$ is positive for
$\lambda<1$, zero for $\lambda=1$ (no contribution from instantaneous
component), and negative for $\lambda>1$. As time evolves the swimmer rotates
to an orthogonal configuration and $V_{z}^{T}$ increases, independent of
$\lambda$. In turn, $V_{z}^{T}(t=0)$ of a pusher nearly aligned to the
interface positive, zero or negative if $\lambda<1$, $\lambda=1$ or
$\lambda\geq 1$ respectively. Therefore, for $t<O(1)$, a swimmer will always
have an attractive component for $\lambda\geq 1$, and
$r_{z_{0}}(t)-r_{z_{0}}(0)<0$.

Figure 12: The relative vertical swimmer trajectory
$r_{z_{0}}(t)-r_{z_{0}}(0)$ of pushers and pullers plotted as a function of
time, initial orientations $\theta(0)=30^{\circ}$ and $60^{\circ}$. In the
insets, the corresponding $V_{z}^{T}(t)$ plotted as a function of time. In all
the plots $\lambda=0.5$, $\kappa=10$, $\Gamma=1$ and
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1$.

We now understand that hydrodynamic interactions with a deformable interface
lead to pushers reorienting into an orthogonal configuration, pullers reorient
into a parallel configuration. Therefore, unlike their response to a rigid
boundary, both swimmer types seek a configuration that drives them away from
the interface. As shown in Fig. 12, pushers and pullers that start at the same
location will eventually segregate spatially, independent of their initial
orientation but to a degree that increases as $\theta(0)$ decreases. Note that
$\bm{p}=0$ is a stable configuration for pushers only the absence of an
imposed flow, which may be present in many microswimmer settings, as we
discuss in Sec. V.

##  IV  Microswimmers confined between a rigid boundary and an interface -
Swimming parallel to boundaries

Here, we generalize our approach to treat the coupled hydrodynamics of a
swimmer confined between a rigid boundary and an underlying deformable
interface (Fig. 13). We note that in this configuration the rotation rate
$\dot{\bm{p}}=0$, similar to a swimmer oriented parallel to a single
deformable interface (see Sec. III.1), or to a single rigid boundary Berke et
al. (2008). Therefore, the hydrodynamics is characterized by the interface
deformation $u_{z}$ and the vertical component of the translation velocity
$V_{z}^{T}$. It is convenient to choose a reference frame moving with the
swimmer but shifted to the plane of the undeformed interface, and hence the
definition of $r_{z_{0}}$ used here is the opposite of that used for a swimmer
near a single interface in Sec. III (see the direction of arrows in Figs. 1
and 13).

The equation for the Fourier transformed interface deformation $\hat{u}_{z}$
is

$\displaystyle\frac{\text{d}\hat{u}_{z}}{\text{d}t}$ $\displaystyle+$
$\displaystyle\frac{\pi
k}{\mathcal{E}}\left(4\pi^{2}k^{2}+\Gamma\right)\big{[}1+\lambda-(1-\lambda)\exp(-8\pi
kH)-2(\lambda+4\pi kH+8\pi^{2}k^{2}H^{2}\lambda)\exp(-4\pi
kH)\big{]}\hat{u}_{z}$ (22)
$\displaystyle=\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\left(\frac{D}{\pi
k\ln\kappa}\right)\sin^{2}\left(\frac{\pi}{2}kp_{l}\right)\frac{1}{\mathcal{E}}\bigg{[}(1+\lambda)r_{z_{0}}\exp(-2\pi
kr_{z_{0}})-(1-\lambda)r_{z_{0}}\exp(-2\pi k(4H-r_{z_{0}}))$
$\displaystyle\quad+\;(1+\lambda)(-r_{z_{0}}-4\pi kHr_{z_{0}}+4\pi
kH^{2})\exp(-2\pi k(2H-r_{z_{0}})$
$\displaystyle\quad+\;(1-\lambda)(r_{z_{0}}-4\pi kHr_{z_{0}}+4\pi
kH^{2})\exp(-2\pi k(2H+r_{z_{0}})\bigg{]},$

and the expression for $V_{z}^{T}$ is

$\displaystyle V_{z}^{T}\equiv\frac{\mathrm{d}r_{z_{0}}}{\mathrm{d}t}$
$\displaystyle=$
$\displaystyle-\frac{2D}{\ln\kappa}\int\mathrm{d}\bm{k}\frac{\sin(\pi
kp_{l})}{\pi
kp_{l}}\frac{1}{\mathcal{E}}\sin^{2}\left(\frac{\pi}{2}kp_{l}\right)\bigg{[}(1-\lambda^{2})r_{z_{0}}^{2}\big{[}\exp(-4\pi
kr_{z_{0}})+\exp(-4\pi k(2H-r_{z_{0}}))\big{]}$
$\displaystyle\qquad+(H-r_{z_{0}})^{2}\big{[}(1+\lambda)^{2}\exp(-4\pi
k(H-r_{z_{0}}))+(1-\lambda)^{2}\exp(-4\pi k(H+r_{z_{0}}))\big{]}$
$\displaystyle\qquad+2H(1-\lambda^{2})(H-2r_{z_{0}})\exp(-4\pi kH)\bigg{]}$
(23)
$\displaystyle+\left(\frac{\kappa_{\beta}}{\eta_{1}V_{s}L^{2}}\right)\int\mathrm{d}\bm{k}\frac{\sin(\pi
kp_{l})}{\pi kp_{l}}\pi
k\left(4\pi^{2}k^{2}+\Gamma\right)\frac{\hat{u}_{z}}{\mathcal{E}}\bigg{[}-(1+\lambda)(1+2\pi
kr_{z_{0}})\exp(-2\pi kr_{z_{0}})$ $\displaystyle\qquad+(1-\lambda)(1-2\pi
kr_{z_{0}})\exp(-2\pi k(4H-r_{z_{0}}))$
$\displaystyle\qquad-(1-\lambda)\big{[}(1-4\pi kH)(1+2\pi
kr_{z_{0}})+8\pi^{2}k^{2}H^{2}\big{]}\exp(-2\pi k(2H+r_{z_{0}}))$
$\displaystyle\qquad+(1+\lambda)\big{[}(1+4\pi kH)(1-2\pi
kr_{z_{0}})+8\pi^{2}k^{2}H^{2}\big{]}\exp(-2\pi k(2H-r_{z_{0}}))\bigg{]},$

where, $\mathcal{E}$ in both Eqs. (22) and (23) is

$\displaystyle\mathcal{E}=(1+\lambda)^{2}+(1-\lambda)^{2}\exp(-8\pi
kH)-2(1-\lambda)^{2}(1+8\pi^{2}k^{2}H^{2})\exp(-4\pi kH).$ (24)

We note that the additional $H$-dependent terms in Eq. (22) for the interface
deformation are at least exponentially smaller by a factor $\exp(-2\pi
k(2H-r_{z_{0}})$. Unless $H\to r_{z_{0}}$ (extreme confinement), the
contributions from these terms remains small, apart from there being an
$H$-dependent amplitude, as the terms with the smaller exponents are
proportional to $1-\lambda$. Hence, varying $\lambda$ about unity does not
yield a qualitatively different deformation than that of a swimmer parallel a
single interface. Therefore, in the following we focus only on the change in
the swimmer translation effected by the coupled hydrodynamics between the
swimmer and the boundaries.

Figure 13: A sectional schematic representation of a slender swimmer
translating with speed $V_{s}$ along its director vector $\bm{p}$ in a channel
with a deformable interface below and a rigid boundary on top. The distance
between the plane of the undeformed interface and the rigid boundary is $H$.
The circles along the axial length of the fore-aft symmetric swimmer represent
a line distribution of stokeslets characterizing the head and tail for
pushers. The disturbance flow field generated by the swimming motion deforms
the interface, and $u_{z}$ is the interface deformation (solid gray line)
relative to its initially flat undeformed (dotted line) state $z_{0}$.

(a)

(b)

Figure 14: (a) The relative swimmer trajectory $r_{z_{0}}(t)-r_{z_{0}}(0)$ of
a pusher plotted as a function of time for different distance $H$ between the
interface and the rigid boundary: $H=$ 2, 2.2, 2.4, 10, when the viscosity
ratio is $\lambda=0.5$. (b) $r_{z_{0}}(t)-r_{z_{0}}(0)$ of a pusher plotted as
a function of time for $\lambda=$ 0.5, 1 and 1.5, and fixing $H$ = 2.3. The
insets contains the corresponding vertical component of the translation
velocity $V_{z}^{T}\equiv\mathrm{d}r_{z_{0}}/\mathrm{d}t$. The “repulsion” and
“attraction” are to be interpreted relative to the interface. In both figures,
the initial swimmer distance to the interface is $\left|r_{z_{0}}\right|=1$
and the ratio of surface tension to bending stress is $\Gamma=$ 1\. The
attraction and repulsion are specified relative to the deformable interface.

In Fig. 14a we plot the relative swimmer trajectory
$r_{z_{0}}(t)-r_{z_{0}}(0)$ as a function of the distance between the
boundaries $H$, for $r_{z_{0}}(0)=1$ and $\lambda=0.5$. In the limit $H\gg 1$,
we recover the trajectory corresponding to that of a swimmer in the vicinity
of a single interface as discussed in Sec. III.1. Upon increased confinement
by reducing $H$ and increasing confinement, the residence time of the swimmer
in the repulsive state increases. When the swimmer begins at the midpoint, it
always moves towards the rigid boundary. This is also shown in the inset of
Fig. 14a, where for $H=2$, $V_{z}^{T}$ does not undergo a zero-crossing in
finite time. Physically, a rigid boundary provides an infinite resistance to
both bending and shear, and thus facilitates a stronger attraction relative to
an interface that only supports finite bending Daddi-Moussa-Ider et al.
(2019). In other words, the flow field induced by the rigid surface dominates
that of the deformable interface.

In Fig. 14b we plot the relative swimmer trajectory
$r_{z_{0}}(t)-r_{z_{0}}(0)$ for three different viscosity ratios $\lambda$ =
0.5, 1, 1.5, with the distance between the interface and the rigid boundary
fixed to $H=2.3$. Interestingly, unlike a swimmer in the vicinity of a
standalone interface, the swimmer experiences a transient attraction and
repulsion for all values of $\lambda$. In Sec. III.1 we explained that for
$\lambda<1$ the swimmer translation had both repulsive and attractive
component due to the difference in the instantaneous contribution, dominant at
short times, and the interface deformation dependent contributions of
$V_{z}^{T}$, dominant at long times (see Eq. (13)). Here, however, the rigid
boundary creates repulsion, through an $H$-dependent positive shift of the
instantaneous component of $V_{z}^{T}$, the largest contribution being
proportional to $(1+\lambda)^{2}$, and thus positive for any $\lambda$. As the
confinement increases, the contribution of this term increases, and eventually
for $H\gtrsim 2r_{z_{0}}$, the swimmer only moves towards the rigid boundary.
Although not shown, the dependence on the ratio of surface tension to bending
stress, $\Gamma$, remains the same as that for a swimmer near a single
deformable interface, and seen in Eqs. (22) and (23).

##  V  Summary and Conclusion

We have investigated the coupled hydrodynamics between finite sized orientable
swimmers and a deformable interface. By treating the swimmers as slender
bodies, we have gone beyond a far-field picture, solving the hydrodynamics
appropriate at $O(1)$ swimmer lengths from the interface. Moreover, we predict
the hydrodynamics for a robust range of the ratio of the viscous stress to
bending stress, $\eta_{1}V_{s}L^{2}/\kappa_{\beta}$, and the capillary number,
$\eta_{1}V_{s}/\gamma$. Our analysis reveals that the swimmer orientation
plays a crucial role in the response of both the swimmer and deformable
interface. Because parallel and perpendicular orientations exhibit differing
overall trajectories, when swimmers take an arbitrary orientations we find a
rich dynamical behavior. Importantly, given that an arbitrarily oriented
pusher (puller) preferentially rotates to a perpendicular (parallel)
orientation, the migration pattern for both swimmer types remains the same.
However, given the difference in the migration speed in the two orientations,
a pusher and puller that start from the same location with any orientation
will spatially segregate.

When confined to swimming between a rigid and a deformable boundary, pushers
closer to the latter interface experience repulsion that is extended to time
scales longer than the elastic response time. Such confinement, between rigid
and soft boundaries, provide new controls on swimmer migration. For example,
the different rotational preferences near deformable versus rigid boundaries
provides one such confinement-dependent control of swimmers.

In general, the dynamics of swimmer translation is extremely sensitive to the
distance from the interface. Even an $O(1)$ increase in this distance results
in a large slowdown in the migration towards or away from the interface. The
interfacial properties play a crucial role in migration, second only to the
distance from the interface. Therefore, it is possible that two scenarios in
which swimmers have different distances from the interface display similar
migration patterns. Owing to the strong sensitivity of swimmer migration to
distance from the interface, even a small change in the latter would require
an $O(1)$ change in the interface properties to obtain a nearly matching
migration pattern. This is particularly relevant from a practical standpoint,
because interfacial properties are not necessarily spatially homogeneous. In
consequence, swimmers in a suspension at different distances from the
interface can, in principle, redistribute in a spatially similar manner.

We find that the viscosity ratio plays a controlling role in whether a swimmer
is attracted or repelled from the interface. Indeed, viscosity variations may
be dictate the migration of planktonic biota in marine ecology, by modifying
both swimming and the nutrient dispersion Guasto et al. (2012); Finnegan et
al. (2015); Simpson (2021). On the other hand, whereas interfacial deformation
is clearly controlled by the relative importance of bending elasticity and
surface tension, it also depends on the swimmer orientation and distance from
the interface. However, independent of the swimmer orientation, the
interfacial deformation in the interfacial plane extends well beyond distances
of order the swimmer size $L$. Moreover, the far-field scaling of the
deformation in the interfacial plane remains $O\left(r_{\|}\right)^{-3}$ for
swimmers oriented parallel or perpendicular to the interface. Even in the
dilute limit, such far-field scaling can have important consequences in
analyzing the deformation due to a suspension of microswimmers near deformable
interfaces. This is because inclusion of the collective contribution from the
swimmers upon the interface involves an integral over the domain volume, which
decays as $O\left(r_{\|}\right)^{-3}$, implying the potential for long-ranged
radial contributions.

It is important to note that time scales short relative to the elastic
response time, $t_{c1}\left(=\eta_{1}L^{3}/\kappa_{\beta}\right)$, are still
usually very long compared to the time scales of swimming. For typical soft
interfaces, the bending modulus $\kappa_{\beta}\approx 10^{-19}$ $J$ Fygenson
et al. (1997); Freund (2014), so if one considers a viscosity approximately
that of water, $\eta_{1}\approx 10^{-3}$ $Ns/m^{2}$, and swimmers of length
$L\approx 10$ $\mu m$ (cell body+flagellar bundle combined), then
$t_{c1}\approx 10$ $s$. In general, microscopic swimmers also respond to
external cues, such as an externally imposed flow Marcos et al. (2012);
Rusconi et al. (2014); Costanzo et al. (2014); Barry et al. (2015); Mathijssen
et al. (2016); Wioland et al. (2016); Lauga (2016); Vennamneni et al. (2020)
or some form of chemical actuation Mesibov and Adler (1972); Pedley et al.
(1988); Vincent and Hill (1996); Cisneros et al. (2010); Painter and Hillen
(2011); Subramanian et al. (2011); Ezhilan et al. (2012); Desai and Ardekani
(2018), that can influence their migration on such long time scales.
Typically, even in the absence of such an external driving force, the
reorientations take seconds, as most microscopic swimmers have an inherent
mechanism allowing them to change their orientation Berg (2008); Koch and
Subramanian (2011); Lauga (2016); Saintillan (2018). We note that when
$\lambda<1$, a swimmer aligned with the interface might only sense either a
repulsion (pusher) or an attraction (puller), before reorienting itself.
Therefore, the nature of the short time dynamics shown in Fig. 5 implies the
possibility of pushers migrating away from the interface on average. A similar
argument can be made for other swimmer orientations.

Although we have considered the motion of a single swimmer near a deformable
interface and confined between a deformable interface and a rigid boundary,
one can readily extend this study to a dilute suspension of non-interacting
swimmers. This requires a more careful consideration of the single-swimmer
statistics developed here, as the expressions for the instantaneous field
variables depend on the swimmer configuration. Moreover, given the complex
dynamics involved, deriving a position-orientation space swimmer distribution
function is challenging. Finally, a host of new dynamical processes will be
revealed when extending this approach beyond the linear regime.

###### Acknowledgements.

The authors gratefully acknowledge support from the Swedish Research Council,
under grant no. 638-2013-9243. Nordita is partially supported by Nordforsk.

##  Appendix A Two-dimensional Fourier transform summary - Swimmers parallel
to interface

Here we outline the two-dimensional Fourier transform technique used in the
main text to solve the governing equations and boundary conditions for a
slender swimmer translating parallel to a deformable interface. For a more
detailed description of the method, the reader is directed to Bickel (2007)
and Daddi-Moussa-Ider et al. (2016).

###  A.1 Formulating the governing equations in Fourier space

For a swimmer translating parallel to the interface, the Stokes equations and
the continuity equation in fluid region 1 given in Eqs. (6a) and (6b) are:

$\displaystyle-2\pi{i\mkern
1.0mu}k_{x}\hat{P}_{1}+\frac{\partial^{2}\hat{v}_{1x}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1x}$ $\displaystyle=$
$\displaystyle\frac{2D\,p_{x}}{\pi{i\mkern
1.0mu}\bm{k}\cdot\bm{p}\ln\kappa}\delta(r_{z})\sin^{2}\left(\frac{\pi}{2}\bm{k}\cdot\bm{p}\right),$
(25a) $\displaystyle-2\pi{i\mkern
1.0mu}k_{y}\hat{P}_{1}+\frac{\partial^{2}\hat{v}_{1y}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1y}$ $\displaystyle=$
$\displaystyle\frac{2D\,p_{y}}{\pi{i\mkern
1.0mu}\bm{k}\cdot\bm{p}\ln\kappa}\delta(r_{z})\sin^{2}\left(\frac{\pi}{2}\bm{k}\cdot\bm{p}\right),$
(25b) $\displaystyle-\frac{\partial\hat{P}_{1}}{\partial
r_{z}}+\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1z}$ $\displaystyle=$ $\displaystyle
0,\;\;\text{and}$ (25c) $\displaystyle 2\pi{i\mkern
1.0mu}(k_{x}\hat{v}_{1x}+k_{y}\hat{v}_{1y})+\frac{\partial\hat{v}_{1z}}{\partial
r_{z}}$ $\displaystyle=$ $\displaystyle 0.$ (25d)

In fluid region 2, the right-hand side of the corresponding Stokes equations
for $\hat{\bm{v}}_{2}$, $P_{2}$ vanishes. It is convenient to solve this
system of equations in Fourier space in a coordinate system aligned with the
wavevector $\bm{k}$ Bickel (2007); Daddi-Moussa-Ider et al. (2016). We define
the longitudinal and transverse coordinate systems $\hat{\bm{l}}$ and
$\hat{\bm{t}}$, as

$\displaystyle\hat{\bm{l}}$ $\displaystyle=$
$\displaystyle\frac{k_{x}}{k}\hat{\bm{x}}+\frac{k_{y}}{k}\hat{\bm{y}},\;\;\text{and}$
(26a) $\displaystyle\hat{\bm{t}}$ $\displaystyle=$
$\displaystyle\frac{k_{y}}{k}\hat{\bm{x}}-\frac{k_{x}}{k}\hat{\bm{y}},$ (26b)

respectively. This choice of coordinates has a few advantages. First, it is
evident from Eq. (25d) that the continuity equation provides a direct relation
between the longitudinal and the $z$-component of the fluid velocity. Second,
as will be described below, the pressure gradient term drops out of the
$\hat{t}$-component of the Fourier transformed Stokes equation. Lastly, one
can obtain a single $4^{th}$-order differential equation for the $z$-component
of the fluid velocity, again without the pressure gradient term, thereby
decoupling the flow-field field variables. In the
($\hat{\bm{l}},\hat{\bm{t}},\bm{z}$) coordinates Eqs. (25a)-(25d) become

$\displaystyle-2\pi{i\mkern
1.0mu}k\hat{P}_{1}+\frac{\partial^{2}\hat{v}_{1l}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1l}$ $\displaystyle=$
$\displaystyle\frac{2D}{\pi{i\mkern
1.0mu}k\ln\kappa}\delta(r_{z})\sin^{2}\left(\frac{\pi}{2}kp_{l}\right),$ (27a)
$\displaystyle\frac{\partial^{2}\hat{v}_{1t}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1t}$ $\displaystyle=$
$\displaystyle\frac{2D\,p_{t}}{\pi{i\mkern
1.0mu}kp_{l}\ln\kappa}\delta(r_{z})\sin^{2}\left(\frac{\pi}{2}kp_{l}\right),$
(27b) $\displaystyle-\frac{\partial\hat{P}_{1}}{\partial
r_{z}}+\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1z}$ $\displaystyle=$ $\displaystyle
0,\;\;\text{and}$ (27c) $\displaystyle 2\pi{i\mkern
1.0mu}k\hat{v}_{1l}+\frac{\partial\hat{v}_{1z}}{\partial r_{z}}$
$\displaystyle=$ $\displaystyle 0,$ (27d)

where we note that the $z$-component is the same as Eq. (25c). In Eqs. (27a)
and (27b), $p_{l}\equiv\bm{p}\cdot\hat{\bm{l}}$ is the component of the
orientation vector parallel to the wavevector and
$p_{t}\equiv\bm{p}\cdot\hat{\bm{t}}$ is that transverse to the wavevector.
Using Eqs. (27a) and (27c) to eliminate the pressure term, and Eq. (27d) to
represent $\hat{v}_{1l}$ in terms of $\hat{v}_{1z}$, one arrives at the
following differential equation to solve for $\hat{v}_{1z}$:

$\frac{\partial^{4}\hat{v}_{1z}}{\partial
r_{z}^{4}}-8\pi^{2}k^{2}\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}+16\pi^{4}k^{4}\hat{v}_{1z}=-\frac{4\,D}{\ln\kappa}\delta^{\prime}(r_{z})\sin^{2}\left(\frac{\pi}{2}kp_{l}\right),$
(28)

where $\delta^{\prime}(r_{z})$ is the derivative of $\delta(r_{z})$. Similar
to the above derivation, in fluid region 2 we have the following set of
equations:

$\displaystyle\frac{\partial^{4}\hat{v}_{2z}}{\partial
r_{z}^{4}}-8\pi^{2}k^{2}\frac{\partial^{2}\hat{v}_{2z}}{\partial
r_{z}^{2}}+16\pi^{4}k^{4}\hat{v}_{2z}$ $\displaystyle=$ $\displaystyle 0,$
(29a) $\displaystyle\frac{\partial^{2}\hat{v}_{2t}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{2t}$ $\displaystyle=$ $\displaystyle
0,\;\;\text{and}$ (29b) $\displaystyle 2\pi{i\mkern
1.0mu}k\hat{v}_{2l}+\frac{\partial\hat{v}_{2z}}{\partial r_{z}}$
$\displaystyle=$ $\displaystyle 0.$ (29c)

Once $\hat{v}_{\alpha z}$ is determined, we can use the continuity equations
(27d) and (29c) to determine the longitudinal component $\hat{v}_{\alpha l}$,
where $\alpha\in[1,2]$ refers to the two fluid regions.

###  A.2 Formulating boundary conditions in Fourier space

The equations for $\hat{v}_{\alpha z}$ and $\hat{v}_{\alpha t}$ are solved
subject to the appropriate boundary conditions; the impenetrability and the
no-slip conditions at the interface, and the continuity of tangential stress
and the normal-stress jump across the interface. The Fourier transformed
velocity boundary conditions are:

$\displaystyle\left.\hat{v}_{1z}\right|_{r_{z_{0}}^{+}}$ $\displaystyle=$
$\displaystyle\left.\hat{v}_{2z}\right|_{r_{z_{0}}^{-}},$ (30a)
$\displaystyle\left.\hat{v}_{1t}\right|_{r_{z_{0}}^{+}}$ $\displaystyle=$
$\displaystyle\left.\hat{v}_{2t}\right|_{r_{z_{0}}^{-}}\;\;\text{and}$ (30b)
$\displaystyle\left.\frac{\partial\hat{v}_{1z}}{\partial
r_{z}}\right|_{r_{z_{0}}^{+}}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\hat{v}_{2z}}{\partial
r_{z}}\right|_{r_{z_{0}}^{-}},$ (30c)

where we have used the continuity equations (27d) and (29c) to arrive at Eq.
(30c) as an alternative to the no-slip boundary condition.

The continuity of tangential stress given by Eqs. (7c) and (7d) in the main
text can be simplified in the $\hat{\bm{l}}$-$\hat{\bm{t}}$ coordinates to
yield the following pair of equations for $\hat{v}_{\alpha z}$ and
$\hat{v}_{\alpha t}$ at the interface

$\displaystyle\left.\left(\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}+4\pi^{2}k^{2}\hat{v}_{1z}\right)\right|_{r_{z_{0}}^{+}}$
$\displaystyle=$
$\displaystyle\lambda\left.\left(\frac{\partial^{2}\hat{v}_{2z}}{\partial
r_{z}^{2}}+4\pi^{2}k^{2}\hat{v}_{2z}\right)\right|_{r_{z_{0}}^{-}}\;\;\text{and}$
(31a) $\displaystyle\left.\frac{\partial\hat{v}_{1t}}{\partial
r_{z}}\right|_{r_{z_{0}}^{+}}$ $\displaystyle=$
$\displaystyle\lambda\left.\frac{\partial\hat{v}_{2t}}{\partial
r_{z}}\right|_{r_{z_{0}}^{-}}.$ (31b)

The Fourier transformed normal stress boundary condition at the interface give
by Eq. (7e) is

$-\left(\left.\hat{P}_{1}\right|_{r_{z_{0}}^{+}}-\left.\hat{P}_{2}\right|_{r_{z_{0}}^{-}}\right)+2\left(\left.\frac{\partial\hat{v}_{1z}}{\partial
r_{z}}\right|_{r_{z_{0}}^{+}}-\lambda\left.\frac{\partial v_{2z}}{\partial
r_{z}}\right|_{r_{z_{0}}^{-}}\right)=4\pi^{2}k^{2}\frac{\gamma}{\eta_{1}V_{s}}\hat{u}_{z}+16\pi^{4}k^{4}\frac{\kappa_{\beta}L^{2}}{\eta_{1}V_{s}}\hat{u}_{z}.$
(32)

Eq. (32) can be further simplified by using Eq. (25c) (and its equivalent
equation for $\hat{\bm{v}}_{2}$ and $P_{2}$) along with the continuity
equations (27d) and (29c) to find:

$\displaystyle\left(\left.\frac{\partial^{3}\hat{v}_{1z}}{\partial
r_{z}^{3}}\right|_{r_{z_{0}}^{+}}-\lambda\left.\frac{\partial^{3}v_{2z}}{\partial
r_{z}^{3}}\right|_{r_{z_{0}}^{-}}\right)$ $\displaystyle+$ $\displaystyle
12\pi^{2}k^{2}\left(\left.\frac{\partial\hat{v}_{1z}}{\partial
r_{z}}\right|_{r_{z_{0}}^{+}}-\lambda\left.\frac{\partial v_{2z}}{\partial
r_{z}}\right|_{r_{z_{0}}^{-}}\right)=16\pi^{4}k^{4}\left(\frac{\kappa_{\beta}L^{2}}{\eta_{1}V_{s}}\right)\left(4\pi^{2}k^{2}+\Gamma\right)\hat{u}_{z},$
(33)

for the boundary condition of the normal component of the stress without the
pressure term; note that $\Gamma\equiv\gamma L^{2}/\kappa_{\beta}$ in Eq.
(33).

###  A.3 Solving for the velocity field in Fourier space

Here, we first solve for the transverse components of the Fourier transformed
velocity fields and then the normal and the longitudinal components. A general
solution for equations of the form of Eqs. (27b) and (29b) is:
$\hat{v}_{\alpha t}=A\exp(\pm 2\pi kr_{z})$ (e.g., Bickel, 2007; Bender and
Orszag, 2013). In the specific problem considered here, the disturbance flow
field must decay in the far field, and hence, the transverse velocity
components are

$\displaystyle\hat{v}_{1t}$ $\displaystyle=$
$\displaystyle\begin{cases}A_{1}\exp(-2\pi kr_{z});r_{z}>0\\\ A_{2}\exp(-2\pi
kr_{z})+A_{3}\exp(2\pi kr_{z});r_{z_{0}}<r_{z}<0,\;\;\text{and}\\\
\end{cases}$ (34a) $\displaystyle\hat{v}_{2t}$ $\displaystyle=$ $\displaystyle
A_{4}\exp(2\pi kr_{z});r_{z}<r_{z_{0}}.$ (34b)

Now, the right-hand side of Eq. (27b) is proportional to $\delta(r_{z})$, and
hence, we seek the Greens function of the differential equation. From Appendix
A.2, we have two boundary conditions for $\hat{v}_{\alpha t}$, namely, the no-
slip condition given by Eq . (30b) and the continuity of tangential stress
given by Eq. (31b), and four unknown constants $A_{1}$-$A_{4}$ to determine.
Therefore, we seek two conditions to uniquely determine the constants. To
extract additional boundary conditions, we use the continuity properties of
the Greens function; (a) about the location of forcing ($r_{z}=0$), and (b)
the first derivative with respect to $r_{z}$ exhibiting a finite-jump
discontinuity about $r_{z}=0$ Bickel (2007); Bender and Orszag (2013). These
additional conditions are:

$\displaystyle\left.\hat{v}_{1t}\right|_{r_{z}=0^{+}}$ $\displaystyle=$
$\displaystyle\left.\hat{v}_{1t}\right|_{r_{z}=0^{-}},\;\;\text{and}$ (35a)
$\displaystyle\left.\frac{\partial\hat{v}_{1t}}{\partial
r_{z}}\right|_{r_{z}=0^{+}}-\left.\frac{\partial\hat{v}_{1t}}{\partial
r_{z}}\right|_{r_{z}=0^{-}}$ $\displaystyle=$
$\displaystyle\frac{2D\,p_{t}}{\pi{i\mkern
1.0mu}kp_{l}\ln\kappa}\sin^{2}\left(\frac{\pi}{2}kp_{l}\right).$ (35b)

Using Eqs. (30)b, (31b) and (35), the four constants are uniquely determined,
and we obtain the following expressions for the transverse velocity components
in the two fluid region:

$\displaystyle\hat{v}_{1t}$ $\displaystyle=$ $\displaystyle\frac{{i\mkern
1.0mu}p_{t}}{2\pi^{2}k^{2}p_{l}\ln\kappa}\sin^{2}\left(\frac{\pi}{2}kp_{l}\right)\left[\exp(-2\pi
k|r_{z}|)+\left(\frac{1-\lambda}{1+\lambda}\right)\exp(4\pi
kr_{z_{0}})\exp(-2\pi kr_{z})\right]$ (36a)
$\displaystyle;r_{z}>r_{z_{0}},\;\;\text{and}$ $\displaystyle\hat{v}_{2t}$
$\displaystyle=$ $\displaystyle\frac{{i\mkern
1.0mu}p_{t}}{\pi^{2}k^{2}p_{l}(1+\lambda)\ln\kappa}\sin^{2}\left(\frac{\pi}{2}kp_{l}\right)\exp(2\pi
kr_{z});r_{z}<r_{z_{0}}.$ (36b)

Next, we consider the normal velocity components in Fourier space
$\hat{v}_{\alpha z}$. A general solution for the $4^{th}$-order differential
equations Eqs. (28) and (29a) can be written as: $\hat{v}_{\alpha
z}=(B+Cr_{z})\exp(\pm 2\pi kr_{z})$ Bickel (2007); Daddi-Moussa-Ider et al.
(2016); Bender and Orszag (2013). Again, noting that the disturbance flow
field must decay in the far-field, we write the solutions as:

$\displaystyle\hat{v}_{1z}$ $\displaystyle=$
$\displaystyle\begin{cases}(B_{1}+B_{2}r_{z})\exp(-2\pi kr_{z});r_{z}>0\\\
(B_{3}+B_{4}r_{z})\exp(2\pi kr_{z})+(B_{5}+B_{6}r_{z})\exp(-2\pi
kr_{z});r_{z_{0}}<r_{z}<0,\;\;\text{and}\\\ \end{cases}$ (37a)
$\displaystyle\hat{v}_{2z}$ $\displaystyle=$
$\displaystyle(B_{7}+B_{8}r_{z})\exp(2\pi kr_{z});r_{z}<r_{z_{0}}.$ (37b)

In Appendix A.2, we have four boundary conditions for the normal velocity
components, namely, the impenetrability of the velocity across the interface
Eq. (30a), the no-slip velocity at the interface Eq. (30c), the continuity of
tangential stress Eq. (31a) and the normal stress jump Eq. (33). Therefore, we
require four relations to uniquely determine the eight unknown constants
$B_{1}-B_{8}$, and we follow the protocol outlined above for $\hat{v}_{\alpha
t}$, using the Greens function properties. In this case, however, the forcing
on the right-hand side of Eq. (28) is proportional to
$\delta^{\prime}(r_{z})$. Therefore, the additional boundary conditions are
the continuity of $\hat{v}_{1z}$, its first and third derivative with respect
$r_{z}$ at $r_{z}=0$, whereas the second derivative undergoes a jump
discontinuity Bickel (2007); Daddi-Moussa-Ider et al. (2016); Bender and
Orszag (2013). Note that if the third derivative exhibited a finite jump
discontinuity, then a $\delta(r_{z})$ would appear on the right-hand side of
Eq. (28). The four additional boundary conditions are

$\displaystyle\left.\hat{v}_{1z}\right|_{r_{z}=0^{+}}$ $\displaystyle=$
$\displaystyle\left.\hat{v}_{1z}\right|_{r_{z}=0^{-}},$ (38a)
$\displaystyle\left.\frac{\partial\hat{v}_{1z}}{\partial
r_{z}}\right|_{r_{z}=0^{+}}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\hat{v}_{1z}}{\partial
r_{z}}\right|_{r_{z}=0^{-}},$ (38b)
$\displaystyle\left.\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}\right|_{r_{z}=0^{+}}-\left.\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}\right|_{r_{z}=0^{-}}$ $\displaystyle=$
$\displaystyle-\frac{4\,D}{\ln\kappa}\sin^{2}\left(\frac{\pi}{2}kp_{l}\right),\;\;\text{and}$
(38c) $\displaystyle\left.\frac{\partial^{3}\hat{v}_{1z}}{\partial
r_{z}^{3}}\right|_{r_{z}=0^{+}}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial^{3}\hat{v}_{1z}}{\partial
r_{z}^{3}}\right|_{r_{z}=0^{-}},$ (38d)

allowing $\hat{v}_{\alpha z}$ to be uniquely determined as:

$\displaystyle\hat{v}_{1z}$ $\displaystyle=$ $\displaystyle\frac{D}{2\pi
k\ln\kappa}\sin^{2}\left(\frac{\pi}{2}kp_{l}\right)$ (39a)
$\displaystyle\left[r_{z}\exp(-2\pi
k|r_{z}|)+\left(\frac{1-\lambda}{1+\lambda}\right)(r_{z}+4\pi
kr_{z_{0}}(r_{z}-r_{z_{0}}))\exp(4\pi kr_{z_{0}})\exp(-2\pi kr_{z})\right]$
$\displaystyle-\left(\frac{\kappa_{\beta}}{\eta_{1}V_{s}L^{2}}\right)\frac{\pi
k\hat{u}_{z}}{(1+\lambda)}\left[1+2\pi
k(r_{z}-r_{z_{0}})\right]\left(4\pi^{2}k^{2}+\Gamma\right)\exp(-2\pi
k(r_{z}-r_{z_{0}}))$ $\displaystyle;r_{z}>r_{z_{0}},\;\;\text{and}$
$\displaystyle\hat{v}_{2z}$ $\displaystyle=$ $\displaystyle\frac{D}{\pi
k(1+\lambda)\ln\kappa}\sin^{2}\left(\frac{\pi}{2}kp_{l}\right)r_{z}\exp(2\pi
kr_{z})$ (39b)
$\displaystyle-\left(\frac{\kappa_{\beta}}{\eta_{1}V_{s}L^{2}}\right)\frac{\pi
k\hat{u}_{z}}{(1+\lambda)}\left[1+2\pi
k(r_{z_{0}}-r_{z})\right]\left(4\pi^{2}k^{2}+\Gamma\right)\exp(-2\pi
k(r_{z_{0}}-r_{z}))$ $\displaystyle;r_{z}<r_{z_{0}}.$

We now use the continuity equations (27d) and (29c) to determine the
longitudinal velocity components $\hat{v}_{\alpha l}$. Note that the
longitudinal component of the fluid velocity field is independent of the
interface deformation at $r_{z}=r_{z_{0}}$.

##  Appendix B Swimmers orthogonal to interface - Two-dimensional Fourier
transform summary

Here, for swimmers oriented orthogonal to the interface we summarize the
expressions for the disturbance flow field in Fourier space. In this case,
$p_{x}=p_{y}=0$, whereas $p_{z}=+1(-1)$, for swimmers oriented away from
(towards) the interface. While the continuity equations given by Eqs. (27d)
and (29c) remain unchanged, the Stokes equations become:

$\displaystyle-2\pi{i\mkern
1.0mu}k_{x}\hat{P}_{1}+\frac{\partial^{2}\hat{v}_{1x}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1x}$ $\displaystyle=$ $\displaystyle 0,$
(40a) $\displaystyle-2\pi{i\mkern
1.0mu}k_{y}\hat{P}_{1}+\frac{\partial^{2}\hat{v}_{1y}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1y}$ $\displaystyle=$ $\displaystyle
0,\;\;\text{and}$ (40b) $\displaystyle-\frac{\partial\hat{P}_{1}}{\partial
r_{z}}+\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1z}$ $\displaystyle=$
$\displaystyle\frac{D}{\ln\kappa}\text{sgn}(r_{z})H\left[\frac{1}{2}-\frac{r_{z}}{p_{z}}\right]H\left[\frac{1}{2}+\frac{r_{z}}{p_{z}}\right],$
(40c)

where $\text{sgn}(x)$ is the sign-function and $H[x]$ is the Heaviside step
function Abramowitz et al. (1988). The corresponding equations for
$\hat{\bm{v}}_{2}$ and $P_{2}$ have a right-hand side equal to zero. Following
the steps described in Appendix A.1, in the
($\hat{\bm{l}},\hat{\bm{t}},\bm{z}$) coordinates, the Stokes equations become:

$\displaystyle\frac{\partial^{4}\hat{v}_{1z}}{\partial
r_{z}^{4}}-8\pi^{2}k^{2}\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}+16\pi^{4}k^{4}\hat{v}_{1z}$ $\displaystyle=$
$\displaystyle-\frac{D}{\ln\kappa}4\pi^{2}k^{2}\text{sgn}(r_{z})H\left[\frac{1}{2}-\frac{r_{z}}{p_{z}}\right]H\left[\frac{1}{2}+\frac{r_{z}}{p_{z}}\right],$
(41a) $\displaystyle\frac{\partial^{4}\hat{v}_{2z}}{\partial
r_{z}^{4}}-8\pi^{2}k^{2}\frac{\partial^{2}\hat{v}_{2z}}{\partial
r_{z}^{2}}+16\pi^{4}k^{4}\hat{v}_{2z}$ $\displaystyle=$ $\displaystyle
0,\;\;\text{and}$ (41b) $\displaystyle\frac{\partial^{2}\hat{v}_{1t}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1t}$ $\displaystyle=$ $\displaystyle 0,$
(41c)

with the corresponding continuity equations being given by Eqs. (27d) and
(29c). We note that unlike swimmers with a non-zero $p_{x}$ and/or $p_{y}$, in
this case the $\hat{v}_{\alpha t}$ satisfy a homogeneous second order
differential equation. Therefore, the no-slip and the tangential stress
boundary conditions are sufficient to determine the transverse component of
the velocity field, and give $\hat{v}_{\alpha t}=0$.

To determine the normal component, we first determine the Green’s function of
the differential operator in fluid region 1, which is

$\frac{\partial^{4}\hat{v}_{1z}}{\partial
r_{z}^{4}}-8\pi^{2}k^{2}\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}+16\pi^{4}k^{4}\hat{v}_{1z}=F_{1}\delta(r_{z}-r_{z}^{\prime}),$ (42)

for an arbitrary forcing $F_{1}$, and then express the solution of Eq. (41a)
as a convolution integral with the Green’s function Bender and Orszag (2013).
In other words, we solve for Eqs. (41b) and (42). The solution form for
$\hat{v}_{\alpha z}$ remains the same as described in Appendix A.3. In this
case, the additional boundary conditions are Bickel (2007); Bender and Orszag
(2013):

$\displaystyle\left.\hat{v}_{1z}\right|_{r_{z}=r_{z}^{\prime+}}$
$\displaystyle=$
$\displaystyle\left.\hat{v}_{1z}\right|_{r_{z}=r_{z}^{\prime-}},$ (43a)
$\displaystyle\left.\frac{\partial\hat{v}_{1z}}{\partial
r_{z}}\right|_{r_{z}=r_{z}^{\prime+}}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\hat{v}_{1z}}{\partial
r_{z}}\right|_{r_{z}=r_{z}^{\prime-}},$ (43b)
$\displaystyle\left.\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}\right|_{r_{z}=r_{z}^{\prime+}}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}\right|_{r_{z}=r_{z}^{\prime-}},\;\;\text{and}$ (43c)
$\displaystyle\left.\frac{\partial^{3}\hat{v}_{1z}}{\partial
r_{z}^{3}}\right|_{r_{z}=r_{z}^{\prime+}}-\left.\frac{\partial^{3}\hat{v}_{1z}}{\partial
r_{z}^{3}}\right|_{r_{z}=r_{z}^{\prime-}}$ $\displaystyle=$ $\displaystyle
F_{1}.$ (43d)

We note that the additional boundary conditions obtained here remain true for
any orientation of the swimmer with $p_{z}\neq 0$, as the
$\bm{\delta}$-function in Eq. (1) integrates out. Therefore, swimmers oriented
parallel to the interface are a special case. Using Eqs. (30a), (30c), (31a),
(33) and (43) yields:

$\displaystyle\hat{v}_{1z-G}(r_{z}|r_{z}^{\prime})$ $\displaystyle=$
$\displaystyle\frac{F_{1}}{32\pi^{3}k^{3}}\left(1+2\pi
k\left|r_{z}-r_{z}^{\prime}\right|\right)\exp(-2\pi
k\left|r_{z}-r_{z}^{\prime}\right|)$ (44a)
$\displaystyle+\frac{F_{1}}{32\pi^{3}k^{3}}\left(\frac{1-\lambda}{1+\lambda}\right)\left[1+2\pi
kr_{z}^{\prime}+2\pi kr_{z}-4\pi
kr_{z_{0}}+8\pi^{2}k^{2}(r_{z}-r_{z_{0}})(r_{z}^{\prime}-r_{z_{0}})\right]$
$\displaystyle\qquad\qquad\exp(4\pi kr_{z_{0}})\exp(-2\pi
k(r_{z}+r_{z}^{\prime}))$
$\displaystyle-\left(\frac{\kappa_{\beta}}{\eta_{1}V_{s}L^{2}}\right)\frac{\pi
k}{(1+\lambda)}\left[1+2\pi
k(r_{z}-r_{z_{0}})\right]\left(4\pi^{2}k^{2}+\Gamma\right)\hat{u}_{z}\exp(-2\pi
k(r_{z}-r_{z_{0}}))$ $\displaystyle;r_{z}>r_{z_{0}},\;\;\text{and}$
$\displaystyle\hat{v}_{2z-G}(r_{z}|r_{z}^{\prime})$ $\displaystyle=$
$\displaystyle\frac{F_{1}}{(1+\lambda)}\frac{1}{16\pi^{2}k^{3}}(1+2\pi
k(r_{z}^{\prime}-r_{z}))\exp(-2\pi k(r_{z}^{\prime}-r_{z}))$ (44b)
$\displaystyle-\left(\frac{\kappa_{\beta}}{\eta_{1}V_{s}L^{2}}\right)\frac{\pi
k}{(1+\lambda)}\left[1+2\pi
k(r_{z_{0}}-r_{z})\right]\left(4\pi^{2}k^{2}+\Gamma\right)\hat{u}_{z}\exp(2\pi
k(r_{z}-r_{z_{0}}))$ $\displaystyle;r_{z}<r_{z_{0}},$

where the subscript $G$ refers to the Green’s function. We can now obtain
$\hat{v}_{\alpha z}$ by convolving $\hat{v}_{\alpha z-G}$ with the original
forcing, that is, replacing $F_{1}$ by the right-hand side of Eq. (41a) and
integrating over $r_{z}^{\prime}$. Thus, we solve

$\hat{v}_{\alpha
z}=-\frac{4\pi^{2}k^{2}D}{\ln\kappa}\int\mathrm{d}r_{z}^{\prime}\hat{v}_{\alpha
z-G}(r_{z}|r_{z}^{\prime})\text{sgn}(r_{z}^{\prime})H\left[\frac{1}{2}-\frac{r_{z}^{\prime}}{p_{z}}\right]H\left[\frac{1}{2}+\frac{r_{z}^{\prime}}{p_{z}}\right],$
(45)

where only the terms involving $F_{1}$ in Eq. (44) are to be integrated over.
This yields

$\displaystyle\hat{v}_{1z}$ $\displaystyle=$
$\displaystyle-\left(\frac{\kappa_{\beta}}{\eta_{1}V_{s}L^{2}}\right)\frac{\pi
k}{(1+\lambda)}\left[1+2\pi
k(r_{z}-r_{z_{0}})\right]\left(4\pi^{2}k^{2}+\Gamma\right)\hat{u}_{z}\exp(-2\pi
k(r_{z}-r_{z_{0}}))$ (46a)
$\displaystyle+\frac{D}{8\pi^{2}k^{2}\ln\kappa}\left(\frac{1-\lambda}{1+\lambda}\right)\exp(4\pi
kr_{z_{0}})\exp(-2\pi kr_{z})$ $\displaystyle\qquad\qquad\bigg{[}2(\cosh(\pi
k)-1)(1+3\pi kr_{z}-4\pi kr_{z_{0}}$
$\displaystyle\qquad\qquad-4\pi^{2}k^{2}r_{z_{0}}(r_{z}-r_{z_{0}}))-\pi
k\sinh(\pi k)(1+4\pi k(r_{z}-r_{z_{0}}))\bigg{]}$
$\displaystyle-\frac{D}{8\pi^{2}k^{2}\ln\kappa}\begin{cases}\exp(2\pi
kr_{z})\Big{[}2(\pi kr_{z}-1)(\cosh(\pi k)-1)+\pi k\sinh(\pi k)\Big{]}\\\
\qquad\qquad;r_{z}<-\left|\frac{p_{z}}{2}\right|\\\ \\\ \left(1+\frac{\pi
k}{2}(1+2r_{z})\right)\exp(-\pi k(1+2r_{z}))\\\ -\left(1+\frac{\pi
k}{2}(1-2r_{z})\right)\exp(-\pi k(1-2r_{z}))+2\text{sgn}(r_{z})\\\
-\left(2\text{sgn}(r_{z})+2\pi kr_{z}\right)\exp\left(-2\pi
k\left|r_{z}\right|\right);r_{z}\in\left[-\left|\frac{p_{z}}{2}\right|,\left|\frac{p_{z}}{2}\right|\right]\\\
\\\ \exp(-2\pi kr_{z})\Big{[}2\pi kr_{z}(\cosh(\pi k)-1)-\pi k\sinh(\pi
k)\Big{]};r_{z}>\left|\frac{p_{z}}{2}\right|,\;\;\text{and}\\\ \end{cases}$
$\displaystyle\hat{v}_{2z}$ $\displaystyle=$
$\displaystyle-\left(\frac{\kappa_{\beta}}{\eta_{1}V_{s}L^{2}}\right)\frac{\pi
k}{(1+\lambda)}\left[1+2\pi
k(r_{z_{0}}-r_{z})\right]\left(4\pi^{2}k^{2}+\Gamma\right)\hat{u}_{z}\exp(2\pi
k(r_{z}-r_{z_{0}}))$ (46b)
$\displaystyle+\frac{D}{4\pi^{2}k^{2}\ln\kappa}\frac{\exp(2\pi
kr_{z})}{(1+\lambda)}\left[2(1-\pi kr_{z})\left(\cosh(\pi k)-1\right)-\pi
k\sinh(\pi k)\right];r_{z}<r_{z_{0}}.$

Again, using $\hat{v}_{\alpha z}$ from Eq. (46) in the continuity equations
(27d) and (29c) readily yield the longitudinal velocity components
$\hat{v}_{\alpha l}$.

##  Appendix C Swimmers orthogonal to interface - Validation of the kinematic
boundary approximation

Here, for swimmers oriented orthogonal to the interface we compare the
kinematic boundary condition from Sec. II.2 given by Eq. (9), which contains
the non-linear terms, with the approximate condition Eq. (10) that we have
used in the main manuscript. The decoupling of the hydrodynamics for this
swimmer configuration enables us to solve both of the equations numerically in
a convenient manner. The expanded form Eq. (9) in Cartesian coordinate system
is:

$\frac{\partial u_{z}}{\partial
t}+\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\bigg{[}v_{x}\frac{\partial
u_{z}}{\partial r_{x}}+v_{y}\frac{\partial u_{z}}{\partial
r_{y}}-u_{z}\frac{\partial v_{z}}{\partial
r_{z}}\bigg{]}=\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\Big{[}\left.v_{z}\right|_{r_{z_{0}}}+p_{z}\Big{]}.$
(47)

Given that $\hat{v}_{\alpha t}$ is zero for this swimmer configuration
($\alpha\in[1,2]$ for the two fluid regions), the Fourier transformed planar
velocity components are proportional only to $\hat{v}_{\alpha l}$, and are
$\hat{v}_{\alpha x}=k_{x}\hat{v}_{\alpha l}/k$ and $\hat{v}_{\alpha
y}=k_{y}\hat{v}_{\alpha l}/k$. We can then readily carry out the inverse
Fourier transform of the planar and normal components of the disturbance
velocity and the gradient of the normal velocity components to get the
following expressions for $v_{x}$, $v_{y}$, $v_{z}$ and $\partial
v_{z}/\partial r_{z}$ at $r_{z_{0}}$:

$\displaystyle\left.v_{x}\right|_{r_{z_{0}}}$ $\displaystyle=$
$\displaystyle-\frac{Dr_{x}}{2\pi(1+\lambda)\ln\kappa}\left[\frac{1}{\left(4r_{\|}^{2}+\left(1-2r_{z_{0}}^{2}\right)^{2}\right)^{\frac{1}{2}}}+\frac{1}{\left(4r_{\|}^{2}+\left(1+2r_{z_{0}}^{2}\right)^{2}\right)^{\frac{1}{2}}}-\frac{1}{\left(r_{\|}^{2}+r_{z_{0}}^{2}\right)^{\frac{1}{2}}}\right],$
$\displaystyle\left.v_{y}\right|_{r_{z_{0}}}$ $\displaystyle=$
$\displaystyle-\frac{Dr_{y}}{2\pi(1+\lambda)\ln\kappa}\left[\frac{1}{\left(4r_{\|}^{2}+\left(1-2r_{z_{0}}^{2}\right)^{2}\right)^{\frac{1}{2}}}+\frac{1}{\left(4r_{\|}^{2}+\left(1+2r_{z_{0}}^{2}\right)^{2}\right)^{\frac{1}{2}}}-\frac{1}{\left(r_{\|}^{2}+r_{z_{0}}^{2}\right)^{\frac{1}{2}}}\right],$
$\displaystyle\left.v_{z}\right|_{r_{z_{0}}}$ $\displaystyle=$
$\displaystyle\frac{D}{2\pi(1+\lambda)\ln\kappa}\int_{0}^{\infty}\mathrm{d}k\frac{1}{k}J_{0}\left(2kr_{\|}\right)\exp\left(2kr_{z_{0}}\right)\Big{[}2\left(1-kr_{z_{0}}\right)\left(\cosh(k)-1\right)-k\sinh(k)\Big{]},\;\;\text{and}$
(50) $\displaystyle\left.\frac{\partial v_{z}}{\partial
r_{z}}\right|_{r_{z_{0}}}$ $\displaystyle=$
$\displaystyle\frac{D}{2\pi(1+\lambda)\ln\kappa}\left[\frac{2\left(2r_{\|}^{2}+\left(1+2r_{z_{0}}^{2}\right)^{2}\right)}{\left[4r_{\|}^{2}+\left(1+2r_{z_{0}}^{2}\right)^{2}\right]^{\frac{3}{2}}}+\frac{2\left(2r_{\|}^{2}+\left(1-2r_{z_{0}}^{2}\right)^{2}\right)}{\left[4r_{\|}^{2}+\left(1-2r_{z_{0}}^{2}\right)^{2}\right]^{\frac{3}{2}}}-\frac{r_{\|}^{2}+2r_{z_{0}}^{2}}{\left[r_{\|}^{2}+r_{z_{0}}^{2}\right]^{\frac{3}{2}}}\right].$

The radial symmetry of $v_{x}$ and $v_{y}$ allows us to solve Eq. (47) in
cylindrical polar coordinates ($\hat{\bm{r}}_{\|}$, $\hat{\bm{\theta}}$,
$\hat{\bm{z}}$). Because $v_{\theta}=0$, there is a dimensional reduction in
the polar coordinate system; $r_{\|}^{2}=r_{x}^{2}+r_{y}^{2}$ and
$\theta=\arctan(r_{y}/r_{x})$. Thus, we solve:

$\frac{\partial u_{z}}{\partial
t}+\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\bigg{[}v_{r_{\|}}\frac{\partial
u_{z}}{\partial r_{\|}}-u_{z}\frac{\partial v_{z}}{\partial
z}\bigg{]}=\left(\frac{\eta_{1}L^{2}V_{s}}{\kappa_{\beta}}\right)\Big{[}\left.v_{z}\right|_{r_{z_{0}}}+p_{z}\Big{]}.$
(52)

where $v_{r_{\|}}=r_{x}v_{x}/r+r_{y}v_{y}/r$. Exploiting radial symmetry, we
calculate $\partial u_{z}/\partial r_{\|}$ as an inverse Hankel transform of
$\hat{u}_{z}$ with the Bessel function $J_{1}(2\pi kr_{\|})$. Here, we refer
to the interface deformation obtained from solving the approximate kinematic
condition Eq. (17) as $u_{z-approx}$.

In the discussion preceeding Eq. (10), we argued that the approximation to the
kinematic boundary condition will hold when the ratio of the viscous stress to
the bending stress $\left(\eta_{1}V_{s}L^{2}/\kappa_{\beta}\right)$ was $O(1)$
or smaller. As this ratio becomes larger, one would expect the largeness in
$\left(\eta_{1}V_{s}L^{2}/\kappa_{\beta}\right)$ to dominate over the weak
perturbation of the inverse logarithm of the slender swimmer aspect ratio
$\kappa$, and hence necessitate the inclusion of the neglected nonlinear
terms. In Fig. 15 we compare $u_{z}$ obtained from Eq. (52) with
$u_{z-approx}$ (both solved in along with Eq. (18) for $V_{z}^{T}$), for
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=1,10$ to assess the relative importance of
the nonlinear terms on the left-hand side of Eq. (52). We focus only on
equations specific to shakers, and thus $p_{z}=0$ on the right-hand side of
both Eqs. (17) and (52). It is clear that even for an initial distance of the
swimmer from the interface of $r_{z_{0}}(0)=1$, $u_{z-approx}$ provides a very
good approximation to the full solution obtained from Eq. (52). In the inset,
a linear plot is shown for $\eta_{1}V_{s}L^{2}/\kappa_{\beta}=10$, showing
that the deviation between the two is small, with a maximum in the region of
$r_{\|}=0$ having a relative error of about $0.24$. Although not shown here,
$u_{z-approx}$ approaches the $u_{z}$ obtained form Eq. (52) as $r_{z_{0}}$
increases. This validates the boundary approximation for a slender swimmer
even at $O(1)$ distances from the interface and for values of
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}\sim O(10)$.

Figure 15: The interface deformation $u_{z}$ due to a pusher plotted as a
function of the distance parallel to the interface $r_{\|}$ on a log-log scale
for the ratio of viscous stress to bending stress:
$\eta_{1}V_{s}L^{2}/\kappa_{\beta}=$ 1 and 10, at $t=1$ when $\lambda=0.5$ and
$r_{z_{0}}(0)=1$. Solid lines represent the full solution obtained from Eq.
(52), and the dashed lines the approximate solution from Eq. (17). The inset
shows a linear plot of $u_{z}$ for $\eta V_{s}L^{2}/\kappa_{\beta}=10$.

##  Appendix D Swimmers arbitrarily oriented to the interface - Two-
dimensional Fourier transform summary

To solve for a swimmer arbitrarily oriented to the interface, we follow the
same approach to that of swimmers oriented orthogonal to the interface. Here,
the equations for $\hat{v}_{z}$ and $\hat{v}_{t}$ are similar to that for
swimmers orthogonal to the interface, as given by Eq. (41) in Appendix B. For
an arbitrarily oriented slender swimmer, we have the following set of
equations to solve for in the two fluid regions:

$\displaystyle\frac{\partial^{4}\hat{v}_{1z}}{\partial
r_{z}^{4}}-8\pi^{2}k^{2}\frac{\partial^{2}\hat{v}_{1z}}{\partial
r_{z}^{2}}+16\pi^{4}k^{4}\hat{v}_{1z}$ $\displaystyle=$
$\displaystyle-2\pi{i\mkern
1.0mu}k\frac{D}{\ln\kappa}\frac{1}{p_{z}^{2}}\exp\left(-2\pi{i\mkern
1.0mu}kr_{z}\frac{p_{l}}{p_{z}}\right)\bigg{\\{}2p_{l}p_{z}\;\delta(r_{z})H\left[\frac{1}{2}-\frac{r_{z}}{p_{z}}\right]H\left[\frac{1}{2}+\frac{r_{z}}{p_{z}}\right]$
(53a)
$\displaystyle\qquad+p_{l}\;\text{sgn}(r_{z})H\left[\frac{1}{2}-\frac{r_{z}}{p_{z}}\right]\delta\left[\frac{1}{2}+\frac{r_{z}}{p_{z}}\right]-p_{l}\;\text{sgn}(r_{z})\delta\left[\frac{1}{2}-\frac{r_{z}}{p_{z}}\right]H\left[\frac{1}{2}+\frac{r_{z}}{p_{z}}\right]$
$\displaystyle\qquad-2\pi{i\mkern
1.0mu}k(p_{z}^{2}+p_{l}^{2})\;\text{sgn}(r_{z})H\left[\frac{1}{2}-\frac{r_{z}}{p_{z}}\right]H\left[\frac{1}{2}+\frac{r_{z}}{p_{z}}\right]\bigg{\\}},$
$\displaystyle\frac{\partial^{4}\hat{v}_{2z}}{\partial
r_{z}^{4}}-8\pi^{2}k^{2}\frac{\partial^{2}\hat{v}_{2z}}{\partial
r_{z}^{2}}+16\pi^{4}k^{4}\hat{v}_{2z}$ $\displaystyle=$ $\displaystyle
0,\;\;\text{and}$ (53b) $\displaystyle\frac{\partial^{2}\hat{v}_{1t}}{\partial
r_{z}^{2}}-4\pi^{2}k^{2}\hat{v}_{1t}$ $\displaystyle=$
$\displaystyle\frac{D}{\ln\kappa}\frac{p_{t}}{p_{z}}\exp\left(-2\pi{i\mkern
1.0mu}kr_{z}\frac{p_{l}}{p_{z}}\right)\text{sgn}(r_{z})H\left[\frac{1}{2}-\frac{r_{z}}{p_{z}}\right]H\left[\frac{1}{2}+\frac{r_{z}}{p_{z}}\right].$
(53c)

After some algebra, we obtain the following expressions for the transverse
component of the disturbance flow field:

$\displaystyle\hat{v}_{1t}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi^{2}k^{2}}\frac{D}{\ln\kappa}\frac{p_{t}\;({i\mkern
1.0mu}p_{l}-p_{z})}{p_{l}^{2}+p_{z}^{2}}\left(\frac{1-\lambda}{1+\lambda}\right)\exp(4\pi
kr_{z_{0}})\exp(-2\pi kr_{z})\sinh^{2}\left(\frac{\pi}{2}k({i\mkern
1.0mu}p_{l}+p_{z})\right)$ (54a)
$\displaystyle-\frac{1}{4\pi^{2}k^{2}}\frac{D}{\ln\kappa}\frac{p_{t}}{p_{l}^{2}+p_{z}^{2}}\begin{cases}2({i\mkern
1.0mu}p_{l}-p_{z})\exp(2\pi kr_{z})\sinh^{2}\left(\frac{\pi}{2}k({i\mkern
1.0mu}p_{l}+p_{z})\right);r_{z}<-\left|\frac{p_{z}}{2}\right|\\\ \\\
\frac{1}{2}\Big{[}({i\mkern 1.0mu}p_{l}-p_{z})\exp(2\pi kr_{z})\exp(-\pi
k\;\text{sgn}(p_{z})({i\mkern 1.0mu}p_{l}+p_{z}))\\\ \quad+({i\mkern
1.0mu}p_{l}+p_{z})\exp(-2\pi kr_{z})\exp(\pi k\;\text{sgn}(p_{z})({i\mkern
1.0mu}p_{l}-p_{z}))\\\ \quad-2({i\mkern
1.0mu}p_{l}+\text{sgn}(r_{z})p_{z})\exp(-2\pi
k|r_{z}|)+2\;\text{sgn}(r_{z})p_{z}\exp\left(-2\pi{i\mkern
1.0mu}kr_{z}\frac{p_{l}}{p_{z}}\right)\Big{]}\\\
;r_{z}\in\left[-\left|\frac{p_{z}}{2}\right|,\left|\frac{p_{z}}{2}\right|\right]\\\
\\\ 2({i\mkern 1.0mu}p_{l}+p_{z})\exp(-2\pi
kr_{z})\sinh^{2}\left(\frac{\pi}{2}k({i\mkern
1.0mu}p_{l}-p_{z})\right);r_{z}>\left|\frac{p_{z}}{2}\right|,\;\;\text{and}\\\
\end{cases}$ $\displaystyle\hat{v}_{2t}$ $\displaystyle=$
$\displaystyle-\frac{1}{\pi^{2}k^{2}}\frac{D}{(1+\lambda)\ln\kappa}\frac{p_{t}\;({i\mkern
1.0mu}p_{l}-p_{z})}{p_{l}^{2}+p_{z}^{2}}\exp(2\pi
kr_{z})\sinh^{2}\left(\frac{\pi}{2}k({i\mkern
1.0mu}p_{l}+p_{z})\right);r_{z}<r_{z_{0}}.$ (54b)

For a swimmer oriented orthogonal to the interface, $p_{t}=0$, implying that
$\hat{v}_{\alpha t}=0$, as pointed out in Appendix B.

The Green’s function for the normal components $\hat{v}_{\alpha z}$ remain the
same as given by Eq. (44), where $\alpha\in[1,2]$ for the two fluid regions.
Convolving the Green’s function with the right-hand side of Eq. (53a) gives
the normal component of the velocity. The longitudinal components
$\hat{v}_{\alpha l}$ can then be obtained from the use of the continuity
equations (27d) and (29c). We omit writing the rather cumbersome expressions
for the normal and longitudinal components of the fluid velocity that emerge
in this swimmer configuration.

##  Appendix E Swimmers confined between a deformable interface and a rigid
boundary - Two-dimensional Fourier transform summary

Here, we briefly describe the formulation for a swimmer confined between a
rigid boundary and an initially undeformed interface, specifically when its
oriented parallel to the boundaries. We follow the development in Appendix A,
and outline only the additional boundary conditions that emerge owing to the
additional rigid boundary. For this swimmer orientation only the $z$-component
of the velocity is relevant, given that $\dot{\bm{p}}=0$ as mentioned in Sec.
IV. We therefore focus solely on the additional boundary conditions for
$\hat{v}_{\alpha z}$, where $\alpha\in[1,2]$ for the two fluid regions. For a
swimmer of arbitrary orientation, however, the rotation rate would be non-
trivial and it would be essential to derive the transverse and longitudinal
velocity components, $\hat{v}_{\alpha t}$ and $\hat{v}_{\alpha l}$.

The equations for $\hat{v}_{\alpha z}$ are solved subject to the appropriate
boundary conditions; the impenetrability and the no-slip conditions at the
deformable interface and the rigid boundary, the continuity of tangential
stress and the normal-stress jump across the deformable interface. The
boundary conditions at the deformable interface are as in Appendix A.2, namely
Eqs. (30a), (30c), (31a) and (33). The additional Fourier transformed velocity
boundary conditions owing to the rigid boundary are:

$\displaystyle\left.\hat{v}_{1z}\right|_{H^{-}}$ $\displaystyle=$
$\displaystyle 0,$ (55a)
$\displaystyle\left.\frac{\partial\hat{v}_{1z}}{\partial
r_{z}}\right|_{H^{-}}$ $\displaystyle=$ $\displaystyle 0.$ (55b)

A general solution for the $4^{th}$-order differential equations Eqs. (28) and
(29a), noting that the disturbance flow field must decay in the far-field, is

$\displaystyle\hat{v}_{1z}$ $\displaystyle=$
$\displaystyle\begin{cases}(B_{1}+B_{2}r_{z})\exp(-2\pi
kr_{z})+(B_{3}+B_{4}r_{z})\exp(-2\pi kr_{z});H>r_{z}>r_{z_{0}}\\\
(B_{5}+B_{6}r_{z})\exp(2\pi kr_{z})+(B_{7}+B_{8}r_{z})\exp(-2\pi
kr_{z});0<r_{z}<r_{z_{0}},\;\;\text{and}\\\ \end{cases}$ (56a)
$\displaystyle\hat{v}_{2z}$ $\displaystyle=$
$\displaystyle(B_{9}+B_{10}r_{z})\exp(2\pi kr_{z});r_{z}<0.$ (56b)

We use the six boundary conditions mentioned above, in addition to the four
conditions emerging from the properties of the Green’s function, which remain
the same as in Eq. (38), to obtain the constants $B_{1}-B_{10}$. We omit for
brevity the cumbersome expressions for $\hat{v}_{\alpha z}$.

## References

  * Rothschild (1963) Rothschild, Nature (London) 198, 1221 (1963).
  * Frymier et al. (1995) P. D. Frymier, R. M. Ford, H. C. Berg, and P. T. Cummings, Proc. Nat. Acad. Sci. U.S.A. 92, 6195 (1995).
  * DiLuzio et al. (2005) W. R. DiLuzio, L. Turner, M. Mayer, P. Garstecki, D. B. Weibel, H. C. Berg, and G. M. Whitesides, Nature 435, 1271 (2005).
  * Berke et al. (2008) A. P. Berke, L. Turner, H. C. Berg, and E. Lauga, Phys. Rev. Lett. 101, 038102 (2008).
  * Di Leonardo et al. (2011) R. Di Leonardo, D. Dell’Arciprete, L. Angelani, and V. Iebba, Phys. Rev. Lett. 106, 038101 (2011).
  * Ferracci et al. (2013) J. Ferracci, H. Ueno, K. Numayama-Tsuruta, Y. Imai, T. Yamaguchi, and T. Ishikawa, PLOS ONE 8, 1 (2013).
  * Molaei et al. (2014) M. Molaei, M. Barry, R. Stocker, and J. Sheng, Phys. Rev. Lett. 113, 068103 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.113.068103.
  * Montecucco and Rappuoli (2001) C. Montecucco and R. Rappuoli, Nat. Rev. Mol. Cell Biol. 2, 457 (2001).
  * Moriarty et al. (2008) T. J. Moriarty, M. U. Norman, P. Colarusso, T. Bankhead, P. Kubes, and G. Chaconas, PLOS Pathogens 4, 1 (2008).
  * Lee et al. (2008) S. Lee, J. W. M. Bush, A. E. Hosoi, and E. Lauga, Phys. Fluids 20, 082106 (2008).
  * Sekirov et al. (2010) I. Sekirov, S. L. Russell, L. C. M. Antunes, and B. B. Finlay, Physiological Reviews 90, 859 (2010).
  * Trouilloud et al. (2008) R. Trouilloud, T. S. Yu, A. E. Hosoi, and E. Lauga, Phys. Rev. Lett. 101, 048102 (2008).
  * Spagnolie and Lauga (2012) S. E. Spagnolie and E. Lauga, J. Fluid Mech. 700, 105 (2012).
  * Dias and Powers (2013) M. A. Dias and T. R. Powers, Phys. Fluids 25, 101901 (2013).
  * Ledesma-Aguilar and Yeomans (2013) R. Ledesma-Aguilar and J. M. Yeomans, Phys. Rev. Lett. 111, 138101 (2013).
  * Bechinger et al. (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 045006 (2016).
  * Shaik and Ardekani (2017) V. A. Shaik and A. M. Ardekani, J. Fluid Mech. 824, 42 (2017).
  * Ishikawa (2019) T. Ishikawa, J. App. Phys. 125, 200901 (2019).
  * Daddi-Moussa-Ider et al. (2019) A. Daddi-Moussa-Ider, C. Kurzthaler, C. Hoell, A. Zöttl, M. Mirzakhanloo, M.-R. Alam, A. M. Menzel, H. Löwen, and S. Gekle, Phys. Rev. E 100, 032610 (2019).
  * Škultéty and Morozov (2020) V. Škultéty and A. Morozov, J. Fluid Mech. 882, A1 (2020).
  * Dalal et al. (2020) S. Dalal, A. Farutin, and C. Misbah, Soft Matter 16, 1599 (2020).
  * Sokolov et al. (2009) A. Sokolov, R. E. Goldstein, F. I. Feldchtein, and I. S. Aranson, Phys. Rev. E 80, 031903 (2009), URL https://link.aps.org/doi/10.1103/PhysRevE.80.031903.
  * Kurtuldu et al. (2011) H. Kurtuldu, J. S. Guasto, K. A. Johnson, and J. P. Gollub, Proc. Nat. Acad. Sci. U.S.A. 108, 10391 (2011).
  * Woodhouse and Goldstein (2012) F. G. Woodhouse and R. E. Goldstein, Phys. Rev. Lett. 109, 168105 (2012).
  * Lushi et al. (2014) E. Lushi, H. Wioland, and R. E. Goldstein, Proc. Nat. Acad. Sci. U.S.A. 111, 9733 (2014).
  * Costanzo et al. (2014) A. Costanzo, J. Elgeti, T. Auth, G. Gompper, and M. Ripoll, EPL (Europhysics Letters) 107, 36003 (2014).
  * Wioland et al. (2016) H. Wioland, E. Lushi, and R. E. Goldstein, New J. Physics 18, 075002 (2016).
  * Mathijssen et al. (2016) A. J. T. M. Mathijssen, T. N. Shendruk, J. M. Yeomans, and A. Doostmohammadi, Phys. Rev. Lett. 116, 028104 (2016).
  * Theillard et al. (2017) M. Theillard, R. Alonso-Matilla, and D. Saintillan, Soft Matter 13, 363 (2017).
  * Li and Tang (2009) G. Li and J. X. Tang, Phys. Rev. Lett. 103, 078101 (2009).
  * Rusconi et al. (2014) R. Rusconi, J. S. Guasto, and R. Stocker, Nature Phys. 10, 212 (2014).
  * Ezhilan et al. (2015) B. Ezhilan, R. Alonso-Matilla, and D. Saintillan, J. Fluid Mech. 781, R4 (2015).
  * Bearon and Hazel (2015) R. Bearon and A. Hazel, J. Fluid Mech. 771, R3 (2015).
  * Barry et al. (2015) M. T. Barry, R. Rusconi, J. S. Guasto, and R. Stocker, J. Royal Soc. Interf. 12, 20150791 (2015).
  * Ezhilan and Saintillan (2015) B. Ezhilan and D. Saintillan, J. Fluid Mech. 777, 482–522 (2015).
  * Manabe et al. (2020) J. Manabe, T. Omori, and T. Ishikawa, J. Fluid Mech. 892, A15 (2020).
  * Vennamneni et al. (2020) L. Vennamneni, S. Nambiar, and G. Subramanian, J. Fluid Mech. 890, A15 (2020).
  * Yan and Brady (2015) W. Yan and J. F. Brady, J. Fluid Mech. 785, R1 (2015).
  * Lee et al. (2017) A. A. Lee, D. Vella, and J. S. Wettlaufer, Proc. Nat. Acad. Sci. U.S.A. 114, 9255 (2017).
  * Hu et al. (2015) J. Hu, A. Wysocki, R. G. Winkler, and G. Gompper, Sci. Rep. 5, 1 (2015).
  * lea (2007) _Advanced transport phenomena: fluid mechanics and convective transport processes_ , Cambridge Series in Chemical Engineering (Cambridge University Press, 2007), ISBN 9781139462068\.
  * Bickel (2007) T. Bickel, Phys. Rev. E 75, 041403 (2007).
  * Pozrikidis (2007) C. Pozrikidis, J. Fluid Mech. 575, 333–357 (2007).
  * Daddi-Moussa-Ider et al. (2016) A. Daddi-Moussa-Ider, A. Guckenberger, and S. Gekle, Phys. Rev. E 93, 012612 (2016).
  * Rallabandi et al. (2018) B. Rallabandi, N. Oppenheimer, M. Y. B. Zion, and H. A. Stone, Nature Phys. 14, 1211 (2018).
  * J. and Amy (1995) F. L. J. and M. Amy, Bull. Math. Biol. 57, 679 (1995), ISSN 0092-8240.
  * Berdan and Leal (1982) C. Berdan and L. Leal, J. Colloid Interface Sci. 87, 62 (1982), ISSN 0021-9797.
  * Blake (1971) J. R. Blake, J. Fluid Mech. 46, 199–208 (1971).
  * Berg (2008) H. Berg, _E. coli in Motion_ , Biological and Medical Physics, Biomedical Engineering (Springer New York, 2008), ISBN 9780387216386.
  * Elgeti et al. (2015) J. Elgeti, R. G. Winkler, and G. Gompper, Rep. Prog. Phys. 78, 056601 (2015).
  * Kasyap et al. (2014) T. V. Kasyap, D. L. Koch, and M. Wu, Phys. Fluids 26, 081901 (2014).
  * Nambiar et al. (2021) S. Nambiar, P. Garg, and G. Subramanian, J. Fluid Mech. 907, A26 (2021).
  * Batchelor (1970) G. Batchelor, J. Fluid Mech. 44, 419 (1970).
  * Faucon et al. (1989) J. F. Faucon, M. D. Mitov, P. Méléard, I. Bivas, and P. l. Bothorel, J. Phys. France 50, 2389 (1989).
  * Fygenson et al. (1997) D. K. Fygenson, J. F. Marko, and A. Libchaber, Phys. Rev. Lett. 79, 4497 (1997).
  * Rosen et al. (2005) M. J. Rosen, H. Wang, P. Shen, and Y. Zhu, Langmuir 21, 3749 (2005).
  * Freund (2014) J. B. Freund, Annu. Rev. Flu. Mech. 46, 67 (2014).
  * Lighthill (1958) M. J. Lighthill, _An introduction to Fourier analysis and generalised functions_ , Cambridge Monographs on Mechanics (Cambridge University Press, 1958).
  * Lauga and Powers (2009) E. Lauga and T. R. Powers, Rep. Prog. Phys 72, 096601 (2009).
  * Subramanian and Nott (2011) G. Subramanian and P. R. Nott, J. Indian Inst. Sc. 91, 283 (2011).
  * Brennen and Winet (1977) C. Brennen and H. Winet, Annu. Rev. Flu. Mech. 9, 339 (1977).
  * Patteson et al. (2016) A. E. Patteson, A. Gopinath, P. K. Purohit, and P. E. Arratia, Soft Matter 12, 2365 (2016).
  * (63) W. Helfrich, Zeitschrift für Naturforschung C 28, 693 (????).