Hydrodynamic interactions between a sedimenting squirmer and a planar wall

Abstract

The hydrodynamic interactions between a sedimenting microswimmer and a solid wall have ubiquitous biological and technological applications. A plethora of gravity-induced swimming dynamics near a planar no-slip wall provide a platform for designing artificial microswimmers that can generate directed propulsion through their translation–rotation coupling near a wall. In this work, we provide exact solutions for a squirmer (a model swimmer of spherical shape with a prescribed slip velocity) facing either towards or away from a planar wall perpendicular to gravity. These exact solutions are used to validate a numerical code based on the boundary integral method with an adaptive mesh for distances from the wall down to 0.1 % of the squirmer radius. This boundary integral code is then used to investigate the rich gravity-induced dynamics near a wall, mapping out the detailed bifurcation structures of the swimming dynamics in terms of orientation and distance to the wall. Simulation results show that a squirmer may traverse the wall, move to a fixed point at a given height with a fixed orientation in a monotonic way or in an oscillatory fashion, or oscillate in a limit cycle in the presence of wall repulsion.

1. Introduction

Microswimmers behave very differently near a wall as their interactions with a solid boundary alter their speed and direction, and how they interact with each other (Shum et al. 2010; Takagi et al. 2014; Elgeti & Gompper 2016), giving rise to many interesting phenomena, such as the swirling of bacteria next to a substrate, and clustering of phoretic Janus particles and bacteria near a wall. While the far-field flow attracts and aligns a pusher (puller) to move along (normal to) a no-slip wall, near-field hydrodynamics, steric interactions and contact dynamics give rise to wall scattering, with the swimmer escaping from the wall at a characteristic angle independent of the initial direction of approach to the wall (Berke et al. 2008; Li & Tang 2009).

The hydrodynamic interactions between a microswimmer and a solid wall are more complex when the swimmer sediments to the wall under gravity (due to the density mismatch between the swimmer and the surrounding fluid). Several types of dynamics of a sedimenting swimmer have been reported: scattering (escaping) from the wall, swimming along the wall at a fixed distance and tilted orientation, and periodic bouncing on the wall (Or & Murray 2009; Crowdy & Or 2010). Under gravity, artificial surface walkers or micro rollers stay close to the wall. At the same time, they rotate under an external force field, exploiting their interactions with a solid surface to generate directed propulsion (Tierno et al. 2008; Sing et al. 2010; Driscoll et al. 2017). These microswimmers are easy to manipulate for directed transport, and offer wide applications in targeted therapeutics and microsurgery (Alapan et al. 2020; Ahmed et al. 2021).

The interactions between a flagellated swimmer and a flat solid wall have been modelled using a multipole approach (Spagnolie & Lauga 2012), which is shown to give good agreement with boundary integral simulations. Through the contribution of each singularity to the effect of a wall on the flagellated swimmer, the reduced model captures the main dynamic features of the wall-induced hydrodynamics of a microswimmer. Alternatively, a microswimmer is often simplified and modelled as a squirming sphere with a slip velocity on the surface to mimic the surrounding flow created by the layer of beating cilia on the microswimmer (Lighthill 1952; Blake 1971b ). Such simplification allows the usage of the Lorentz reciprocal theorem to derive an exact solution for a squirming sphere close to a no-slip surface (Papavassiliou & Alexander 2017). Exact solutions for a sphere moving towards or away from a flat wall have been derived by Brenner (1961). Cox & Brenner (1967) later derived the near-field solution that shows the asymptotic divergence of the viscous drag coefficient as the sphere approaches the solid wall. For a squirmer interacting with a solid wall through the hydrodynamic interactions, Théry et al. (2023) combined the far-field flow of a sedimenting sphere (Kim & Karrila 2013) with the squirming flow to illustrate the different swimming dynamics near a wall, and noted that the far-field approximation may not be uniformly valid across different types of squirmer dynamics. For example, oscillatory (bouncing) dynamics of a squirmer may involve motion both near and far from the wall, and the sliding squirmer can also occur in the near field (Li & Ardekani 2014; Rühle et al. 2018; Kuhr et al. 2019).

In this work, we seek to elucidate the detailed swimming dynamics of a single squirmer sedimenting towards a flat wall, using a boundary integral code validated for both the near-field and far-field hydrodynamics. In particular, we seek to quantify how the swimming dynamics of a squirmer sedimenting to a wall depends on $\alpha$ (the ratio of swimming speed to sedimenting velocity) and $\beta$ (the ratio of the first two squirming mode amplitudes). We first derive an exact solution for a squirming sphere sedimenting towards a flat no-slip wall, using the approach in Brenner (1961). We use this analytic solution to validate a boundary integral code, which is highly efficient and accurate for us to examine the dynamics of a sedimenting squirmer near a flat wall over a wide range of parameters.

This paper is organised as follows. In § 2, we present the formulation for a squirmer under gravity in the presence of a planar bottom wall. We assume that the squirmer is immersed in a viscous Stokes flow, and there may be a steric repulsion between the solid wall and the squirmer when close to the wall. We summarise the boundary integral formulation for the numerical implementation in § 2.1. We present the exact solution for a squirmer perpendicular to a flat no-slip wall in § 3. This exact solution allows us to derive an extended far-field formula, which we compare against the numerical results to validate the boundary integral code, and also examine the range of validity for the far-field approximation in § 3.1. We further study the near-field approximation to the flow, and compare between the exact solution, boundary integral simulation results, and the asymptotic results in the literature in § 3.2. In § 4, we classify the swimming dynamics of a squirmer interacting with a no-slip planar wall under gravity, and show the detailed bifurcation structures for mixed squirming modes (§ 4.1) and pure squirming modes (§ 4.2). In § 5, we provide discussions of our results and implications for future directions.

2. Problem formulation

We consider a three-dimensional incompressible viscous fluid governed by the equations of Stokes flow,

(2.1) \begin{equation} -{\boldsymbol\nabla} p + \mu _f\,{\boldsymbol\nabla} ^2\boldsymbol {u} = \boldsymbol {0}, \quad {\boldsymbol\nabla} {\boldsymbol\cdot}\, \boldsymbol {u} = 0, \end{equation}

for $\textbf { x}\in \Omega$ , the space between squirmers and a planar wall in figure 1. The squirmer is a sphere of radius $R$ located at a height $h$ above the planar wall, and has an orientation vector $\hat {\boldsymbol {e}}$ at angle $\theta$ with respect to the wall: $\theta =0$ when the squirmer is parallel to the wall, and $\theta =\unicode{x03C0} /2$ when the squirmer is upright. Furthermore, we assume that there is a density mismatch $\Delta \rho$ between the viscous fluid and the squirmer, which sediments under gravity. For spherical squirmers of constant excess density $\Delta \rho$ relative to the surrounding fluid, the excess gravitational force on the squirmer is given by

