Spontaneous locomotion of a symmetric squirmer

Abstract

The squirmer is a popular model to analyse the fluid mechanics of a self-propelled object, such as a micro-organism. We demonstrate that some fore–aft symmetric squirmers can spontaneously self-propel above a critical Reynolds number. Specifically, we numerically study the effects of inertia on spherical squirmers characterised by an axially and fore–aft symmetric ‘quadrupolar’ distribution of surface-slip velocity; under creeping-flow conditions, such squirmers generate a pure stresslet flow, the stresslet sign classifying the squirmer as either a ‘pusher’ or ‘puller’. Assuming axial symmetry, and over the examined range of the Reynolds number $Re$ (defined based upon the magnitude of the quadrupolar squirming), we find that spontaneous symmetry breaking occurs in the puller case above $Re \approx 14.3$, with steady swimming emerging from that threshold consistently with a supercritical pitchfork bifurcation and with the swimming speed growing monotonically with $Re$.

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 1981; Degen 2014; Alapan et al. 2018). Many cellular organisms exhibit some form of self-propulsion (Bray 2000; Lauga 2020), which is usually achieved by flagella or cilia acting on the surrounding fluid (Brennen & Winet 1977). The motion of microscopic organisms has been widely studied (Lauga & Powers 2009; Marchetti et al. 2013). 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. 2015; Maertens, Gao & Triantafyllou 2017). However, the swimming of mesoscale organisms, at $Re$ of order unity, is relatively unexplored (Klotsa 2019). 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 (1952) and Blake (1971), 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:

(1.1)\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 2010), locomotion in viscoelastic fluids (Zhu et al. 2011; Zhu, Lauga & Brandt 2012) and nutrient uptake (Magar, Goto & Pedley 2003; Magar & Pedley 2005). 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 1971). The second ‘quadrupolar’ (fore–aft symmetric) mode contributes a stresslet flow, whose sign distinguishes between ‘puller’ and ‘pusher’ swimmers (Ishikawa & Pedley 2007): 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 (2012) developed an asymptotic expansion through $O(Re)$ for the swimming speed of a two-mode squirmer at small $Re$, which Khair & Chisholm (2014) later extended to $O(Re^2)$. (In these works, $Re$ is defined based upon the magnitude of the dipolar squirming mode.) Chisholm et al. (2016) 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. (2016) 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 (2012) and Khair & Chisholm (2014). 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 2013; Izri et al. 2014) and Leidenfrost drops (Bouillant et al. 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 2004; Alben & Shelley 2005; Vandenberghe, Childress & Zhang 2006).

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

(2.1a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$

(2.1b)$$\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

(2.2)\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,

(2.3)\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

(2.4)\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$:

(3.1)\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

(3.2)\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

(3.3)\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

(3.4a)$$\begin{gather} \psi=0\quad \text{and}\quad \frac{\partial\psi}{\partial r} ={-}v_s\sin\theta \quad \text{at }r=1, \end{gather}$$

(3.4b)$$\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

(3.5)\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 (2014). 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:

(3.6)\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. (2016), which employs the Galerkin method of weighted residuals (Karniadakis & Sherwin 2005). 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 2009), 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 (2022) and Cobos & Khair (2023).

In the second step, we update the squirmer speed according to (2.4),

(3.7)\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 (1993) 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. 2013). 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 2019; Saha, Yariv & Schnitzer 2021; Kailasham & Khair 2022; Schnitzer 2023). 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 (2023). 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.

1. (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 2009), would facilitate a comparison between spontaneous and conventional swimming at moderate Reynolds numbers.

2. (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. 2013; Izri et al. 2014), and Leidenfrost drops (Bouillant et al. 2018).

3. (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?

4. (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. 2021; Kailasham & Khair 2022; Li & Koch 2022; Schnitzer 2023; Peng & Schnitzer 2023; Zhu & Zhu 2023) 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.