1. Introduction
Recently there has been a great deal of work on studying locomotion in complex fluids, because, indeed, most biological fluids exhibit non-Newtonian behaviour owing to the presence of large biomolecules embedded within them (Spagnolie Reference Spagnolie2015; Patteson, Gopinath & Arratia Reference Patteson, Gopinath and Arratia2016). Notably, these biofluids are often viscoelastic, meaning they exhibit both a fluid-like and solid-like response to deformation (D'Avino & Maffettone Reference D'Avino and Maffettone2015). This fluid elasticity can have a profound impact on the motility of swimming micro-organisms, either hindering (Fu, Powers & Wolgemuth Reference Fu, Powers and Wolgemuth2007; Lauga Reference Lauga2007; Fu, Wolgemuth & Powers Reference Fu, Wolgemuth and Powers2009; Shen & Arratia Reference Shen and Arratia2011; Zhu et al. Reference Zhu, Do-Quang, Lauga and Brandt2011; Zhu, Lauga & Brandt Reference Zhu, Lauga and Brandt2012; Thomases & Guy Reference Thomases and Guy2014; Li et al. 2017; Binagia, Guido & Shaqfeh Reference Binagia, Guido and Shaqfeh2019) or enhancing (Teran, Fauci & Shelley Reference Teran, Fauci and Shelley2010; Liu, Powers & Breuer Reference Liu, Powers and Breuer2011; Spagnolie, Liu & Powers Reference Spagnolie, Liu and Powers2013; Riley & Lauga Reference Riley and Lauga2014; Thomases & Guy Reference Thomases and Guy2014, Reference Thomases and Guy2017; Patteson et al. Reference Patteson, Gopinath, Goulian and Arratia2015; Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020) propulsion depending upon the swimmer's gait, the structural properties of the immersed body and, of course, the rheology of the fluid (Dasgupta et al. Reference Dasgupta, Liu, Fu, Berhanu, Breuer, Powers and Kudrolli2013). Likewise, in regard to collective motion, viscoelasticity has been shown to have a substantial effect, engendering both aggregation and alignment of swimming cells (Li & Ardekani Reference Li and Ardekani2016; Tung et al. Reference Tung, Lin, Harvey, Fiore, Ardon, Wu and Suarez2017; Ishimoto & Gaffney Reference Ishimoto and Gaffney2018).
A workhorse model to study the effect of fluid elasticity on micro-organism motility is the squirmer model (Lighthill Reference Lighthill1952; Blake Reference Blake1971), which, in short, describes an active particle as being a nearly spherical body with a prescribed slip velocity on its surface (Pedley Reference Pedley2016). Using this mathematical description of a micro-swimmer, researchers have leveraged both numerical simulations (Zhu et al. Reference Zhu, Do-Quang, Lauga and Brandt2011, Reference Zhu, Lauga and Brandt2012; Li, Karimi & Ardekani Reference Li, Karimi and Ardekani2014Reference Li, Qin, Gopinath, Arratia, Thomases and Guy; De Corato & D'Avino Reference De Corato and D'Avino2017; Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020) as well as asymptotic analysis (Lauga Reference Lauga2009; Yazdi, Ardekani & Borhan Reference Yazdi, Ardekani and Borhan2014, Reference Yazdi, Ardekani and Borhan2015; De Corato, Greco & Maffettone Reference De Corato, Greco and Maffettone2015; Datt et al. Reference Datt, Natale, Hatzikiriakos and Elfring2017; Yazdi & Borhan Reference Yazdi and Borhan2017; Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020) to better understand the hydrodynamics of swimming in elastic fluids. In particular, though, no one to date has yet considered if the squirmer model itself is mathematically well-behaved in the context of the nonlinear constitutive models commonly used to model viscoelastic fluids such as the upper-convected Maxwell (UCM), Oldroyd-B, Giesekus and finite-extensibility nonlinear elastic model with the Peterlin approximation (FENE-P) models. Indeed, for the UCM and Oldroyd-B models, we demonstrate the presence of a singularity that is only removed for models that include regularization of the elastic stress, such as the Giesekus and FENE-P models. Furthermore, we show that the radius of convergence of asymptotic, perturbation, solutions in terms of the Weissenberg number, Wi, applied to squirming in elastic fluids is very small, necessitating the use of techniques that accelerate the convergence of series (Housiadas Reference Housiadas2017), where Wi, defined in § 2, gives a measure of the elasticity of the ambient fluid. The singularity at small Wi and the even smaller radius of convergence of the relevant perturbation solutions fully explain the observations made by other researchers that the higher-order corrections significantly add to the results, leading to qualitatively different speeds (Datt et al. Reference Datt, Natale, Hatzikiriakos and Elfring2017; Datt & Elfring Reference Datt and Elfring2019, Reference Datt and Elfring2020). They also provide the range of Wi at which the perturbation expansions can render accurate predictions.
While the spherical squirmer model has been described in detail in the literature and is summarized in our latest article (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020), here we provide a brief description for completeness. We assume a spherical, axisymmetric, neutrally buoyant, micro-swimmer with radius R in a highly viscous matrix fluid. The fluid is considered Newtonian with constant mass density ρ and constant shear viscosity 
${\eta _s}$. The micro-swimmer is translated with velocity V and rotates with rotation rate 
$\omega$. The motion of the body takes place under the influence of the force exerted on its surface by the ambient fluid, in the absence of any external forces and torques, and under isothermal and steady-state conditions. In this case, the relevant Reynolds number for the fluid, 
$Re = \rho VR/{\eta _s}$, is negligible, 
$Re \ll 1$, and it can be safely considered zero. Because of the axisymmetric shape of the body, the directions of translation and rotation are the same, which implies that the swimmer can only swim in a straight line (Pak & Lauga Reference Pak and Lauga2014).
We choose an inertialess (co-moving) frame of reference and a spherical coordinate system with origin located at the centre of the swimmer. The spherical coordinates are 
$\tilde{r},\;\theta ,\;\phi$ (a tilde above a variable indicates a dimensional quantity) and the unit vectors are 
$\{ {\boldsymbol{e}_r},{\boldsymbol{e}_\theta },{\boldsymbol{e}_\phi }\}$, where 
$\tilde{r}$ is the distance from the centre of the swimmer, 
$\theta$ is the polar angle measured from the axis of symmetry 
$(0 \le \theta \le {\rm \pi})$ and 
$\phi$ is the azimuthal angle 
$(0 \le \phi \lt 2{\rm \pi} )$; see figure 1 for more details. The velocity vector and the isotropic pressure are denoted by 
$\tilde{\boldsymbol{v}} = {\tilde{v}_r}{\boldsymbol{e}_r} + {\tilde{v}_\theta }{\boldsymbol{e}_\theta } + {\tilde{v}_\phi }{\boldsymbol{e}_\phi }$ and 
$\tilde{p}$, respectively. Under the assumptions mentioned above, the governing equations for the flow around the micro-swimmer are
\begin{equation}\boldsymbol{\nabla} \boldsymbol{\cdot}\tilde{\boldsymbol{v}} = 0,\quad - {\boldsymbol{\nabla}}\tilde{p} + {\eta _s}{\nabla ^2}\tilde{\boldsymbol{v}} = {\boldsymbol{0}}.\end{equation}Geometry, coordinate systems and dimensionless boundary conditions. (a) Cartesian: xyz with the x-axis being normal to the yz plane. (b) Spherical: 
$r\theta \phi $, where 
$\theta $ is the polar angle 
$(0 \le \theta \le {\rm \pi})$ and 
$\phi $ is the azimuthal angle 
$(0 \le \phi \le 2{\rm \pi} )$.
We are interested in velocity fields characterized by no radial velocity on the surface of the body, i.e. 
${\tilde{v}_r}(\tilde{r} = R,\theta ,\phi ) = 0$, which corresponds to the no-penetration boundary condition on a non-deformable object. Equations (1.1a,b) are linear and can be solved analytically to find the general solution for the field variables in the domain exterior to the swimmer 
$(\tilde{r} \gt R)$. Then, by evaluating the solution on the surface of the body, 
$\tilde{r} = R$, one can impose the so-called ‘slip velocity’, which denotes a purely tangential motion 
${\tilde{\boldsymbol{v}}^{(S)}} \equiv \tilde{\boldsymbol{v}}(\tilde{r} = R,\theta ,\phi ) = v_\theta ^{(S)}{\boldsymbol{e}_\theta } + v_\phi ^{(S)}{\boldsymbol{e}_\phi }$, where hereafter the superscript (S) indicates a value on the surface of the body. Considering only axisymmetric and creeping flow, 
${\tilde{\boldsymbol{v}}^{(S)}}$ is generally given by the expression (Pak & Lauga Reference Pak and Lauga2014)
\begin{equation}{\tilde{\boldsymbol{v}}^{(S)}}(\theta ) = \sin (\theta )\left\{ {\left( {2\sum\limits_{n = 1}^\infty {{B_n}\frac{{{{P^{\prime}}_n}(\mu )}}{{n(n + 1)}}} } \right){\boldsymbol{e}_\theta } + \left( {\sum\limits_{n = 1}^\infty {{C_n}\frac{{{{P^{\prime}}_n}(\mu )}}{{{R^{n + 1}}}}} } \right){\boldsymbol{e}_\phi }} \right\},\end{equation}
where 
${B_n}$, 
${C_n}$, 
$n = 1,2,3, \ldots$, are constants, 
$\mu = \cos (\theta )$ and 
${P_n} = {P_n}(\mu )$ is the Legendre orthogonal polynomial in [−1, +1] of degree n. We focus on flows driven only by the first two polar modes, the first two azimuthal modes and a rigid-body motion in the z-direction with velocity 
$\omega R\sin (\theta )$. Note, however, that the first azimuthal mode, 
${C_1}$, and the rigid-body rotation rate, 
$\omega$, can be merged into a single parameter. Thus, the velocity at the surface of a micro-swimmer is expressed as
\begin{equation}{\tilde{\boldsymbol{v}}^{(S)}} = \left( {{B_1}\sin (\theta ) + \frac{{{B_2}}}{2}\sin (2\theta )} \right){\boldsymbol{e}_\theta } + \left( {\left( {\frac{{{C_1}}}{{{R^2}}} + \omega R} \right)\sin (\theta ) + \frac{{3{C_2}}}{{2{R^3}}}\sin (2\theta )} \right){\boldsymbol{e}_\phi }.\end{equation}
This simplification is made since, in a Newtonian fluid, the swimming speed is determined solely by the value of 
${B_1}$, while 
${B_2}$ is the only coefficient that appears in the force that the fluid exerts on the particle. Furthermore, these are the only two modes necessary to distinguish between pushers (swimmers that generate thrust from behind their body) and pullers (swimmers that generate thrust from the front of their body). Regarding the first two swirling modes considered here, 
${C_1}$ is closely related to 
$\omega$, and for a Newtonian fluid one determines the other (see below); while 
${C_2}$ corresponds to a rotlet dipole and is the cause of the slowest-decaying flow in the azimuthal direction (Pak & Lauga Reference Pak and Lauga2014). Thus, the surface velocity as described by (1.3) contains the most important information for this type of motion of the micro-swimmer.
We choose the characteristic velocity as the coefficient of the first polar mode, i.e. 
${B_1}$, and we define the dimensionless quantities
\begin{equation}\xi \equiv \frac{{{B_2}}}{{{B_1}}},\quad \zeta \equiv \frac{{{C_2}}}{{{R^3}{B_1}}},\quad U \equiv \frac{V}{{{B_1}}},\quad \varOmega \equiv \frac{{{C_1}}}{{{R^2}{B_1}}} + \frac{{\omega R}}{{{B_1}}}.\end{equation}
We will be referring to 
$\xi$ as the slip parameter, to 
$\zeta$ as the swirl parameter, to U as the speed of the swimmer (or translation velocity) and to 
$\varOmega$ as the dimensionless net rotation rate of the swimmer. Thus, the dimensionless boundary conditions on the surface of the body become
\begin{equation}v_r^{(S)} = 0,\quad v_\theta ^{(S)} = \sin (\theta ) + \frac{\xi }{2}\sin (2\theta ),\quad v_\phi ^{(S)} = \varOmega \sin (\theta ) + \frac{{3\zeta }}{2}\sin (2\theta ).\end{equation}Evaluating the continuity equations (1.1a,b) at r = 1 and using (1.5a–c) gives
\begin{equation}{\left. {\frac{{\partial {v_r}}}{{\partial r}}} \right|_{r = 1}} ={-} \frac{\xi }{2} - 2\cos (\theta ) - \frac{{3\xi }}{2}\cos (2\theta ).\end{equation}
For ξ < 0 the micro-swimmer is defined as a pusher, for ξ > 0 it is a puller and for ξ = 0 the swimmer is considered neutral. Also, all the components of the velocity gradient tensor, 
${\boldsymbol{\nabla}}\boldsymbol{v}$, at r = 1 can be found using (1.5a–c) and (1.6) and the fact that the flow field is axisymmetric (i.e. 
$\partial \boldsymbol{v}/\partial \phi = 0$), with the exception of the components 
${({\boldsymbol{\nabla}}\boldsymbol{v})_{r\theta }} = \partial {v_\theta }/\partial r$ and 
${({\boldsymbol{\nabla}}\boldsymbol{v})_{r\phi }} = \partial {v_\phi }/\partial r$. To determine the latter two components, full knowledge of the velocity field around the swimmer is needed.
Finally, by using the superscript 
$(\infty )$ to denote the far flow field, the velocity 
$- U{\boldsymbol{e}_z}$ and pressure become
\begin{equation}v_r^{(\infty )} ={-} U\cos (\theta ),\quad v_\theta ^{(\infty )} = U\sin (\theta ),\quad v_\phi ^{(\infty )} = 0,\quad {p^{(\infty )}} = 0.\end{equation}Since the swimmer is force- and torque-free, two additional auxiliary conditions can be applied, the force-free and torque-free conditions (FFC and TFC, respectively). In the spherical coordinate system, these conditions are given as
\begin{equation}\int_{r = 1} {\boldsymbol{\mathsf{T}}\boldsymbol{\cdot }{\boldsymbol{e}_r}\,\textrm{d}S} = {\boldsymbol{0}},\quad \int_{r = 1} {{\boldsymbol{e}_r} \times (\boldsymbol{\mathsf{T}}\boldsymbol{\cdot }{\boldsymbol{e}_r})\,\textrm{d}S} = {\boldsymbol{0}},\end{equation}
where 
$\boldsymbol{\mathsf{T}} ={-} p\boldsymbol{\mathsf{I}} + \dot{\boldsymbol{\gamma }}$ is the dimensionless total stress tensor of the fluid, 
$\dot{\boldsymbol{\gamma }} = {\boldsymbol{\nabla}}\boldsymbol{v} + {({\boldsymbol{\nabla}}\boldsymbol{v})^\textrm{T}}$ is the rate-of-deformation tensor, the superscript ‘T’ denotes the transpose and 
$\boldsymbol{\mathsf{I}}$ is the unit tensor. By applying (1.8), the unknown swimming velocity U and net rotation rate 
$\Omega$ of the micro-swimmer are determined, leaving the remaining two dimensionless quantities ζ and ξ, reported in (1.4a–d), as free parameters. Consequently, we do not investigate either the effect of 
${C_1}$ or the effect of ω, but we discuss only the net rotation rate, 
$\Omega$, as this results from the TFC.
The analytical solution of (1.1a,b) with the associated boundary conditions (1.5a–c)–(1.7a–d) is
\begin{gather}{v_r} = \frac{\xi }{{4{r^2}}}\left( { - 1 + \frac{1}{{{r^2}}}} \right)(1 + 3\cos (2\theta )) + \left( { - U + \frac{{3U - 2}}{{2r}} + \frac{{2 - U}}{{2{r^3}}}} \right)\cos (\theta ),\end{gather}
\begin{gather}{v_\theta } = \left( {U + \frac{{2 - 3U}}{{4r}} + \frac{{2 - U}}{{4{r^3}}}} \right)\sin (\theta ) + \frac{\xi }{{2{r^4}}}\sin (2\theta ),\end{gather}
\begin{gather}{v_\phi } = \frac{\varOmega }{{{r^2}}}\sin (\theta ) + \frac{{3\zeta }}{{2{r^3}}}\sin (2\theta ),\end{gather}
\begin{gather}p ={-} \frac{\xi }{{2{r^3}}}(1 + 3\cos (2\theta )) + \frac{1}{{{r^2}}}\left( {\frac{{3U}}{2} - 1} \right)\cos (\theta ).\end{gather}Moreover, by applying (1.8), we determine the translation velocity and the rotation rate of the swimmer as
\begin{equation}U = 2/3,\quad \varOmega = 0.\end{equation}
In this case, the 
$O(1/r)$ contribution in the velocity components and the 
$O(1/{r^2})$ contribution in the pressure field are eliminated. Thus, as is well known, the velocity of the micro-swimmer in a simple Newtonian fluid under creeping conditions is two-thirds of the coefficient of the first polar mode (Lighthill Reference Lighthill1952; Blake Reference Blake1971), i.e. 
$V = 2{B_1}/3$ in dimensional units, and the rotation rate is 
$\omega ={-} {C_1}/{R^3}$ (obviously, 
$\omega = 0$ for 
${C_1} = 0$). The linearity of the governing equations, (1.1a,b), has the consequence that the translational and rotational velocities of the body are independent. It is also clear that the slip and swirl parameters, ξ and ζ, do not affect the velocity of the body, U. Finally, the analytical solution reveals that at the poles, 
$\theta = 0,{\rm \pi}$, and at any radial position 
$r \ge 1$, 
${\dot{\gamma }_{ii}} \ne 0$ and 
${\dot{\gamma }_{ij}} = 0,\;i \ne j$, where 
$i,j = r,\theta ,\phi$. Therefore, at the poles of the spherical coordinate system, the flow is purely extensional in character.
2. The squirmer model in viscoelastic fluids
In our previous publication (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020), we described the fact that, for real swimming organisms in elastic fluids, for which the squirmer model is a reasonable model, e.g. Escherichia coli, the Weissenberg numbers based on their swimming speed and characteristic dimension are quite modest (Wi < 0.25). However, the Weissenberg number based on their swirl may be much larger, and, moreover, such swirling flow can make the swimming velocity, for the same swimming gait, faster than that in a Newtonian fluid of the same viscosity. Since this result was only demonstrated over a relatively small parameter range, and primarily for the Giesekus fluid, we revisit this problem focusing on the low -Weissenberg-number regime where significant analytic progress can be made.
Thus, we consider a viscoelastic highly viscous ambient fluid that is composed of a pure Newtonian solvent with viscosity 
${\eta _s}$ and a viscoelastic solute with zero -shear-rate viscosity 
${\eta _p}$ and single relaxation time 
$\lambda$. We also assume that the spherical micro-swimmer changes direction of swimming with a frequency 
$\nu$. We make dimensionless the governing equations by scaling all lengths with R, the velocity components with 
${B_1}$, the time with 
$1/\nu$, the pressure with 
$({\eta _s} + {\eta _p}){B_1}/R$ and the viscoelastic extra-stress components with 
${\eta _p}{B_1}/R$. Thus, the dimensionless total stress tensor can be expressed as
\begin{equation}\boldsymbol{\mathsf{T}} ={-} p\boldsymbol{\mathsf{I}} + \beta \dot{\boldsymbol{\gamma }} + (1 - \beta )\boldsymbol{\tau }.\end{equation}
In the above equation, 
$\boldsymbol{\tau }$ is the symmetric elastic extra-stress tensor and 
$\beta$ is the solvent viscosity ratio, 
$\beta \equiv {\eta _s}/({\eta _s} + {\eta _p})$, which by definition varies in the range 
$0 < \beta < 1$. For 
$\beta = 1$ the fluid is purely Newtonian, and for 
$\beta = 0$ it is purely polymeric. Then, the dimensionless mass and momentum balances are
\begin{equation}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{v} = 0,\quad \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\mathsf{T}} ={-} {\boldsymbol{\nabla}}p + \beta {\nabla ^2}\boldsymbol{v} + (1 - \beta )\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\tau } = S\frac{{\partial \boldsymbol{v}}}{{\partial t}} + Re\,\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{v},\end{equation}
where 
$Re \equiv \rho R{B_1}/({\eta _s} + {\eta _p})$ and 
$S \equiv \rho {R^2}\nu /({\eta _s} + {\eta _p})$ are the Reynolds and frequency- related Reynolds number, respectively. However, due to the micrometre size of the swimmer, the high viscosity of the ambient fluid and the small magnitude of the slip velocity, the Reynolds number in (2.2) is vanishingly small, 
$Re{\kern 1pt} = 0$. Moreover, for typical biological systems, the frequency 
$\nu$ is low; hence 
$S \ll 1$ and the transient term in (2.2) is negligible compared to the remaining terms (Eastham & Shoele Reference Eastham and Shoele2020).
The polymer extra-stress tensor is determined via the following generalized dimensionless equation:
\begin{equation}\left. {\begin{array}{* {20}{c}} {f(\boldsymbol{\mathsf{c}})\boldsymbol{\mathsf{c}} + De\dfrac{{\partial \boldsymbol{\mathsf{c}}}}{{\partial t}} + Wi(\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{\mathsf{c}} - \boldsymbol{\mathsf{c}}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{v} - {{({\boldsymbol{\nabla}}\boldsymbol{v})}^\textrm{T}}\boldsymbol{\cdot }\boldsymbol{\mathsf{c}}) + {a_m}{{(\boldsymbol{\mathsf{c}} - \boldsymbol{\mathsf{I}})}^2} = \boldsymbol{\mathsf{I}},}\\ {\boldsymbol{\tau } = \dfrac{{f(\boldsymbol{\mathsf{c}})\boldsymbol{\mathsf{c}} - \boldsymbol{\mathsf{I}}}}{{Wi}},} \end{array}} \right\}\end{equation}
where 
$\boldsymbol{\mathsf{c}}$ is the symmetric and positive definite conformation tensor, 
$f = f(\boldsymbol{\mathsf{c}})$ is a strictly positive scalar function of 
$\boldsymbol{\mathsf{c}}$ defined below, 
${a_m}$ is a rheological parameter, Wi is the Weissenberg number, 
$Wi \equiv \lambda {\kern 1pt} {B_1}/R$, which is the well-known measure of the elasticity of the fluid, and 
$De \equiv \lambda \nu$ is the frequency-related Deborah number. Equation (2.3) shows that the transient term 
$\partial \boldsymbol{\mathsf{c}}/\partial t$ can be neglected provided that 
$De \ll 1$, namely, as long as the characteristic time 
$1/\nu$ is much longer than the relaxation time 
$\lambda$ of the fluid. In summary, we utilize the spherical squirmer model for viscoelastic ambient fluids, namely (2.1)–(2.3) with 
$Re = S = De = 0$. We also mention here that the suitability of the steady-state squirmer model for the study of real micro-swimmers in biological Newtonian systems has been investigated by Ito, Omori & Takuji (Reference Ito, Omori and Takuji2019). These authors concluded that the steady squirmer model can predict very well the swimming speed of the body but drastically underestimates the power consumption, i.e. the energy cost of swimming.
As far as the boundary and auxiliary conditions are concerned, these are the same as in the Newtonian case, namely the boundary conditions, (1.5a–c)–(1.7a–d), and the FFC and TFC, (1.9), are used (see also below in § 6). We summarize specific constitutive models that are special cases of the generalized equation (2.3); for more details about the models see Bird, Armstrong & Hassager (Reference Bird, Armstrong and Hassager1987a).
UCM/Oldroyd-B models
The UCM/Oldroyd-B models are the most fundamental nonlinear differential models in viscoelasticity. Both models are recovered from (2.3) with 
${a_m} = 0$ and 
$f(\boldsymbol{\mathsf{c}}) = 1$. When 
$\beta = 0$, one gets the UCM model, when 
$0 < \beta < 1$ the Oldroyd-B model is recovered, and for 
$\beta = 1$ the model reduces to the simple Newtonian fluid.
Giesekus model
The Giesekus model (Giesekus Reference Giesekus1982) is recovered from (2.3) for 
${a_m} > 0$ and 
$f(\boldsymbol{\mathsf{c}}) = 1$. Thermodynamically admissible values for the Giesekus mobility parameter are 
$0 < {a_m} < 1/2$.
FENE-P model
The FENE-P model (Peterlin Reference Peterlin1966; Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987b) is recovered from (2.3) with 
${a_m} = 0$ and with the scalar function 
$f$, known as the ‘Peterlin function’, given by the expression
\begin{equation}f(\boldsymbol{\mathsf{c}})= \frac{{1 - 3\chi }}{{1 - \chi \,\textrm{tr}(\boldsymbol{\mathsf{c}})}},\,\quad \chi \equiv \frac{1}{{{L^2}}},\end{equation}
where 
$L$ is the maximum extensibility parameter of the FENE-P model (L corresponds to the maximum ensemble average length of the polymer molecules).
In the analysis that follows, it is important to emphasize that the momentum balance does not play any role, and thus the results are purely kinematic in nature. First, evaluating all the components of the generalized constitutive model, (2.3), on the surface of the swimmer (at r = 1) in a spherical coordinate system, and using the applied boundary conditions (namely the no-penetration and slip velocity as given by (1.5a–c)–(1.7a–d)), we obtain the following:
rr-component
\begin{align}& f{{\mathsf{c}}_{rr}} + (1 -
2{{\mathsf{c}}_{rr}} + {\mathsf{c}}_{rr}^2 + {\mathsf{c}}_{r\theta }^2 + {\mathsf{c}}_{r\phi }^2){a_m}\nonumber\\
&\quad + Wi\left\{ {v_\theta^{(S)}\frac{{\partial
{{\mathsf{c}}_{rr}}}}{{\partial \theta }} + (\xi + 4\cos (\theta ) +
3\xi \cos (2\theta )){{\mathsf{c}}_{rr}}} \right\} =
1,\end{align}
$r\theta $-component
\begin{align}
& f{{\mathsf{c}}_{r\theta }} + ( - 2{{\mathsf{c}}_{r\theta }} + {{\mathsf{c}}_{rr}}{{\mathsf{c}}_{r\theta
}} + {{\mathsf{c}}_{r\theta }}{{\mathsf{c}}_{\theta \theta }} + {{\mathsf{c}}_{r\phi
}}{{\mathsf{c}}_{\theta \phi }}){a_m}\nonumber\\ & \quad + Wi\left\{
{v_\theta^{(S)}\dfrac{{\partial {{\mathsf{c}}_{r\theta }}}}{{\partial
\theta }} + \left( {\sin (\theta ) - \dfrac{{\partial
{v_\theta }}}{{\partial r}} + \dfrac{\xi }{2}\sin (2\theta
)} \right){{\mathsf{c}}_{rr}}} \right.\nonumber\\ & \quad \left. { + \left(
{\dfrac{\xi }{2} + \cos (\theta ) + \dfrac{\xi }{2}\cos
(2\theta )} \right){{\mathsf{c}}_{r\theta }}} \right\} = 0,
\end{align}
$r\phi $-component
\begin{align}
& f{{\mathsf{c}}_{r\phi }} + ( - 2{{\mathsf{c}}_{r\varphi }} + {{\mathsf{c}}_{rr}}{{\mathsf{c}}_{r\phi }}
+ {{\mathsf{c}}_{r\theta }}{{\mathsf{c}}_{\theta \phi }} + {{\mathsf{c}}_{r\phi }}{{\mathsf{c}}_{\phi
\phi }}){a_m}\nonumber\\ & \quad + Wi\left\{
{v_\theta^{(S)}\dfrac{{\partial {{\mathsf{c}}_{r\phi }}}}{{\partial
\theta }} - {{\mathsf{c}}_{rr}}\dfrac{{\partial {v_\phi }}}{{\partial
r}} + (\cos (\theta ) + \xi \cos (2\theta )){{\mathsf{c}}_{r\phi }}} + \varOmega\sin(\theta){\mathsf{c}}_{rr}
\right.\nonumber\\ & \quad \left. +\,{\dfrac{{3\zeta }}{2}(\sin
(2\theta ){{\mathsf{c}}_{rr}} + {{\mathsf{c}}_{r\theta }} - \cos (2\theta
){{\mathsf{c}}_{r\theta }})} \vphantom{\left\{ {v_\theta^{(S)}\dfrac{{\partial {{\mathsf{c}}_{\theta
\theta }}}}{{\partial \theta }} + 2\left( {\sin (\theta ) -
\dfrac{{\partial {v_\theta }}}{{\partial r}}}
\right){{\mathsf{c}}_{r\theta }}} \right.}\right\} = 0,
\end{align}θθ-component
\begin{align} &
f{{\mathsf{c}}_{\theta \theta }} + (1 + {\mathsf{c}}_{r\theta }^2 - 2{{\mathsf{c}}_{\theta
\theta }} + {\mathsf{c}}_{\theta \theta }^2 + {\mathsf{c}}_{\theta \phi }^2){a_m} \nonumber\\
&\quad + Wi\left\{ {v_\theta^{(S)}\dfrac{{\partial {{\mathsf{c}}_{\theta
\theta }}}}{{\partial \theta }} + 2\left( {\sin (\theta ) -
\dfrac{{\partial {v_\theta }}}{{\partial r}}}
\right){{\mathsf{c}}_{r\theta }}} \right.\nonumber\\ & \quad \left.{ -\, 2\cos (\theta
){{\mathsf{c}}_{\theta \theta }} + \xi (\sin (2\theta ){{\mathsf{c}}_{r\theta }}
- 2\cos (2\theta ){{\mathsf{c}}_{\theta \theta }})} \vphantom{\left\{ {v_\theta^{(S)}\dfrac{{\partial {{\mathsf{c}}_{\theta
\theta }}}}{{\partial \theta }} + 2\left( {\sin (\theta ) -
\dfrac{{\partial {v_\theta }}}{{\partial r}}}
\right){{\mathsf{c}}_{r\theta }}} \right.}\right\}= 1,
\end{align}
$\theta \phi$-component
\begin{align} &
f{{\mathsf{c}}_{\theta \varphi }} \!+\! ({{\mathsf{c}}_{r\theta }}{{\mathsf{c}}_{r\phi }} -
2{{\mathsf{c}}_{\theta \phi }} + {{\mathsf{c}}_{\theta \theta }}{{\mathsf{c}}_{\theta \phi
}} + {{\mathsf{c}}_{\theta \phi }}{{\mathsf{c}}_{\phi \phi }}){a_m}\nonumber\\
&\quad + Wi\,\left\{ {v_\theta^{(S)}\dfrac{{\partial {{\mathsf{c}}_{\theta \phi
}}}}{{\partial \theta }} + \left( {\sin (\theta ) -
\dfrac{{\partial {v_\theta }}}{{\partial r}}}
\right){{\mathsf{c}}_{r\phi }}}+ \varOmega\sin(\theta){\mathsf{c}}_{r\theta} \right.\nonumber\\ & \quad - 2\cos (\theta
){{\mathsf{c}}_{\theta \phi }} - {{\mathsf{c}}_{r\theta }}\dfrac{{\partial
{v_\phi }}}{{\partial r}} + \dfrac{\xi }{2}(\sin (2\theta
){{\mathsf{c}}_{r\phi }} - {{\mathsf{c}}_{\theta \phi }} - 3\cos (2\theta
){{\mathsf{c}}_{\theta \phi }})\nonumber\\ & \quad \left. { +\, \dfrac{{3\zeta
}}{2}(\sin (2\theta ){{\mathsf{c}}_{r\theta }} + {{\mathsf{c}}_{\theta \theta }}
- \cos (2\theta ){{\mathsf{c}}_{\theta \theta }})} \right\} = 0,
\end{align}
$\phi \phi $-component
\begin{align}
& f{{\mathsf{c}}_{\phi \phi }} + (1 + {\mathsf{c}}_{r\phi }^2 + {\mathsf{c}}_{\theta \phi }^2
- 2{{\mathsf{c}}_{\phi \phi }} + {\mathsf{c}}_{\phi \phi }^2){a_m}\nonumber\\
&\quad + Wi\left\{{v_\theta^{(S)}\dfrac{{\partial {{\mathsf{c}}_{\varphi \varphi}}}}{{\partial \theta }} - 2{{\mathsf{c}}_{r\phi }}\dfrac{{\partial
{v_\phi }}}{{\partial r}}} +2\varOmega\sin(\theta){\mathsf{c}}_{r\phi}\right.\nonumber\\ &\quad - \left. (2\cos (\theta ) + \xi (1 + \cos (2\theta ))){{\mathsf{c}}_{\phi \phi }} + 3\zeta (\sin (2\theta ){{\mathsf{c}}_{r\phi }} + {{\mathsf{c}}_{\theta \phi }} -
\cos (2\theta ){{\mathsf{c}}_{\theta \phi }}) \vphantom{\left\{ {v_\theta^{(S)}\dfrac{{\partial {{\mathsf{c}}_{\theta
\theta }}}}{{\partial \theta }} + 2\left( {\sin (\theta ) -
\dfrac{{\partial {v_\theta }}}{{\partial r}}}
\right){{\mathsf{c}}_{r\theta }}} \right.}\right\}= 1.
\end{align} Next, we proceed with each individual model separately, and we focus on the singular points of the coordinate system, i.e. at 
$\theta =0$ (the north pole) and 
$\theta= {\rm \pi}$ (the south pole). One is reminded here that the poles are not physical singular points of the coordinate system but are a consequence of the transformation from Cartesian to spherical coordinates. Indeed, the Jacobian determinant of the transformation is 
$|J| = {r^2} \sin (\theta )$, which becomes zero at 
$r = 0$ or at 
$\theta = 0,{\rm \pi}$. However, since we are interested only in the flow exterior to the spherical swimmer, the centre of the spherical coordinate system 
$(r = 0)$ does not belong to the domain of definition of the governing equations. On the contrary, the points 
$\theta = 0,{\rm \pi}$ belong to the domain both interior and exterior to the swimmer. Notice also that in the spherical coordinate system the domain 
$D = \{ (r,\theta )\mathrm{\mid }r \ge 0,\theta \in (0,{\rm \pi} )\}$ actually describes the z-axis of the Cartesian coordinate system, along which the body moves; see figure 1. Owing to the singular mathematical nature of the poles, the solution of the equations must satisfy certain conditions, and thus special attention to the poles needs to be given. For instance, and since we are interested in flow fields that do not depend on the azimuthal angle, the governing equations are well defined provided that the polar and azimuthal components of the velocity vector are zero at both poles. In addition, the off-diagonal components of the total stress tensor must be zero, while the associated diagonal components in the angles of the spherical coordinate system must be equal, viz.
\begin{equation}{{\mathsf{T}}_{\theta \theta }} = {{\mathsf{T}}_{\phi \phi }},\quad {{\mathsf{T}}_{\theta \phi }} = {{\mathsf{T}}_{r\theta }} = {{\mathsf{T}}_{r\phi }} = {v_\theta } = {v_\phi } = 0\;\textrm{at}\;r \ge 1,\quad \theta = 0\;\textrm{or}\;{\rm \pi} .\end{equation}Because of these properties, the flow along the z-axis is purely extensional. Finally, one can trivially confirm that the analytical solution for the Newtonian fluid given by (1.9)–(1.12) satisfies all the conditions given in (2.11a,b).
3. Analysis of squirming for the UCM/Oldroyd-B models
3.1. Analysis at the poles
For the UCM/Oldroyd-B models 
$({a_m} = 0,f(\boldsymbol{\mathsf{c}}) = 1)$, we evaluate (2.5)–(2.10) at 
$(r = 1,\textrm{ }\theta = 0)$ and 
$(r = 1, \theta = {\rm \pi})$. The resulting equations can be solved easily, showing that all the off-diagonal components of the conformation tensor are zero for both poles, 
${{\mathsf{c}}_{r\theta }} = {{\mathsf{c}}_{r\phi }} = {{\mathsf{c}}_{r\theta }} = 0$, while the diagonal components are
\begin{equation}{{\mathsf{c}}_{rr}} = \frac{1}{1 + 4\tilde{W}},\quad {{\mathsf{c}}_{\theta \theta}} = {{\mathsf{c}}_{\phi \phi }} = \frac{1}{1 - 2\tilde{W}},\end{equation}
where 
$\tilde{W}$ is a modified Weissenberg number, which is used hereafter for all the constitutive models and is defined as
\begin{equation}\tilde{W}: = \left\{ \begin{array}{@{}l@{}} (\xi + 1)Wi,\quad \theta = 0,\\ (\xi - 1)Wi,\quad \theta = {\rm \pi}. \end{array} \right.\end{equation}Equation (3.1a,b) reveals the following features for the UCM/Oldroyd-B models:
(a) As
$Wi \to 0$ the diagonal components of the conformation tensor go to unity, i.e. ${{\mathsf{c}}_{rr}} = {{\mathsf{c}}_{\theta \theta }} = {{\mathsf{c}}_{\phi \phi }} = 1$ (as they should). Unexpectedly, and for any $Wi > 0$, the same limit is also achieved for $\xi = -1$ at the north pole, and for $\xi = 1$ at the south pole.(b) The flow at the poles is purely extensional, as previously identified for a Newtonian fluid. However, the modified Weissenberg number,
$\tilde{W}$, which determines the degree of extensional character of the flow, is different between the poles. Note that, for the neutral squirmer, i.e. for ξ = 0, the modified Weissenberg number has the same magnitude but opposite signs at the poles ($\tilde{W} = Wi > 0$ on the north pole, and $\tilde{W} = - Wi < 0$ on the south). This implies a uniaxial extensional flow on the north pole, with z being the main axis of elongation, and a biaxial stretching in the x- and y-directions (i.e. on the z = −1 plane) at the south pole.(c) The solution for the conformation tensor is physically valid only in a specific window for the Weissenberg number. Based on the limit as Wi goes to zero, as well as the positive definite property of the conformation tensor, which requires that its diagonal components must be strictly positive, the window of validity of the UCM/Oldroyd-B models depends on the slip parameter ξ:
(3.3)\begin{equation}0 \le Wi \lt \left\{ \begin{array}{@{}l} \dfrac{1}{{4(1 - \xi )}},\quad \xi \le \dfrac{1}{3},\\ \dfrac{1}{{2(1 + \xi )}},\quad \xi \gt \dfrac{1}{3}. \end{array} \right.\end{equation}Hereafter, we will be using the symbol
$W{i_u}$ to denote the upper limit of validity for Wi, as given in (3.3). For $\xi = - \!1$, (3.3) gives $0 \le Wi < 1/8$, and for $\xi = 0$ and 1, it gives $0 \le Wi < 1/4$. The maximum window of validity is observed for $\xi = 1/3$; in this case one gets $0 \le Wi < 3/8$. We emphasize that these values correspond to the range of validity of the rheological models in steady uniaxial or biaxial elongational flow if the local rate of extension is correctly factored into the Weissenberg number at both poles.(d) The solution given by (3.1a,b) can be expanded in terms of
$\tilde{W}$ (or $Wi$) as follows:
(3.4)\begin{equation}\left. {\begin{array}{c@{}} {{{\mathsf{c}}_{rr}} \approx 1 - 4\tilde{W} + {{(4\tilde{W})}^2} - {{(4\tilde{W})}^3} + {{(4\tilde{W})}^4} + \cdots ,}\\ {{{\mathsf{c}}_{\theta \theta }} = {{\mathsf{c}}_{\phi \phi }} \approx 1 + 2\tilde{W} + {{(2\tilde{W})}^2} + {{(2\tilde{W})}^3} + {{(2\tilde{W})}^4} + \cdots .} \end{array}} \right\}\end{equation}The radius of convergence,
$W{i_\rho }$, of these series expansions is of great interest, for instance when perturbation methods are used to solve the full problem. It is determined as the minimum distance from the closest singularity in the complex plane. From (3.1a,b), one can see that for $\xi \ne \pm 1$ the singular points of the solution with respect to Wi are
(3.5a–d)\begin{align} W{i_{S,1}} &={-} \dfrac{1}{{4(1 + \xi )}},\quad W{i_{S,2}} ={-} \dfrac{1}{{4( - 1 + \xi )}},\nonumber\\ W{i_{S,3}} &= \dfrac{1}{{2(1 + \xi )}},\quad W{i_{S,4}} = \dfrac{1}{{2( - 1 + \xi )}}. \end{align}Their signs depend on the slip parameter, ξ. A negative
$W{i_{S,j}},\textrm{ }j = 1,2,3,4,$ reveals a non-physical singularity, while a positive one reveals a physical singularity at a finite Weissenberg number. Based on the singular points, the radius of convergence of the series solutions given in (3.4) is
(3.6)\begin{equation}W{i_\rho } = \left\{ \begin{array}{l@{}} \dfrac{1}{8},\quad \xi ={\pm} 1,\\ \dfrac{1}{4}min\left\{ {\dfrac{1}{{|1 + \xi |}},\dfrac{1}{{|1 - \xi |}}} \right\},\quad \xi \ne \pm 1. \end{array} \right.\end{equation}Since only positive values of Wi are of interest, the series solutions converge in the range
$0 \le Wi < W{i_\rho }$. The upper limit of validity of the exact solution, $W{i_u}$, along with the radius of convergence, $W{i_\rho }$, given by (3.3) and (3.6), respectively, are shown as a function of the slip parameter ξ, in figure 2. It is seen that, for $\xi \le 0$, the curves coincide, while, for $\xi > 0$, $W{i_\rho }$ is less than $W{i_u}$. As revealed by (3.6), $W{i_\rho }$ is symmetric with respect to ξ = 0, namely $W{i_\rho }(\xi ) = W{i_\rho }( - \xi )$, while no symmetry is predicted for $W{i_u}$. It is also seen that the magnitude of $W{i_\rho }$ is very small and decreases monotonically with the magnitude of the slip parameter; for instance, for $\xi = 0$, 1 and 3 one finds $W{i_\rho } = 1/4$, 1/8 and 1/16, respectively, while for $\xi = 1/3$, for which (3.3) predicts the largest window of validity of the exact solution $(0 \le Wi < 3/8)$, $W{i_\rho } = 3/16$. The magnitude of $W{i_u}$ also decreases monotonically as $\xi$ increases or decreases from $\xi = 1/3$.
The upper limit of validity of the solution for the UCM/Oldroyd-B models, Wiu (black, solid line), and the radius of convergence of the series solution, Wip (red, dashed line), as functions of the slip parameter, ξ.
Based on the analysis above, it is clear that, although the UCM/Oldroyd-B are fundamental models for studying the effect of elasticity in the flow, they can be utilized only for very small values of Wi, i.e. for weakly elastic fluids.
3.2. Analysis at the surface of the swimmer
We can further the analysis presented above by focusing on the rr-component of the constitutive model on the surface of the particle, i.e. in (2.5) with 
${a_m} = 0,f(\boldsymbol{\mathsf{c}}) = 1$,
\begin{equation}Wi(1 + \xi \cos (\theta ))\sin (\theta )G^{\prime}(\theta ) + (1 + Wi\xi + 4Wi\cos (\theta ) + 3Wi\xi \cos (2\theta ))G(\theta ) = 1,\end{equation}
where 
$G(\theta ) \equiv {{\mathsf{c}}_{rr}}(1,\theta )$ has been used for brevity. Equation (3.7) is a linear, inhomogeneous, first-order ordinary differential equation with non-constant coefficients. We emphasize that the swirling parameter, ζ, does not enter in the equation, and, most importantly, that (3.7) is decoupled from the other components of the constitutive model and therefore it can be studied independently. Thus, 
${{\mathsf{c}}_{rr}}(1,\theta )$ is fully determined through (3.7). However, a singular (or discontinuous) solution of (3.7) will affect the spatial evolution of the conformation components over the surface of the body, which in turn will affect all the field variables exterior to the body.
To examine this possibility, we will first change the independent variable and use 
$x = \cos (\theta )$ instead of the polar angle 
$\theta$. This change of variable transforms equation (3.7) into the following:
\begin{equation}-\! Wi(1 - {x^2})(1 + \xi x)G^{\prime}(x) + (1 + 2Wi( - \xi + 2x + 3\xi {x^2}))G(x) = 1,\quad - 1 \le x \le 1.\end{equation}Equation (3.8) shows that there are singular spatial points that depend on the magnitude of the slip parameter ξ:
(i) For
$|\xi |\lt 1$, there are two regular singular points, i.e. $x ={-} 1$ and $x = 1$.(ii) For
$|\xi |= 1$, there is one irregular and one regular singular point. In particular, for $\xi = 1$, $x ={-} 1$ is an irregular singular point and $x = 1$ is a regular singular point, while for $\xi ={-} 1$, $x = 1$ and $x ={-} 1$ are irregular and regular singular points, respectively.(iii) For
$|\xi |\gt 1$, there are three regular singular points, i.e. $x ={-} 1$, 1 and $- 1/\xi $.
At the singular points, G can be directly determined from (3.8) as
\begin{equation}\begin{array}{c} G( - 1) = \dfrac{1}{{1 + 4Wi( - 1 + \xi )}},\quad G(1) = \dfrac{1}{{1 + 4Wi(1 + \xi )}},\\ G( - {\xi ^{ - 1}}) = \dfrac{1}{{1 + 2Wi({\xi ^{ - 1}} - \xi )}}, \end{array}\end{equation}
where the last expression is valid only when 
$|\xi |\gt 1$, which implies that 
$- {\xi ^{ - 1}} \in ( - 1,1)$. The exact analytical solution of (3.8) has been found for any value of ξ and is reported below.
3.2.1. The case $|\xi |\lt 1$
In this case, the solution for 
$G(x) \equiv {{\mathsf{c}}_{rr}}(1,x)$ is
\begin{equation}\left. {\begin{array}{c@{}} {G(x) = \dfrac{1}{{Wi(1 + x)(1 + \xi x)}}\displaystyle\int_0^1 {{{(1 - t)}^a}{{\left( {1 + \dfrac{{1 - x}}{{1 + x}}t} \right)}^\delta }{{\left( {1 + \dfrac{{\xi (1 - x)}}{{1 + x\xi }}t} \right)}^\gamma }} \,\textrm{d}t,}\\ {a = 1 + \dfrac{1}{{2Wi(1 + \xi )}} \gt 0,\quad \delta = 1 - \dfrac{1}{{2Wi(1 - \xi )}} \lt - 1,\quad \gamma = 1 + \dfrac{\xi }{{Wi(1 - {\xi^2})}}.} \end{array}} \right\}\end{equation}
Equation (3.10) is given such that the solution is finite for any 
$x \in ( - 1,1]$, while for 
$x ={-} 1$, G is found as a limit. Also, one can confirm that 
$G(x = 1) = 1/(1 + 4Wi(1 + \xi ))$ and, as 
$x \to - {1^ + }$, 
$G \to 1/(1 + 4Wi( - 1 + \xi ))$; these values are in full agreement with (3.1a,b).
The solution for a neutral swimmer can be trivially recovered by setting ξ = 0 in (3.10). The resulting expression can be simplified as
\begin{equation}\left. {\begin{array}{c@{}} {G(x) = \dfrac{8}{{Wi{{(1 - x)}^{1 + a}}{{(1 + x)}^{1 + \delta }}}}\int_0^{(1 - x)/2} {{t^a}{{(1 - t)}^\delta }\,\textrm{d}t,} }\\ {a = 1 + \dfrac{1}{{2Wi}} \gt 0,\quad \delta = 1 - \dfrac{1}{{2Wi}} \lt - 1,} \end{array}} \right\}\end{equation}
where the second integral in (3.11) is simply the incomplete Beta function, 
$B((1 - x)/2,1 + a,1 + \delta )$. Again, from the above, 
$G(1) = 1/(1 + 4Wi)$, and 
$G \to 1/(1 - 4Wi)$ as 
$x \to - {1^ + }$ (valid only for 
$0 \le Wi \lt 1/4$) can be confirmed; otherwise the solution diverges. This is in full agreement with (3.1a,b) and (3.2), which for ξ = 0 give 
$0 \le Wi \lt 1/4$.
3.2.2. The case $|\xi |= 1$
For 
$\xi = 1$, the south pole (x = −1) becomes an irregular singular point of (3.8), while the north pole (x = 1) is a regular singular point. The solution is
\begin{equation}\left. {\begin{array}{c@{}} {G(x) = \dfrac{1}{{Wi}}\,{\textrm{e}^{ - 1/(2Wi(1 + x))}}\dfrac{1}{{(1 - {x^2})(1 + x)}}\int_x^1 {{\textrm{e}^{1/(2Wi(1 + x))}}{{\left( {\dfrac{{1 - t}}{{1 - x}}} \right)}^a}{{\left( {\dfrac{{1 + x}}{{1 + t}}} \right)}^\delta }\,\textrm{d}t,} }\\ {a = 1 + \dfrac{1}{{4Wi}}>0,\quad \delta ={-} 2 + \dfrac{1}{{4Wi}}>-1.} \end{array}} \right\}\end{equation}
For 
$\xi ={-} 1$, the north pole (x = 1) becomes an irregular singular point of (3.8), while the south pole (x = −1) is a regular singular point. The solution is as follows:
\begin{equation}\left. {\begin{array}{c@{}} {G(x) = \dfrac{1}{{Wi}}{\textrm{e}^{1/(2Wi(1 - x))}}\dfrac{1}{{(1 - {x^2})(1 - x)}}\int_x^1 {{\textrm{e}^{ - 1/(2Wi(1 - x))}}{{\left( {\dfrac{{1 - t}}{{1 - x}}} \right)}^a}{{\left( {\dfrac{{1 + x}}{{1 + t}}} \right)}^\delta }\,\textrm{d}t,} }\\ {a = 2 + \dfrac{1}{{4Wi}}>0,\quad \delta ={-} 1 + \dfrac{1}{{4Wi}}>1.} \end{array}} \right\}\end{equation}3.2.3. The case $|\xi |\gt 1$
In this case, (3.10) has three regular spatial singular points, the north pole (x = 1), the south pole 
$(x ={-} 1)$ and the point 
$x ={-} 1/\xi$. First, we define the auxiliary integral function 
$I = I(x;q,a,\delta ,\gamma )$
\begin{equation}I(x;q,a,\delta ,\gamma ): = \int_x^q {{{\left( {\frac{{1 + t}}{{1 + x}}} \right)}^a}{{\left( {\frac{{1 + \xi x}}{{1 + \xi t}}} \right)}^\delta }{{\left( {\frac{{1 - t}}{{1 - x}}} \right)}^\gamma }\,\textrm{d}t} ,\quad - 1 \lt x \lt 1,\end{equation}
where 
$a,\delta ,\gamma$ are constants. For 
$\xi > 1$, the analytical solution of (3.10) is
\begin{equation}G(x) = \frac{1}{{Wi(1 - {x^2})(1 + \xi x)}} \times \left\{ \begin{array}{l} I(x; - 1,a,\delta ,\gamma ),\quad - 1 \lt x \lt - 1/\xi ,\\ I(x;1,a,\delta ,\gamma ),\quad - 1/\xi \lt x \lt 1. \end{array} \right.\end{equation}
For 
$\xi < - \!1$, the analytical solution of (3.10) is
\begin{equation}G(x) = \frac{{I(x; - 1/\xi ,a,\delta ,\gamma )}}{{Wi(1 - {x^2})(1 + \xi x)}},\quad - 1 \lt x \lt 1,\;x \ne - 1/\xi .\end{equation}In (3.15) and (3.16) the following constants appear:
\begin{equation}a = 1 + \frac{1}{{2Wi(\xi - 1)}},\quad \delta ={-} 1 + \frac{\xi }{{Wi({\xi ^2} - 1)}},\quad \gamma = 1 + \frac{1}{{2Wi(\xi + 1)}}.\end{equation} Results for 
$G(x) \equiv {{\mathsf{c}}_{rr}}(1,x)$ as a function of 
$x$, for 
$\xi = 2,1,0, - 1$ and −2, are shown in figure 3 with solid lines. In all cases the Weissenberg number is 
$Wi = 1/9$, except for 
$\xi = - 2$, for which 
$Wi = 1/13$. Equation (3.8) is also solved numerically with a second-order finite difference method and the results are shown with dashed lines. Excellent agreement between numerical and analytical results is demonstrated in all cases. Significant differences can be witnessed in the surface stresses between pushers, neutral swimmers and pullers. Indeed, pullers (ξ = 2 and 1) induce large stretching of the polymer molecules on the surface of the micro-swimmer, with a maximum stretch which is observed at the southern hemisphere. Neutral swimmers show a monotonic decrease as we move away from the south and approach the north pole. Finally, pushers (ξ = −2 and −1) show an exponentially decreasing stretch close to the south pole, maintaining relatively small stretch values over most of the surface.
Plots of crr(1,x) as a function of x = cos(θ) for the UCM/Oldroyd-B models: (a) ξ = 2, Wi = 1/9, (b) ξ = 1, Wi = 1/9, (c) ξ = −1, Wi = 1/9, (d) ξ = −2, Wi = 1/13 and (e) ξ = 0, Wi = 1/9. Solid line, analytical solution; dots, numerical solution.
4. Analysis for squirming for regularized constitutive models
4.1. The Giesekus model
As before, we evaluate (2.5)–(2.10) at the poles. For 
$0 < {a_m} < 1/2$, we find that the off-diagonal components of the conformation tensor are zero at both poles, 
${{\mathsf{c}}_{r\theta }} = {{\mathsf{c}}_{r\phi }} = {{\mathsf{c}}_{r\theta }} = 0$, and the relevant (physical) solutions of interest for the normal components are
\begin{equation}\left. {\begin{array}{c@{}} {{{\mathsf{c}}_{rr}} = 1 - \dfrac{{4\tilde{W} + 1}}{{2{a_m}}} + \dfrac{{\sqrt {1 + 8\tilde{W}(1 - 2{a_m} + 2\tilde{W})} }}{{2{a_m}}},}\\ {{{\mathsf{c}}_{\theta \theta }} = {{\mathsf{c}}_{\phi \phi }} = 1 + \dfrac{{2\tilde{W} - 1}}{{2{a_m}}} + \dfrac{{\sqrt {1 + 4\tilde{W}( - 1 + 2{a_m} + \tilde{W})} }}{{2{a_m}}},} \end{array}} \right\}\end{equation}
where 
$\tilde{W}$ is given by (3.2). Equation (4.1) reveals interesting features of the Giesekus model.
(a) The limit of the above solution as
$Wi \to 0$ is ${{\mathsf{c}}_{rr}} = {{\mathsf{c}}_{\theta \theta }} = {{\mathsf{c}}_{\phi \phi }} = 1$. The same limit is achieved for $\xi = - \!1$ at the north pole, and for $\xi = 1$ at the south pole (i.e. for $\tilde{W} = 0$). Thus, in this limit, the behaviour of the solution is the same as that of the UCM/Oldroyd-B models.(b) The flow at the poles is purely extensional, as previously identified for Newtonian and UCM/Oldroyd-B fluids. The modified Weissenberg number,
$\tilde{W}$ (see (3.2)), determines the extensional character of the flow, which is different at the two poles.(c) In contrast to the UCM/Oldroyd-B models, the solution is valid for any value of the Weissenberg number,
$0 \le Wi < \infty$, and thus the singularity predicted for the UCM/Oldroyd-B models with respect to the Weissenberg number is fully removed.(d) The solution given by (4.1) can be expanded in terms of
$Wi$ (or $\tilde{W}$) as follows:
(4.2)\begin{equation}\left. \begin{array}{l} {\mathsf{c}}_{rr} \approx 1 - 4\tilde{W} +
{(4\tilde{W})^2}(1 - {a_m}) - {(4\tilde{W})^3}(1 - {a_m})(1 - 2{a_m})\\ \qquad +\, {(4\tilde{W})^4}(1 - {a_m})(1 - 5{a_m} + 5a_m^2) + \cdots, \\
{{\mathsf{c}}_{\theta \theta }} = {{\mathsf{c}}_{\phi \phi }} \approx 1 + 2\tilde{W} + {(2\tilde{W})^2}(1 - {a_m}) + {(2\tilde{W})^3}(1 - {a_m})(1 - 2{a_m})\\
\qquad +\, {(2\tilde{W})^4}(1 - {a_m})(1 - 5{a_m} + 5a_m^2) + \cdots . \end{array}\right\}\end{equation}As
${a_m} \to {0^ + }$, (4.2) reduces to (3.4), namely to the perturbation solution for the UCM/Oldroyd-B models. Since no singularities with respect to Wi exist for the Giesekus model, the region of convergence of these series expansions is determined by solving the inequalities $|8\tilde{W}(1 - 2{a_m} + 2\tilde{W})| < 1$ and $|4\tilde{W}( - 1 + 2{a_m} + \tilde{W})| < 1$, and finding the intersection of the resulting solutions. We find two branches, which we present in terms of $\tilde{W}$ as follows:
(4.3)and\begin{equation}\begin{array}{ccccc} & - \dfrac{1}{2} - {a_m} - \sqrt {a_m^2 - {a_m} + \dfrac{1}{2}} \lt \tilde{W} \lt - \dfrac{1}{4} + \dfrac{{{a_m}}}{2} + \dfrac{1}{2}\sqrt {a_m^2 - {a_m} + \dfrac{1}{2}} ,\\ & \quad \textrm{for}\;0 \le {a_m} \lt \dfrac{1}{2} - \dfrac{{\sqrt 2 }}{8}, \end{array}\end{equation}(4.4)\begin{equation}\begin{array}{ccccc} & - \dfrac{1}{4} + \dfrac{{{a_m}}}{2} - \dfrac{1}{2}\sqrt {a_m^2 - {a_m} + \dfrac{1}{2}} \lt \tilde{W} \lt - \dfrac{1}{4} + \dfrac{{{a_m}}}{2} + \dfrac{1}{2}\sqrt {a_m^2 - {a_m} + \dfrac{1}{2}} ,\\ & \quad \textrm{for}\;\dfrac{1}{2} - \dfrac{{\sqrt 2 }}{8} \lt {a_m} \lt \dfrac{1}{2}. \end{array}\end{equation}The ranges given by (4.3) and (4.4) are narrow and become even narrower when the admissible regions are given in terms of the original Weissenberg number; recall that
$\tilde{W}$ takes different values at the poles. In order to make this clear, we present the region of convergence of the perturbation series at the north and south poles, in figure 4(a). The regions are shown for both poles, in terms of the slip parameter, $\xi $, and the Weissenberg number, $Wi$, for a Giesekus mobility parameter ${a_m} = 0.2$. The intersection of these regions determines the final window of validity of the series solutions. The results reveal that, although the singularity predicted by the UCM/Oldroyd-B models has been removed, the corresponding final region of convergence is still very small, making (once again) the perturbation solution of quite limited value, while the perturbation results for $Wi \ge W{i_\rho }$ are divergent and completely erroneous. For comparison, the corresponding regions of convergence for the UCM/Oldroyd-B models are shown in figure 4(b). This is merely the region $0 \le |Wi| < W{i_\rho }$, where $W{i_\rho }$ is given by (3.6). Clearly, the difference in the region of convergence for the series between the constitutive models is very small. Thus, for all practical reasons, we can safely claim that the perturbation solutions for the UCM/Oldroyd-B and Giesekus models exhibit approximately the same radius of convergence, and they can be used only for weakly viscoelastic fluids.
The domain of convergence of the perturbation solution, (4.2), denoted by the common, shadowed area: (a) Giesekus model with αm = 0.2 and (b) UCM/Oldroyd-B models (αm = 0).
The diagonal components of the conformation tensor, 
${{\mathsf{c}}_{rr}}$ and 
${{\mathsf{c}}_{\phi \phi }} = {{\mathsf{c}}_{\theta \theta }}$, as well as the trace of the conformation tensor, 
$\textrm{tr}(\boldsymbol{\mathsf{c}}) = {{\mathsf{c}}_{rr}} + 2{{\mathsf{c}}_{\theta \theta }}$, are presented in figure 5 as functions of the slip coefficient, 
$\xi$. The results are shown with black solid lines for the north pole and with red dashed lines for the south pole; the Weissenberg number is Wi = 0.25, 0.5 and 1.0. For the individual components, a value larger than one denotes extension of the elastic molecules, while a value less than one denotes compression; to better distinguish between these two situations, the constant lines 
${{\mathsf{c}}_{rr}} = 1$ and 
${{\mathsf{c}}_{\theta \theta }} = 1$ have also been drawn. Interestingly, although 
$\xi = 0$ (a neutral squirmer) is the limit at which the micro-swimmer changes type, the values 
$\xi = \pm 1$ define the transition from compression to extension for the individual conformation components. For 
$\xi = 0$, 
${{\mathsf{c}}_{rr}}$ is less than one at the north pole, and larger than one at the south pole. Notice that the opposite is observed for 
${{\mathsf{c}}_{\theta \theta }}$. We can also see in figure 5(c) that the trace of the conformation tensor is always larger than three, at both poles, and achieves very large values throughout the whole range of the slip coefficient.
Plots of (a) c rr(r = 1), (b) 
${{\mathsf{c}}_{\theta \theta }}(r = 1)$ and (c) tr(c)(r = 1) as functions of the slip parameter ξ, for the Giesekus model with αm = 0.2, and Wi = 0.25, 0.5 and 1. The arrow shows the direction of increasing Wi. The results at the north pole are shown with black solid lines, and at the south pole with red dashed lines.
4.2. The FENE-P model
Similarly, we evaluate (2.5)–(2.10) at the poles for the FENE-P model. We find that the acceptable solution for the off-diagonal components of 
$\boldsymbol{\mathsf{c}}$ are again zero, at both poles, 
${{\mathsf{c}}_{r\theta }} = {{\mathsf{c}}_{r\phi }} = {{\mathsf{c}}_{r\theta }} = 0$. The diagonal components are found in terms of the Peterlin function, 
$f$,
\begin{equation}{{\mathsf{c}}_{rr}} = \frac{1}{{f + 4\tilde{W}}},\quad {{\mathsf{c}}_{\theta \theta }} = {{\mathsf{c}}_{\phi \phi }} = \frac{1}{{f - 2\tilde{W}}},\end{equation}
where 
$\tilde{W}$ is given by (3.2) and 
$f$ is given by
\begin{equation}f = \frac{{1 - 3\chi }}{{1 - ({{\mathsf{c}}_{rr}} + 2{{\mathsf{c}}_{\theta \theta }})\chi }}.\end{equation}Substituting equations (4.5a,b) in (4.6) and rearranging gives the following cubic equation, which determines the Peterlin function:
\begin{equation}{f^3} + {f^2}( - 1 + 2\tilde{W}) - f2\tilde{W}(1 + 4\tilde{W}) + 8{\tilde{W}^2}(1 - 3\chi ) = 0.\end{equation}
One can easily show that 
$f \ge 1 - 3\chi$. When the FENE-P parameter goes to zero 
$(\chi \equiv {L^{ - 2}} \to 0)$, i.e. for the UCM/Oldroyd-B models, the solutions of (4.7) are 
$f = 1,2\tilde{W}, - 4\tilde{W}$. In this case, only the 
$f = 1$ root is acceptable in the window 
$0 < \tilde{W} < 1/2$ (for 
$\tilde{W} > 1/2$ one gets the unphysical solution 
${{\mathsf{c}}_{\theta \theta }},{{\mathsf{c}}_{\phi \phi }} < 0$), while the other two roots are singular points for the diagonal components of the conformation tensor (see (4.5a,b)). For 
$0 < \chi < 1/3 \Rightarrow L > \sqrt 3$, there are three real roots, only one of which is acceptable:
\begin{equation}\begin{array}{ccccc} & f= \dfrac{{1 - 2\tilde{W}}}{3} + \dfrac{{2\sqrt {1 + 2\tilde{W}(1 + 14{{\tilde{W}}^2})} }}{3}\\ & \quad \times \cos \left( {\dfrac{{\rm \pi} }{6} + \dfrac{1}{3}{{\sin }^{ - 1}}\left( {\dfrac{{80{{\tilde{W}}^3} + 2{{\tilde{W}}^2}(39 - 162\chi ) - 3\tilde{W} - 1}}{{{{(1 + 2\tilde{W} + 28{{\tilde{W}}^2})}^{3/2}}}}} \right)} \right). \end{array}\end{equation}
The above analysis shows that there are no singular points with respect to the Weissenberg number. As previously revealed for the Giesekus model, the analytical solution has the correct limit as 
$\tilde{W} \to 0$, the flow at the poles is purely extensional, and the solution is valid for any value of 
$\tilde{W}$, or, alternatively, for any positive value of the Weissenberg number 
$(Wi > 0)$. Equations (4.5a,b) can be expressed as a series solution about 
$\tilde{W} = 0$, or 
$Wi = 0$, as follows:
\begin{equation}\left. {\begin{array}{c@{}} \begin{array}{ccccc} {{\mathsf{c}}_{rr}} \!\!\!\!\!\!& \approx 1 - 4\tilde{W} + {(4\tilde{W})^2}\left( {1 - \dfrac{{3\chi }}{2}} \right) - {(4\tilde{W})^3}\left( {1 - \dfrac{{15\chi }}{4}} \right)\\ & \quad +\, {(4\tilde{W})^4}\left( {1 - \dfrac{{57\chi }}{8} + \dfrac{{54{\chi^2}}}{8}} \right) + \cdots , \end{array}\\ \begin{array}{ccccc} {{\mathsf{c}}_{\theta \theta }} \!\!\!\!\!\! & = {{\mathsf{c}}_{\phi \phi }} \approx 1 + 2\tilde{W} + {(2\tilde{W})^2}(1 - 6\chi ) + {(2\tilde{W})^3}(1 - 6\chi )\\ & \quad +\, {(2\tilde{W})^4}(1 - 6\chi )(1 - 18\chi ) + \cdots . \end{array} \end{array}} \right\}\end{equation}
As 
$\chi \to {0^ + }$, (4.9) reduces to (3.4), namely to the perturbation solution for the UCM/Oldroyd-B models.
The domain of convergence of the series expressions, (4.9), is found by first substituting equation (4.8) in (4.7). Then, one needs to solve the inequalities 
$|\,f - 2/3 + 4\tilde{W}| < 1$ and 
$|\,f - 2/3 - 2\tilde{W}| < 1$ and determine the intersection of the solutions. Note that this procedure is required for both poles. Although the exact domain of convergence has not been found analytically, we can confirm the rather unexpected result that the domain of convergence for the FENE-P model is even smaller than that of the UCM/Oldroyd-B models. As 
$\chi \to 0$, (4.9) converges only for 
${-} 1/4 < \tilde{W} < 1/4$, while as 
$\chi$ increases the region of convergence decreases very slightly. Therefore, although the FENE-P removes the singularity of the UCM/Oldroyd-B models and provides a regular solution for the conformation tensor for any 
$\tilde{W}$, the corresponding series solutions exhibit the opposite trend, and the effect of 
$\chi$ (i.e. of the finite extensibility of the elastic molecules) is minor. As we mentioned for the Giesekus model, the perturbation solutions exhibit approximately the same radius of convergence, and they can be used only for weakly viscoelastic fluids; for 
$Wi \ge W{i_\rho }$ the perturbation solutions are divergent.
The diagonal components of the conformation tensor and the trace of the conformation tensor, 
${{\mathsf{c}}_{rr}}$, 
${{\mathsf{c}}_{\phi \phi }} = {{\mathsf{c}}_{\theta \theta }}$ and 
$\textrm{tr}(\boldsymbol{\mathsf{c}}) = {{\mathsf{c}}_{rr}} + 2{{\mathsf{c}}_{\theta \theta }}$, respectively, are presented in figure 6 as functions of the slip coefficient, 
$\xi$, for a maximum extensibility parameter 
$L = 10$ (χ = 0.01). The results are shown with black solid lines for the north pole and with red dashed lines for the south pole; the Weissenberg number is the same as in figure 5. For all quantities, we see quite similar results to those for the Giesekus model. It is also worth mentioning the fact that 
$\textrm{tr}(\boldsymbol{\mathsf{c}})$ achieves almost its maximum value (recall that for the FENE-P model 
$0 < \textrm{tr}(\boldsymbol{\mathsf{c}}) < {L^2}$) for a high enough absolute value of the slip parameter, even at a small Weissenberg number. Therefore, it is not Wi itself that really affects the degree of extension of the elastic molecules but its product with the slip parameter.
Plots of (a) crr(r = 1), (b) 
${{\mathsf{c}}_{\theta \theta }}(r = 1)$ and (c) tr(c)(r = 1) as functions of the slip parameter ξ, for the FENE-P model with L = 10 (χ = 0.01), and Wi = 0.25, 0.5 and 1. The arrow shows the direction of increasing Wi. The results at the north pole are shown with black solid lines, and at the south pole with red dashed lines.
5. Approximate solutions
All the results presented in the previous sections are analytical and exact. They contribute to the understanding of the flow around the spherical micro-swimmer and they can be used to compare with approximate solutions of the full governing equations obtained with large-scale numerical simulations or asymptotic methods. Both these methods have been developed and used in our previous publication (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020). For completeness we provide a short description of these methods along with additional information regarding the nonlinear processing of the high-order asymptotic solutions that we have derived.
5.1. Numerical solution
To solve the governing equations (2.2) and (2.3), and auxiliary conditions (1.5a–c)–(1.8), mentioned in § 2, numerically, we utilize the same method as described in detail in Binagia et al. (Reference Binagia, Phoa, Housiadas and Shaqfeh2020). In short, the overall methodology is as follows. We consider the co-moving frame of reference such that a body-fitted mesh may be used. In particular, we consider an unstructured mesh of tetrahedral elements with increasing resolution near the squirmer. The computational domain defined by this cylindrical mesh has length 40R and radius 20R. The governing equations and boundary conditions are solved using a third-order -accurate finite volume flow solver developed in a series of previous publications (Ham, Mattsson & Iaccarino Reference Ham, Mattsson and Iaccarino2006; Richter, Iaccarino & Shaqfeh Reference Richter, Iaccarino and Shaqfeh2010; Padhy et al. Reference Padhy, Shaqfeh, Iaccarino, Morris and Tonmukayakul2013; Yang, Krishnan & Shaqfeh Reference Yang, Krishnan and Shaqfeh2016; Castillo et al. Reference Castillo, Murch, Einarsson, Mena, Shaqfeh and Zenit2019; Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020; Zhang et al. Reference Zhang, Murch, Einarsson and Shaqfeh2020). The components of the conformation tensor are simultaneously determined by solving the polymer constitutive equation using the log-conformation method (Fattal & Kupferman Reference Fattal and Kupferman2004). The swimmer's translational and rotational velocity are found iteratively based on the measured net force and torque of the swimmer, respectively, until the micro-swimmer is found to be force- and torque-free. Extensive mesh refinement has been performed in order to achieve convergence and results of high accuracy. The exact analytical solutions found and reported in the previous sections are also used to confirm the accuracy of our numerical results; more comments and results follow further below.
5.2. Asymptotic solution
The asymptotic solution is derived for all the constitutive models employed here by following a regular perturbation scheme valid for weakly viscoelastic fluids. We mention that the first step of the solution procedure is to introduce an auxiliary symmetric tensor, 
$\boldsymbol{\sigma }$, according to the expression 
$\boldsymbol{\tau } = \dot{\boldsymbol{\gamma }} - Wi\boldsymbol{\sigma }$. Substituting in (2.1) gives the total stress tensor as 
$\boldsymbol{\mathsf{T}} = - p\boldsymbol{\mathsf{I}} + \dot{\boldsymbol{\gamma }} - (1 - \beta )Wi\boldsymbol{\sigma }$. Likewise, the new form of the balance equation (2.2) with 
$Re = S = 0$ is
\begin{equation}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{v} = 0,\quad - {\boldsymbol{\nabla}}p + {\nabla ^2}\boldsymbol{v} - (1 - \beta )Wi \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\sigma } = {\boldsymbol{0}}.\end{equation}
The Giesekus equation with 
$De = 0$ reduces to
\begin{equation}\boldsymbol{\sigma } + Wi\frac{{\delta \boldsymbol{\sigma }}}{{\delta t}} = \frac{{\delta \dot{\boldsymbol{\gamma }}}}{{\delta t}} + {a_m}(\dot{\boldsymbol{\gamma }}\boldsymbol{\cdot }\dot{\boldsymbol{\gamma }} - Wi(\boldsymbol{\sigma }\boldsymbol{\cdot }\dot{\boldsymbol{\gamma }} + \dot{\boldsymbol{\gamma }}\boldsymbol{\cdot }\boldsymbol{\sigma }) + W{i^2}\boldsymbol{\sigma }\boldsymbol{\cdot }\boldsymbol{\sigma }).\end{equation}
Similarly, the FENE-P model with 
$De = 0$ reduces to
\begin{equation}(\dot{\boldsymbol{\gamma }} - Wi\boldsymbol{\sigma })f(\boldsymbol{\sigma }) + Wi\frac{{\delta \dot{\boldsymbol{\gamma }}}}{{\delta t}} - W{i^2}\frac{{\delta \boldsymbol{\sigma }}}{{\delta t}} = \boldsymbol{v}\cdot \boldsymbol{\nabla}(\ln(f(\sigma)))({\boldsymbol{\mathsf{I}}} + Wi\dot{\boldsymbol{\gamma }} - W{i^2}\boldsymbol{\sigma }) + \dot{\boldsymbol{\gamma }},\end{equation}
where 
$f(\boldsymbol{\sigma }) = 1 - \chi W{i^2}tr (\boldsymbol{\sigma })$. In (5.2) and (5.3) 
$\delta ({\bullet} )/\delta t: = (\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla})({\bullet} ) - ({\bullet} )\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{v} - {({\boldsymbol{\nabla}}\boldsymbol{v})^\textrm{T}}\boldsymbol{\cdot }({\bullet} )$ is the steady-state upper-convected (or Maxwell) derivative of a second-order real and symmetric tensor.
For all the models employed here, the conformation tensor is given in general as
\begin{equation}\boldsymbol{c} = \frac{{\boldsymbol{\mathsf{I}} + Wi\dot{\boldsymbol{\gamma }} - W{i^2}\boldsymbol{\sigma }}}{{1 - \chi W{i^2}tr (\boldsymbol{\sigma })}},\end{equation}
where for the UCM/Oldroyd-B and Giesekus models, χ = 0. Although the constitutive equation in terms of 
$\boldsymbol{\sigma }$ is more complicated than its original form (see (2.3)), the solution procedure using regular perturbation methods is facilitated by this construction; for more details the interested reader is referred to Housiadas (Reference Housiadas2019). We also note that the new auxiliary tensor represents the deviation of the polymer extra stress 
$\boldsymbol{\tau }$ from the pure viscous tensor 
$\dot{\boldsymbol{\gamma }}$ suitably scaled with the Weissenberg number. Obviously, in the limit of vanishing Weissenberg number 
$\boldsymbol{\sigma } = {\lim _{Wi \to 0}}(\dot{\boldsymbol{\gamma}}-\boldsymbol{\tau})/Wi$ and 
$\boldsymbol{\sigma }$ reduces to the polymer extra-stress tensor of the second-order fluid model (Bird et al. Reference Bird, Armstrong and Hassager1987a). This implies that 
$\boldsymbol{\sigma } = O(1)$ as 
$Wi \to 0$, as can be seen from (5.2) and (5.3). It is worth mentioning however, that all formulations of the governing equations and the accompanying auxiliary conditions in terms of 
$\boldsymbol{c}$, 
$\boldsymbol{\tau }$ or 
$\boldsymbol{\sigma }$ are equivalent.
The details of the solution procedure are described in Housiadas (Reference Housiadas2019) and Binagia et al. (Reference Binagia, Phoa, Housiadas and Shaqfeh2020). In particular, we apply a regular perturbation scheme with the small parameter being the Weissenberg number, while the remaining dimensionless parameters 
$\zeta ,\xi ,\eta$ and 
${a_m}$ are considered 
$O(1)$ quantities. According to the perturbation method, each dependent flow variable is expanded asymptotically in a standard power series in terms of 
$Wi$:
\begin{equation}X \approx {X_0} + Wi{X_1} + W{i^2}{X_2} + \cdots \quad \textrm{as}\;Wi \to {0^ + }\textrm{,}\quad X = U,\varOmega ,p,\boldsymbol{v},\boldsymbol{\sigma }.\end{equation}
Equation (5.5) is substituted in (5.1)–(5.3) and in boundary conditions (1.5a–c)–(1.8) and thus a sequence of partial differential equations and boundary/auxiliary conditions at 
$O(W{i^j}),\;j = 0,1,2, \ldots$, are determined. The leading-order equations correspond to the Stokes equations, the solution of which is simply given by (1.9)–(1.13a,b). At higher orders in Wi, the equations are solved analytically with the aid of the ‘Mathematica’ software (Wolfram Research Inc. 2019) up to fourth order for all variables except for the azimuthal component of the velocity vector, 
${v_\phi }$, and the net rotation rate of the swimmer, 
$\varOmega$, which are found up to fifth order. The solutions are too involved to be presented in print, but they are available upon request. We do mention, however, that, for a neutral swimmer, i.e. for 
$\xi = 0$, the odd coefficients in the solution for U and the even coefficients in the solution for 
$\Omega$ are zero,
\begin{equation}U(\xi = 0) \approx {\textstyle{2 \over 3}} + (1 - \beta ){\tilde{U}_2}W{i^2} + (1 - \beta ){\tilde{U}_4}W{i^4},\end{equation}
\begin{equation}\varOmega (\xi = 0) \approx (1 - \beta )\zeta Wi({\tilde{\varOmega }_1} + W{i^2}{\tilde{\varOmega }_3} + W{i^4}{\tilde{\varOmega }_5}),\end{equation}
where in the expression for 
$\Omega$ the quantity 
$(1 - \beta )\zeta Wi$ is a common factor for all 
${\varOmega _j}$, and to avoid confusion with the original 
${\varOmega _j}$ symbols, we have used the scaled components 
${\tilde{\varOmega }_j} = {\varOmega _j}/((1 - \beta )\zeta Wi),\;j = 1,2,3, \ldots$. For the UCM/Oldroyd-B models and 
$\xi = 0$, the quantities that appear in (5.6) and (5.7) are provided in the appendix.
The correctness of our high-order perturbation solution has been confirmed by performing a variety of tests on the properties that the solution must satisfy but are not used or enforced during the solution procedure (Housiadas Reference Housiadas2019). Furthermore, we investigate and confirm that our analytical solution satisfies the properties given by (2.11a,b) at any order of the Weissenberg number. In particular, the polar and azimuthal components of the velocity vector are found as sine-Fourier series, which implies that, at both poles and regardless of the radial coordinate 
$r$, these velocity components are zero. Hence, the velocity gradients 
$\partial {v_\theta }/\partial r$ and 
$\partial {v_\phi }/\partial r$ are sine-Fourier series too, which become zero at the poles. On the contrary, the radial component of the velocity vector, 
${v_r}$, and its radial derivative 
$\partial {v_r}/\partial r$, are found as cosine-Fourier series, resulting in non-zero values for any 
$r > 1$. Therefore, only the diagonal components of the velocity gradient tensor are different from zero, revealing the purely extensional character of the flow along the z-axis, as previously found for a simple Newtonian fluid. Moreover, the off-diagonal components of the auxiliary tensor 
$\boldsymbol{\sigma }$, the polymer extra-stress tensor 
$\boldsymbol{\tau }$ and the conformation tensor 
$\boldsymbol{c}$ are zero at both poles for any 
$r \ge 1$. Lastly, it is confirmed that, at the poles, 
${\dot{\gamma }_{\theta \theta }} = {\dot{\gamma }_{\phi \phi }}$ and 
${\sigma _{\theta \theta }} = {\sigma _{\phi \phi }}$.
Finally, if we evaluate at the poles the asymptotic solution for the components of the conformation tensor and compare with (4.2) for the Giesekus model and with (4.9) for the FENE-P model, we confirm that the same expressions are derived. Recall that (4.2) and (4.9) have been derived quite independently from the exact analytical solution at the poles, i.e. from (4.1).
5.3. ‘Accelerated’ asymptotic solution
It is an almost impossible task to determine the domain of convergence of an asymptotic power series representation of the field variables of a nonlinear fluid mechanics problem. However, and based on the analysis at the poles presented in §§ 3 and 4, we conclude that the asymptotic solution, (5.5), is divergent for 
$Wi \ge W{i_\rho }$. In this region, the inclusion of the higher-order terms in the asymptotic solution leads to a fast divergence from the exact solution; indeed, this was observed in our previous paper (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020). On the contrary, we anticipate a good or reasonable agreement with the simulation results in the region 
$0 < Wi < W{i_\rho }$. Furthermore, high-order perturbation solutions, usually with three or more terms such as those given in (5.6) and (5.7), can be processed further in order to increase their accuracy and extend their domain of convergence. This has been done successfully in a variety of viscoelastic fundamental flows (Housiadas Reference Housiadas2017), Poiseuille-type flows with singularities (Housiadas Reference Housiadas2020) and viscoelastic flows around spherical bodies (Housiadas Reference Housiadas2019; Zhang et al. Reference Zhang, Murch, Einarsson and Shaqfeh2020). For the problem under consideration, two nonlinear techniques have been implemented; the Shanks transform (Shanks Reference Shanks1955) and the diagonal Padé [2/2] approximant (Padé Reference Padé1892); in order to distinguish between an original and a transformed solution, we will be referring to the latter as the ‘accelerated asymptotic solution’, and we will be using the subscript ‘acc’. These techniques are applied to constant and/or integral quantities of interest only, such as the speed of the micro-swimmer, U, its net rotation rate, 
$\Omega$, and the force contributions on the swimmer, which are presented and discussed in the subsequent section. We also emphasize that the accuracy and efficiency of these formulae were tested through extensive comparison with the results obtained from our large-scale numerical simulations and the exact analytical solutions derived in §§ 3 and 4.
Both the Shanks transform and the diagonal Padé [2/2] approximant yield the same expressions when applied on 
$U(\xi = 0)$ and 
$\varOmega (\xi = 0)/((1 - \beta )\zeta Wi)$, i.e.
\begin{equation}{U_{acc}}(\xi \!=\! 0) \!=\! \frac{2}{3} + \frac{{(1 - \beta )W{i^2}\tilde{U}_2^2}}{{{{\tilde{U}}_2} - W{i^2}{{\tilde{U}}_4}}},\quad {\varOmega _{acc}}(\xi \!=\! 0) = (1 - \beta )\zeta Wi\left( {{{\tilde{\varOmega }}_1} + \frac{{W{i^2}\tilde{\varOmega }_3^2}}{{{{\tilde{\varOmega }}_3} \!-\! W{i^2}{{\tilde{\varOmega }}_5}}}} \right).\end{equation}
When 
${\tilde{U}_2},{\tilde{U}_4}$ (and/or 
${\tilde{\varOmega }_3},{\tilde{\varOmega }_5}$) have the same sign, the expressions in (5.8a,b) indicate a singular point of the solution 
$W{i_s} = \sqrt {{{\tilde{U}}_3}/{{\tilde{U}}_5}}$ (and/or 
$W{i_s} = \sqrt {{{\tilde{\varOmega }}_3}/{{\tilde{\varOmega }}_5}}$). However, one cannot easily determine whether this singular point is not spurious. Nevertheless, when a singular point exists, the solutions are valid in the range 
$0 \le Wi < W{i_s}$. Note that, for a neutral swimmer and the range of parameters that are of interest, no singularities for the Giesekus and FENE-P models have been detected.
For a non-neutral swimmer, 
$\xi \ne 0$, the speed of the micro-swimmer and its rotation rate are constructed using the first five available terms in the perturbation solution, i.e. up to 
$O(W{i^4})$ and 
$O(W{i^5})$, respectively. In this case, the Shanks transform implemented with two iterations and the Padé [2/2] approximant result in different transformed solutions. Although, in general, the Shanks transform gives better results compared to Padé-type approximants (Hinch Reference Hinch1991; Housiadas Reference Housiadas2017, Reference Housiadas2020), the opposite trend has been observed for the current problem because the former technique predicts spurious singular points. Thus, hereafter we will be using only the diagonal Padé [2/2] for all types of swimmers.
We close this section, first, by emphasizing that the original perturbation solutions, derived for all the employed constitutive models, by no means should be used for 
$Wi \ge W{i_\rho }$, where 
$W{i_\rho }$ is given by (3.6). Moreover, detailed comparison with the simulation results revealed that the series solutions, (5.6), appear to predict the right trends only in the range 
$0 \le Wi < 0.1$ approximately. Even for the Giesekus and FENE-P models, (5.6) give completely erroneous results at 
$Wi \approx W{i_\rho }$. On the contrary, the transformed accelerated solutions restore the right trends for all the constant and integral quantities mentioned above in the range 
$0 \le Wi < W{i_\rho }$.
6. Results and mechanisms for speed and rotation rate enhancement
In our recent paper (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020) we predicted that the swirl parameter induces a speed enhancement in an elastic fluid over and above the speed in a Newtonian fluid of the same viscosity. We also analysed the force contributions acting on the body for the Giesekus model and discussed how these changed as the model parameters are varied. This led to a brief discussion of the mechanism of speed enhancement for a squirmer swimming in a Giesekus fluid. From our discussion above, especially in regard to (2.5)–(2.10) for the stress on the surface of the squirmer, as well as the pure extensional character of the flow at the poles, we can anticipate that these surface force contributions will be sensitive to the rheological model. Thus, we present new aspects of the force contributions and the mechanism behind the speed enhancement for the different rheological models examined in this paper.
We reiterate that the FFC and TFC, as given in (1.8), are required to determine the velocity and rotation rate of the swimmer. Owing to its spherical shape, the symmetry with respect to the z-axis and the form of the total stress tensor given by (2.1), there is only one non-trivial component of the FFC which can be expressed in terms of the components of the polymer extra stress 
$\boldsymbol{\tau }$, the conformation tensor 
$\boldsymbol{c}$ or the auxiliary tensor 
$\boldsymbol{\sigma }$. All formulations are equivalent, and in terms of 
$\boldsymbol{\sigma }$ we have
\begin{equation}\int_0^{\rm \pi} {\{ ( - p + {{\dot{\gamma }}_{rr}} - (1 - \beta )Wi\,{\sigma _{rr}})\cos (\theta ) - ({{\dot{\gamma }}_{r\theta }} - (1 - \beta )Wi\,{\sigma _{r\theta }}){\rm sin}(\theta )\} {\rm sin}(\theta )\,\textrm{d}\theta } = 0.\end{equation}Similarly, there is only one non-trivial component of the TFC:
\begin{equation}\int_0^{\rm \pi} {({{\dot{\gamma }}_{r\phi }} - (1 - \beta )\,Wi\,{\sigma _{r\phi }}){{\sin }^2}(\theta )\,\textrm{d}\theta } = 0.\end{equation}
In (6.1) and (6.2), the required integration along the azimuthal angle 
$\phi$ has been ignored since the flow is axisymmetric, and all the field variables are independent of 
$\phi$. Although (6.1) and (6.2) are applied at r = 1, i.e. at the surface of the swimmer, one can prove, by exploiting the integral form of the momentum balance, and separately the 
$\phi$-component of the momentum balance, that these are valid at any radial position 
$r \ge 1$. Consequently, a consistent solution must satisfy these equations and, indeed, this has been confirmed with the perturbation solution up to 
$O(W{i^4})$ for (6.1) and up to 
$O(W{i^5})$ for (6.2).
Furthermore, by taking into account that 
${\dot{\gamma }_{r\phi }} = r\partial ({v_\phi }/r)/\partial r$, one can integrate (6.2) along the radial coordinate r, and use the boundary conditions for 
${v_\phi }$ at the surface and far from the swimmer, (1.5a–c) and (1.7a–d), respectively, to obtain
\begin{equation}\varOmega ={-} \frac{3}{4}(1 - \beta )Wi\int_1^\infty {\int_0^{\rm \pi} {\frac{{{\sigma _{r\phi }}}}{r}{{\sin }^2}(\theta )\,\textrm{d}\theta \,\textrm{d}r} } .\end{equation}
Equation (6.2), evaluated at r = 1, can be used to determine 
$\varOmega$. Alternatively, 
$\Omega$ can be calculated using (6.3). The latter, however, can also be used as an independent check of the accuracy and correctness of an approximate solution of the governing equations (either a numerical or an asymptotic solution).
We proceed with (6.1) by distinguishing between the different sources of force and taking into account that 
${\dot{\gamma }_{rr}} = 2\partial {v_r}/\partial r$ and 
${\dot{\gamma }_{r\theta }} = \partial {v_\theta }/\partial r - {v_\theta }/r$, so that
\begin{gather}{F_p} ={-} \frac{1}{2}\int_0^{\rm \pi} {p\sin (2\theta )\,\textrm{d}\theta } ,\end{gather}
\begin{gather}{F_V} = \int_0^{\rm \pi} {\left( {\frac{{\partial {v_r}}}{{\partial r}}{\rm sin}(2\theta ) - \left( {\frac{{\partial {v_\theta }}}{{\partial r}} - \frac{{{v_\theta }}}{r}} \right){\rm sin}{^2}(\theta )} \right)} \,\textrm{d}\theta ,\end{gather}
\begin{gather}{F_E} ={-} (1 - \beta )Wi\int_0^{\rm \pi} {\left( {\frac{1}{2}{\sigma_{rr}}\sin (2\theta ) - {\sigma_{r\theta }}{{\sin }^2}(\theta )} \right)\,\textrm{d}\theta } ,\end{gather}
where 
${F_p}$ is the contribution due to the isotropic pressure, 
${F_V}$ is the contribution due to the viscous stresses, and 
${F_E}$ is the pure elastic contribution due to the polymer extra stress. In principle, all components are functions of the radial coordinate, but due to (6.1) their sum is zero, i.e. 
${F_p} + {F_V} + {F_E} = 0$. Also, although 
${F_p}$, 
${F_V}$ and 
${F_E}$ are given in terms of spherical components and coordinates, one can show that these quantities are equivalent to those reported in cylindrical coordinates in our recent publication (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020).
Furthermore, both the viscous and elastic contributions can be decomposed into a normal and a shear component, i.e. 
${F_V} = {F_{V,n}} + {F_{V,s}}$, where
\begin{equation}{F_{V,n}} = \int_0^{\rm \pi} {\frac{{\partial {v_r}}}{{\partial r}}{\rm sin}(2\theta )\,\textrm{d}\theta } ,\quad {F_{V,s}} ={-} \int_0^{\rm \pi} {\left( {\frac{{\partial {v_\theta }}}{{\partial r}} - \frac{{{v_\theta }}}{r}} \right)\,{\rm sin}{^2}(\theta )\,\textrm{d}\theta } ,\end{equation}
and 
${F_E} = {F_{E,n}} + {F_{E,s}}$, where
\begin{equation}{F_{E,n}} ={-} \frac{{(1 - \beta )}}{2}Wi\int_0^{\rm \pi} {{\sigma _{rr}}\sin (2\theta )\,\textrm{d}\theta } ,\quad {F_{E,s}} = (1 - \beta )Wi\int_0^{\rm \pi} {{\sigma _{r\theta }}{{\sin }^2}(\theta )\,\textrm{d}\theta } .\end{equation}Note, however, that the individual normal and shear contributions as shown in (6.8a,b) are different from those reported in Binagia et al. (Reference Binagia, Phoa, Housiadas and Shaqfeh2020), since the latter were defined in a cylindrical system.
Evaluating all force contributions on the surface of the swimmer, we find
\begin{gather}{F_P} ={-} \frac{1}{2}\int_0^{\rm \pi} {p{|_{r = 1}}\sin (2\theta )\,\textrm{d}\theta } ,\end{gather}
\begin{gather}{F_{V,n}} ={-} \frac{8}{3},\quad {F_{V,s}} = \frac{4}{3} - \int_0^{\rm \pi} {{{\left. {\frac{{\partial {v_\theta }}}{{\partial r}}} \right|}_{r = 1}}{\rm sin}{^2}(\theta )\,\textrm{d}\theta } ,\end{gather}
\begin{gather}{F_{E,n}} ={-} \frac{1}{2}(1 - \beta )Wi\int_0^{\rm \pi} {{\sigma _{rr}}{|_{r = 1}}\sin (2\theta )\,\textrm{d}\theta } ,\quad {F_{E,s}} = (1 - \beta )Wi\int_0^{\rm \pi} {{\sigma _{r\theta }}{|_{r = 1}}{{\sin }^2}(\theta )\,\textrm{d}\theta } ,\end{gather}
where (1.6) has been used in the evaluation of 
${F_{V,n}}$.
For a Newtonian fluid Wi = 0 and thus 
${F_{E,n}} = {F_{E,s}} = 0$. Also, substituting the pressure and velocity fields given by (1.9)–(1.12) in (6.9) and (6.10a,b) yields
\begin{equation}{F_P} = {\textstyle{2 \over 3}} - {U_N},\quad {F_{V,s}} = 4 - 2{U_N},\end{equation}
where the subscript ‘N’ is used to denote the Newtonian speed of the swimmer. Owing to the FFC, the balance of all contributions must be zero, i.e. 
${F_P} + {F_{V,n}} + {F_{V,s}} + {F_{E,n}} + {F_{E,s}} = 0$. The latter allows for the evaluation of the swimmer speed along the axial direction, 
${U_N} = 2/3$, as reported in the Introduction, (1.13a,b). Therefore, for a Newtonian fluid, (6.9)–(6.11a,b) give 
${F_p} = 0$, 
${F_{V,n}} = - 8/3,\;{F_{V,s}} = 8/3$ and 
${F_{E,n}} = {F_{E,s}} = 0$. Notice that, despite the fact that the slip and swirl parameters enter into the solution for the pressure and the velocity fields, they do not affect 
${U_N}$, as well as the force contributions 
${F_p}$ and 
${F_V}$.
Finally, we mention that, if we impose 
${U_N} = 0$ in (1.9)–(1.12) instead of the FFC, (6.9) and (6.10a,b) give 
${F_P} = 2/3,\;{F_{V,s}} = 4$ and 
${F_{V,n}} = - 8/3$, and therefore the net force on the body along the axial direction is 
$F = 2$. In other words, if the body is kept stationary, the slip boundary conditions on its surface will generate a flow field that will exert a positive, i.e. propulsive, force on the body. This type of analysis also reveals that both the pressure contribution, 
${F_P}$, and the shear component of the viscous forces, 
${F_{V,s}}$, generate thrust (with the shear force being six times larger than the pressure force), while the normal component of the viscous forces, 
${F_{V,n}}$, is purely resistive. Also, note that the FFC condition gives a zero rotation rate of the swimmer, 
$\varOmega = 0$, which is a consequence of the fact that rotation and translation in a simple Newtonian fluid are uncoupled and do not affect each other. A similar analysis for an elastic fluid has not been undertaken here but will follow in a future publication.
6.1. UCM and Oldroyd-B models
For a viscoelastic fluid 
$Wi > 0$, 
${F_p}$, 
${F_V}$ and 
${F_E}$ are 
$O(Wi)$ quantities and contribute to the FFC affecting the final speed of the swimmer. In order to comment on the force contributions and how the slip and swirl parameters affect U, the solution for the field variables on the surface of the swimmer is required. However, (6.9)–(6.11a,b) reveal an unexpected feature of the flow. In particular, for the UCM/Oldroyd-B models, and due to (1.6) and (5.4), one finds
\begin{equation}{{\mathsf{c}}_{rr}}{|_{r = 1}} = 1 - Wi(\xi + 4\cos (\theta ) + 3\xi \cos (2\theta )) - W{i^2}{\sigma _{rr}}{|_{r = 1}}.\end{equation}
As discussed extensively in § 3, 
${{\mathsf{c}}_{rr}}{|_{r = 1}}$ depends exclusively on the slip parameter ξ, the Weissenberg number Wi and the polar angle θ. Owing to (6.13), the same holds for 
${\sigma _{rr}}{|_{r = 1}}$, and therefore 
${F_{V,n}}$ and 
${F_{E,n}}$ do not depend on the swirl parameter ζ, which implies that they are invariants of the swirling flow. Consequently, the variation of the swirl parameter causes a redistribution of 
${F_P},{F_{V,s}}$ and 
${F_{E,s}}$ but keeps their sum unaltered and equal to 
${-} ({F_{V,n}} + {F_{E,n}})$. This suggests that the variation of U for the UCM/Oldroyd-B models is caused by the development and redistribution of the shear components of the rate-of-deformation and conformation tensors, 
${\dot{\gamma }_{r\theta }}$ and 
${{\mathsf{c}}_{r\theta }}$, respectively, due to the viscoelasticity of the ambient fluid.
In figure 7(a), we present results for the normalized velocity U/UN of a neutral squirmer as a function of Wi up to Wi = 0.225, where UN is the Newtonian speed given by (1.13a,b); in this case the Oldroyd-B model is defined up to Wi = 0.25 (see (3.3)). The results are shown for the no-swirl case (ζ = 0) and a swirl case with ζ = 3; the solid lines are theoretical predictions according to the accelerated technique, (5.8a,b), while the red lines with symbols are results from the large-scale simulations. For the no-swirl case one can see a small decrease of the velocity of the swimmer with increasing fluid elasticity. Also, the agreement between theory and simulations is excellent. As the swirl parameter increases to ζ = 3, both theory and simulations predict a monotonic increase of the squirmer velocity; note that the increase at Wi = 0.225 is similar to that previously calculated at a much larger Wi number for the Giesekus fluid (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020).
(a) Normalized squirmer velocity and (b) net rotation rate for a neutral squirmer (ξ = 0) using the Oldroyd-B model and viscosity ratio β = 0.5. The theoretical results are denoted with black solid lines, and the numerical results with red lines with symbols.
Moreover, we observe that the swirl generates a rigid-body rotation of the micro-swimmer, as figure 7(b) shows. The rotation rate of the swimmer is determined using the TFC, (6.2), evaluated at r = 1. Owing to the conservation of momentum along the azimuthal angle, 
$\Omega$ can also be determined from (6.3) as the volume integral of the shear elastic stress 
${\sigma _{r\phi }}$ weighted with the positive spatial function 
${\sin ^2}(\theta )/r$. The latter suggests that the major contribution of shear elastic forces to the development of rigid-body motion comes from the region close to the body. For the no-swirl case, the asymptotic solution shows that 
${\sigma _{r\phi }} = 0$ everywhere in the flow domain and thus 
$\varOmega = 0$. The same has been confirmed by the numerical simulations within the accuracy of the calculations. When the swirl parameter increases to 
$\zeta = 3$, at Wi = 0.225 the rotation rate has increased to nearly 0.3, showing that weak viscoelasticity in the presence of swirl significantly affects the rotation rate. Even at 
$Wi = 0.2$, the elasticity of the fluid coupled with the slip velocity at the surface of the body produce a tangential stress field, 
${\sigma _{r\phi }}$, everywhere negative in sign, the contours of which are presented in figure 8. We see that 
${\sigma _{r\phi }}$ reduces very quickly as the distance from the body increases. Near the swimmer, very large values in magnitude appear at the wake of the body. On the contrary, much lower absolute values are observed away from the south pole. This asymmetry of the tangential elastic stress with respect to the equator, i.e. to 
$\theta = {\rm \pi}/2$, contributes to the volume integral in (6.3), and causes the rotation of the swimmer. Finally, we observe that, precisely at the poles, 
${\sigma _{r\phi }}(r,\theta = 0) = {\sigma _{r\phi }}(r,\theta = {\rm \pi}) = 0$, confirming equation (2.11a,b) and the theoretical predictions reported in § 4, regarding the purely extensional character of the flow at the poles.
The σrφ contours for a neutral squirmer (ξ = 0) using the Oldroyd-B model with Wi = 0.2, β = 0.5 and swirl parameter ζ = 3.
In figure 9, we present the force contributions 
${F_p},{F_V}$ and 
${F_E}$ on a neutral swimmer as functions of Wi. First, we observe an excellent agreement between the accelerated asymptotic theory and the numerical results. Recalling that 
${F_{V,n}}$ and 
${F_{E,n}}$ are invariant quantities for the Oldroyd-B model, the differences between the no-swirl and swirl cases is due to the variation of 
${F_p}$, 
${F_{V,s}}$ and 
${F_{E,s}}$. Looking closely at the magnitudes and signs of these quantities, we can extract useful information about the type of forces acting on the swimmer. Indeed, the pressure contribution generates a major thrust compared to the Newtonian fluid and becomes even more propulsive in the swirl case. The total viscous contribution goes from propulsive to slightly resistive; this implies that the shear viscous contribution changes sign from positive (
${F_{V,s}} > 0$ for 
$\zeta = 0$) to negative (
${F_{V,s}} < 0$ for 
$\zeta = 3$). Lastly, the total elastic contribution, which is absent for a Newtonian fluid, is purely resistive in nature, and increases substantially in magnitude in the swirl case. Both the normal and shear components, 
${F_{E,n}}$ and 
${F_{E,s}}$, respectively, are negative, but, due to the fact that 
${F_{E,n}}$ does not depend on the swirl parameter, this implies that the shear contribution of the elastic stress 
${F_{E,s}}$ is affected a lot by the increase of swirl. Based on these observations, we conclude that, for an elastic, UCM/Oldroyd-B type fluid, the pressure contribution is the only propulsive force on the swimmer with swirl, adequate to overcome the resistance caused by the viscous and elastic forces, and is the major mechanism that generates thrust and increases the speed of the swimmer.
Force contributions for a neutral squirmer (ξ = 0) using the Oldroyd-B model and viscosity ratio β = 0.5. The theoretical results are denoted with solid lines, and the numerical results with dotted lines with symbols. (a) No-swirl case (ζ = 0) and (b) swirl case with ζ = 3.
Further evidence on the propulsive role of the isotropic pressure can be seen in figure 10, where the contours of the pressure are drawn. Comparing the results between ζ = 0 and ζ = 3, two major observations can be made. First, the magnitude of the contours is much larger in the swirl case; and second, the pressure difference between the south and north poles is greatly intensified in the swirl case. Note that the Weissenberg number is really very small (Wi = 0.2), revealing again the great effect that even weak viscoelasticity can produce on the flow. We also notice that 
$p(r,{\rm \pi} ) \gg p(r,0)$, in contrast to the no-swirl case, where 
$p(r,{\rm \pi} )$ is only slighter larger than 
$p(r,0)$, and in contrast to the Newtonian pressure, for which 
$p(r,{\rm \pi} ) = p(r,0)$ (see (1.12) with U = 2/3). The large pressure difference between the poles generates the major thrust along the z-direction which dominates the resistive forces on the body.
Pressure field for a neutral squirmer (ξ = 0) using the Oldroyd-B model with Wi = 0.2 and β = 0.5. The contours are shown in the YZ plane. (a) No-swirl case (ζ = 0) and (b) swirl case with ζ = 3.
Finally, the resistive nature of the elastic forces can be realized from the results shown in figure 11 in conjunction with the definitions of 
${F_{E,n}}$ and 
${F_{E,s}}$, given by (6.11a,b). In particular, 
${\sigma _{rr}}(r = 1,\theta )$ and 
${\sigma _{r\theta }}(r = 1,\theta )$ are presented as functions of 
$x = \cos (\theta )$, for both the no-swirl (black, solid lines) and swirl cases (red, dotted lines). It is seen that very large values are developed close to 
$x\, = - \!1$ 
$(\theta = {\rm \pi})$, and hence both 
${\sigma _{rr}}$ and 
${\sigma _{r\theta }}$ produce resistance to the body, most of which comes from the stress close to the south pole.
Numerical results for the (a) normal elastic stress σrr and (b) tangential elastic stress σr
$_{\theta}$ at the surface of a neutral squirmer (ξ = 0) for an Oldroyd-B fluid with viscosity ratio β = 0.5.
6.2. Giesekus and FENE-P models
As has already been discussed in the previous sections, for the Giesekus and the FENE-P models, (2.5)–(2.10) show that all the components of the conformation tensor on the surface of the swimmer are coupled and affect each other. As the mobility parameter 
${a_m}$ of the Giesekus model, or the maximum extensibility parameter 
$L$ of the FENE-P model, increases, these models deviate significantly from the UCM/Oldroyd-B models, resulting in a substantial redistribution of the force contributions to the FFC. This actually leads to different mechanisms of speed enhancement even in the low-Weissenberg-number regime.
For the Giesekus model, U/UN and 
$\Omega$ are shown in figure 12 as functions of Wi, for mobility parameter αm = 0.1 and viscosity ratio β = 0.5. In the range of the Weissenberg number shown, the Giesekus model exhibits negligible shear thinning. The no-swirl case and a swirl case with ζ = 3 have been chosen. The accelerated asymptotic theory is presented with solid black lines, while the numerical results are denoted with red lines with symbols. The agreement between theory and simulations is excellent for the no-swirl case, while for the swirl case the agreement is very good up to 
$Wi \approx 0.1$ for the normalized velocity, and up to 
$Wi \approx 0.15$ for the rotation rate. As Wi increases, the theory starts to deviate from the numerical results, although the degree of the deviation is reasonable, and the right trends are predicted. On the contrary, the original perturbation solutions (not shown in figure 12) give completely erroneous results for 
$Wi$ larger than approximately 0.1. We observe that the speed enhancement is small, and the rotation rate is significant. Compared to the results for the Oldroyd-B model, shown in figure 7, the speed enhancement is reduced, and the same holds for the rotation rate. Obviously, the choice of the constitutive model plays an important role in understanding swimming (with or without swirl) in an elastic fluid.
(a) Normalized squirmer velocity and (b) net rotation rate for a neutral squirmer (ξ = 0) using the Giesekus model with αm = 0.1 and β = 0.5. Black solid lines, accelerated asymptotic results; red lines with symbols, numerical results.
As explained before for the UCM/Oldroyd-B models, the rotation rate of the swimmer under creeping flow conditions for a viscoelastic fluid is produced exclusively by the development of the shear elastic stress 
${\sigma _{r\phi }}$. Since the rotation rate for the Giesekus model is smaller than that of the Oldroyd-B model, one expects the same to be true for the magnitude of 
${\sigma _{r\phi }}$. Indeed, if we compare the 
${\sigma _{r\phi }}$ contours, shown in figure 13, with those shown in figure 8, we see that the Giesekus model predicts smaller stresses tangential to the body than the Oldroyd-B model.
Tangential stress σrϕ contours for the Giesekus model with Wi = 0.2, αm = 0.1 and β = 0.5. The contours are shown in the YZ plane for a neutral squirmer (ξ = 0) with ζ = 3.
In figure 14, we present the force contributions 
${F_p},{F_V}$ and 
${F_E}$ on a neutral swimmer as functions of Wi; all parameters are the same as in the previous two figures. The agreement between the accelerated asymptotic theory and the numerical results for the no-swirl case is worth noting, while for 
$\zeta = 3$ reasonable differences are observed. In all cases, however, the accelerated asymptotic theory predicts the right trends for all quantities. As mentioned above as well, the additional quadratic term involving the conformation tensor in the original form of the constitutive equation, (2.3), has the consequence that the normal contribution to the force due to the elastic stress is not invariant to the swirl in the flow; only the normal viscous contribution 
${F_{V,n}}\, = - 8/3$ remains constant. Again, by looking closely at the details, we can extract information about the type of forces acting on the swimmer and compare this with the results for the Oldroyd-B model. First, we observe that the pressure contribution does not generate a major thrust compared to UCM/Oldroyd-B fluids, as well as the fact that the introduction of swirl has a negligible effect on the pressure contribution. Second, the total viscous contribution goes from propulsive to resistive, namely, the shear viscous contribution changes sign from positive (
${F_{V,s}} > 0$ for 
$\zeta = 0$) to negative (
${F_{V,s}} < 0$ for 
$\zeta = 3$). Lastly, the total elastic contribution, 
${F_E} = {F_{E,n}} + {F_{E,s}}$, changes type and becomes propulsive. Since, 
${F_{E,n}}$ is negative as the mobility parameter goes to zero, we conclude that the non-isotropic term of the Giesekus model is responsible for the development of substantial normal stresses that vary significantly with the increase of the mobility parameter. It also appears that the elevated normal elastic stresses hinder the development of the isotropic pressure. Consequently, the major propulsive character of pressure, previously identified for the UCM/Oldroyd-B fluids, diminishes and the speed enhancement is minimized.
Force contributions for a neutral squirmer (ξ = 0) using the Giesekus model with αm = 0.1 and β = 0.5. The theoretical results are denoted with solid lines, and the numerical results with dotted lines with symbols. (a) No-swirl case (ζ = 0) and (b) swirl case with ζ = 3.
The contours of the isotropic pressure can be seen in figure 15. First, we mention that for the no-swirl case the magnitude of the pressure, and its spatial variations, are even smaller than those in the Oldroyd-B model (shown in figure 10a). When swirl is introduced, ζ = 3, slightly larger values of the pressure and its gradient are observed compared to the no-swirl case. However, these are much smaller than in the Oldroyd-B model (shown in figure 10b) and thus the contribution of the pressure to the speed of the swirling swimmer is inconsequential.
Pressure field for a neutral squirmer (ξ = 0) using the Giesekus model with Wi = 0.2, αm = 0.1, β = 0.5 and swirl parameter ζ = 3.
In figure 16, we show 
${\sigma _{rr}}(r = 1,\theta )$ and 
${\sigma _{r\theta }}(r = 1,\theta )$ as functions of 
$x = \cos (\theta )$, for both the no-swirl (black, solid lines) and swirl cases (red, dotted lines). Very steep gradients are observed at the south pole, which cause the normal component of the elastic stress to switch from negative to positive values. Recalling the definitions of 
${F_{E,n}}$ and 
${F_{E,s}}$, given by (6.11a,b), we can easily confirm that the tangential elastic stress 
${\sigma _{r\theta }}$ remains resistive in character for both the no-swirl and swirl cases. On the contrary, the normal elastic stress 
${\sigma _{rr}}$ is resistive for the no-swirl case, but changes to propulsive for swirling swimmers. Finally, it is worth noting that, although both the pressure contribution and the total elastic contribution generate thrust to the swirling body, their magnitude is small, resulting in a minor speed enhancement compared to the simple Newtonian fluid with the same viscosity.
Numerical results for the (a) normal elastic stress σrr and (b) tangential elastic stress σrθ as functions of x = cos(θ) at the surface of the body. A neutral squirmer (ξ = 0) in a Giesekus fluid with Wi = 0.2, αm = 0.1 and β = 0.5 is shown.
Finally, in figures 17 and 18, we present results for the FENE-P model with ξ = 0, L = 10 (χ = 1/102), β = 0.5 and ζ = 0 and ζ = 3. In particular, U/UN and 
$\Omega$ are shown in figure 17, and the force contributions 
${F_p},{F_V}$ and 
${F_E}$, in figure 18; all quantities are presented as functions of Wi. In the range of the Weissenberg number shown, 
$0 \le Wi \le 1/4$, the FENE-P model exhibits negligible shear thinning. One can see that the results are very similar to those for the Oldroyd-B model. Although the agreement between theory and simulations is not as good as for the Oldroyd-B model, the general predictions and trends remain the same. Therefore, the UCM/Oldroyd and FENE-P models exhibit the same mechanism for speed enhancement and rotation rate of the body. These results indicate, for one more time, that the choice of the constitutive model plays an important role in understanding swimming (with or without swirl) in an elastic fluid.
(a) Normalized squirmer velocity and (b) net rotation rate for a neutral squirmer (ξ = 0) using the FENE-P model with L = 10 and viscosity ratio β = 0.5. The theoretical results are denoted with black solid lines, and the numerical results with red solid lines with symbols.
Force contributions for a neutral squirmer (ξ = 0) using the FENE-P model with L = 10 and β = 0.5. The theoretical results are denoted with solid lines, and the numerical results with dotted lines with symbols. (a) No-swirl case (ζ = 0) and (b) swirl case with ζ = 3.
7. Conclusions
In our previous publication (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020), we discussed the fact that the low-Weissenberg-number region is appropriate for real micro-organism swimming, because of the characteristic swimming speeds, the size of the organisms and the fluids in which they swim. Here, we thoroughly investigated this region, both analytically and numerically. In summary, we can conclude as follows.
(a) We proved that the UCM/Oldroyd-B models become singular at a small critical Weissenberg number which depends on the slip parameter ξ. For a puller with ξ = 1/3, one finds that the critical Wi is 3/8 and decreases for any other value of the slip parameter. For a neutral swimmer, the critical Wi is 1/4. The singularity is revealed based on the analysis at the surface of the body, and it is the first time that such a singularity for the viscoelastic flow past a spherical body has been predicted. The same analysis also shows that the flow at the poles is purely extensional in character. Unexpectedly, the swirl parameter does not affect the solution for the conformation tensor at the poles.
(b) We proved that the singularity predicted for the UCM/Oldroyd-B fluids is completely removed when regularized and more realistic constitutive equations are employed such as the Giesekus and FENE-P models. However, the purely extensional character of the flow at the poles still holds.
(c) For the UCM/Oldroyd-B models, and based on the exact analytical solution for the conformation tensor at the poles, we found analytically the radius of convergence of the series solutions with respect to the Weissenberg number. We showed that the radius of convergence is even smaller than the critical Wi at which the UCM/Oldroyd-B models become singular. Surprisingly, and although the Giesekus and FENE-P models are valid for any value of Wi, the radius of convergence of the corresponding series solutions is slightly smaller than in the UCM/Oldroyd-B fluids. Therefore, the perturbation solution(s) can be used as a good approximation of the exact solution only for very small values of Wi (
$Wi\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel<\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }0.1$). This information in conjunction with that reported in points (a) and (b) above fully explains the fact that other researchers in the field (Datt et al. Reference Datt, Natale, Hatzikiriakos and Elfring2017; Datt & Elfring Reference Datt and Elfring2019, Reference Datt and Elfring2020) have reported large changes to their results, even at relatively small values of Wi, when higher-order corrections are taken into account. It also explains the large differences between the simulation and asymptotic results seen in our recent publication (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020) in the intermediate- and high-Weissenberg -number regions.(d) When techniques that accelerate the convergence of series are applied, the transformed solutions are more accurate than the original perturbation solutions. Even when the agreement with the simulation results is not excellent, the transformed solutions predict the right trends for all the integral and/or constant quantities of interest in the range
$0 \le Wi < W{i_\rho }$, where $W{i_\rho }$ is given by (3.6).(e) Both theory and simulations reveal that all viscoelastic models, at a sufficiently large swirl parameter ζ, predict a speed enhancement compared to the speed of the swimmer in a Newtonian fluid with the same viscosity, UN = 2/3. Use of the Oldroyd-B model results in a prediction of the largest increase in U. Also, all viscoelastic models predict a rigid-body rotation of the micro-swimmer, and again the Oldroyd-B model, among all models employed, predicts the largest
$\Omega$.(f) We showed that, for an Oldroyd-B fluid, the speed enhancement for a swimmer with swirl is driven primarily by an increase in the force due to pressure. The same was shown for the FENE-P model, namely the FENE-P model is only a small modification to the Oldroyd-B model. In contrast, while speed enhancement is also seen for a Giesekus fluid, the propulsion is instead created from a large increase in the normal elastic force. For this model, the role of pressure remains propulsive, but is much smaller than that predicted using the Oldroyd-B and FENE-P models.
We emphasize that, although the Oldroyd-B model is the most basic and fundamental differential fluid mechanics model for viscoelastic fluids, these results reveal phenomena and mechanisms, for both non-swirling and swirling swimmers, not previously identified by other researchers in the field. From our analysis, we concluded that the choice of the constitutive model is absolutely crucial for understanding the locomotion of micro-swimmers in viscoelastic fluids, and therefore guidance from experiments is necessary to determine which model is best suited to predict the major features of these types of flows. We also revealed that the low-Weissenberg-number region (i.e. of weak viscoelasticity) is much more important in understanding swimming in an elastic fluid than generally believed thus far. Finally, in a future publication we will present results on the time-dependent spherical squirmer model, addressing thoroughly the effect of the transient, frequency-related, terms in the governing equations; that work is under way.
Acknowledgements
J.P.B. is supported by a National Science Foundation (NSF) Graduate Research Fellowship (Grant No. DGE - 1656518). This work is also supported partially by NSF Grant No. CBET 1803765.
Declaration of interests
The authors declare no conflict of interest.
Appendix A
For the Oldroyd-B model, and using the symbol 
$\eta \equiv 1 - \beta$ for brevity, the high-order terms given in (5.6) and (5.7) are
\begin{gather}{\tilde{U}_2} ={-} \frac{{772}}{{2145}} + {\zeta ^2}\left( {\frac{{468}}{{385}} - \frac{{4392}}{{5005}}\eta } \right),\end{gather}
\begin{gather}{{\tilde{U}}_4} ={-} \dfrac{{314744}}{{980343}} + \dfrac{{267001496}}{{233648415}}\eta - \dfrac{{381076}}{{855855}}{\eta ^2}\nonumber\\ \quad +\, {\zeta ^2}\left( { - \dfrac{{135928}}{{75411}} - \dfrac{{463142193656}}{{83412484155}}\eta + \dfrac{{2681778613412}}{{417062420775}}{\eta^2} - \dfrac{{2216920922}}{{1168242075}}{\eta^3}} \right)\nonumber\\
\quad +\, {\zeta ^4}\eta \left( { - \dfrac{{922566056}}{{237642405}} + \dfrac{{1570635736}}{{468083525}}\eta - \dfrac{{118941418678}}{{112540653225}}{\eta^2}} \right),\nonumber \end{gather}
\begin{gather}{\tilde{\varOmega }_1} = \frac{6}{5},\end{gather}
\begin{gather}{\tilde{\varOmega }_3} ={-} \frac{{2052}}{{455}} + \frac{{534}}{{455}}\eta + \frac{{48}}{{65}}{\eta ^2} + {\zeta ^2}\eta \left( { - \frac{{202914}}{{35035}} + \frac{{135528}}{{35035}}\eta } \right),\end{gather}
\begin{gather} {{\tilde{\varOmega }}_5} = \dfrac{{553568}}{{19845}} + \dfrac{{95095649288}}{{14424715305}}\eta - \dfrac{{3296321288824}}{{216370729575}}{\eta ^2} + \dfrac{{115117929392}}{{30910104225}}{\eta ^3}\nonumber\\ \quad +\, \dfrac{{937088}}{{1117935}}{\eta ^4} + {\zeta ^2}\eta \left( { - \dfrac{{172245369776}}{{213165237285}} + \dfrac{{1195093253533414}}{{19711928261025}}\eta } \right.\nonumber\\ \quad \left. { -\, \dfrac{{117780004568598902}}{{1793785471753275}}{\eta^2} + \dfrac{{181723992108056408}}{{8968927358766375}}{\eta^3}} \right)\nonumber\\ \quad +\, {\zeta ^4}{\eta ^2}\left( {\dfrac{{179145233734466}}{{9118735150525}} - \dfrac{{9796607184217891}}{{465055492676775}}\eta + \dfrac{{53424844660340698}}{{6975832390151625}}{\eta^2}} \right). \end{gather}
Open access