(2.2) \begin{equation} \boldsymbol {F}^{\mathit {grav}} = -F^{\mathit {grav}}\, \hat {\boldsymbol {e}}_z = -\frac {4}{3}\unicode{x03C0} R^3 g\, \Delta \rho \, \hat {\boldsymbol {e}}_z, \end{equation}

where $g$ is the constant gravitational acceleration in the $-\hat {\boldsymbol {e}}_z$ direction; see figure 1.

Figure 1. Schematic of a spherical squirmer of radius $R$ at height $h$ above a no-slip planar boundary with gravity pointing towards the wall.

We focus on fluid flow generated by the activity on the squirmer surface $S$ , and assume that the fluid flow vanishes in the far field. Furthermore, we consider the spherical squirmer with up to the first two squirming modes, prescribing the tangential velocity distribution. The prescribed surface tangential velocity distribution is purely in the polar ( $\hat {\boldsymbol{e}}_\psi$ ) direction and is given by

(2.3) \begin{equation} \boldsymbol {u}^{\mathit {s}} = u_\psi \hat {\boldsymbol {e}}_\psi , \quad u_\psi = B_1 \sin (\psi ) + B_2 \sin (\psi )\cos (\psi ), \end{equation}

where $\psi$ is the angle of the radial vector at a point on the surface to the orientation vector $\hat {\boldsymbol {e}}$ of the squirmer; see figure 1. Here, $B_1$ is the neutral, self-propelling mode with swimming speed $V = 2B_1/3$ in free space ( $V=2/3$ when $B_1 = 1$ ), and $B_2$ is the stresslet mode: $B_2 \gt (\lt )\ 0$ for contractile puller (extensile pusher) squirmers. The velocity continuity and stress balance at the squirmer boundary and the planar wall provide the boundary conditions that close the system of equations.

In the simulations, we apply a repulsive force on the squirmer at a distance $r$ (the bottom of the sphere to the wall) to the no-slip wall (Brady & Bossis 1985; Ishikawa et al. 2006):

(2.4) \begin{equation} F^{\mathit {rep}} = \frac {C^{\mathit {rep}} \exp (-a^{\mathit {rep}} r)}{1-\exp (- a^{\mathit {rep}} r)}, \end{equation}

pointing away from the wall. We found that using numerical values $C^{\mathit {rep}} = 10^3,\ a^{\mathit {rep}} = 100$ was sufficient to maintain equilibrium separation $0.04R$ for a sphere of radius $R=1$ and free space sedimentation speed $V_g=1$ with no active squirming. In most of the results presented, we used $C^{\mathit {rep}} = 10^4$ so that adequate separations could be maintained over the ranges of $B_1$ , $B_2$ and $V_g$ values considered.

2.1. Numerical algorithm

The incompressible velocity field $\boldsymbol {u}$ that satisfies the Stokes equations (2.1) can be expressed in terms of integrals of force and/or stress densities on the surfaces (Pozrikidis 1992). In particular, in the absence of a background flow, the $i$ th component of the fluid velocity at a point $\boldsymbol {x}$ outside or on the surface $S$ of a particle can be represented by a single-layer potential as

(2.5) \begin{equation} u_i(\boldsymbol {x}) = \int _S \mathbf {G}_{ij}(\boldsymbol {x},\boldsymbol {y})\,q_j(\boldsymbol {y})\,\textrm{d}S(\boldsymbol {y}), \end{equation}

where $\mathbf {G}_{ij}$ are the $ij$ th components of the Green’s function tensor $\mathbf {G}$ for Stokes flow. In free space, the Green’s function has the formula

(2.6) \begin{equation} \mathbf {G}^{\mathit {FS}}_{ij}(\boldsymbol {x},\boldsymbol {y}) = \frac {\delta _{ij}}{r} + \frac {r_i r_j}{r^3}, \end{equation}

where $\boldsymbol {r} = \boldsymbol {x}-\boldsymbol {y}$ and $r = |{\boldsymbol {r}}|$ . For simulations near a no-slip plane boundary, we use the modified Green’s function $\mathbf {G} =\mathbf {G}^{\mathit {FS}}+\mathbf {G}^{\mathit {im}}$ so that the velocity field satisfies the no-slip boundary condition on the plane $z=0$ by including image terms given by (Blake 1971a )

(2.7) \begin{equation} \mathbf {G}^{\mathit {im}}_{ij}(\boldsymbol {x},\boldsymbol {y}) = -\frac {\delta _{ij}}{\tilde {r}} - \frac {\tilde {r}_i \tilde {r}_j}{\tilde {r}^3} + 2y_3\left (\delta _{j {\nu }} \delta _{{\nu } l}-\delta _{j 3}\delta _{3 l}\right )\frac {\partial }{\partial \tilde {r}_l}\left [\frac {y_3 \tilde {r}_i}{\tilde {r}^3}-\left (\frac {\delta _{i3}}{\tilde {r}} + \frac {\tilde {r}_i\tilde {r}_3}{\tilde {r}^3}\right )\right ], \end{equation}

where $\tilde {\boldsymbol {r}} = (x_1-y_1,x_2-y_2,x_3+y_3)$ , $\tilde {r}=|\tilde {\boldsymbol {r}}|$ , and summations are implied over ${\nu } = 1,2$ and $l=1,2,3$ . In the case of a rigid body motion of the particle, the density $\boldsymbol {q}$ of the single-layer potential is proportional to the traction vector $\boldsymbol {f}$ :

(2.8) \begin{equation} \boldsymbol {q} = -\frac {1}{8\unicode{x03C0} \mu _f}\boldsymbol {f}. \end{equation}

For a squirmer with tangential surface velocity distribution $\boldsymbol {u}^{\mathit {s}}$ moving with translational velocity $\boldsymbol {U}$ and rotational velocity $\boldsymbol {\Omega }$ about its centre $\boldsymbol {x}_0$ , we have the boundary condition

(2.9) \begin{equation} \boldsymbol {U} + \boldsymbol {\Omega }\times (\boldsymbol {x}-\boldsymbol {x}_0) + \boldsymbol {u}^{\mathit {s}} = \int _S \mathbf{G}(\boldsymbol {x},\boldsymbol {y})\,\boldsymbol {q}(\boldsymbol {y})\,\textrm {d}S(\boldsymbol {y}) \end{equation}

for $\boldsymbol {x}\in S$ .

The total hydrodynamic force acting on the squirmer is given by

(2.10) \begin{equation} \boldsymbol {F}^{\mathit {hydro}} = -{8\unicode{x03C0} \mu _f}\int _S \boldsymbol {q}(\boldsymbol {y}) \,\textrm {d}S(\boldsymbol {y}), \end{equation}

and the total hydrodynamic torque is

(2.11) \begin{equation} \boldsymbol {L}^{\mathit {hydro}} = -{8\unicode{x03C0} \mu _f}\int _S \boldsymbol {y}\times \boldsymbol {q}(\boldsymbol {y}) \,\textrm {d}S(\boldsymbol {y}). \end{equation}

