%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
Swimming microorganisms live in a variety of fluid environments. While locomotion in Newtonian environments at low Reynolds number is relatively well established (Lauga Reference Lauga2016), swimming in complex fluids is an area of active research due to the relevance of biological fluids, which often have highly non-Newtonian properties. Understanding the role of the non-Newtonian fluid properties on swimming performance is important for the design of artificial swimmers (Bunea & Taboryski Reference Bunea and Taboryski2020; Tsang et al. Reference Tsang, Demir, Ding and Pak2020; Wu et al. Reference Wu, Chen, Mukasa, Pak and Gao2020).
Helicobacter pylori is an infectious microorganism that moves through dense, gel-like gastric mucus by actively changing its local environment to enable locomotion (Bansil et al. Reference Bansil, Celli, Hardcastle and Turner2013). A gastric mucus and other gel-like materials, unlike Newtonian fluids, can have solvent–gel transition and, in general, do not follow a single canonical fluid model. Because the fluid region near H. pylori is understood to be easier to swim through than the outer gel, one model of this situation combines an inner Stokes fluid with either an outer Stokes fluid with a different viscosity (Reigh & Lauga Reference Reigh and Lauga2017) or an outer Brinkman fluid (Nganguia et al. Reference Nganguia, Zhu, Palaniappan and Pak2020). Similar to many others, these studies used an idealized model of a swimmer as a squirmer (Lighthill Reference Lighthill1952; Blake Reference Blake1971; Pedley Reference Pedley2016). The squirmer model has been employed to study the locomotion in various complex fluids (Lauga Reference Lauga2009; Zhu, Lauga & Brandt Reference Zhu, Lauga and Brandt2012; Li, Karimi & Ardekani Reference Li, Karimi and Ardekani2014; Datt et al. Reference Datt, Natale, Hatzikiriakos and Elfring2017; Pietrzyk et al. Reference Pietrzyk, Nganguia, Datt, Zhu, Elfring and Pak2019; Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020; van Gogh et al. Reference van Gogh, Demir, Palaniappan and Pak2022). When both the swimmer and the inter-fluid boundary are spherical or ellipsoids, these models are analytically tractable and consequently can be used to produce exact results for a wide range of parameters. However, these models suffer from the deficiency that the confinement is a fixed, independent parameter. A more realistic model would allow the confinement radius to be induced by the swimming body itself, through either surface kinematics or chemical reaction. In this paper, we ignore any chemistry-induced changes to the fluid and focus purely on hydrodynamic effects. The role of chemistry-induced fluid changes on locomotion is discussed in Mirbagheri & Fu (Reference Mirbagheri and Fu2016).
To study the above effects, a specific model is needed that intrinsically links the fluidity in the near field to the fluid boundary conditions, i.e. the degree of confinement is not an independent parameter. This will lead to more natural results as would be encountered by a typical swimmer such as H. pylori. One canonical model that achieves this goal is a Bingham viscoplastic fluid subject to yield stress. This means that the material behaves as a plastic solid – capable of rigid translation but not deformation – when the interior stress is below a certain value, while it is allowed to deform when the stress exceeds this value. Common examples include hair gel and toothpaste, where the material can maintain its shape indefinitely under low shear but flows easily once adequate shear is applied. While in their most complex state, viscoplastic models can include the effect of shear-thinning or -thickening, i.e. the Herschel–Bulkley model; here, we use a Bingham fluid model (Saramito Reference Saramito2016) for simplicity while retaining the key viscoplastic component relevant to our biological motivation. Previous studies on swimming in Bingham fluids have been confined to two dimensions (Hewitt & Balmforth Reference Hewitt and Balmforth2017; Supekar, Hewitt & Balmforth Reference Supekar, Hewitt and Balmforth2020) or for slender bodies (Hewitt & Balmforth Reference Hewitt and Balmforth2018). Here, we extend these studies and consider an idealized steady spherical squirmer that employs tangential motions on its boundary.
A closely related problem is the locomotion in shear-thinning and viscoelastic fluids (Datt et al. Reference Datt, Zhu, Elfring and Pak2015; Elfring & Lauga Reference Elfring and Lauga2015; Li, Lauga & Ardekani Reference Li, Lauga and Ardekani2021; Wu et al. Reference Wu, Solano, Shoele and Mohammadigoushki2022). By employing the squirmer model, it is shown that while a swimmer slows down in a shear-thinning fluid, it could have higher swimming efficiency, and its maximum swimming efficiency occurs at a particular boundary actuation rate (Nganguia, Pietrzyk & Pak Reference Nganguia, Pietrzyk and Pak2017). A similar conclusion is reached for a wide range of spheroidal shape swimmers, and it is found that spheroidal swimmers compared to spherical squirmers can reach higher swimming velocity and energetic efficiency in shear-thinning fluids (van Gogh et al. Reference van Gogh, Demir, Palaniappan and Pak2022). In addition, it is shown that the nonlinear rheology of the fluid can modify the efficient swimming gaits in non-Newtonian fluids compared to Newtonian fluids, and allows non-motile surface actuation modes in Newtonian fluids to generate net motion in shear-thinning fluids (Pietrzyk et al. Reference Pietrzyk, Nganguia, Datt, Zhu, Elfring and Pak2019). The changes in the swimming performance are associated with the relative importance of local and non-local effects in shear-thinning fluids. Such effects are more complex in yield fluids as the flow properties undergo rapid transition across the yield surface, which itself is formed based on the swimmer geometry and mode of actuation. To date, it is still unclear how the hydrodynamics of the squirmer model and its swimming kinematics are modified in yield fluids. This constitutes the main objective of this work.
The paper is organized as follows. The governing equations of a squirmer in a Bingham fluid are described in § 2. The computational model is described in § 3. Section 4 contains results for treadmill and higher-order modes. We conclude with § 5 by describing what these results suggest about the nature of swimming in viscoplastic materials, particularly commenting on the possibility of different optimal swimming strategies that employ higher-order modes.
2. Governing equations
We employ the regularized Bingham viscoplastic Stokes flow model. The regularization is adequate to model the dynamics of a small squirmer, on the length scale of microns, swimming through a yield stress fluid. At these scales, diffusion results in a smooth transition between yielded and unyielded regions. Yet, for a sufficiently small smoothing parameter, as will be explained later, the yield surface approaches the sharp interface encountered in the singular Bingham formulation.
Mathematically, we consider an exterior flow problem of a steady squirmer swimming in an infinite viscoplastic fluid given by the regularized Bingham model (Saramito Reference Saramito2016; Saramito & Wachs Reference Saramito and Wachs2017), written in a non-dimensional form as
\begin{equation} -\text{div}\,\boldsymbol{\mathsf{\tau}} + \boldsymbol{\nabla} p = \boldsymbol{0},\quad \boldsymbol{\mathsf{\tau}} = 2\eta(\dot{\gamma})\,\boldsymbol{\mathsf{D}}(\boldsymbol{u}),\quad \text{div}\,\boldsymbol{u}=0, \end{equation}where
\begin{equation} \eta(\dot{\gamma}) = 1 + \textit{Bi}\,(\dot{\gamma}^2+\epsilon^2)^{{-1}/{2}} \end{equation}
is the effective viscosity, which depends on the shear rate 
$\dot {\gamma }=|2\boldsymbol{\mathsf{D}}|$. Here, 
$p$ is the pressure, 
$\boldsymbol{\mathsf{\tau}}$ and 
$\boldsymbol{\mathsf{D}}$ are the deviatoric stress and rate-of-strain tensors, respectively, 
$\boldsymbol {u}$ is the fluid velocity, and 
$\epsilon$ is a small regularization parameter. The tensor norm is defined as 
$|\boldsymbol {\delta }|^2=(\boldsymbol {\delta }\boldsymbol {:}\boldsymbol {\delta })/2$ for any tensor 
$\boldsymbol {\delta } \in \mathbb {R}^{3\times 3}$.
The Bingham number 
$\textit {Bi}{}= \tau _y R / (\mu _0 V)$ is employed as the key characteristic parameter of the non-dimensional equations, where 
$R$ is the length scale, taken to be the radius of the squirmer, 
$\tau _y$ is the yield stress, 
$\mu _0$ is the characteristic (Newtonian) fluid viscosity, and 
$V$ is some characteristic velocity scale. Because we consider swimming problems where the swimmer might not achieve motion, we let 
$V=\sqrt {\mathcal {P}_N / \mu _0 R}$, where 
$\mathcal {P}_N$ is the power expenditure of the squirmer in the Newtonian case (Blake Reference Blake1971; Michelin & Lauga Reference Michelin and Lauga2011). The yield surface is defined to be associated with 
$|\boldsymbol{\mathsf{\tau}} | = \textit {Bi}{}$. We assume perfect viscoplasticity, and the rigid region is placed at 
$|\boldsymbol{\mathsf{\tau}} | < \textit {Bi}{}$ (figure 1).
Schematic of an axisymmetric squirmer in a Bingham fluid environment. The diffuse yield surface (
$|\boldsymbol{\mathsf{\tau}}| = \textit {Bi}$) is given by the dashed line. The fluid region exterior to this line has negligible rate of deformation, i.e. 
$\dot {\gamma }\ll 1$.
The non-dimensional power 
$\mathcal {P}$ expended by the squirmer into the fluid is defined to be 
${\mathcal {P} = -\int _\mathcal {S}\boldsymbol {u}^S\boldsymbol {\cdot }( \hat {\boldsymbol {n}} \boldsymbol {\cdot }\boldsymbol {\sigma })\, {\rm d}S}$, where 
$\boldsymbol {\sigma } = -p\boldsymbol{\mathsf{I}} + \boldsymbol{\mathsf{\tau}} $ is the total stress on the body 
$\mathcal {S}$, and 
$\hat {\boldsymbol {n}}$ is the unit normal vector, and 
$\boldsymbol{\mathsf{I}}$ is the identity matrix. The boundary conditions in the body-attached coordinate are 
$\boldsymbol {u}(|\boldsymbol {x}|=1) = \boldsymbol {u}_s$ and 
$\boldsymbol {u}(\boldsymbol {x}\to \infty ) = -U \boldsymbol {e}_x$, where 
$\boldsymbol {u}_s$ is the tangential squirming surface motion, and 
$U$ is the translational swimming velocity of the squirmer in the 
$\boldsymbol {e}_x$ direction.
The surface motion 
$\boldsymbol {u}_s$ consists of orthogonal ‘stroke modes’ with the surface velocity of 
$\boldsymbol {u}_s = u_\theta \boldsymbol {e}_\theta$, where 
$\boldsymbol {e}_\theta$ is the unit tangent vector on the surface of the sphere. The surface velocity magnitude 
$u_\theta$ is represented by (Blake Reference Blake1971)
\begin{equation} u_\theta(\nu) = \sum_{n=1}^\infty \beta_n\,K_n(\nu),\quad K_n(\nu) = \sqrt{\frac{3}{n(n+1)}}\sqrt{1-\nu^2}\,P_n'(\nu), \end{equation}
where 
$\nu = \cos \theta$. Here, 
$\theta$ is is the polar angle with respect to the swimming direction. For each orthogonal mode 
$n$, 
$P_n'(\nu )$ is the first derivative of the 
$n{\text {th}}$-order Legendre polynomial. The formulation permits us to fully characterize the swimming stroke using the coefficients 
$\beta _n$, and in the case of a swimmer in the Newtonian fluid, to correlate the swimming velocity to the first squirming stroke mode as 
$U = \sqrt {2/3}\,\beta _1$. The choice of non-dimensionalization based on the rate of energy dissipation of Newtonian swimmers leads to a constraint relation between 
$\beta _n$ coefficients (Blake Reference Blake1971; Michelin & Lauga Reference Michelin and Lauga2011). In particular, the stroke modes 
$\beta _n$ are subject to the constraint
\begin{equation} \frac{2}{3}\,\beta_1^2 + \sum_{n=2}^\infty\beta_n^2 = 1. \end{equation}While in the case of a classic Stokes fluid, this fixes the power expenditure of the squirmer (Blake Reference Blake1971), we will see that this is not the case for swimmers in a Bingham fluid. Nevertheless, this choice is made to enable us to cross-compare different swimmers based on their energy expenditure and to study both swimming and non-swimming strokes.
3. Numerical method
A continuous Galerkin finite element method (FEM) is used to calculate numerically the flow around a free-swimming steady squirmer (Eastham Reference Eastham2019). The quadrilateral mesh is used, and the simulations with 
$\textit {Bi}{}=0$ are done with outer radius 20 000. All simulations with 
$\textit {Bi}{}>0$ are performed with outer radius 20; the presence of the yield surface means that velocity field deformations decay to negligible amounts rapidly, so more elements can be included near the swimmer to capture the yield surface. We use a Picard iteration scheme to solve the nonlinear system (2.1a–c). The procedure starts with the assumption 
$\textit {Bi}=0$ and with the choice of larger regularization parameter. At each iteration, the fluid equation (2.1a–c) can be interpreted as a variable-viscosity Stokes equation. The system is solved iteratively until convergence with absolute difference tolerance 
$10^{-6}$. Following the convergence, the yield surface is identified and the mesh is resolved near the interface. Then the solutions with higher 
$\textit {Bi}{}$ and smaller regularization parameter are recalculated, and the procedure continues until the final convergence criterion of 
$10^{-6}$ is achieved. The axisymmetric weak formulation at the 
$n{\text {th}}$ iteration is
\begin{gather} \int_D2\eta\left(\dot{\gamma}^{n-1}\right)\left[\boldsymbol{\mathsf{D}}(\boldsymbol{u}^n):\boldsymbol{\mathsf{D}}(\boldsymbol{v}) + \frac{u_r^n v_r}{r^2}\right]r\,{\rm d}S - \int_D{p^n(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}) r}\,{\rm d}S= 0, \end{gather}
\begin{gather}\int_D q(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}^n) r \,{\rm d}S = 0, \end{gather}
where 
$D$ is the axisymmetric domain, 
$r$ is the radial variable, 
$(\boldsymbol {u}^n, p^n)$ refers to unknowns to be solved for, and 
$\dot {\gamma }^{n-1}$ is assumed to be known from the previous iteration's solution. The test functions 
$(\boldsymbol {v}, q)$ are chosen to be Taylor–Hood 
$Q_2$–
$Q_1$ shape functions. More details on the axisymmetric weak form for the Stokes equation and the validation of the model can be found in Eastham & Shoele (Reference Eastham and Shoele2020). The yield surface is identified by applying a bisection method radially to the scalar field 
$|\boldsymbol{\mathsf{\tau}}| - \textit {Bi}{}$ for different 
$\theta$ values.
4. Results
We investigate a treadmill squirmer in a Bingham fluid, and compare the results to those for an equivalent squirmer in a Newtonian fluid. Then we explore the combinations of surface modes, and analyse the qualitative differences between even and odd modes. Our goal is to understand the performance of the squirmer relative to its performance in a Newtonian fluid (
$\textit {Bi}{}=0$). Therefore, we define a normalized power as 
$\mathcal {P}^* = \mathcal {P} / \mathcal {P}_N$, where the 
$N$ subscript indicates the power expended by the same squirmer in a Newtonian fluid environment. Because we consider ‘purely mixing’ stroke modes with 
$\beta _1=0$, where 
$U_N=0$, we do not normalize the swimming speed 
$U$.
4.1. Treadmill squirmer
We first simulate a treadmill squirmer, i.e. 
$\beta _i = \sqrt {3/2}\,\delta _{i1}$ where 
$\delta _{ij}$ is the Kronecker delta. In Newtonian fluids, this is the only mode that generates locomotion, while all other modes, generally termed ‘mixing modes’, do not. We compute both the swimming speed and power expenditure for Bingham number 
$\textit {Bi}{}$. The value of 
$\epsilon$ is selected to be sufficiently small such that it does not modify the trends observed here. In particular, it is found that the results are insensitive to the choice of 
$\epsilon$ if 
$\epsilon <0.01$. This can also be seen from figures 2 and 3. In all cases, we found that the viscoplastic squirmer requires more power and swims slower than its Newtonian counterpart (figures 2a,b). The decrease of the swimming velocity is much faster at the intermediate range 
$\textit {Bi}{} \approx 1$. For a higher 
$\textit {Bi}{}$ range, the power expenditure increases exponentially with 
$\textit {Bi}{}$. To compare the efficiency of the swimmer in the viscoplastic and Newtonian fluids, we use the standard definition of swimming efficiency proposed by Lighthill (Reference Lighthill1975), and define the efficiency 
$\mathcal {E}$ as
\begin{equation} \mathcal{E}=\frac{U\mathcal{D}}{\mathcal{P}}, \end{equation}
where the power needed to translate a rigid sphere at the same swimming speed, 
$U$, is compared to the power dissipation 
$\mathcal {P}$ in the viscoplastic fluid case. Here, 
$\mathcal {D}$ is the force required to drag a rigid sphere at the swimming speed in the same fluid medium, calculated from separate simulations with the imposed swimming velocity. This definition of swimming efficiency has been adopted widely to characterize the efficiency of different swimmers (Stone & Samuel Reference Stone and Samuel1996; Tam & Hosoi Reference Tam and Hosoi2007; Michelin & Lauga Reference Michelin and Lauga2010; Nganguia & Pak Reference Nganguia and Pak2018) and is adopted here. A similar expression can be employed to find the efficiency of a swimmer in Newtonian fluid, 
$\mathcal {E}_N$. It is possible to show mathematically that for spherical squirmers with the pure treadmill motion, 
$\mathcal {E}_N=\frac {1}{2}$ (Lighthill Reference Lighthill1952; Blake Reference Blake1971).
(a) Swimming speed 
$U$, (b) normalized power expenditure 
$\mathcal {P}^*$, and (c) ratio of the swimming efficiency compared with Newtonian efficiency 
$\mathcal {E}_N$, of a squirmer with 
$\beta _1=\sqrt {3/2}$ for three different regularization parameters 
$\epsilon =0.001,0.005,0.01$.
The velocity (streamlines) and pressure fields (colours) for (a,b) 
$\textit {Bi}{}=0.1$, (c,d) 
$\textit {Bi}=1.0$ and (e,f) 
$\textit {Bi}{}=10$, and two regularization values, 
$\epsilon =0.001$ and 0.005. The red dashed line represents the yield surface.
Figure 2(c) shows the efficiency of a squirmer as a function of 
$\textit {Bi}{}$. It is observed that for 
$\textit {Bi}{}<10$, the efficiency of the squirmer in the yield flow is higher than that of the Newtonian swimmer, attaining its maximum value 
$\mathcal {E}/\mathcal {E}_N=1.4$ at 
$\textit {Bi}{}=1.0$. Only for very large yield limits, where the power increases exponentially, do we see lower efficiency than for the Newtonian swimmer. The results reveal that the swimmer is most efficient at 
$\textit {Bi}{}=O(1)$. The increase of the hydrodynamic efficiency with 
$\textit {Bi}{}$ has been observed recently for the slender wavy cylindrical filament in yield fluids (Hewitt & Balmforth Reference Hewitt and Balmforth2022). Although different from the current prediction, it was found that the efficiency increases monotonically with an increase of 
$\textit {Bi}{}$ for wavy filaments.
The flow fields are visualized in figures 3(c,d), where the yield surface is shown with the red line for three yield limits of 
$\textit {Bi}{}=0.1$ (figures 3a,b), 
$\textit {Bi}{}=1$ (figures 3c,d) and 
$\textit {Bi}{}=10$ (figures 3e,f). The streamlines outside the yield region are marked in grey. Figures 3(a,c,e) use one regularization parameter 
$\epsilon$, and figures 3(b,d,f) use another. First, it is found that the shape of the yield surface is independent of the choice of 
$\epsilon$ for sufficiently small values. As 
$\textit {Bi}{}$ increases, the yield surface shrinks closer to the surface and the plasticity becomes more dominant. The decreased swimming speed of treadmill squirmers with increased confinement is similar to results from other models that involve confinement (Reigh & Lauga Reference Reigh and Lauga2017; Nganguia et al. Reference Nganguia, Zhu, Palaniappan and Pak2020). The most drastic change is related to the pressure field. The pressure becomes larger with increasing 
$\textit {Bi}{}$, mainly because of the confinement effect from the yield surface. Mathematically, this is related to rapidly varying effective viscosity 
$\eta (\dot {\gamma })$ near the yield surface in (2.1a–c). This effect, where pressure compensates for spatially varying viscosity, is seen in other complex fluid models with non-constant viscosity (Eastham & Shoele Reference Eastham and Shoele2020; Nganguia et al. Reference Nganguia, Pietrzyk and Pak2017).
The results are different from the observed behaviour in shear-thinning fluids (Nganguia et al. Reference Nganguia, Pietrzyk and Pak2017). Contrary to shear-thinning fluids and how the swimming behaviour changes with Carreau number (equivalently, the ratio of strain rate and relaxation rate of the fluid medium), the velocity now decreases continuously with 
$\textit {Bi}{}$. Moreover, unlike shear-thinning fluids, where the efficiency is always larger than in the Newtonian fluid, the swimmer efficiency in viscoplastic fluids can be higher or lower than that for the swimmer in Newtonian fluid due to a rapid increase of power expenditure at large 
$\textit {Bi}{}$ cases. Nevertheless, both viscoplastic and shear-thinning fluids show local maxima in the swimming efficiency at intermediate values of 
$\textit {Bi}{}$ or Carreau number.
Although the swimming speed and power are clearly affected by viscoplasticity, the qualitative features of the velocity field are consistent between Newtonian and non-Newtonian cases. This will change when we include certain higher-order modes. As shown in the figures, the results do not change with the choice of the regularization parameter 
$\epsilon$. Following the sensitivity analysis, 
$\epsilon =0.005$ is fixed for the remainder of this paper.
4.1.1. Surface perturbation analysis
To form a conceptual bridge between the treadmill squirmers in Newtonian and viscoplastic fluids, we approximate the nonlinear viscoplastic environment as a confined linearly viscous environment around the swimmer. We follow Reigh & Lauga (Reference Reigh and Lauga2017) and derive an analytic solution to the problem of a spherical squirmer of radius 
$R$ swimming inside a Stokes fluid confined by a concentric viscoplastic yield surface, which is a perturbed spherical surface. In the following, we investigate how small surface perturbations of the outer surface affect the swimming velocity.
Mathematically, this problem is defined by considering the Stokes squirmer with boundary conditions in the reference frame
\begin{equation} \left.u_\theta\right|_R = \left.u_\theta^0\right|_R - \epsilon U_1 V_1, \quad \left. u_r\right|_R = \left. u_r^0\right|_R + \epsilon U_1 P_1, \quad \left. \boldsymbol{u}\right|_{R'_o(\xi)} = 0, \end{equation}
where 
$V_1$ is the first associated Legendre polynomial, 
$\boldsymbol {u}^0$ is the solution to the pure treadmill problem with swimming speed 
$U_0$ provided in Reigh & Lauga (Reference Reigh and Lauga2017), and 
$R'_o(\xi ) = R_o[1-\epsilon \sum _{k} b_k\,P_k(\xi )]$ is the perturbed outer surface expanded based on the Legendre polynomial, with 
$b_k$ representing the amplitude of mode 
$k\geq 0$. This is an alternative formulation of (2.3a,b) with 
$V_n(\nu )=2/\sqrt {3n(n+1)}\,K_n(\nu )$. The new swimming speed is defined by 
$U = U_0 + \epsilon U_1$. We calculate how 
$U_1$ depends on the mode shapes of the outer surface perturbation.
The first-order solution 
$\boldsymbol {u}^1$ can be found in the same manner as in the original problem, with a modified boundary condition at 
$r=R_o$ obtained from the Taylor series approximation at the bounding surface 
$R'_o(\xi )$. The 
${O}(\epsilon )$ problem has updated boundary conditions
\begin{gather} \left. u_\theta^1\right|_R ={-} U_1V_1,\quad \left. u_r^1\right|_R = U_1P_1,\quad \left. u_\theta^1\right|_{R_o} = \left.\frac{\partial u_\theta^0}{\partial r}\right|_{R_o} P_k,\quad \left. u_r^1\right|_{R_o} = \left. \frac{\partial u_r^0}{\partial r}\right|_{R_o} P_k. \end{gather}
Because we know 
$\boldsymbol {u}^0$ analytically (Reigh & Lauga Reference Reigh and Lauga2017), we can compute these derivatives exactly. Furthermore, we can use the expansion 
$u_\theta ^1|_{R_o} = \sum _{n=1}^\infty \gamma _n V_n$, and solve for the resulting coefficients using orthogonality, which results in
\begin{equation} U_1 = \gamma_1\left[1- \frac{5}{3}\,\frac{\displaystyle\sum\limits_{i=0}^{2} \lambda^i}{\displaystyle\sum\limits_{i=0}^{4} \lambda^i}\right], \end{equation}
where 
$\lambda =R_o/R$. The term in brackets is positive for all 
$\lambda$, so the adjustment to swimming speed depends entirely on the coefficient 
$\gamma _1$. We can use the orthogonality of associated Legendre polynomials to obtain the expression for 
$\gamma _1$ as
\begin{equation} \gamma_1 = \frac{B_1 15 {\rm \pi}^2\lambda}{8\,R(\lambda^5-1)} \left[-\left(\frac{{\rm \pi}+4}{\rm \pi}\right)^2 b_0 -\frac{5}{64}\,b_2 - \frac{9}{4096}\,b_4 + \cdots\right]. \end{equation} The presence of the 
$b_0$ mode corresponds to a simple decrease in the confining sphere's radius. This portion of the asymptotic result is consistent with a result from Reigh & Lauga (Reference Reigh and Lauga2017) that showed that increased confinement has a purely negative effect on squirmer swimming speed. The increased swimming speed 
$U_1$ will have the same sign as 
$\gamma _1$. This result shows that if the sphere is perturbed with the 
$0\text {th}$ or even higher order Legendre mode, then the swimming speed decreases. Interestingly, this asymptotic result shows that the swimming speed can also result from the presence of the even Legendre modes in the perturbed outer surface. Depending on the sign of the modes 
$b_i$, these results show that perturbed surfaces can even lead to positive 
$U_1$.
A Legendre decomposition of perturbations to the yield surface of figure 3 shows that the surface perturbations consist mostly of the second Legendre mode regardless of 
$\textit {Bi}$. For example, at 
$\textit {Bi}{}=1$, approximately 
$87\,\%$ of the surface Legendre mode energy is in the (negative) second mode. Therefore, the surface perturbation contributes to the increase of the swimming performance, compared to the prediction based on the concentric spheres. This suggests that modifications to swimming speed come through two competing effects: the size and shape of the yield surface.
4.2. Multimode squirmer
As discussed above, in a Newtonian fluid, all modes after the first fail to produce propulsion. This is because Newtonian fluids are homogeneous and the boundary modes 
$\beta _n$, 
$n\geq 2$, are symmetric in a sense that ‘forward’ components of the surface stroke motion exactly cancel out the ‘backward’ components. It is reasonable to wonder whether the introduction of viscoplasticity could break this symmetry and, consequently, higher-order modes could contribute to locomotion. To see if this is the case in a Bingham fluid, we examine cases whose surface stroke motion comprises the first three modes. We visualize the results using a unit circle where each point corresponds to the 
$\beta _2$ and 
$\beta _3$ modes, with all higher modes set to zero; this sets the 
$\beta _1$ mode automatically via the mode normalization.
The results for swimming velocity and normalized power expenditure for 
$\textit {Bi}{}=1$ and 
$10$ are shown in figure 4. They are largely consistent with those from the treadmill case in that the squirmer uses more energy in viscoplastic fluids and swims more slowly. However, an exception to this trend is when there are pure ‘mixing modes’. In particular, when 
$\beta _3\neq 0$, there is a non-zero swimming speed, a situation impossible in a Newtonian fluid. The fact that 
$\beta _2=1$ still results in no locomotion suggests that not all higher-order modes can lead to locomotion – this is explored further in the next subsection. Figures 4(b,d) show that the power expenditure behaviour changes substantially depending on 
$\textit {Bi}{}$ value. In comparison, the power expenditure is mostly affected by the second mode 
$\beta _2$ for 
$\textit {Bi}{}=1$, while at a higher yield limit of 
$\textit {Bi}{}=10$, it mostly modifies with the 
$\beta _3$ mode. Furthermore, interestingly, the highest swimming velocity is observed at finite 
$\beta _3$ values rather than in pure treadmill swimming mode. In figure 6, the optimal 
$\beta _3$ value associated with the maximum swimming velocity is shown versus 
$\textit {Bi}{}$. The optimization is performed for the combined 
$\beta _1$, 
$\beta _2$ and 
$\beta _3$ multimode squirming motion, while the value of 
$\beta _1$ is specified through the constraint relation (2.4). For all the cases, the maximum swimming velocity always occurs when 
$\beta _2=0$, as expected. For low 
$\textit {Bi}{}<1$, the pure treadmill motion results in the highest swimming velocity, while for larger 
$\textit {Bi}{}$, the third mode of surface actuation enhances the locomotion of the swimmer. For very high yield limits 
$\textit {Bi}{}>10$, the maximum swimming velocity always occurs at 
$\beta _3=0.17$.
The effects of combinations of the first three modes of a Bingham fluid on (a,c) swimming speed 
$U$, and (b,d) normalized power expenditure 
$\mathcal {P}^*$, for (a,b) 
$\textit {Bi}{}=1$ and (c,d) 
$\textit {Bi}{}=10$. Flow fields at 
$(\beta _1, \beta _2, \beta _3) = (0, 0.79, 0.61)$ are examined in (e) 
$\textit {Bi}{}=1$, (f) 
$\textit {Bi}{}=10$ and (g) Newtonian fluid. Streamlines represent the velocity field, and colour represents pressure. The yield surface is shown with a red dashed line. The Bingham case can swim while the Newtonian case is immobile.
In figures 4(e)–4(g), we examine the flow fields of the purely mixing mode 
$(\beta _1, \beta _2, \beta _3) = (0, 0.79,0.61)$ in two Bingham environments with 
$\textit {Bi}{}=1$ and 10, as well as Newtonian fluid. The Bingham flow profiles are quite different from that for Newtonian fluids; the Bingham fluid displays a stagnation point away from the surface in the aft body, which does not exist in the equivalent Newtonian fluid. This stagnation point is caused by the confinement effect of the yield surface. The next subsection examines pure higher-order modes to determine why the third mode can lead to locomotion.
In Newtonian fluids, the addition of higher modes results in lower swimming efficiency as they do not contribute to the swimming speed and only add to the power expenditure by the swimmer (Blake Reference Blake1971; Stone & Samuel Reference Stone and Samuel1996). In non-Newtonian fluids, however, previous studies observed that higher modes could be effective in swimming depending on the rheology of the fluid medium (Datt et al. Reference Datt, Zhu, Elfring and Pak2015; Nganguia et al. Reference Nganguia, Pietrzyk and Pak2017). Here, we examine how the yield limit of viscoplastic fluid affects the swimmer's efficiency with the surface stroke made up of the first three modes.
Figure 5 shows the swimming efficiency of a squirmer for 
$\textit {Bi}{}=1$ and 10, and for Newtonian fluid. It is found that the swimmer has higher efficiency in yield fluids when 
$\textit {Bi}{}<30$, and the most efficient cases are associated with the combination of 
$\beta _1$ and 
$\beta _3$ modes. The region with high efficiency is shifted to higher 
$\beta _3$ modes when 
$\textit {Bi}{}$ increases. To explore the level of efficiency enhancement, in figure 6, the maximum attainable efficiency is shown as a function of the yield limit of Bingham fluid. Also, the corresponding 
$\beta _3$ magnitude that results in maximum efficiency is shown as a blue dashed line. It is also found (not shown here) that the inclusion of higher modes does not alter the results. The maximum efficiency is associated with pure treadmill motion for 
$\textit {Bi}{}<0.6$, while for a wide range 
$0.6<\textit {Bi}{}<30$, the optimal third mode increases logarithmically with 
$\textit {Bi}{}$ according to 
$\beta _3 \approx 0.6+0.21 \log _{10}\,(\textit {Bi}{}-0.6)$. The value of optimal 
$\beta _3$ reduces slightly for much larger 
$\textit {Bi}{}$ cases. The effect of the third mode on efficiency and its relation to the yield strength is reminiscent of previous observations in shear-thinning fluids where the efficiency change is interconnected to the relative magnitude of 
$\beta _2$ and 
$\beta _3$ modes (Nganguia et al. Reference Nganguia, Pietrzyk and Pak2017).
The swimming efficiency for the combinations of the first three modes in a Bingham fluid with (a) 
$\textit {Bi}{}=1$ and (b) 
$\textit {Bi}=10$, as well as (c) a Newtonian fluid.
The changes of maximum attainable swimming efficiency and optimal third mode actuation amplitude with 
$\textit {Bi}{}$. The solid line shows the maximum efficiency, the dashed line is the value of 
$\beta _3$ associated with maximum efficiency, and the dash-dotted line is the value of 
$\beta _3$ where the maximum swimming velocity is observed. The optimization is performed for a multimode squirmer with 
$\beta _1$, 
$\beta _2$ and 
$\beta _3$, while 
$\beta _1$ is calculated using the constraint relation (2.4). The optimal performance is always associated with 
$\beta _2=0$.
4.3. Higher pure modes
In the Newtonian fluid, modes with 
$\beta _1=0$ produce zero swimming speed. As we saw in the previous subsection, this is not the case in a Bingham fluid so long as 
$\beta _3\neq 0$. To determine the extent of this phenomenon, we conduct simulations for 
$\textit {Bi}{}=0.1,1,10$ for all pure higher-order modes 
$\beta _n$, 
$n=1,\dots,12$. It is noticed that odd actuation modes can lead to swimming, while even modes do not. Furthermore, this effect decreases with the mode number, e.g. 
$\beta _3$ has a significantly higher swimming speed than 
$\beta _5$, which has a larger swimming speed than 
$\beta _7$, and so on. For all even modes, the swimming speed is always zero.
We examine the flow fields for odd (figure 7a) and even (figure 7b) modes in both Bingham and Newtonian environments. The flow profiles for even modes are qualitatively similar between Bingham and Newtonian environments, while odd modes (
$n\geq 3$) display qualitatively different flow profiles. Specifically, the only difference for even modes is that the vortices generated near the surface due to the boundary conditions shrink more tightly towards the body as the yield surface gets closer to the swimmer; this manifests as a confinement effect. On the other hand, odd-numbered modes, in particular mode 
$\beta _3$, have substantially different flow fields wherein a vortex presented in the Newtonian fluid disappears and a non-periodic flow path takes its place in the viscoplastic environment. We hypothesize that the confinement effect of the yield surface breaks up this vortex and, importantly, allows the top slip boundary condition to ‘reach out’ to the far field, thus generating non-zero propulsion. This observation is aligned with prior results in shear-thinning fluids that the nonlinear rheology of the fluid medium can modify the near-body flow and pressure field, and results in new forms of non-Newtonian local and non-local forces on the body. More information about the contribution of the third stroke mode of the flow field around a squirmer in shear-thinning can be found in Pietrzyk et al. (Reference Pietrzyk, Nganguia, Datt, Zhu, Elfring and Pak2019).
Flow fields for (a) odd modes and (b) even modes. The odd modes display qualitative differences and maintain front–back antisymmetry while the even modes do not display qualitative differences and maintain front–back symmetry. Streamlines represent kinematic velocity field, and colour represents pressure field. The yield surface is shown with a red dashed line. All squirmers shown in (a) are mobile, and those shown in (b) are stationary.
Examining the effect of 
$\textit {Bi}{}$ on swimming speed at higher-order odd modes, we visualize the swimming speed (normalized by the swimming speed of a treadmill squirmer in Newtonian fluid) and power expenditure of pure 
$\beta _3$ and 
$\beta _5$ modes in figure 8. It can be seen that as 
$\textit {Bi}{}$ increases, both modes can achieve higher swimming speeds, although 
$\beta _5$ is growing significantly more slowly than 
$\beta _3$, the trend observed for higher-order odd modes. Higher modes attain their optimal swimming speed at particular 
$\textit {Bi}{}$ values that increase with the mode number. This is unlike the pure treadmill motion case discussed in § 4.1, where the velocity continuously decays with 
$\textit {Bi}{}$. For example, the maximum swimming speed of 
$\beta _3$ is observed at 
$\textit {Bi}{}=2.2$, and it shifts to 
$\textit {Bi}{}=7$ for 
$\beta _5$. Additionally, the power expenditure increases as 
$\textit {Bi}{}$ increases. This is related to the increasing confinement effect, which means that the same kinematic slip velocity boundary conditions ‘push harder’ against the surrounding fluid due to the close presence of plastic confinement.
Modes 
$\beta _3$ and 
$\beta _5$ generate non-zero swimming motion. An increasing power trend is also observed, similar to the treadmill case. Each odd-numbered mixing mode can generate its maximum swimming speed at a particular 
$\textit {Bi}{}$ range, different from the pure treadmill case shown in figure 2.
5. Discussion
We have explored an axisymmetric squirmer in a Bingham fluid using numerical simulation. It is shown that swimming in a Bingham fluid requires more power with the same kinematic boundary conditions for all modes. We argue that this effect is due to the increasing confinement effect from the yield surface. For the treadmill squirmer, we found that swimming speed decreased as the Bingham number 
$\textit {Bi}{}$ increased. Nevertheless, the shape of the yield surface can alleviate some of the decreasing effects. The efficiency of the swimmer, on the other hand, is higher at moderate 
$\textit {Bi}{}$ range before worsening at much higher Bingham numbers. We found that specific pure higher-order modes, specifically odd-numbered modes, could achieve locomotion on their own, a result impossible in a Newtonian fluid. Also, the 
$\beta _3$ and 
$\beta _5$ mode swimming speeds increase with increasing 
$\textit {Bi}{}$ for moderate Bingham numbers – an opposite trend compared to the treadmill mode. The swimming velocity and efficiency of the swimmer can be improved by including higher odd modes into the pure treadmill stroke (in particular, the third mode). The combination of these results suggests that there is a ‘crossover’ value of 
$\textit {Bi}{}$ at which the 
$\beta _3$ mode can be employed to enhance the efficiency and swimming speed compared to the treadmill mode. This suggests the existence of different optimal strategies in stroke selection in viscoplastic fluids, similar to other types of non-Newtonian fluids, especially shear-thinning fluids (Nganguia et al. Reference Nganguia, Pietrzyk and Pak2017; Pietrzyk et al. Reference Pietrzyk, Nganguia, Datt, Zhu, Elfring and Pak2019).
The results here and in previous studies on non-Newtonian modelled swimmers reveal new physics that we invite the community to investigate further. Many authors consider only the first two modes, 
$\beta _1$ and 
$\beta _2$, as it is believed that these approximate sufficiently the far-field flow dynamics of real microswimmers. Our results suggest that investigating higher-order modes is important for discovering possible alternative optimal swimming strategies in complex fluids, especially viscoplastic and elastoviscoplastic fluids. We believe that it is worth examining existing analytic results in viscoplastic fluid involving self-induced confinement with a focus on higher-order modes; e.g. see Reigh & Lauga (Reference Reigh and Lauga2017) and Nganguia et al. (Reference Nganguia, Zhu, Palaniappan and Pak2020). The nonlinear interaction between the modes could also be an interesting future extension of this study.
Finally, we address the direct applicability of these results to H. pylori. The Bingham fluid model ignores an important property of mucus, namely the fact that it is thixotropic. Therefore, our steady analysis presented here should be considered only an approximation, and we have possibly missed temporal effects present in real stomach mucus. Additionally, typical scales for H. pylori include length 
$R=1\ \mathrm {\mu }{\rm m}$, swimming speed 
$V=10\ {\mathrm {\mu }}{\rm m}\ {\rm s}^{-1}$, and characteristic viscosity of gastric mucus 
$\mu _0=0.05\ {\rm Pa}\ {\rm s}$ (Curt & Pringle Reference Curt and Pringle1969). Unfortunately, we could not find an estimate for the yield stress of gastric mucus, but using a range of estimates for common viscoplastic materials, from ketchup (
${15}\ {\rm Pa}$) to hair gel (
${135}\ {\rm Pa}$), provides an estimated range for 
$\textit {Bi}{}$ of 30–270. As our results suggest that higher-order odd modes become more effective at driving microswimmer propulsion at higher 
$\textit {Bi}{}$, this insinuates the possibility that H. pylori incorporates a modified mode of locomotion to swim more effectively through the stomach mucus. Related to this, results for two-dimensional squirmers showed that coupling treadmill and mixing modes achieved maximum propulsion at an intermediate 
$\textit {Bi}{}$ value (Supekar et al. Reference Supekar, Hewitt and Balmforth2020). Finally, real gastric mucus is chemically changed by the bacteria itself. Therefore, gastric mucus should not be treated as a homogeneous viscoplastic material regarding H. pylori locomotion – a more realistic model could include spatially varying yield stress 
$\tau _y$ subject to the spatial distribution of an emanating chemical.
Acknowledgements
The authors would also like to acknowledge the Florida State University Research Computing Center at which these simulations were carried out.
Declaration of interests
The authors report no conflict of interest.
Author contributions
P.S.E. and K.S. derived the theory and performed the simulations. All authors contributed equally to analysing data and reaching conclusions, and in writing the paper.
