1. Introduction
The study of the motion of living matter in fluids is a cornerstone of biological fluid mechanics, and important to the design of synthetic active matter (Childress Reference Childress1981; Degen Reference Degen2014; Alapan et al. Reference Alapan, Yasa, Schauer, Giltinan, Tabak, Sourjik and Sitti2018). Many cellular organisms exhibit some form of self-propulsion (Bray Reference Bray2000; Lauga Reference Lauga2020), which is usually achieved by flagella or cilia acting on the surrounding fluid (Brennen & Winet Reference Brennen and Winet1977). The motion of microscopic organisms has been widely studied (Lauga & Powers Reference Lauga and Powers2009; Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013). At this scale, inertial forces are negligible, i.e. the Reynolds number 
$Re$ is small. Swimming at large 
$Re$, where inertia dominates, has also been extensively investigated (Becker et al. Reference Becker, Masoud, Newbolt, Shelley and Ristroph2015; Maertens, Gao & Triantafyllou Reference Maertens, Gao and Triantafyllou2017). However, the swimming of mesoscale organisms, at 
$Re$ of order unity, is relatively unexplored (Klotsa Reference Klotsa2019). This is due to the simplifications that can be made in Stokes flows 
$(Re\ll 1)$ and in Euler flows 
$(Re\gg 1)$, by neglecting inertial and viscous forces, respectively, being invalid at intermediate 
$Re$.
Simplified, or reduced-order, models have been proposed to analyse the locomotion of swimming organisms. A popular example is the squirmer model developed by Lighthill (Reference Lighthill1952) and Blake (Reference Blake1971), wherein self-propulsion is achieved by prescription of a surface velocity, or swimming gait, at the instantaneous surface of the squirmer. Most studies have focused on axisymmetric and impenetrable spherical squirmers, for which the fluid velocity at the surface, and relative to it, can be represented by a modal expansion:
\begin{equation} \text{relative surface velocity} ={-}2\sum_{n=1}^\infty\frac{B_n}{n(n+1)}P_n^1(\cos\theta)\hat{\boldsymbol{e}}_\theta, \end{equation}
in which 
$n$ denotes mode number, 
$B_n$ the corresponding mode amplitude, 
$\theta$ the polar angle (and 
$\hat {\boldsymbol {e}}_{\theta }$ its associated unit vector) measured from an arbitrarily chosen ‘forward’ direction along the symmetry axis of the squirmer, and 
$P_n^m(\cos \theta )$ the associated Legendre polynomials. Note that 
$n$ odd (even) implies fore–aft antisymmetric (symmetric) squirming, e.g. 
$P_1^1(\cos \theta )=-\sin \theta$, 
$P_2^1(\cos \theta )=-3\cos \theta \sin \theta$.
The squirmer model has been instrumental in examining various aspects of swimming at zero Reynolds number, including hydrodynamic interactions (Llopis & Pagonabarraga Reference Llopis and Pagonabarraga2010), locomotion in viscoelastic fluids (Zhu et al. Reference Zhu, Do-Quang, Lauga and Brandt2011; Zhu, Lauga & Brandt Reference Zhu, Lauga and Brandt2012) and nutrient uptake (Magar, Goto & Pedley Reference Magar, Goto and Pedley2003; Magar & Pedley Reference Magar and Pedley2005). In that regime, the flow field and squirmer swimming velocity induced by the swimming-gait modal expansion (1.1) can be obtained by superposing those motions induced by each mode separately. Only the first ‘dipolar’ (fore–aft antisymmetric) mode contributes to a non-zero swimming velocity (Blake Reference Blake1971). The second ‘quadrupolar’ (fore–aft symmetric) mode contributes a stresslet flow, whose sign distinguishes between ‘puller’ and ‘pusher’ swimmers (Ishikawa & Pedley Reference Ishikawa and Pedley2007): pusher corresponds to negative 
$B_2$ (relative surface velocity from the equator to the symmetry axis, or poles), and puller to positive 
$B_2$ (relative surface velocity from the poles to the equator).
Beyond the Stokes-flow regime, the nonlinear nature of inertia implies that the above superposition principle no longer holds. Previous studies of squirmers at non-zero 
$Re$ have focused on squirmers whose swimming gait involves only the first two, dipolar and quadrupolar modes in (1.1). Wang & Ardekani (Reference Wang and Ardekani2012) developed an asymptotic expansion through 
$O(Re)$ for the swimming speed of a two-mode squirmer at small 
$Re$, which Khair & Chisholm (Reference Khair and Chisholm2014) later extended to 
$O(Re^2)$. (In these works, 
$Re$ is defined based upon the magnitude of the dipolar squirming mode.) Chisholm et al. (Reference Chisholm, Legendre, Lauga and Khair2016) performed numerical computations of such a two-mode squirmer for 
$0< Re< 1000$, bridging the gap between Stokes and Euler flows. They found that, in contrast to the Stokes-flow regime, the swimming speed for a given non-zero value of 
$B_1$ is affected by the value of 
$B_2$. For 
$B_2<0$ (pusher at zero 
$Re$), increasing 
$Re$ leads to a monotonic increase in the swimming speed and the axisymmetric flow remains stable to at least 
$Re=1000$. For 
$B_2>0$ (puller at zero 
$Re$), the swimming speed has a non-monotonic dependence on 
$Re$ and the axisymmetric flow becomes unstable at sufficiently large 
$Re$.
The scenario considered by Chisholm et al. (Reference Chisholm, Legendre, Lauga and Khair2016) was the effect of inertia on a squirmer that is motile at 
$Re=0$, namely a fore–aft asymmetric squirmer with 
$B_1\ne 0$. It is intuitive that fore–aft asymmetric squirmers that are non-motile at 
$Re=0$ (i.e. squirmers with 
$B_1=0$ but 
$B_n\ne 0$ for at least one odd, non-unity value of 
$n$), generally become motile for 
$Re>0$, though with their swimming speed vanishing as 
$Re\searrow 0$; this is readily demonstrable by adapting the small-
$Re$ analyses of Wang & Ardekani (Reference Wang and Ardekani2012) and Khair & Chisholm (Reference Khair and Chisholm2014). Here, we ask whether inertia can also enable fore–aft symmetric squirmers to swim via nonlinear symmetry breaking (figure 1). To this end, we shall numerically study the effects of inertia on quadrupolar squirmers.
Figure 1. Schematic of spontaneous symmetry breaking of a fore–aft symmetric squirmer. (a) Symmetric steady state, wherein the squirmer is stationary. (b) Symmetry-broken steady state, wherein the squirmer swims. Blue arrows: symmetrically prescribed surface-slip velocity (we show poles-to-equator squirming as in the case of a quadrupolar puller). Green arrows: induced flow in a frame comoving with the squirmer.
We were led to this question by analogy with recent discoveries of symmetry-breaking locomotion of isotropically active droplets and particles (Michelin, Lauga & Bartolo Reference Michelin, Lauga and Bartolo2013; Izri et al. Reference Izri, Van Der Linden, Michelin and Dauchot2014) and Leidenfrost drops (Bouillant et al. Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and Quéré2018). While those examples rely on nonlinear coupling between hydrodynamics and other physics, inertia alone is well known to result in symmetry breaking in many flow scenarios (e.g. the axial asymmetry of wakes downstream of a sufficiently fast-moving blunt body). In fact, spontaneous locomotion enabled by inertial symmetry breaking has already been demonstrated for flapping bodies at sufficiently high Reynolds numbers (Vandenberghe, Zhang & Childress Reference Vandenberghe, Zhang and Childress2004; Alben & Shelley Reference Alben and Shelley2005; Vandenberghe, Childress & Zhang Reference Vandenberghe, Childress and Zhang2006).
2. Problem formulation
We consider a spherical, axisymmetric squirmer of radius 
$a$ that is suspended in a Newtonian, incompressible fluid of density 
$\rho$ and viscosity 
$\eta$. Our goal is to study the effects of inertia on a quadrupolar squirmer, whose swimming gait is described by just the second mode in (1.1), with 
$B_2$ constant. For the sake of demonstration, we include in the formulation below the possibility to perturb the fore–aft symmetric squirmer forwards/backwards by means of a time-localised dipolar-mode (
$n=1$) contribution in (1.1). For simplicity, we assume axial symmetry, and that the squirmer is density matched to the suspending fluid.
Henceforth, velocities are normalised by 
$|B_2|$, the pressure and hydrodynamic stress tensor by 
$\eta |B_2|/a$ and time by 
$\rho a^2/\eta$. These scales are associated with the Reynolds number 
$Re=\rho a|B_2|/\eta$. We work in a frame moving with the particle, and employ spherical coordinates 
$(r,\theta,\phi )$ and associated unit vectors 
$(\hat {\boldsymbol {e}}_r,\hat {\boldsymbol {e}}_{\theta },\hat {\boldsymbol {e}}_{\phi })$, with the origin at the particle centre and 
$\theta$ the polar angle from the ‘forward’ 
$\hat {\boldsymbol {\imath }}$ direction along the axis of symmetry.
With these conventions, the fluid flow is governed by the continuity and momentum equations
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$
$$\begin{gather}\frac{\partial\boldsymbol{u}}{\partial t}+\hat{\boldsymbol{\imath}}\frac{{\rm d}U}{{\rm d}t}+Re \,\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}p- \nabla^2\boldsymbol{u}=\boldsymbol{0}, \end{gather}$$
in which 
$\boldsymbol {u}$ and 
$p$ are the fluid-velocity and pressure fields, 
$\hat {\boldsymbol {\imath }}U$ is the squirmer velocity in the laboratory frame of reference and 
$t$ denotes time. The second term on the left-hand side of (2.1b) represents the fictitious force due to the reference frame's acceleration. With reference to (1.1), the velocity on the squirmer's boundary is
\begin{equation} \boldsymbol{u}=v_s\hat{\boldsymbol{e}}_{\theta} \quad \text{at }r=1, \quad \text{where} \ v_s={\pm}\sin\theta\cos\theta+\lambda(t)\sin\theta. \end{equation}
Here, the sign is that of 
$B_2$ and thus indicates a pusher or puller, and 
$\lambda (t)$ corresponds to the time-localised dipolar perturbation. (The pusher–puller terminology is based on Stokes-flow theory, where the sign of 
$B_2$ determines the directionality of the induced force dipole. We have numerically checked that in all cases presented herein inertia does not affect that directionality.) Far from the squirmer,
\begin{equation} \boldsymbol{u}\to-\hat{\boldsymbol{\imath}}U\quad \text{as}\ r\to\infty. \end{equation}Lastly, the squirmer velocity is coupled to the induced flow via Newton's second law
\begin{equation} \frac{4{\rm \pi}}{3}\frac{{\rm d}U}{{\rm d}t}=\hat{\boldsymbol{\imath}}\boldsymbol{\cdot}\oint_{r=1}\boldsymbol{\mathsf{N}}\boldsymbol{\cdot}\hat{\boldsymbol{e}}_r\,{\rm d}S, \end{equation}
wherein 
$\boldsymbol{\mathsf{N}}=-p\boldsymbol{\mathsf{I}}+\boldsymbol {\nabla } \boldsymbol {u}+(\boldsymbol {\nabla } \boldsymbol {u})^{\rm T}$ is the hydrodynamic stress tensor, 
$\boldsymbol{\mathsf{I}}$ being the identity tensor and the superscript T the tensor transpose, and 
$dS$ is a dimensionless areal element.
3. Methodology
We next overview our numerical approach to solving the problem formulated in § 2. We shall perform both time-dependent simulations – initial conditions will be specified later – and steady-state calculations. Some readers may wish to skip to § 4, where we present and discuss our results.
3.1. Stream function–vorticity formulation
It is convenient to represent the incompressible flow field 
$\boldsymbol {u}$ in terms of its stream function 
$\psi$:
\begin{equation} \boldsymbol{u}=\frac{1}{r^2\sin\theta}\frac{\partial \psi}{\partial \theta}\hat{\boldsymbol{e}}_r -\frac{1}{r\sin\theta}\frac{\partial \psi}{\partial r}\hat{\boldsymbol{e}}_\theta, \end{equation}whereby the continuity equation (2.1a) is trivially satisfied. The momentum equation (2.1b) can then be written as
\begin{equation} \frac{\partial\boldsymbol{\omega}}{\partial t}+Re[\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\omega}-\boldsymbol{\omega} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}]-\nabla^2\boldsymbol{\omega}=\boldsymbol{0}, \end{equation}
wherein the vorticity field 
$\boldsymbol {\omega }=\boldsymbol {\nabla }\times \boldsymbol {u}$ is azimuthal, i.e. 
$\boldsymbol {\omega }=\hat {\boldsymbol {e}}_{\phi }\omega$, the azimuthal component 
$\omega$ being given in terms of the stream function as
\begin{equation} \omega={-}\frac{1}{r\sin\theta}\left\{\frac{\partial^2\psi}{\partial r^2}+\frac{\sin\theta}{r^2}\frac{\partial}{\partial\theta}\left(\frac{1}{\sin\theta}\frac{\partial\psi}{\partial \theta}\right)\right\}. \end{equation}The boundary and far-field conditions, (2.2) and (2.3), become
$$\begin{gather} \psi=0\quad \text{and}\quad \frac{\partial\psi}{\partial r} ={-}v_s\sin\theta \quad \text{at }r=1, \end{gather}$$
$$\begin{gather}\frac{\partial\psi}{\partial r}\sim{-}Ur\sin^2\theta\quad \text{and}\quad \frac{\partial\omega}{\partial r}\to 0\quad \text{as}\ r\to\infty. \end{gather}$$Lastly, the hydrodynamic-force integral in (2.4) can be expressed as
\begin{align} \hat{\boldsymbol{\imath}}\boldsymbol{\cdot}\oint_{r=1}\boldsymbol{\mathsf{N}}\boldsymbol{\cdot} \hat{\boldsymbol{e}}_r\,{\rm d}S={-}\frac{4{\rm \pi}}{3}\frac{{\rm d}U}{{\rm d}t}+ {\rm \pi}\int_0^{\rm \pi}\left[\frac{Re}{2}v_s^2\sin(2\theta)+\left(\frac{\partial(r\omega)}{\partial r}-2\omega\right)\sin^2\theta\right]{\rm d}\theta.\end{align}
The drag formula (3.5) has been adapted from Khair & Chisholm (Reference Khair and Chisholm2014). The first term takes into account the fictitious force affecting the pressure due to a non-inertial frame of reference. This term acts as an added mass in (2.4), effectively doubling the squirmer's inertia. The contribution proportional to 
$v_s^2$ vanishes for a fore–aft symmetric squirmer, thus in the absence of the dipolar perturbation in the present formulation.
3.2. Numerical scheme
Our method for solving the above formulation as an initial-value problem involves the following two steps. The first consists of solving for the stream function and vorticity fields at a given time, given the squirmer velocity and the vorticity field at a previous time. The time derivative in (3.2) is discretised by the backward Euler method:
\begin{equation} \boldsymbol{\omega}^{(t)}=\boldsymbol{\omega}^{(t-\Delta t)}+\Delta t[\nabla^2\boldsymbol{\omega}^{(t)}-Re(\boldsymbol{u}^{(t)}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{\omega}^{(t)}-\boldsymbol{\omega}^{(t)}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^{(t)})], \end{equation}
where 
$\Delta t$ is a time step. Equations (3.3) and (3.4) are written at the present time, except that the previous-time squirmer velocity is used in (3.4b). The resulting nonlinear flow problem is solved by a spectral-element method adapted from Chisholm et al. (Reference Chisholm, Legendre, Lauga and Khair2016), which employs the Galerkin method of weighted residuals (Karniadakis & Sherwin Reference Karniadakis and Sherwin2005). The two-dimensional basis set is obtained from a tensor product of one-dimensional Lagrange polynomials of order 
$N=8$, leading to 
$(N+1)^2$ degrees of freedom per node, and we make use of the Gauss–Lobatto quadrature to integrate over each parametric subdomain. The mesh is generated using Gmsh (Geuzaine & Remacle Reference Geuzaine and Remacle2009), employing a half-ring configuration of inner radius 
$R_i=1$ and outer radius 
$R_o=200$. The number of nodes is 
$28^2$, with a radial geometric progression outwards with a factor of 
$1.25$, while in the polar direction, a factor of 
$1.1$ is used and the progression direction is towards 
${\rm \pi} /2$. The radial and polar progressions are implemented to handle sharp changes near the squirmer and the wake that occur near the symmetry axis. The resulting set of nonlinear algebraic equations is solved using the Newton–Raphson algorithm. For further details about the discretisation, and validation of the method in different scenarios, the reader is referred to Kailasham & Khair (Reference Kailasham and Khair2022) and Cobos & Khair (Reference Cobos and Khair2023).
In the second step, we update the squirmer speed according to (2.4),
\begin{align} U^{( t)}=U^{(t-\Delta t)}+\frac{3}{8}\Delta t\int_0^{\rm \pi}\left[\frac{Re}{2}(v_s^{(t)})^2\sin(2\theta)+\left(\frac{\partial}{\partial r}(r\omega^{(t)})-2\omega^{(t)}\right)\sin^2\theta\right]{\rm d}\theta, \end{align}
where, consistently with the first step, we discretise the time derivative by the backward Euler method. We choose this method since it is unconditionally stable, while maintaining consistency with the first step mentioned above. The two steps are applied iteratively until a given final time is reached (
$t=100$ in our simulations). This scheme was validated against Lovalenti & Brady (Reference Lovalenti and Brady1993) for the time dependent velocity of a sphere subject to a step force at non-zero 
$Re$. The steady-state problem is solved similarly: the first step remains the same, only that the time-derivative term is dropped from (3.2), while the second step is replaced by a secant method to find the value of 
$U$ that makes the hydrodynamic force vanish.
4. Results and discussion
We have performed time-dependent simulations over the range 
$0\leqslant Re\leqslant 50$ of pusher and puller quadrupolar squirmers. The squirmers are initially at rest, and begin to move in response to a time-localised dipolar perturbation. The flow field at the initial time, 
$t=0$, is that obtained by the steady-state solver in the absence of the dipolar perturbation and with fore–aft symmetry enforced; we identify that flow as a fore–aft symmetric steady base state of the quadrupolar swimmer (see figure 1a), which constitutes a continuation of the stresslet flow at 
$Re=0$ to 
$Re>0$. The dipolar perturbation is represented by the function 
$\lambda (t)$ in (2.2), which is chosen to be a Gaussian centred about 
$t=0.5$, with amplitude 
$0.1$ and standard deviation 
$0.1$. (Since the function 
$\lambda (t)$ is exceedingly small at the initial time, the incompatibility between the perturbed surface velocity and the base-state flow is negligible.)
Figure 2(a) shows the resulting time evolution of the squirmer velocity in the pusher case, for 
$Re=0,10$ and 
$50$. The swimming induced by the time-localised dipolar perturbation attenuates for all examined 
$Re$, more slowly with increasing 
$Re$. For 
$Re=0$, the velocity nearly traces the perturbation function 
$\lambda (t)$, as would be expected from Stokes-flow theory (though not precisely owing to a linear inertial effect associated with the rapid variation of the perturbation). For 
$Re=10$ and 
$Re=50$, the attenuation of the swimming speed to zero exhibits overshoot; for 
$Re=50$, the maximum swimming speed in the initial forward-motion phase is actually smaller than that in the later backwards-motion phase. The time evolution of the streamlines is presented in figure 3(a,c,e,g), for 
$Re=20$. Note the fore–aft symmetry at the initial and last times (corresponding to the base state), versus the upstream recirculation in both the forward- and backward-motion phases.
Figure 2. Dimensionless squirmer velocity from time-dependent simulations of initially stationary quadrupolar (a) pushers and (b) pullers, for the indicated 
$Re$ values. The squirmers are set into motion by means of a time-localised dipolar perturbation represented by the Gaussian function 
$\lambda (t)$ (cf. (2.2)), depicted by the dashed curves.
Figure 3. Time evolution of streamlines corresponding to time-dependent simulations as in figure 2, for 
${Re=20}$. The dipolar perturbation is maximum at 
$t=0.5$ and negligible at the other times. The streamlines at 
$t=100$ are indicative of the steady-state flow patterns. (a) Pusher, 
$t = 0$; (b) puller, 
$t = 0$; (c) pusher, 
$t = 0.5$; (d) puller, 
$t = 0.5$; (e) pusher, 
$t = 2$; ( f) puller, 
$t = 2$; (g) pusher, 
$t = 100$; (h) puller, 
$t = 100$.
We conclude that for a quadrupolar pusher, the symmetric steady base state is stable, at least up to 
$Re=50$ and under the dipolar perturbations considered.
Figure 2(b) shows the time evolution of the squirmer velocity in the puller case, for 
$Re=0,10,15,20$ and 
$25$. The response of pullers to the dipolar perturbation is seen to dramatically differ from that of pushers. Immediately following the attenuation of the perturbation, the swimming velocity attenuates as well, more slowly with increasing 
$Re$ and without reversing direction as in the pusher case. For 
$Re\leqslant 10$, this attenuation persists such that the swimming speed vanishes at long times. In contrast, for 
$Re\geq 15$, the swimming speed approaches a non-zero value, which increases with 
$Re$; for 
$Re\geq 20$, the approach to steady-state swimming is non-monotonic. The time evolution of the streamlines is presented in figure 3(b,d, f,h), for 
$Re=20$. The downstream recirculation generated by the squirmer's motion, which contrasts the upstream recirculation observed in the pusher case, is qualitatively similar to that generated by a translating no-slip sphere at moderate Reynolds number.
We conclude that the symmetric steady base state of a quadrupolar puller is unstable beyond a critical Reynolds number. Following a time-localised dipolar perturbation, the dynamics are seen to approach a symmetry-broken steady state where the squirmer exhibits self-sustained locomotion and the flow around the squirmer is fore–aft asymmetric (figure 1b). While the spontaneous locomotion is ‘forward’ in our time-dependent simulations, it is clear from the symmetry of the problem that a mirror-reflected swimming state could be excited by flipping the sign of the dipolar perturbation.
From figure 2(b), we conclude that the critical Reynolds number for a quadrupolar puller (at least under dipolar perturbations), say 
$Re_c$, lies between 
$10$ and 
$15$. However, pinning down the specific value of 
$Re_c$ with time-dependent computations is prohibitively expensive because the steady swimming speed is asymptotically small as 
$Re\searrow Re_c$. Therefore, we turn to steady-state calculations. The resulting bifurcation diagram of 
$U$ as a function of 
$Re$ is shown in figure 4. The critical Reynolds number is found to be 
$Re_c\approx 14.3$, with the steady swimming speed beyond the bifurcation monotonically growing with 
$Re$. As shown in the inset, 
$|U|\propto \sqrt {Re-Re_c}$ for 
$Re$ near 
$Re_c$, in agreement with the canonical scaling expected for a pitchfork bifurcation. The figure also depicts sample streamlines of the symmetric steady state at 
$Re=10$ and 
$Re=20$, and the symmetry-broken swimming state at 
$Re=20$.
Figure 4. Steady swimming velocity 
$U$ versus 
$Re$ for a quadrupolar-puller squirmer. Blue curves: steady-state computations employing a fore–aft asymmetric (solid) and symmetric (dashed) initial guess. Red squares: final velocity in the time-dependent simulations. The insets show the streamlines at the indicated 
$Re$ and confirm the 
$|U|\propto (Re-Re_c)^{1/2}$ behaviour near the swimming threshold, which is canonical of a pitchfork bifurcation.
It is interesting to contrast the observed scaling for 
$U$ near the bifurcation threshold with that for an isotropic autophoretic particle at 
$Re=0$ (Michelin et al. Reference Michelin, Lauga and Bartolo2013). In the autophoretic problem, an instability leading to spontaneous locomotion is observed at a sufficiently large dimensionless rate of solute emission, or intrinsic Péclet number 
$Pe$. Instead of the canonical scaling for a pitchfork bifurcation, as found herein, the particle speed is found to obey 
$|U|\propto Pe-Pe_c$, for 
$Pe$ near its critical value 
$Pe_c$ (Morozov & Michelin Reference Morozov and Michelin2019; Saha, Yariv & Schnitzer Reference Saha, Yariv and Schnitzer2021; Kailasham & Khair Reference Kailasham and Khair2022; Schnitzer Reference Schnitzer2023). This linear scaling can be traced to the fact that the base state for a spherical autophoretic particle involves no flow, hence the effective Péclet number (indicating the true ratio of advection to diffusion) near the bifurcation is actually small, such that nonlinear advection is to leading-order negligible except at large distances from the particle. Nonetheless, an analogy can be drawn with the case of a fore–aft symmetric, yet non-isotropic autophoretic particle, such as the homogeneous elliptical particles recently studied by Zhu & Zhu (Reference Zhu and Zhu2023). In such cases, the base state involves a non-trivial fore–aft symmetric flow, as herein, and indeed the bifurcation is canonical.
5. Concluding remarks
We have shown via axisymmetric numerical simulations of the Navier–Stokes equations that quadrupolar-puller squirmers, which possess axial and fore–aft symmetry, are capable of self-sustained locomotion above a moderate critical Reynolds number, 
$Re_c\approx 14.3$. Beyond that threshold, the steady swimming speed monotonically increases with 
$Re$ over the range of 
$Re$ examined, initially like 
$(Re-Re_c)^{1/2}$. Our simulations have further demonstrated that the symmetric base state, in which the squirmer is stationary, becomes unstable above the swimming threshold; in particular, when that state is disturbed by a time-localised dipolar perturbation, the squirmer relaxes towards steady swimming. These results together suggest that the spontaneous swimming emerges through a supercritical pitchfork bifurcation.
As far as we are aware, this paper is the first to suggest, let alone demonstrate, the possibility of spontaneous squirmer locomotion arising from an inertial symmetry breaking. As such, our findings give rise to many intriguing questions, which call for more extensive numerical investigations, as well as theoretical analyses.
(i) Besides quadrupolar pullers, it is clear that other, more general fore–aft symmetric squirmers are also capable of spontaneous locomotion. How does the critical Reynolds number and swimming speed depend on the combination of even modes in the modal expansion (1.1)? Computationally optimising the swimming gait with respect to some appropriate cost function, such as the swimming efficiency used at low Reynolds numbers (Lauga & Powers Reference Lauga and Powers2009), would facilitate a comparison between spontaneous and conventional swimming at moderate Reynolds numbers.
(ii) Our simulations were limited to moderate
$Re$. Up to what $Re$ can a quadrupolar puller sustain stable locomotion, with a speed monotonically increasing with $Re$? A maximum with respect to the bifurcation parameter and arrest of the swimming for values of that parameter sufficiently away from the onset of symmetry breaking were found for active droplets and particles (Michelin et al. Reference Michelin, Lauga and Bartolo2013; Izri et al. Reference Izri, Van Der Linden, Michelin and Dauchot2014), and Leidenfrost drops (Bouillant et al. Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and Quéré2018).(iii) In our simulations, the squirmer motion is collinear and the flow field is axisymmetric. For what range of
$Re$, if at all, do the spontaneous-swimming states identified here remain stable under general three-dimensional perturbations? Even if the axisymmetric swimming states are unstable, they may have stable variants featuring axial-symmetry breaking, in addition to the fore–aft-symmetry breaking. Furthermore, can some fore–aft symmetric squirmers also sustain swimming normal to their axis of symmetry?(iv) How would the spontaneous swimming of a symmetric squirmer be affected by weak asymmetric perturbations, such as owing to a constant fore–aft asymmetric swimming gait perturbation, gravity or interactions with other swimmers and with boundaries? Given the nature of spontaneous swimming, we generally expect such perturbations to significantly affect the squirmer's dynamics – not directly, but rather by slowly redirecting the motion. For example, studies of perturbed active colloids (Saha et al. Reference Saha, Yariv and Schnitzer2021; Kailasham & Khair Reference Kailasham and Khair2022; Li & Koch Reference Li and Koch2022; Schnitzer Reference Schnitzer2023; Peng & Schnitzer Reference Peng and Schnitzer2023; Zhu & Zhu Reference Zhu and Zhu2023) have demonstrated that such perturbations can introduce asymmetry between forward and backward locomotion manifested in both speed and stability, and in some cases also promote meandering or curvilinear spontaneous motion.
We conclude by noting two potential broader implications of this work. First, given the biological context of the squirmer model, the uncovered inertial symmetry breaking presents a possible mechanism for locomotion of moderate-Reynolds-number organisms. Second, our findings may lead to new design strategies for robotic swimmers that spontaneously undergo locomotion depending on the flow conditions they encounter, for example the viscosity of the surrounding fluid.
Funding
This work was supported by the National Science Foundation through grant number CBET 2002120. O.S. acknowledges the support of the Leverhulme Trust through research project grant RPG-2021-161.
Declaration of interests
The authors report no conflict of interest.
Open access