For squirmers that experience forces due to gravity and short-range repulsion, we impose the force balance equation

(2.12) \begin{equation} \boldsymbol {F}^{\mathit {hydro}} + \boldsymbol {F}^{\mathit {rep}} + \boldsymbol {F}^{\mathit {grav}} = \boldsymbol {0}. \end{equation}

By symmetry, gravity and short-range repulsion do not exert torques on the spherical squirmers, so the torque balance equation is

(2.13) \begin{equation} \boldsymbol {L}^{\mathit {hydro}} = \boldsymbol {0}. \end{equation}

The full system of equations to be solved at a given time consists of (2.9), (2.12) and (2.13). We use a boundary element method (BEM) to solve this system numerically (Pozrikidis 2002). The surface of the squirmer is discretised into curved triangular elements that are described geometrically by quadratic interpolations through six nodes consisting of the three vertices and midpoints along each of the three edges. The single-layer density $\boldsymbol {q}$ is likewise approximated by a quadratic interpolation of the values at the six nodes of each triangular element. In this work, we prescribe the tangential velocities $\boldsymbol {u}^s$ in (2.9) at each of the $N$ nodes on the surface of the squirmer, and the single-layer density $\boldsymbol {q}$ at the nodes are unknowns, as are the translational and rotational velocity vectors. This yields $3N$ equations in $(3N + 6)$ unknowns. Six further equations arise from the force and torque balance constraints.

3. Exact solutions for a squirmer perpendicular to a planar wall under gravity

We next assume that the squirmer orientation vector $\hat {\boldsymbol {e}}$ is normal to the planar wall, pointing towards the wall ( $\theta =-\unicode{x03C0} /2$ ). Under such axial symmetry, we compute the axisymmetric flow generated by the squirmer interacting with a planar wall using bipolar spherical coordinates defined as (see figure 2)

(3.1) \begin{align} z+i\varrho &= i c \cot \frac {1}{2}(\eta + i\xi), \end{align}

where $\varrho$ and $z$ are the cylindrical coordinates that can be expressed explicitly in terms of $\eta$ and $\xi$ :

(3.2) \begin{align} \varrho = \frac {c \sin \eta }{\cosh \xi -\cos \eta }, \;\;\; {z} = \frac {c \sinh \xi }{\cosh \xi - \cos \eta }, \end{align}

with $2c$ the distance between the two poles of the bispherical coordinates. In the present application, it is necessary to consider only the situation $\varrho \gt 0$ , which corresponds to $0\leqslant \xi \lt \infty$ and $0\leqslant \eta \leqslant \unicode{x03C0}$ . Here, $\xi =0$ corresponds to a plane, and $\xi =\zeta = \text {arcosh}(h/R)$ corresponds to a spherical surface of radius $c/\sinh \zeta$ centred at $c\coth \xi$ . In axisymmetric flow, one can write the velocity as a curl of a vector, $\boldsymbol{u} = {\boldsymbol\nabla} \times (\varphi (\xi ,\eta )\, \hat {\boldsymbol{e}}_{\phi } )$ , with $\hat {\boldsymbol{e}}_{\phi }$ the unit vector in the azimuthal direction, and $\varphi (\xi ,\eta )$ being the Stokes stream function. Accordingly, the governing differential equation for the axisymmetric viscous flow reduces to a fourth-order linear differential equation for $\varphi$ :

(3.3) \begin{align} D^4\varphi &= (D^2)^2\varphi = 0, \end{align}

where the differential operator $D^2$ in bispherical coordinates is given by

(3.4) \begin{align} D^2 \equiv \frac {\cosh \xi - \mu }{c^2}\left \{ \frac {\partial }{\partial \xi }\left [\left (\cosh \xi -\mu \right )\frac {\partial }{\partial \xi }\right ] + (1-\mu ^2)\frac {\partial }{\partial \mu }\left [\left (\cosh \xi -\mu \right )\frac {\partial }{\partial \mu }\right ]\right \}, \end{align}

with $\mu (\eta ) = \cos \eta$ . The general expression for the stream function $\varphi$ can be given in the bispherical coordinates (Stimson & Jeffery 1926)

(3.5) \begin{align} \varphi (\xi ,\eta ) &= \left (\cosh \xi -\cos \eta \right )^{-3/2}\sum ^{\infty }_{n=1}\chi _n\left (\xi \right )\,V_n\left (\mu \left (\eta \right )\right ), \end{align}

(3.6) \begin{align} \chi _n (\xi ) &= A_n\cosh \left (n-\frac {1}{2}\right )\xi + B_n\sinh \left (n-\frac {1}{2}\right )\xi \nonumber\\ &\quad + C_n \cosh \left (n+\frac {3}{2}\right )\xi + D_n\sinh \left (n+\frac {3}{2}\right )\xi , \end{align}

(3.7) \begin{align} V_n(\mu \left (\eta \right )) &= P_{n-1}(\cos \eta ) - P_{n+1}(\cos \eta ), \\[12pt] \nonumber \end{align}

Figure 2. Schematic of bispherical coordinates $(\xi , \eta )$ and the cylindrical coordinates $(\varrho , z)$ , adapted from Prakash et al. (2013).

where $P_n$ is the $n$ th-order Legendre polynomial. We can express the components of velocity in terms of the stream function as

(3.8) \begin{align} u_{\xi } &= \frac {\left (\cosh \xi -\cos \eta \right )^2}{c\sin \eta }\frac {\partial \varphi }{\partial \eta }\\ &=-\frac {3}{2c\sqrt {\cosh \xi -\cos \eta }}\sum ^{\infty }_{n=1}\chi _n\left (\xi \right )\,V_n\left (\mu \right ) + \frac {\sqrt {\cosh \xi -\cos \mu }}{c \sin \eta }\sum ^{\infty }_{n=1}\chi _n\left (\xi \right )\frac {\partial V_n(\mu )}{\partial \eta }, \nonumber \end{align}

(3.9) \begin{align} u_{\eta } &= -\frac {\left (\cosh \xi -\cos \eta \right )^2}{c\sin \eta }\frac {\partial \varphi }{\partial \xi } \\ &=\frac {3\sinh \xi }{2c \sin \eta \sqrt {\cosh \xi -\cos \eta }}\sum ^{\infty }_{n=1}\chi _n\left (\xi \right ) V_n\left (\mu \right ) - \frac {\sqrt {\cosh \xi -\cos \mu }}{c\sin \eta }\sum ^{\infty }_{n=1}\chi '_n\left (\xi \right )V_n\left (\mu \right ). \nonumber \end{align}

For each ( $n$ th) term of the stream function expansion $\varphi$ , there are four coefficients $A_n$ , $B_n$ , $C_n$ and $D_n$ to be determined as functions of the squirmer speed $U$ according to the boundary conditions of the rigid body motion and velocity continuity on the squirmer surface and the planar wall. Assuming axisymmetry (see figure 1), we use the force-free condition to compute the squirmer velocity $U$ in the downward vertical direction as

(3.10) \begin{align} U &= 2\frac {\mathcal {B}_1(\zeta )}{\lambda _B(\zeta )}B_1 - \frac {8}{3}\frac {\mathcal { B}_2(\zeta )}{\lambda _B(\zeta )}B_2 + \frac {F^{\mathit {grav}}-F^{\mathit {rep}}}{8\unicode{x03C0} \mu _f R\,\lambda _B(\zeta )}, \end{align}

(3.11) \begin{align} \lambda _B(\zeta ) &= \sinh \zeta \sum ^{\infty }_{n=1}\frac {n(n+1)}{(2n-1)(2n+3)\varDelta_n}\left [2\sinh (2n+1)\zeta + (2n+1)\sinh 2\zeta -\varDelta_n\right ], \end{align}

(3.12) \begin{align} \mathcal {B}_1\left (\zeta \right ) &= \sinh ^3\zeta \sum ^{\infty }_{n=1}\frac {n(n+1)}{\varDelta_n}\left [1-\cosh (2n+1)\zeta +\sinh (2n+1)\zeta \right ], \end{align}

(3.13) \begin{align} \mathcal {B}_2\left (\zeta \right ) &= \sinh ^2\zeta \sum ^{\infty }_{n=1}\frac {n(n+1)}{(2n-1)(2n+3)}\sinh \left (\left (n+\frac {1}{2}\right )\zeta \right )\frac {M_2(n,\zeta )}{\varDelta_n}, \end{align}

(3.14) \begin{align} M_2\left (n,\zeta \right ) &= \sinh \zeta (n-1)(n+2)\left [(2n+3)\textrm{e}^{-\left (n-{1}/{2}\right )\zeta }-(2n-1)\textrm{e}^{-\left (n+ {3}/{2}\right )\zeta }\right ] \nonumber\\ &\quad - {5}/{4}\bigg [(n-1)(2n+3)\textrm{e}^{-\left (n- {3}/{2}\right )\zeta }+(2n+1)\textrm{e}^{-\left (n+ {1}/{2}\right )\zeta } \nonumber\\ &\quad -(n+2)(2n-1)\textrm{e}^{-\left (n+ {5}/{2}\right )\zeta }\bigg ], \end{align}

(3.15) \begin{align} \varDelta_n &= 4\sinh ^2\left (n+\frac {1}{2}\right )\zeta - (2n+1)^2\sinh ^2\zeta . \end{align}

The drag coefficient for vertical motion in the presence of the no-slip wall is $c_D(\zeta )=8\unicode{x03C0} \mu _f\,\lambda (\zeta )$ , and the individual contributions to the net velocity from $B_1$ , $B_2$ , gravity and repulsion are identified as

(3.16) \begin{align} U_1 &= 2\frac {\mathcal {B}_1(\zeta )}{\lambda _B(\zeta )}B_1, \\[-10pt] \nonumber \end{align}

(3.17) \begin{align} U_2 &= - \frac {8}{3}\frac {\mathcal {B}_2(\zeta )}{\lambda _B(\zeta )}B_2,\\[-10pt] \nonumber \end{align}

(3.18) \begin{align} U^{\mathit {grav}} &= \frac {F^{\mathit {grav}}}{8\unicode{x03C0} \mu _f R\,\lambda _B(\zeta )} = V_g\, \frac {3}{4\lambda _B(\zeta )}, \\[-10pt] \nonumber \end{align}

(3.19) \begin{align} U^{\mathit {rep}} &= -\frac {F^{\mathit {rep}}}{8\unicode{x03C0} \mu _f R\,\lambda _B(\zeta )}, \end{align}

respectively, where $V_g=2 g\,\Delta\rho\, R^2/9\mu _f$ is the terminal velocity of the passive sphere in free space. We remark that in the situation where the squirmer is pointing vertically upwards ( $\theta = +\unicode{x03C0} /2$ ), the velocities (still expressed in the downward direction) are unchanged from the corresponding formulas in (3.16)–(3.19), apart from a change of sign in (3.16).

3.1. Far-field expansion of $U$ for a squirmer near a wall

The expression for the squirmer speed in (3.10) allows for a far-field expansion of $U$ in the limit $\zeta \rightarrow \infty$ . In the absence of squirming activity ( $B_1 = B_2 = 0$ ), the speed of a sphere under gravity is

(3.20) \begin{align} U^{\mathit {grav}} &=V_g \left (1-\frac {9}{8}\frac {R}{h} + \frac {1}{2} \left (\frac {R}{h}\right )^3-\frac {135}{256}\left (\frac {R}{h}\right )^4 - \frac {1}{8}\left (\frac {R}{h}\right )^5 \right.\nonumber\\ &\quad\left. + \frac {401}{512}\left (\frac {R}{h}\right )^6-\frac {675}{1024}\left (\frac {R}{h}\right )^7+\cdots \right ). \end{align}

Equation (3.20) is identical, up to ${{\mathcal{O}}} ((R/h)^3 )$ , to the often-used expression derived from the method of images with one image. The activity on the squirmer surface ( $B_1$ and $B_2$ ) contributes to the far-field squirmer speed as

(3.21) \begin{align} U &=U^{\mathit {grav}} + B_1\left (\frac {2}{3} - \frac {1}{3}\left (\frac {R}{h}\right )^3 + \frac {1}{6}\left (\frac {R}{h}\right )^5 -\frac {45}{128}\left (\frac {R}{h}\right )^6+\cdots \right ) \\ &\quad {}+ B_2 \left (\frac {3}{8}\left (\frac {R}{h}\right )^2 - \frac {1}{2}\left (\frac {R}{h}\right )^4 + \frac {15}{32}\left (\frac {R}{h}\right )^5 + \frac {5}{24}\left (\frac {R}{h}\right )^6 - \frac {441}{512}\left (\frac {R}{h}\right )^7 + \cdots \right ). \nonumber \end{align}

3.2. Near-field $U$ for a squirmer near a wall

Cox & Brenner (1967) and Cooley & O’Neill (1969) provided an expression for the near-field velocity of a sedimenting rigid sphere ( $B_1 = B_2 = 0$ ) at distance $h$ above a planar wall, with $\varepsilon = h/R-1 \ll 1$ and $\zeta = \text {arcosh} (1+\varepsilon )$ :

(3.22) \begin{align} U^{\mathit {grav}} &= \frac {\frac {4}{3}\unicode{x03C0} R^3\,\Delta \rho\, g}{6 \unicode{x03C0} \mu _f R}\frac {1}{\varepsilon ^{-1}-\left (0.2\ln \varepsilon - 0.971280\right )} \end{align}

(3.23) \begin{align} &=V_g \frac {1}{\varepsilon ^{-1}-\left (0.2\ln \varepsilon - 0.971280\right )}. \end{align}

This result shows that the infinite series $\lambda _B$ in the denominators on the right-hand-side of (3.10) has the following asymptotic behaviour as $\varepsilon \rightarrow 0$ :

(3.24) \begin{align} \lambda _B\left (\zeta (\varepsilon )\right ) &\rightarrow \overline {\lambda _B} \equiv \frac {3}{4}\left (\varepsilon ^{-1}-0.2\ln {\varepsilon }+0.971280\right ). \end{align}

This is confirmed in figure 3, showing that the truncated series solution and the boundary integral simulation results are in good agreement with the near-field asymptotic expression in (3.22).

Figure 3. Comparison of the series expression (3.10) with series truncated to 5000 terms, boundary integral simulation solutions, the near-field formula (3.22), and the far-field formula (3.20) for the vertical speed $v=U^{\mathit {grav}}/V_g$ of a passive sphere ( $B_1 = B_2 = 0$ ) sedimenting under gravity near a no-slip wall without repulsion. (a,b) Plots of the same quantities focusing over different ranges of separation from the wall. (c) Errors with respect to the series solution.

Figure 4. Comparison of the series expression (3.10) with series truncated to 5000 terms, boundary integral simulation solutions, the near-field formula (3.27), and the far-field formula (3.21) for the vertical speed of a neutral squirmer ( $B_1 = 1,\ B_2 = 0,\ V_g=0$ ) perpendicular to and near a no-slip wall without repulsion. For the near-field expansions, we used values of $a_1$ and $b_1$ in (3.28) from least squares fitting, and we used $a_1 = -0.75$ and $b_1= 0.8913$ for the asymptotic results. (a,b) Plots of the same quantities focusing over different ranges of separation from the wall. (c) Errors with respect to the series solution.

Figure 5. Comparison of the series expression (3.10) with series truncated to 5000 terms, boundary integral simulation solutions, the near-field formula (3.27), and the far-field formula (3.21) for the vertical speed of a contractile squirmer ( $B_1 = 0,\ B_2 = 1,\ V_g=0$ ) perpendicular to and near a no-slip wall without repulsion. (a,b) Plots of the same quantities focusing over different ranges of separation from the wall. (c) Errors with respect to the series solution.

For a squirmer with a prescribed slip velocity $B_1\ne 0$ and $B_2=0$ moving towards the wall ( $\theta =-\unicode{x03C0} /2$ ), Yariv (2016) provided an asymptotic expression for the near-field velocity (adapted to our variable definitions):

(3.25) \begin{align} U_1&\approx -2B_1 \varepsilon \left (\ln \varepsilon - 0.1087\right ), \end{align}

which corresponds to

(3.26) \begin{align} \mathcal {B}_1(\zeta (\varepsilon )) &\approx -\frac {3}{4}(\ln \varepsilon - 0.1087). \end{align}

Würger (2016) considered a simplifying assumption and derived an asymptotic result where the constant $-0.1087$ in (3.25) is replaced by $+2.25$ . Based on the results in (3.24)–(3.26), we construct an ansatz for the near-field behaviour of the squirmer speed in (3.10), assuming identical forms for the approximations to $\mathcal {B}_1$ and $\mathcal {B}_2$ , with $\varepsilon \ll 1$ :

(3.27) \begin{align} U &\approx \frac {1}{\overline {\lambda _B}}\left [ 2B_1a_1\left (\ln \varepsilon + b_1\right )-\frac {8}{3}B_2 a_2\left (\ln \varepsilon +b_2\right ) + \frac {3}{4}V_g\right ]. \end{align}

We then compute the coefficients $a_1, b_1, a_2, b_2$ by least squares fitting of the series solution truncated at $5000$ terms over the range $10^{-5}\leqslant \varepsilon \leqslant 10^{-3}$ to (3.27), and obtain

(3.28) \begin{equation} a_1 = -0.7476,\quad b_1 = 0.8633, \quad a_2 = 1.006, \quad b_2 = 1.825. \end{equation}

As an independent check (see Appendix A), we apply the methods in Cox & Brenner (1967) to compute the near-field asymptotic expansion for a squirmer with $B_1=1$ and $B_2=0$ , and find $a_1=-0.75$ and $b_1=0.8913$ , in good agreement with the values from the least squares fitting.

We next compare the asymptotic near-field results from both Yariv (2016) and Würger (2016) with our boundary integral simulation results and our near-field expression in (3.27) in figure 4 with $B_1=1$ and $B_2=0$ . We find that our near-field expression in (3.27)–(3.28) is close to the truncated series solution for $10^{-6}\leqslant \varepsilon \leqslant 0.05$ , and close to the boundary integral simulation for $10^{-3}\leqslant \varepsilon \leqslant 0.05$ (the boundary integral solution is not accurate for $\varepsilon \lt 10^{-3}$ ), while the asymptotic expressions from both Yariv (2016) and Würger (2016) diverge significantly for $\varepsilon \geqslant 0.02$ . For $B_1=0$ and $B_2=1$ , we compare our near-field expression with the truncated series and boundary integral simulation results in figure 5.

4. Swimming dynamics of a squirmer under gravity near a no-slip wall

4.1. Classification of squirmer dynamics in the $\alpha{-}\beta$ plane

Under gravity, various types of swimming dynamics arise from the hydrodynamic interaction between a squirmer and a no-slip flat wall (Li & Ardekani 2014; Lintuvuori et al. 2016; Rühle et al. 2018; Théry et al. 2023). We first define two parameters to quantify such diverse swimming dynamics: (i) $\alpha =V/V_g$ , the ratio of self-propulsion speed due to $B_1$ mode to the gravity-induced speed, and (ii) $\beta =B_2/B_1$ , the ratio of the two squirming mode magnitudes.

For $\alpha \gt 1$ , the squirmer is prone to escape from the wall in the long-time limit, even though its initial height and orientation may lead to a transient contact with the wall. We consider a trajectory to have escaped the wall if $h/R \gt 100$ at any positive time. For a squirmer that is bound to the wall, at least three types of swimming dynamics have been reported. (i) The squirmer is pinned close to the wall at a fixed height, pointing either towards or away from the wall, depending on the values of $(\alpha ,\beta )$ . (ii) The squirmer slides along the wall at a fixed height with a tilted orientation. (iii) The squirmer oscillates (bounces) in both height and orientation. While these near-wall swimming dynamics under gravity have been reported in the literature (Li & Ardekani 2014; Rühle et al. 2018), no detailed investigation on how these states may bifurcate in terms of $(\alpha ,\beta )$ is available to our knowledge.

From both numerical simulations and far-field analysis, we find that the squirmer is less likely to be bound to the wall if its initial orientation points away from the wall. Thus we first focus on the various dynamics of a squirmer initially pointing towards the flat wall. We use the efficient and accurate boundary integral codes to evaluate the velocities due to the two squirming modes, gravity, and repulsion (with $a^{\mathit {rep}} = 100$ ) separately on a grid in $(\theta ,h)$ configuration space. Linear combinations are then taken to construct velocity data on the grid for various values of $B_1$ , $B_2$ , $V_g$ and $C^{\mathit {rep}}$ . We then obtain trajectories from specific initial conditions using the MATLAB function stream2, which evaluates velocities by linear interpolation over the grid. Phase portraits presented below are generated similarly using the function streamslice. We first map out the various squirming dynamics in the $\alpha{-}\beta$ plane for a squirmer initially pointing nearly vertically down (specifically, the initial orientation angle is $\theta _0=-0.99\unicode{x03C0} /2$ ) at height $h_0=10R$ . The various swimming dynamics are summarised in figure 6. For pinned and sliding dynamics, we colour code their regions using the orientation angle, with $\theta ^*=\unicode{x03C0} /2$ pointing up (red) and $\theta ^*=-\unicode{x03C0} /2$ pointing down (blue). The region for oscillatory dynamics is colour coded by the oscillation amplitude ( $h^{\mathit {amp}}$ ) of squirmer height (green–yellow). We note that while oscillations can be obtained over a large region of parameter space, the amplitudes are small except near the boundary with escaping behaviour.

As we will show later, there may be multiple equilibrium states in $(\theta ,h)$ space; the state that is reached depends not only on the parameters $\alpha$ and $\beta$ , but also on the initial conditions. Figure 7 shows the distribution of swimming dynamics for a squirmer initially parallel to the wall. For a squirmer initially pointing away from the wall (vertically up), figure 8 shows that the squirmer can stay bound to the wall for a wider range of $\alpha$ for $\beta \gg 1$ , when the upward pointing squirmer moves towards the wall and stays bound to the wall for $|\alpha | \gg 1$ due to the strong puller mode $B_2 \gg B_1$ .

Figure 6. Long-time behaviours of a squirmer under gravity next to a flat wall, with the squirmer initially pointing nearly vertically downwards ( $\theta (t=0)=-0.99\unicode{x03C0} /2$ ) at starting height $h/R=10$ . The squirmer is not bound to the wall in the ‘Escape’ region. For squirmers bound to the wall under gravity, either they settle to a steady state at a fixed height with a steady tilt angle sliding along the wall (red–blue colour bar), or they oscillate in the ‘Oscillations’ region where the squirmer–wall distance oscillates with amplitudes in height indicated by the green–yellow colour bar. Negative values of $\alpha$ signify that gravity acts vertically away from the wall. The inset shows the detailed distribution of swimming dynamics for $0\lt \alpha \lt 0.15$ and $-10\lt \beta \lt -2$ .

Figure 7. Long-time behaviours of a squirmer under gravity next to a flat wall, with the squirmer initially pointing horizontally (parallel to the wall) at starting height $h/R=10$ . The squirmer is not bound to the wall in the ‘Escape’ region. For squirmers bound to the wall under gravity, either they settle to a steady state at a fixed height with a steady tilt angle (red–blue colour bar), or they oscillate in the ‘Oscillations’ region where the squirmer–wall distance oscillates with amplitudes in height indicated by the green–yellow colour bar. Negative values of $\alpha$ signify that gravity acts vertically away from the wall.

Figure 8. Long-time behaviours of a squirmer under gravity next to a flat wall, with the squirmer initially pointing nearly vertically upwards ( $\theta (t=0)=0.99\unicode{x03C0} /2$ ) at starting height $h/R=10$ . The squirmer is not bound to the wall in the ‘Escape’ region. For squirmers bound to the wall under gravity, they settle to a steady state at a fixed height with a steady tilt angle indicated by the colour bar. In regions labelled as ‘At wall’, the squirmer approaches the minimum wall separation for which velocities were computed, so trajectories could not be continued further in time. Negative values of $\alpha$ signify that gravity acts vertically away from the wall. An inset shows more detail for the range $-0.5\lt \alpha \lt 1.5$ , $-15\lt \beta \lt 15$ .

The squirmer dynamics summarised in figures 6–8 offer some general observations. (i) A squirmer stays bound to the bottom wall under gravity as long as $\alpha$ is in the range $\alpha \in [0,1]$ for all values of $\beta$ . (ii) The oscillatory dynamics of a squirmer under gravity is generally associated with negative $\beta$ (extensile squirmers). When gravity is directed away from the wall (e.g. if the wall is at the top of a chamber), however, it is possible to observe oscillatory dynamics for both contractile and extensile squirmers. (iii) For a squirmer initially pointed towards the wall, the squirmer can be pinned facing down, albeit in a very small region (inset of figure 6). (iv) The squirmer does not bounce (oscillate) around the wall if it initially points away from the wall. (v) For sufficiently large $\beta$ , we expect that a squirmer can be bound to the wall, independent of squirmer’s initial configuration and the value of $\alpha \gt 0$ .

We next examine the detailed bifurcation structures of squirmer swimming dynamics as a function of $\beta$ for various fixed values of $\alpha$ . An example is shown in figure 9, where $\alpha = 2/3$ and we vary $\beta =B_2/B_1$ with $B_1=1$ and $C^{\mathit {rep}}=10^4$ . We quantify the swimming dynamics via the equilibrium states (the fixed points) in the flow map of the squirmer in the plane of height versus angle. In the top row of figure 9, the red curves are for stable spirals, and blue curves are for stable nodes in the flow map. A stable spiral and a stable node are found for $\beta =-10$ in figure 9(a), and $\beta =7.5$ in figure 9(c). In figure 9(b), only one stable node is found for $\beta =-1$ . The bottom row of figure 9 summarises the bifurcation structures for $\beta =B_2/B_1 \in [-15,15]$ , with the dotted curves denoting the unstable branches.

Figure 9. Bifurcation structures of swimming dynamics of a squirmer with $\alpha =2/3$ . Top row: (a) $\beta =-10$ , with three saddles (green crosses), a stable node (blue filled circle) and a stable spiral; (b) $\beta =-1$ , with two saddles and a stable node; and (c) $\beta =7.5$ , with three saddles, a stable node and a stable spiral. Bottom row: dependence of (d) stable height, (e) orientation and (f) wall-parallel speed on $\beta$ . Red curves indicate stable spirals, and blue curves indicate stable nodes. Green dashed curves indicate saddle nodes. Different curve thicknesses are used to visually distinguish branches of stationary points.

Figure 10. Bifurcation structures of swimming dynamics of a squirmer with $\alpha = 5$ . Top row: (a) $\beta =-12$ , with two saddles (green crosses), a stable node (blue filled circle) and an unstable spiral (star) enclosed by a limit cycle; (b) $\beta =-2$ , with a saddle and a stable spiral (for clarity, only a portion of a trajectory attracted to the spiral is shown); and (c) $\beta =5$ , with four saddles, a stable node (blue filled circle) and a stable spiral. Bottom row: dependence of (d) stable height, (e) orientation and (f) wall-parallel speed on $\beta$ . Red curves indicate stable spirals, and blue curves indicate stable nodes. Green dashed curves indicate saddle nodes. Different curve thicknesses are used to visually distinguish branches of stationary points.

Figure 11. Oscillations in height as $\beta$ is varied for $\alpha = 5$ . The bold dashed and solid curves are respectively the unstable and stable spiral branches in figure 10(d). The shaded region around the unstable branch indicates the extent of oscillations that emerge near the unstable spiral point. Double-headed arrows indicate the minimum and maximum heights over the oscillation, the length of the arrows representing $2h^{\mathit {amp}}$ . For $\beta$ between the dashed lines, the squirmer escapes from the wall and we do not compute periodic oscillations. The insets show trajectories for particular initial configurations in the cases $\beta = -1$ , which converges to the stable spiral point, and $\beta = -10$ , which converges to a stable limit cycle around the unstable spiral point.

Figure 12. Effect of repulsion on the swimming dynamics of a squirmer with $\beta =-10$ $(B_1=1,\ B_2 = -10)$ and $\alpha = 5$ . Phase plane dynamics near an unstable spiral point (a) with no wall repulsion, and (b) with short-range wall repulsion given in (2.4) using repulsion strength $C^{\mathit {rep}}=10^3$ , and (c) similarly with $C^{\mathit {rep}}=10^4$ . Red curves show trajectories starting close to the spiral point. The red dot in (a) (with no repulsion) indicates where the trajectory terminates due to close proximity with the wall. With repulsion in (b) and (c), the trajectory approaches a limit cycle, in which the distance from the wall and the orientation of the squirmer oscillate periodically. Trajectories converging to the limit cycle from outside are also shown (green) in (b) and (c).

At a higher squirmer propulsion speed (higher value of $\alpha$ ), we expect the region for stable nodes to shrink. An example with $\alpha =5$ is shown in figure 10. In the top row, the phase plane diagrams are shown for $\beta =-12$ in figure 10(a), $-2$ in figure 10(b), and $5$ in figure 10(c). For $\beta =-12$ , squirmers that approach the wall are attracted either to a stable node with a negative tilt or to a limit cycle with positive tilt; both of these attractors are very close to the wall. For $\beta =-2$ , most initial conditions lead to escape from the wall, but some are attracted to a stable spiral. For $\beta =5$ in figure 10(c), the squirmer either spirals into a fixed point close to the wall or escapes, depending on the initial angle. There is an additional stable node close to the wall that cannot be reached by trajectories that begin far from the wall. The bifurcations with respect to $\beta$ are summarised in the bottom row of figure 10. We find that the branch of stable nodes that was present for $\alpha =2/3$ (blue curves in figure 9) at intermediate values of $\beta$ disappears at $\alpha =5$ . In fact, we find no stable equilibrium points for $\beta \in (-9.5,-2.07)\cup (0,3.5)$ . For $\beta \lt -2.07$ , there is an unstable spiral around which we expect a limit cycle (as shown in figure 10 a), and for $\beta \in (0,3.5)$ , the generic behaviour is to escape from the wall.

We illustrate the swimming dynamics of a squirmer corresponding to the unstable spiral branch (dashed red curve) of figure 10(d) in figure 11, where squirmer oscillations in height are plotted against $\beta$ . An example is shown for a spiral point ( $\beta =-1$ ) and for an unstable spiral point ( $\beta = -10$ ). When a squirmer starts near a stable spiral point, oscillations in height and orientation decay, and the squirmer approaches a steady state. Near an unstable spiral point, the squirmer tends to a stable limit cycle with periodic oscillations in height and orientation. Furthermore, we observe that the oscillation amplitude in height increases as $\beta$ approaches the values at the dashed lines, and eventually the squirmer exceeds a distance $100R$ from the wall, at which point we characterise the outcome as ‘escape from the wall’ (for $\beta$ between the vertical dashed lines in figure 11). We include two supplementary movies to illustrate the oscillation of squirmers with $(\alpha ,\beta )=(5, -2.25)$ (movie 1) and $(\alpha ,\beta )=(5, -10)$ (movie 2).

Figure 12 shows the effect of wall repulsion on the spiralling squirmer dynamics, which transitions from (a) an unstable spiral that eventually becomes too close to the wall for numerical solutions to continue, into (b) a limit cycle, giving rise to oscillatory dynamics of a squirmer in both height and orientation due to a wall repulsion that is strong enough to maintain a minimum separation between the wall and the squirmer.

4.2. Effects of gravity (varying $V_g$ ) for a pure squirmer

Here, we investigate the swimming dynamics of a pure squirmer (either a ‘shaker’ with $B_1=0$ , $B_2 \neq 0$ or a ‘neutral swimmer’ with $B_1\neq 0$ , $B_2=0$ ) under gravity, focusing on three combinations of $(B_1, B_2)$ – $(1,0)$ for a neutral swimmer in figure 13, $(0,1)$ for a contractile shaker (puller) in figure 14, and $(0,-1)$ for an extensile shaker (pusher) in figure 15 – and examine the swimming dynamics as a function of $V_g$ . We exclude short-range repulsion to examine purely hydrodynamic interactions with the wall.

For a neutral swimmer with $B_1=1$ and $B_2=0$ , we find that it can be bound to the wall and stay at an equilibrium height, pointing towards the wall with zero velocity ( $v^*=0$ ) as long as the gravity is sufficiently large so that $V_g \gt 2/3$ ; see figure 13.

Slightly more complicated swimming dynamics is found for a contractile puller shaker ( $B_1=0$ and $B_2=1$ ): we find a branch of saddle nodes where the squirmer points to the wall with zero velocity at a finite height. Co-existent with this saddle node is a stable spiral with a tilt angle and a sliding velocity for intermediate gravity. For sufficiently large $V_g$ , the purely extensile squirmer is bound to the wall at a fixed height, pointing to the wall with zero velocity. Furthermore, we find that for a contractile puller shaker under a sufficiently large $V_g$ , the squirmer can reach an equilibrium height with its director parallel to the wall and a zero sliding velocity; see figure 14(e,f), where $\theta ^*=0$ and $v^*=0$ for $V_g \geqslant 1$ .

The bifurcation structure of the swimming dynamics of an extensile pusher shaker is summarised in figure 15, where we find that the steady state ( $v^*=0$ for large $V_g$ ) is a squirmer at a steady equilibrium height while pointing upright; see figure 15(e,f), where $\theta ^*=\unicode{x03C0} /2$ and $v^*=0$ for $V_g \gtrsim 1.6$ .

Figure 13. Threshold in $V_g$ for a steady equilibrium height and angle for a squirmer with $B_1=1,\ B_2=0$ (neutral swimmer) without wall repulsion ( $C^{\mathit {rep}}=0$ ). (a) Escaping dynamics for a squirmer with $V_g=0.5$ . (b) A squirmer with $V_g=1$ reaches a steady equilibrium height and points upwards. (c) Dependence of stable height in perpendicular-up orientation on the free space sedimentation speed $V_g$ . Dashed vertical lines indicate the values of $V_g$ used in (a) and (b), respectively.

Figure 14. Bifurcation structures of swimming dynamics of a squirmer with $B_1=0,\ B_2=1$ (contractile shaker) without wall repulsion ( $C^{\mathit {rep}}=0$ ). Top row: (a) $V_g=0$ , (b) $V_g=0.5$ , and (c) $V_g=1$ . Bottom row: dependence of (d) stable height, (e) orientation and (f) wall-parallel speed on the free space sedimentation speed $V_g$ . Red curves indicate stable spirals, and blue curves indicate stable nodes. Green dashed curves indicate saddle nodes.

Figure 15. Bifurcation structures of swimming dynamics for a squirmer with $B_1=0,\ B_2=-1$ (extensile shaker) without wall repulsion ( $C^{\mathit {rep}}=0$ ). Top row: (a) $V_g=0.5$ , (b) $V_g=1.5$ , and (c) $V_g=3$ . Bottom row: dependence of (d) stable height, (e) orientation and (f) wall-parallel speed on the free space sedimentation speed $V_g$ . Red curves indicate stable spirals, and blue curves indicate stable nodes. Green dashed curves indicate saddle nodes.

5. Discussion and conclusion

In this work, we provide an exact solution for a spherical squirmer sedimenting to a flat solid wall. We provide both far-field and near-field approximations to the swimming velocity of a squirmer under gravity, and show that our near-field approximations, different from both Yariv (2016) and Würger (2016), are valid over a wider range of squirmer distances to the wall. We next use boundary integral simulations to map out its various swimming dynamics in the $\alpha{-}\beta$ plane, and find that the squirmer may escape from the wall, slide along the wall at a fixed height and orientation, stay at a fixed height pointing towards or away from the wall, or oscillate in both height and orientation.

We further examine the bifurcations in the steady-state configurations of the squirmer interacting with a solid wall as the parameters $\alpha$ and $\beta$ are varied, identifying branches of stable and unstable spirals, stable nodes, and saddle nodes in $(\theta,h)$ phase space. In particular, we find that there are parameter regions where different swimming dynamics co-exist (overlaying solid branches for stable spiral and stable node in figures 9 and 10). Such identification allows us to characterise when the squirmer will escape from the wall, remain bound to the wall at a fixed location, slide along the wall at a fixed height, or bounce along the wall, making it possible to design robotic microswimmers that can adjust their gaits (by changing the values of $B_1$ and $B_2$ ) to navigate along the solid wall in the presence of obstacles.

We note that previous work by Ishimoto & Gaffney (2013) predicted that hydrodynamic interactions with the wall do not allow for stable near-surface swimming of spherical pushers ( $\beta \lt 0$ ) in the absence of gravity, because they found a single fixed point that is unstable for pushers, and stable for pullers. In contrast, we have found that pullers with sufficiently large $\beta$ have an unstable spiral point and a stable one, even without gravity and repulsion. By time-reversal symmetry, a pusher with the same $|\beta |$ has the same fixed points but with opposite stability properties. In particular, pushers do have stable configurations near the wall. We further show that this is also possible in the presence of gravity.

Lintuvuori et al. (2016) showed that by including a short-ranged repulsion between the pusher and the wall, stable oscillations could be achieved by pushers. Our results shown in figure 12 provide further explanation for the emergence of this stable, oscillatory bound state: it occurs when there is an unstable spiral point around which trajectories spiral outwards until the squirmer collides with the wall in the absence of wall repulsion, as observed by Ishimoto & Gaffney (2013). Repulsion stabilises the trajectory by preventing the squirmer from approaching too closely to the wall; see figure 12. While limit cycles are observed only when wall repulsion is included, we note that periodic oscillations (that are not limit cycles) are possible without repulsion in the special case of pure extensile or contractile squirmers ( $B_1=0,\ B_2\ne 0$ ) without gravity ( $V_g=0$ ).

Overall, we find that the diverse swimming states depicted in figures 6–8 generally align with prior observations in Rühle et al. (2018) and Théry et al. (2023). Specifically, puller squirmers can stably bind at a tilted angle for large $\beta$ , while maintaining an upright orientation for small $\beta$ . Similarly, pusher squirmers can bind at a tilted angle for small to intermediate $\beta \lt 0$ , and adopt a downward orientation for highly negative $\beta$ values. The emergence of oscillations was predicted by Théry et al. (2023), based on the characteristics of the far- and near-field behaviours, for sufficiently strong pushers and for sufficiently weak gravity. We have verified that this is, indeed, the case, and numerically determine the parameter regime for this to occur. Interestingly, we identified a new regime in which pullers exhibit stable oscillations – when gravity acts away from the wall and $|\beta |$ is sufficiently large, both pushers and pullers have periodic orbits.

Another difference that we observe compared with the near-field results in Rühle et al. (2018) and Théry et al. (2023) is the existence of multiple stable equilibrium states for a given set of squirmer parameters. These equilibria are typically close to the wall. Additional tests (not shown) confirm that in some cases, the multiple stable equilibria persist without wall repulsion and are therefore hydrodynamic in origin, while in other cases, a stable equilibrium disappears when repulsion is removed; the squirmer approaches the wall until its distance to the wall is below the minimum separation and the solutions cannot be continued.

Interactions between particles and walls in complex biological fluids are significantly influenced by non-Newtonian rheological properties. Biological fluids such as blood and mucus, commonly exhibiting shear-thinning behaviour, markedly affect the locomotion of microswimmers. Recent studies have explored novel particle–wall interactions in shear-thinning, non-Newtonian fluids. For instance, Chen et al. (2021) investigated the dynamics of a rotating sphere rolling near a solid wall, demonstrating that the shear-thinning rheology can induce wall-driven translation counter to the direction of frictional forces along the wall. Similarly, Li & Ardekani (2017) examined how non-Newtonian effects alter the interaction of an undulatory swimmer adjacent to a wall. Exploring the near-field flow dynamics of a squirmer in such fluids remains a compelling avenue for future research. It remains an interesting subject to examine the near-field flow of a squirmer in a non-Newtonian fluid.

In porous media, the interactions between active swimmers and complex boundaries play a crucial role in determining the diffusive transport of active suspensions, such as bacteria undergoing run-and-tumble dynamics (Datta et al. 2024). Our findings reveal that under gravity, a squirmer exhibits multiple trajectory states, either remaining bound to a wall or escaping from it, depending on variations in $\alpha$ (the ratio of propulsion velocity to sedimentation velocity) and $\beta$ (the relative strength of swimming to contractile or extensile modes). For a squirmer sedimenting in a porous medium, each encounter with an obstacle can induce a state transition analogous to the run-and-tumble behaviour of bacteria in similar environments (Mattingly 2023). A quantitative investigation into how the effective diffusivity of a squirmer sedimenting in porous media varies with $\alpha$ and $\beta$ remains an intriguing direction for future research.