1. Introduction
A small motile organism swimming in a marine environment experiences a hydrodynamic force. How does this force depend on the mechanism of propulsion, upon the swimming gait of the organism? How does the shape of the swimmer affect the hydrodynamic force, and how does it depend on the centre-of-mass velocity, and upon the angular velocity of the swimmer? Closely related questions concern the disturbance produced by the organism as it moves through the fluid. How does the disturbance decay away from the swimmer, how does it change as a function of time for a time-dependent swimming gait, and how do sudden accelerations affect the disturbance? These are important questions, because the disturbance flow determines hydrodynamic interactions between organisms, and it is known that motile microorganisms employ hydrodynamic signals to localise prey, or to escape their predators, or to find mates (Guasto, Rusconi & Stocker Reference Guasto, Rusconi and Stocker2012).
Answering these questions in general is obviously a very challenging task. To simplify the problem, the dynamics of motile microorganisms in water is often described using the creeping-flow approximation – the Stokes approximation – which neglects possible effects of unsteady and convective fluid inertia (Lighthill Reference Lighthill1952; Happel & Brenner Reference Happel and Brenner1965; Blake Reference Blake1971; Yates Reference Yates1986; Fauci & Dillon Reference Fauci and Dillon2006; Lauga & Powers Reference Lauga and Powers2009; Visser Reference Visser2011; Pedley Reference Pedley2016). This approach works very well for the majority of microorganisms. But for larger ones, in particular for organisms that swim with highly unsteady time-dependent gaits, fluid inertia cannot be neglected. The flow field of the organism Chlamydomonas reinhardtii, for instance, cannot be described by the Stokes approximation (Wei et al. Reference Wei, Dehnavi, Aubin-Tam and Tam2021). As a second example, consider Mesodinium rubrum, which swims with short jumps interrupted by longer rest periods (Fenchel & Hansen Reference Fenchel and Hansen2006; Jiang Reference Jiang2011). When the organism stops suddenly, fluid inertia cannot be neglected: the Stokes approximation predicts that the surrounding fluid arrests instantly. In reality, the disturbance flow around the organism takes some time to vanish. For a passive sphere, this inertia effect induced by the unsteadiness of the disturbance flow gives rise to a memory term in the expression for the hydrodynamic force that depends on its past accelerations, the Boussinesq–Basset–Oseen (BBO) history force (Boussinesq Reference Boussinesq1885; Basset Reference Basset1888; Oseen Reference Oseen1927). Wang & Ardekani (Reference Wang and Ardekani2012b) used the unsteady Stokes equation to model how the velocity of a spherical squirmer decays after a sudden jump to escape a predator. They determined how the history force affects the velocity decay, and found good agreement with measurements performed by Jiang & Kiorboe (Reference Jiang and Kiorboe2011) on copepods, despite the fact that possible effects of convective fluid inertia were not considered. Spelman & Lauga (Reference Spelman and Lauga2017) calculated the hydrodynamic force and the disturbance flow produced by a deforming swimmer, accounting for unsteady fluid inertia.
At least for passive particles, it is known that convective fluid inertia, induced by a non-zero slip velocity between the organism and the surrounding fluid, can change the long-time behaviour of the history force. For times larger than the Oseen time, the kernel of the history force decays more rapidly than the BBO kernel, namely as 
$t^{-2}$ for a sudden start (Sano Reference Sano1981; Lovalenti & Brady Reference Lovalenti and Brady1993) instead of the 
$t^{-1/2}$ decay of the BBO history force. Sano (Reference Sano1981) and Lovalenti & Brady (Reference Lovalenti and Brady1993) derived their results for weak convective fluid inertia, for small but non-zero particle Reynolds numbers 
$Re_{{p}}$. In the experiments carried out by Jiang & Kiorboe (Reference Jiang and Kiorboe2011), 
$Re_{{p}}$ was quite high, however. The small-
$Re_{{p}}$ theory is expected to fail in this case.
It is not known how the results of Sano (Reference Sano1981) and Lovalenti & Brady (Reference Lovalenti and Brady1993) generalise to motile organisms, yet fluid-inertia forces are thought to play a significant role for swimmers in many situations. Wang & Ardekani (Reference Wang and Ardekani2012a), for instance, determined convective-inertia corrections for a steady, spherical squirmer, modelling the swimming gait as a steady surface velocity (see also Khair & Chisholm Reference Khair and Chisholm2014; Chisholm et al. Reference Chisholm, Legendre, Lauga and Khair2016). However, the slip velocity of motile organisms is usually time-dependent, and the quasi-steady approximation may fail if the slip velocity varies much faster than the disturbance flow.
In order to understand how convective inertia modifies the unsteady dynamics of a small motile organism, we calculated the hydrodynamic force on a small spherical squirmer in an unsteady, spatially homogeneous flow, thus generalising the results of Lovalenti & Brady (Reference Lovalenti and Brady1993) to a small motile organism. While Lovalenti & Brady (Reference Lovalenti and Brady1993) used the reciprocal theorem to obtain their results, we employed asymptotic matching of perturbation expansions (Hinch Reference Hinch1995) in the parameter 
$\varepsilon = \sqrt {{Re}_{{p}}\,{Sl}}$, where 
${Sl}$ is the Strouhal number. We note that Sennitskii (Reference Sennitskii1990) computed the disturbance velocity for a swimmer, far from its surface, using asymptotic matching of perturbation expansions (in a different parameter), but not the hydrodynamic force. Our solution is related to that of Lovalenti & Brady (Reference Lovalenti and Brady1993). In particular, we find that the hydrodynamic force on the squirmer involves memory forces, and that their kernels are closely related to those of Lovalenti & Brady (Reference Lovalenti and Brady1993). For all examples that we checked, the kernels were in fact numerically identical.
We identified two major differences to the passive case. First, for the squirmer, the memory forces involve a source term that contains an active part. Second, the inhomogeneous part of the solution to the inner problem of order 
$\varepsilon$ contributes to the hydrodynamic force on the squirmer. For a passive spherical particle, by contrast, spherical symmetry ensures that such contributions vanish. We discuss the significance of the inhomogeneous contribution for the squirmer. Last, but not least, asymptotic matching allows us to determine the disturbance flow. We show how the disturbance develops for a sudden start and for a sudden stop, and discuss the implications of our findings for the biology of small motile organisms in a marine environment.
2. Formulation of the problem
The spherical squirmer is an idealised model for a microswimmer, introduced by Lighthill (Reference Lighthill1952) and Blake (Reference Blake1971). It is used widely to investigate the dynamics of motile microorganisms, their interactions, and collective behaviours (Lauga & Powers Reference Lauga and Powers2009; Visser Reference Visser2011; Pak & Lauga Reference Pak and Lauga2014; Pedley Reference Pedley2016). In this model, the swimming gait is represented by an axisymmetric surface-velocity field in the frame translating with the squirmer:
\begin{equation} \boldsymbol{v} (\theta,t) = \sum_l A_l(t)\,P_l(\cos\theta)\,\hat{\boldsymbol{e}}_r+ B_l(t)\,{V_l}(\cos\theta)\, \hat{\boldsymbol{e}}_\theta. \end{equation}
Here, 
$\theta$ parametrises points on the surface, and 
$\hat {\boldsymbol {e}}_r$ and 
$\hat {\boldsymbol {e}}_\theta$ are the corresponding basis vectors in the body frame (figure 1). Using the notation of Lighthill (Reference Lighthill1952), 
$P_l$ is the Legendre polynomial of order 
$l$, and 
$V_l$ is related to the first associated Legendre polynomial of order 
$l$, namely 
$V_l = -2 P_l^1/[l (l+1)]$. In general, the coefficients 
$A_l(t)$ and 
$B_l(t)$ are allowed to depend on time, but the simplest version of the model is the steady, spherical squirmer with a time-independent, tangential surface-velocity field (Pedley Reference Pedley2016)
\begin{equation} \boldsymbol{v}(\theta)= \left(B_1\sin \theta + B_2\sin \theta \cos \theta \right) \hat{\boldsymbol{e}}_\theta. \end{equation}When inertial effects are negligible, one solves the steady Stokes equation with the boundary condition (2.2) to determine the hydrodynamic force acting on such a spherical squirmer. The result is
\begin{equation} \boldsymbol{f}'^{(0)} ={-}6 {\rm \pi}\mu a \left( \dot{\boldsymbol{x}} -\tfrac{2}{3}B_1 \boldsymbol n\right), \end{equation}
where 
$\mu$ is the dynamic viscosity of the fluid, 
$\boldsymbol n$ is the unit vector along the symmetry axis of the surface-velocity field (figure 1), 
$a$ is the radius of the spherical squirmer, and 
$\boldsymbol {x}$ is the position vector of its centre of mass in the laboratory frame. The steady swimming velocity in the Stokes approximation is obtained by setting the hydrodynamic force to zero. This gives 
$\dot {\boldsymbol {x}}= (2/3) B_1 \boldsymbol {n}$.
Schematic illustration of the squirmer model described in § 2. The squirmer swims head first, along the direction 
$\boldsymbol {n}$ in the 
$\hat {\boldsymbol {e}}_1$–
$\hat {\boldsymbol {e}}_2$ plane. Points on the surface in this plane are parametrised by the angle 
$\theta$. We show the lab-coordinate system 
$\hat {\boldsymbol {e}}_i$, and the coordinate system 
$\hat {\boldsymbol {e}}_r$ and 
$\hat {\boldsymbol {e}}_\theta$ in the body frame.
We consider an unsteady squirmer, with time-dependent coefficients 
$B_1(t)$ and 
$B_2(t)$ in a quiescent fluid. The surface-velocity field produces a time-dependent disturbance flow that causes the organism to move with a time-dependent centre-of-mass velocity 
$\dot {\boldsymbol {x}}(t)$. In order to determine the hydrodynamic force on the squirmer, one has to solve the continuity and momentum equations for the disturbance flow 
$\boldsymbol {w}$ and the disturbance pressure 
$p$. In non-dimensional variables, these equations take the standard form
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{w} = 0, \end{gather}$$
$$\begin{gather}\left.{Re}_{{p}}\,{Sl}\,\frac{\partial \boldsymbol{w}}{\partial t} \right|_{\boldsymbol{r}} -{Re}_{{p}}\,\dot{\boldsymbol{x}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{w}+{Re}_{{p}} \, \boldsymbol{w} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{w} ={-} \boldsymbol{\nabla} p + \Delta \boldsymbol{w} . \end{gather}$$When the angular velocity of the swimmer is negligible, the boundary conditions for (2.4b) read
\begin{equation} \boldsymbol{w} = \dot{\boldsymbol{x}} + \boldsymbol{v}\ \mbox{for}\ |\boldsymbol r|=1 \quad \mbox{and} \quad \boldsymbol{w} \to \boldsymbol{0}\ \mbox{as} |\boldsymbol r|\to \infty. \end{equation}
All vectors in (2.4) are expressed in the laboratory frame, but using a system of moving coordinates that translates with the squirmer. In particular, (2.4b) is obtained using the following relation that links the partial time derivative at a fixed point 
$\boldsymbol {x}$ in the laboratory frame to the partial time derivative at a fixed point 
$\boldsymbol {r}$ in the moving coordinate system: 
$(\partial /\partial t)|_{\boldsymbol {x}} = (\partial /\partial t)|_{\boldsymbol {r}}-\dot {\boldsymbol {x}} \boldsymbol {\cdot } \boldsymbol {\nabla }$. Here, 
$\dot {\boldsymbol {x}}$ is the centre-of-mass velocity of the squirmer. Furthermore, we non-dimensionalised lengths by dividing with the radius 
$a$ of the spherical squirmer, velocities with a typical velocity 
$u_{c}$, and time with a typical time scale 
$\tau _{c}$. The resulting non-dimensional parameters are the particle Reynolds number and the Strouhal number:
\begin{equation} {Re}_{{p}} = \frac{a u_{c}}{\nu} \quad \mbox{and} \quad {Sl} = \frac{a}{u_{c} \tau_{c}}\:. \end{equation}
The particle Reynolds number determines the importance of convective terms in (2.4b), and 
${Re}_{{p}}\,{Sl}=a^2/(\nu \tau _{c})$ is the ratio of two time scales, the viscous time 
$\tau _\nu =a^2/\nu$, and 
$\tau _{c}$. Equations (2.4a) and (2.4b) are precisely those used by Lovalenti & Brady (Reference Lovalenti and Brady1993) to study the effect of convective inertia on the hydrodynamic force on a passive sphere in a spatially uniform flow, their (2.3a), (2.3b). Only the boundary conditions (2.4c) differ, by the surface-velocity field 
$\boldsymbol {v}$. This additional term makes the difference between an active and a passive particle. In the following, we describe how this term modifies the hydrodynamic force.
3. Matched asymptotic expansions
We solve (2.4) under the assumption that
\begin{equation} {Re}_{p} \ll1 \quad \mbox{and} \quad {Re}_{p}{Sl} \ll 1. \end{equation}
In this limit, the disturbance near the spherical squirmer is well approximated by the quasi-steady Stokes solution, which decays as 
$|\boldsymbol {r}|^{-1}$ for large values of 
$|\boldsymbol {r}|$. As a consequence, the unsteady term 
${Re}_{{p}}\,{Sl}\,\partial _t \boldsymbol {w}$ in (2.4b) becomes of the same order of magnitude as the viscous term 
$\Delta \boldsymbol {w}$, at a distance 
$|\boldsymbol {r}| \sim \ell _p\equiv a /\sqrt {{Re}_{{p}}\,{Sl}}$ called the ‘penetration length’. The convective term 
${Re}_{{p}}$ 
$\dot {\boldsymbol {x}} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {w}$ in (2.4b) is of the same order of magnitude as the viscous term at the Oseen length, 
$|\boldsymbol {r}| \sim \ell _{O} \equiv a /{Re}_{{p}}$.
The problem has two asymptotic limits. In the limit 
${Sl}\to 0$, the problem becomes steady, leading to an Oseen correction to the hydrodynamic force (Wang & Ardekani Reference Wang and Ardekani2012a; Khair & Chisholm Reference Khair and Chisholm2014; Chisholm et al. Reference Chisholm, Legendre, Lauga and Khair2016). Conversely, when 
${Re}_{{p}}\to 0$ and 
${Sl}\to \infty$, keeping 
$\varepsilon ^2\equiv Re_{{p}}\,{Sl}$ constant, one obtains an unsteady Stokes problem. Its solution yields the BBO history force obtained by Wang & Ardekani (Reference Wang and Ardekani2012b).
Since inertial corrections are singular perturbations, we use asymptotic matching of perturbation expansions (Hinch Reference Hinch1995) in the parameter 
$\varepsilon =\sqrt {{Re}_{{p}}\,{Sl}}$ to compute the inertial corrections to the hydrodynamic force on the squirmer. It is natural to take the independent non-dimensional parameters as
\begin{equation} \varepsilon = \sqrt{{Re}_{{p}}\,{Sl}} \quad \mbox{and}\quad {\ell_p}/{\ell_{O}} =\sqrt{ {Re}_{{p}}/{Sl}}. \end{equation}
We require 
$\varepsilon$ to be small – in keeping with (3.1a,b) – and 
${\ell _p}/{\ell _{O}}$ to remain finite as 
$\varepsilon$ becomes small. The length scale ratio 
${\ell _p}/{\ell _{O}}$ characterises the competition between convective and unsteady inertia upon the disturbance flow. At first sight, it might appear that one cannot treat large unsteadiness in this fashion, but we show below that the hydrodynamic force obtained in this way is uniformly valid. In terms of the parameters (3.1a,b), the equations of motion (2.4) take the form
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{w}}= \boldsymbol{0}, \end{gather}$$
$$\begin{gather}\left.\varepsilon^2\,\frac{\partial \boldsymbol{w}}{\partial t} \right|_{\boldsymbol{r}} - \varepsilon \left(\frac{\ell_p}{\ell_{O}} \right)\dot{\boldsymbol{x}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{w}+\varepsilon \left(\frac{\ell_p}{\ell_{O}} \right) \boldsymbol{w} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{w} ={-} \boldsymbol{\nabla} p + \Delta \boldsymbol{w} . \end{gather}$$
To find the disturbance flow 
$\boldsymbol {w}$, configuration space is divided in two different regions: the inner region, 
$|\boldsymbol {r}| \sim 1$, and the outer region, 
$|\boldsymbol {r}| \gg 1$. In these two regions, disturbance flow and pressure are written in the form of two different series expansions in 
$\varepsilon$. The inner and outer expansions are matched at 
$|\boldsymbol {r}| \sim 1/\varepsilon$. This yields the necessary boundary conditions for the inner problem that can then be solved to determine the hydrodynamic force (Hinch Reference Hinch1995).
Before going further, let us comment briefly on the different terms on the left-hand side of (3.3b). Near the squirmer, all terms are small when 
$\varepsilon$ is small. But while the unsteady term scales as 
$\varepsilon ^2$, the convective terms are proportional to 
$\varepsilon$. These terms cause the key difficulty. When the perturbation is of order 
$\varepsilon ^n$ with 
$n>1$, the hydrodynamic force is given by the outer solution alone (Legendre & Magnaudet Reference Legendre and Magnaudet1998; Redaelli et al. Reference Redaelli, Candelier, Mehaddi and Mehlig2022). In the present case, the perturbation is larger – of order 
$\varepsilon$ – therefore the method of Redaelli et al. (Reference Redaelli, Candelier, Mehaddi and Mehlig2022) cannot be applied here. As a consequence, we need to consider the details of the inner solution. We discuss the outer solution first, however, because it yields the boundary conditions needed for solving the inner problem.
3.1. Outer solution
Far from the squirmer, the boundary conditions on the surface of the organism can be replaced by a Dirac delta function with amplitude 
$\boldsymbol {f}^{(0)}$ (Childress Reference Childress1964):
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{w}}_{{out}}= \boldsymbol{0}, \end{gather}$$
$$\begin{gather}\left.\varepsilon^2\,\frac{\partial \boldsymbol{w}_{{out}}}{ \partial t} \right|_{\boldsymbol{r}} -\varepsilon \left(\frac{\ell_p}{\ell_{O}}\right)\dot{\boldsymbol{x}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{w}}_{{out}} ={-} \boldsymbol{\nabla} {p}_{{out}} + \Delta{\boldsymbol{w}}_{{out}} + \boldsymbol{f}^{(0)}\,\delta(\boldsymbol{r}) . \end{gather}$$
In the matching region, the quadratic term 
$\varepsilon (\ell _p/\ell _{O}) \, \boldsymbol {w} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {w}$ in (3.3) is negligible compared to the other terms in the equation, because its magnitude scales as 
$\varepsilon (\ell _p/\ell _{O})\,|\boldsymbol {r}|^{-3}\sim \varepsilon ^4$. The amplitude of the source term is the opposite of the Stokes force (2.3):
\begin{equation} \boldsymbol{f}^{(0)}={-}{6{\rm \pi}}\left({\tfrac{2}{3} B_1 \boldsymbol{n}} - \dot{\boldsymbol{x}}\right). \end{equation}
The first term on the right-hand side is an active part, related to the surface-velocity field 
$\boldsymbol {v}$.
The solution of the outer problem (3.4) is derived in Appendix A. It reads
\begin{align} {\boldsymbol{w}}_{{out}} &= \underbrace{\frac{1}{8{\rm \pi}}\left(\frac{\unicode{x1D7D9}}{r} + \frac{\boldsymbol{r}\otimes \boldsymbol{r}}{r^3}\right) \boldsymbol{\cdot} \boldsymbol{f}^{(0)}}_{\equiv\boldsymbol{\mathcal{T}}_{reg}^{(0)}(\boldsymbol{r},t)} \underbrace{{}-{Re}_{{p}}\,\dot{\boldsymbol{x}}\boldsymbol{\cdot} \boldsymbol{\nabla}\left[\frac{3 r}{32{\rm \pi}} \left({\unicode{x1D7D9}} - \frac{1}{3}\,\frac{\boldsymbol{r} \otimes \boldsymbol{r}}{r^2}\right) \boldsymbol{\cdot} \boldsymbol{f}^{(0)} \right]}_{\equiv\varepsilon \,\boldsymbol{\mathcal{T}}_{reg}^{(1)}(\boldsymbol{r},t)}\nonumber\\ &\quad \underbrace{{}- \varepsilon \int_{0}^{t}\frac{\mbox{d}\tau}{6{\rm \pi}}\, {{\mathbb{K}}^{(1)}(t,\tau)} \boldsymbol{\cdot} \frac{{\rm d}}{{ \rm d} \tau} \boldsymbol{f}^{(0)}(\tau) - {Re}_{{p}} \, \int_{0}^{t}\frac{\mbox{d}\tau}{6{\rm \pi}} \,{\mathbb{K}}^{(2)}(t,\tau) \boldsymbol{\cdot} \boldsymbol{f}^{(0)}(\tau)}_{{{}\equiv\varepsilon \boldsymbol{\mathcal{U}}(t)}}. \end{align}
The integral kernels 
${\mathbb{K}}_1$ and 
${\mathbb{K}}_2$ in (3.5a) have elements
\begin{align} {[}{\mathbb{K}}^{(1)}(t,\tau)]_{ij}&={-} \frac{3}{8}\,\Bigg[({-}2+A(t,\tau)^{{-}2}) \frac{\mbox{erf}(A(t,\tau))}{A(t,\tau)}- 2\,\frac{\exp({-}A(t,\tau)^2)}{A(t,\tau)^2\sqrt{\rm \pi}} \Bigg] \frac{\delta_{ij}}{\sqrt{t-\tau}} \nonumber\\ &\quad - \frac{3}{8}\, \Bigg[\!\left(1-\frac{3}{2}\,A(t,\tau)^{{-}2}\right) \frac{\mbox{erf}(A(t,\tau))}{A(t,\tau)}+ 3\,\frac{\exp({-}A(t,\tau)^2)}{A(t,\tau)^2\sqrt{\rm \pi}} \Bigg] \nonumber\\ &\quad \times\frac{\delta_{ij} - q_i(t,\tau)\,q_j(t,\tau)}{\sqrt{t-\tau}}, \end{align}
\begin{align} {[}{\mathbb{K}}^{(2)}(t,\tau)]_{ij} &={-}\frac{3 }{8}\,\frac{1}{2\,A(t,\tau)}\left[ \left(1-\frac{3}{2}\,A(t,\tau)^{{-}2}\right) \frac{\mbox{erf}(A(t,\tau))}{A(t,\tau)}+ 3\,\frac{\exp({-}A(t,\tau)^2)}{A(t,\tau)^2\sqrt{\rm \pi}} \right] \nonumber\\ &\quad \times \frac{q_i(t,\tau)\,{\dot{x}}_j(\tau) - 3\,q_i(t,\tau)\,q_j(t,\tau) {\sum_{k=1}^3}q_k(t,\tau)\,{\dot{x}}_{k}(\tau) }{t-\tau}\nonumber\\ &\quad - \frac{3 }{8} \left( \frac{1}{2\,A(t,\tau)}\right)\left[ 3\,\frac{\mbox{erf}(A(t,\tau))}{A(t,\tau)^3}-(4+6\, A(t,\tau)^{{-}2})\frac{\exp({-}A(t,\tau)^2)}{\sqrt{\rm \pi}} \right]\nonumber\\ &\quad \times \frac{ [\delta_{ij} - q_i(t,\tau)\, q_j(t,\tau)]{\sum_{k=1}^3}q_k(t,\tau)\,{\dot{x}}_{k}(\tau)}{t-\tau}. \end{align}
Here, 
$A(t,\tau )$ is the norm of the pseudo-displacement vector 
$\boldsymbol {a}(t,\tau )$ introduced by Lovalenti & Brady (Reference Lovalenti and Brady1993) in (6.7d) of their paper:
\begin{equation} \boldsymbol{a} (t,\tau) \equiv \left(\frac{\ell_p}{\ell_{O}}\right) \frac{1}{2 \sqrt{t-\tau}} \int_\tau^t \,\mbox{d}t' \,\dot{\boldsymbol{x}}(t'), \end{equation}
and 
$q_i$ are the components of the (unit) direction vector 
$\hat {\boldsymbol {q}}= \boldsymbol {a}(t,\tau )/A(t,\tau )$. We note that the outer solution derived here differs from that obtained by Shu & Chwang (Reference Shu and Chwang2001); they considered a Dirac delta function force, instead of a continuously time-dependent force.
Let us discuss briefly the different terms in the outer solution (3.5a). Following Meibohm et al. (Reference Meibohm, Candelier, Rosén, Einarsson, Lundell and Mehlig2016), we denote the first term on the right-hand side of (3.5a) by 
$\boldsymbol {\mathcal {T}}_{reg}^{(0)}(\boldsymbol {r},t)$. This is the well-known Stokeslet describing the far-field disturbance flow produced by an active particle in the Stokes limit. The second term on the right-hand side of (3.5a) is denoted by 
$\varepsilon \,\boldsymbol {\mathcal {T}}_{reg}^{(1)}(\boldsymbol {r},t)$. This term represents a regular perturbation to the Stokeslet 
$\boldsymbol {\mathcal {T}}_{reg}^{(0)}(\boldsymbol {r},t)$. The two integral terms, finally, combine to give the history force. It can be written in the form 
$\varepsilon \,\boldsymbol {\mathcal {U}}(t)$, where 
$\boldsymbol {\mathcal {U}}(t)$ is a spatially uniform flow at infinity. A passive sphere has 
$B_1=B_2=0$. In this case, spherical symmetry ensures that all first-order inertia corrections to the hydrodynamic force are due to the outer flow alone, and one obtains 
$\boldsymbol {f}'= \boldsymbol {f}'^{(0)} + 6{\rm \pi} \varepsilon \,\boldsymbol {\mathcal {U}}(t)$ to order 
$\varepsilon$. This expression for the hydrodynamic force is equivalent to that obtained by Lovalenti & Brady (Reference Lovalenti and Brady1993) (their (6.15)) with three minor differences. First, Lovalenti & Brady (Reference Lovalenti and Brady1993) combined the terms multiplying 
$\boldsymbol f^{(0)}$ and the time derivative of 
$\boldsymbol f^{(0)}$. We kept them separate here, because they have different behaviours at short and long times (§ 4). Second, we assumed that 
$\boldsymbol {f}^{(0)}=\boldsymbol 0$ for 
$t \le 0$. As a consequence, the integration domain is 
$[0,t]$ instead of 
$[-\infty,t]$. Third, our expression for the kernel looks slightly different from the kernel in (6.15) of Lovalenti & Brady (Reference Lovalenti and Brady1993). However, we found that the two expressions are numerically equivalent for all centre-of-mass motions that we examined: sudden start (Sano Reference Sano1981), linearly increasing velocity 
$\dot {\boldsymbol {x}}=\boldsymbol v_0 t/t_0$, and sinusoidally varying centre-of-mass velocity 
$\dot {\boldsymbol {x}} =\boldsymbol v_0 \sin (\omega _0 t)$, with coefficients 
$\boldsymbol v_0$, 
$t_0$ and 
$\omega _0$. It is quite common that the two methods, reciprocal theorem and asymptotic matching, yield expressions for the hydrodynamic force and torque that look different but are equivalent (Meibohm et al. Reference Meibohm, Candelier, Rosén, Einarsson, Lundell and Mehlig2016).
3.2. Inner solution
Near the squirmer, for 
$| \boldsymbol {r}| = {O}(1)$, the disturbance flow is expanded as
\begin{equation} {\boldsymbol{w}}_{{in}} = {\boldsymbol{w}}_{{in}}^{(0)} + \varepsilon \, {\boldsymbol{w}}_{{in}}^{(1)} + O (\varepsilon^2) \quad \mbox{and} \quad {p}_{{in}} = {p}_{{in}}^{(0)} + \varepsilon \, {p}_{{in}}^{(1)} + O (\varepsilon^2) . \end{equation}
These inner expansions must be matched, order by order, to the outer expansion (3.5a). At order 
$\varepsilon ^0$, the inner problem to solve is
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{w}}_{{in}}^{(0)}= \boldsymbol{0},\quad - \boldsymbol{\nabla} {p}_{{in}}^{(0)} + \Delta {\boldsymbol{w}}_{{in}}^{(0)} = \boldsymbol{0} , \end{gather}$$
$$\begin{gather}{\boldsymbol{w}}_{{in}}^{(0)} = \dot{\boldsymbol{x}} (t) + \boldsymbol{v}(t) \ \mbox{for}\ {|\boldsymbol r|=1} \quad \mbox{and} \quad {\boldsymbol{w}}_{{in}}^{(0)} \sim \boldsymbol{\mathcal{T}}_{reg}^{(0)} \ \mbox{for}\ |\boldsymbol{r}| \to \infty. \end{gather}$$
This is a homogeneous Stokes problem whose solution is well known (Blake Reference Blake1971). At order 
$\varepsilon$, an inhomogeneous Stokes problem must be solved:
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{w}}_{{in}}^{(1)}= 0,\quad - \boldsymbol{\nabla} {p}_{{in}}^{(1)} + \Delta {\boldsymbol{w}}_{{in}}^{(1)} ={-} \left(\frac{\ell_p}{\ell_{O}}\right) \dot{\boldsymbol{x}}\boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{w}}_{{in}}^{(0)} + \left(\frac{\ell_p}{\ell_{O}}\right) {\boldsymbol{w}}_{{in}}^{(0)}\boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{w}}_{{in}}^{(0)}, \end{gather}$$
$$\begin{gather}{\boldsymbol{w}}_{{in}}^{(1)} = \boldsymbol{0} \ \mbox{for}\ {|\boldsymbol{r}|=1}\quad \mbox{and} \quad \quad {\boldsymbol{w}}_{{in}}^{(1)} \sim \boldsymbol{\mathcal{T}}_{reg}^{(1)} + \boldsymbol{\mathcal{U}}(t) \ \mbox{for}\ |\boldsymbol r|\to \infty. \end{gather}$$Since the problem is inhomogeneous, we seek its solution in the form of a sum of a particular solution and the homogeneous solution:
\begin{equation} {\boldsymbol{w}}_{{in}}^{(1)} = [ {\boldsymbol{w}}_{{p}}^{(1)} + \boldsymbol{\mathcal{U}}(t) ] + {\boldsymbol{w}}_{{h}}^{(1)} \quad \mbox{and} \quad p_{in}^{(1)} = p_{p}^{(1)}+p_{h}^{(1)}. \end{equation}
Consider first the particular solution. The velocity 
${\boldsymbol {w}}_{{p}}^{(1)}$ and the pressure 
$p_{p}^{(1)}$ must satisfy
\begin{align} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{w}}_{{p}}^{(1)}= 0,\quad - \boldsymbol{\nabla} {p}_{{p}}^{(1)} + \Delta {\boldsymbol{w}}_{{p}}^{(1)} ={-} \left(\frac{\ell_p}{\ell_{O}}\right) \dot{\boldsymbol{x}}\boldsymbol{\cdot}\boldsymbol{\nabla} {\boldsymbol{w}}_{{in}}^{(0)} + \left(\frac{\ell_p}{\ell_{O}}\right) {\boldsymbol{w}}_{{in}}^{(0)} \boldsymbol{\cdot}\boldsymbol{\nabla} {\boldsymbol{w}}_{{in}}^{(0)}. \end{align}
Since we added the uniform term 
$\boldsymbol {\mathcal {U}}(t)$ to 
${\boldsymbol {w}}_{{p}}^{(1)}$, the boundary condition at infinity is 
${\boldsymbol {w}}_{{p}}^{(1)} \sim \boldsymbol {\mathcal {T}}_{reg}^{(1)}$ as 
$|\boldsymbol {r}|\to \infty$. The homogeneous part of the full solution satisfies
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{w}}_{{h}}^{(1)}= 0,\quad - \boldsymbol{\nabla} {p}_{{h}}^{(1)} + \Delta {\boldsymbol{w}}_{{h}}^{(1)} = \boldsymbol{0}, \end{gather}$$
$$\begin{gather}{\boldsymbol{w}}_{{h}}^{(1)} = \left.- {\boldsymbol{w}}_{{p}}^{(1)}\right|_{|\boldsymbol{r}|=1} - \boldsymbol{\mathcal{U}}(t) , \quad |\boldsymbol{r}|=1 \quad \mbox{and} \quad {\boldsymbol{w}}_{{in}}^{(1)} \sim \boldsymbol{0}\ \text{as}\ r \to \infty. \end{gather}$$The solutions of (3.7) and (3.8) obtained in this way are quite lengthy; we do not reproduce them here. The inhomogeneous equation (3.9a,b) is solved using a Fourier transform (Candelier et al. Reference Candelier, Mehaddi, Mehlig and Magnaudet2023). Equation (3.11a) is a Stokes problem, solved by Lamb's general solution (Happel & Brenner Reference Happel and Brenner1965). Integrating the corresponding stress tensor over the surface of the swimmer yields the hydrodynamic force. For a passive spherical particle, the contribution of the inhomogeneous part of the inner solution to the hydrodynamic force must vanish due to spherical symmetry. In this case, the inertial corrections to the hydrodynamic force are determined by the outer solution alone. For the swimmer, the inhomogeneous solution contributes to the hydrodynamic force. This explains why the method of Redaelli et al. (Reference Redaelli, Candelier, Mehaddi and Mehlig2022) fails to determine the entire corrections to the hydrodynamic force due to convective-inertia effects.
4. Results
4.1. Hydrodynamic force
For the unsteady squirmer, the calculations outlined in the previous section yield the hydrodynamic force
\begin{align} \boldsymbol f' &={-} 6 {\rm \pi}\left[\dot{\boldsymbol{x}} -\frac{2}{3}\,B_1(t)\,\boldsymbol{n}(t) \right]\nonumber\\ &\quad - \varepsilon \, 6{\rm \pi} \int_{0}^{t} \mbox{d}\tau \, {\mathbb{K}}^{(1)}(t,\tau) \boldsymbol{\cdot} \frac{\mbox{d}}{\mbox{d}\tau} \left[ \dot{\boldsymbol{x}}(\tau) -\frac{2}{3}\,B_1(\tau)\,\boldsymbol{n}(\tau)\right]\nonumber\\ &\quad - {{Re}_{{p}} } \,6{\rm \pi} \int_{0}^{t} {\mbox{d} \tau\,{\mathbb{K}}^{(2)}(t,\tau) } \boldsymbol{\cdot} \left[\dot{\boldsymbol{x}}(\tau) -{\frac{2}{3}}\,B_1(\tau)\,\boldsymbol n(\tau)\right]\nonumber\\ &\quad - {{Re}_{{p}} } \left[ \frac{2 {\rm \pi}}{5}\,\dot{x}(t)\,B_2(t) + \frac{2{\rm \pi}}{15}\,B_1(t)\,B_2(t) \right]\boldsymbol{n}(t) . \end{align}This is our main result, the hydrodynamic force on an unsteady, spherical squirmer in a time-dependent homogeneous flow at small but non-zero particle Reynolds number. Equation (4.1) generalises (6.15) of Lovalenti & Brady (Reference Lovalenti and Brady1993) from a passive to an active sphere. As discussed already, the numerical calculations described above indicate that the kernels in (4.1) combine to those of Lovalenti & Brady (Reference Lovalenti and Brady1993).
Equation (4.1) simplifies to known solutions when unsteadiness is either very weak or very strong. To see this, consider the asymptotic behaviours of the integrals
\begin{equation} {\boldsymbol I_1={-} \varepsilon \int_{0}^{t}\mbox{d}\tau \, {\mathbb{K}}^{(1)}(t,\tau) \boldsymbol{\cdot} \frac{{{ \rm d}}}{{ \rm d} \tau}}\boldsymbol{f}^{(0)}(\tau) {\quad\mbox{and}\quad \boldsymbol I_2={-} {Re}_{{p}} \, \int_{0}^{t} {\rm d}\tau\,{\mathbb{K}}^{(2)}(t,\tau) \boldsymbol{\cdot} \boldsymbol{f}^{(0)}(\tau)} . \end{equation}
In the steady limit, 
$\ell _p / \ell _{O} \to \infty$, we find
\begin{equation} \boldsymbol{I}_1 \to \boldsymbol{0} \quad \mbox{and} \quad \boldsymbol{I}_2 \to - \tfrac{3}{8} {Re_{{p}}}\,|\dot{\boldsymbol{x}}|\,\boldsymbol{f}^{(0)}. \end{equation}Only the convective Oseen correction survives. It combines with the contribution of the inhomogeneous inner solution to give the steady hydrodynamic force obtained by Wang & Ardekani (Reference Wang and Ardekani2012a) and Khair & Chisholm (Reference Khair and Chisholm2014).
Conversely, when unsteadiness is high, 
$\ell _p / \ell _{O} \to 0$. In this limit, 
$A(t,\tau ) \to 0$. As a result, the integrals converge as follows:
\begin{equation} \boldsymbol{I}_1 \to - \varepsilon\int_0^t {\rm{d}} \tau\,\frac{{\mathbb{I}}}{\sqrt{{\rm \pi}(t-\tau)}}\boldsymbol{\cdot} \frac{{\rm{d}} \boldsymbol{f}^{(0)}(\tau)}{{\rm{d}} \tau}\quad \mbox{and} \quad \boldsymbol{I}_2 \to \boldsymbol{0}, \end{equation}
where 
${\mathbb{I}}$ is the identity tensor. So the kernel 
${\mathbb{K}}^{(1)}$ in 
$\boldsymbol {I}_1$ simplifies to the BBO kernel, which describes how the disturbance-velocity gradients relax due to viscous diffusion. In this limit, the 
$\varepsilon$ terms dominate over the 
${Re}_{{p}}$ terms, and the outer flow (3.5a) converges to the outer flow obtained from the unsteady Stokes equation. This is the limit considered by Wang & Ardekani (Reference Wang and Ardekani2012b) and Redaelli et al. (Reference Redaelli, Candelier, Mehaddi and Mehlig2022). The fact that the BBO history force is contained in (4.1) shows that the perturbation theory is valid uniformly in 
$\varepsilon$ (as in Mehaddi, Candelier & Mehlig Reference Mehaddi, Candelier and Mehlig2018). In other words, (4.1) is accurate if the particle Reynolds number is small, regardless of how large the unsteadiness of the problem. Note, however, that when the unsteadiness is very large, for 
$\varepsilon \gg 1$, there is an additional contribution to the hydrodynamic force: the added-mass force (Wang & Ardekani Reference Wang and Ardekani2012b; Candelier et al. Reference Candelier, Mehaddi, Mehlig and Magnaudet2023), which is of order 
$\varepsilon ^2$.
In summary, there are two fundamental differences between the hydrodynamic forces on an active compared to a passive sphere in a time-dependent uniform flow. First, while the history kernels are likely to be the same, the other factors in the integrals are different because they involve an active part for the active particle. Second, for an active particle, there are additional instantaneous contributions to the hydrodynamic force that stem from an inhomogeneous solution of the inner problem. Spherical symmetry causes these contributions to vanish for the passive sphere. This does not happen for the spherical swimmer, because swimming breaks the spherical symmetry.
4.2. Effect of convective and unsteady fluid inertia on the dynamics
As an example, consider the motion of a small neutrally buoyant swimmer that jumps. Its motion starts suddenly, followed by a sudden stop. We mimic the sudden start/stop by imposing a corresponding time dependence on the tangential surface velocity:
\begin{equation} B_1(t) = [1 - \mbox{erf}\left( {2t}/{\tau_\nu}- 10 \right)]\mbox{erf}\left({2t}/{\tau_\nu}\right)^2 (\mbox{mm}\ \mbox{s}^{{-}1}) \quad \mbox{and} \quad B_2(t)=\tfrac{3}{2} B_1(t) \end{equation}
in dimensional variables. The error function squared in (4.2a,b) ensures that 
$\boldsymbol {f}^{(0)}$ and its time derivative both vanish at 
$t=0$, as required by (4.1).
Consider how to non-dimensionalise the problem. The maximal value of 
$B_1(t)$ produces swimming speed 
$\tfrac {4}{3}\ {\rm mm}\ {\rm s}^{-1}$ in the Stokes limit, and we chose to non-dimensionalise with 
$u_c=\tfrac {4}{3}\ {\rm mm}\ {\rm s}^{-1}$. For a particle of radius 
$a = 150\ \mathrm {\mu }{\rm m}$ in water (
$\nu = 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$), this gives 
${Re}_{{p}}=0.2$. How to choose the time scale 
$\tau _c$ is perhaps less obvious. The typical time of variation of the boundary conditions on the surface of the swimmer can be estimated by 
$B_1(t)/\dot {B}_1(t)$. At short times, 
$B_1(t)/\dot {B}_1(t) \approx \tfrac {1}{2} t$. This shows that the characteristic time tends to zero at very short times. After a while, however, the parameter 
$B_1(t)$ reaches a plateau, the disturbance flow becomes steady, and the characteristic time tends to infinity. As discussed at the end of § 4.1, (4.1) is valid regardless of how unsteady the motion is. It is required only that the particle Reynolds number is small. As a consequence, it does not matter precisely how the time scale 
$\tau _c$ is chosen. We took 
$\tau _c= \tau _\nu$.
To determine the trajectory, one needs to solve Newton's equation of motion. In non-dimensional variables, it reads
\begin{equation} \frac{4 {\rm \pi}}{3}\,\varepsilon^2 \ddot{{\boldsymbol{x}}} = \boldsymbol{f}' . \end{equation}
Both fluid and the particle inertia are accounted for in 
$\boldsymbol f^\prime$. Considering the expression (4.1) for the hydrodynamic force, we see that (4.6) is an integro-differential equation. We solved this equation numerically using the method described by Daitche (Reference Daitche2013). When only unsteady fluid inertia is considered, one can solve the equation of motion by Laplace transform (Wang & Ardekani Reference Wang and Ardekani2012b; Ishimoto Reference Ishimoto2013; Fouxon & Or Reference Fouxon and Or2019). In our case, this is not possible, because the equation is nonlinear.
Figure 2 illustrates how the dynamics and the disturbance flow develop after a sudden start followed by a sudden stop, obtained by solving (4.1) and (4.6) together with (4.2a,b). Figure 2(a) shows how the centre-of-mass speed 
$\dot {x}$ changes when both unsteady and convective fluid inertia are taken into account (solid black line). At very short times, the centre-of-mass motion of the squirmer is well described by the BBO equation, which neglects the effect of convective fluid inertia (solid grey line). This is expected, because the dynamics is dominated by the unsteadiness of the disturbance flow at short times. Also shown in figure 2(a) is the Stokes swimming speed 
$(2/3)\,B_1(t)$ (dashed grey line). We see that the actual swimming velocity differs significantly from the Stokes approximation: during the acceleration phase, 
$\dot {x}$ is lower than 
$(2/3)\,B_1$. This produces a negative value of the component of the force 
$\boldsymbol {f}^{(0)}$ along the swimming direction 
$\boldsymbol {n}$ (figure 2b), and this affects the uniform contribution 
$\boldsymbol {\mathcal {U}}$ (figure 2c).
Sudden start followed by sudden stop (see (4.2a,b)). (a) Time dependence of the resulting centre-of-mass speed 
$\dot {x}$ non-dimensionalised with 
$u_c$ and 
$t_\nu$ (solid black line). Also shown is the BBO approximation obtained using (2.12) of Wang & Ardekani (Reference Wang and Ardekani2012b) (solid grey line). The dashed grey line shows the Stokes approximation for the centre-of-mass speed. The horizontal dashed black line corresponds to the steady-state speed (Wang & Ardekani Reference Wang and Ardekani2012a; Khair & Chisholm Reference Khair and Chisholm2014). (b) Time dependence of 
$\boldsymbol {f}^{(0)} \boldsymbol {\cdot } \boldsymbol {n}$ (see (3.4c)). (c) Time dependence of 
$\boldsymbol {\mathcal {U}}\boldsymbol {\cdot } \boldsymbol {n}$. (d–g) Disturbance flow produced by the squirmer at the different non-dimensional times 
$t=0.5$, 3, 5 and 
$5.6$, with the streamlines as well as the magnitude.
At larger times, during the plateau of 
$B_1(t)$, differences between the short-time BBO approximation and the present theory appear, caused by the instantaneous contributions from the inner solution to the hydrodynamic force (the last two terms on the right-hand side of (4.1)). Although the inertial parameters are quite small (as they must be for the theory to be valid), these terms nevertheless have a noticeable effect upon the centre-of-mass speed. Khair & Chisholm (Reference Khair and Chisholm2014) and Wang & Ardekani (Reference Wang and Ardekani2012a) derived an approximation for the centre-of-mass speed in the steady limit, where 
$B_1$ and 
$B_2$ are independent of time:
\begin{equation} \dot{\boldsymbol{x}} \approx \frac{2}{3}\,B_1\left(1 - \frac{3 \beta}{20}\,{Re}_{{p}}\right)\boldsymbol{n}. \end{equation}
This result is shown as a horizontal dashed black line in figure 2(a). Comparing with our theory, we observe that 
$\dot {x}$ reaches its steady limit after a few viscous times. As a consequence, the magnitude of the hydrodynamic force decreases to 
$\boldsymbol {f}^{(0)} \approx -({6 {\rm \pi}}/{10})B_2\, {Re}_{{p}}\,\boldsymbol {n}$ (where 
$B_2 = \beta B_1$ was used), and 
$\boldsymbol {\mathcal {U}}$ decreases to 
$\tfrac {1}{40} B_1 B_2 \,{Re}_{{p}}^2\,\boldsymbol {n}$. Comparing the two expressions, we see that convective-inertia effects due to the singular term 
$\boldsymbol {\mathcal {U}}$ are negligible at small 
$Re_p$, because it contributes only at order 
${O}({Re}_{{p}}^2)$. Now consider the deceleration phase. As 
$B_1 \to 0$, the Stokes velocity decreases to zero. When inertia effects are taken into account, however, 
$\dot {x}$ relaxes more slowly. The present theory and the BBO approximation yield quite similar results in this phase. This suggests that the decrease of the velocity of the organism for a sudden stop is governed mainly by unsteady-inertia effects.
Figures 2(d–g) illustrate how the disturbance near the swimmer changes as a function of time. Streamlines and contours are shown of the magnitude 
$|{\boldsymbol {w}}_{{in}}^{(0)} + \varepsilon {\boldsymbol {w}}_{{in}}^{(1)}|$ of the disturbance-flow velocity at four different times, marked by arrows in figure 2(a). Figure 2(d) shows the flow during the acceleration phase. The term 
$\boldsymbol {\mathcal {U}}\boldsymbol {\cdot } \boldsymbol {n}$ contributes a Stokeslet to the disturbance flow that decays slowly, as 
$r^{-1}$, and thus dominates far from the swimmer. As a consequence, the disturbance flow extends far from the swimmer. Near its surface, on the other hand, the disturbance must match the imposed tangential surface velocity; in the upper half-plane, this velocity opposes the swimming velocity. Therefore, a growing flow cell forms in front of the swimmer. Figure 2(e) shows the disturbance at a later stage, when the motion is quasi-steady, so that 
$\boldsymbol {\mathcal {U}}$ vanishes. The disturbance flow is essentially a stresslet that decays as 
$r^{-2}$. The disturbance is therefore localised closer to the swimmer. We checked that the flow field is well approximated by the steady solution obtained by Chisholm et al. (Reference Chisholm, Legendre, Lauga and Khair2016). In figure 2( f), we plot the disturbance during the deceleration phase. This case is similar to that shown in figure 2(d), except that the sign of 
$\boldsymbol {\mathcal {U}}$ is reversed. Near the swimmer, the flow cell now forms behind the swimmer. Figure 2(g) shows the disturbance flow at large times. Here, the parameter 
$B_1$ is almost zero. The disturbance flow is therefore similar to that produced by a passive sphere, a Stokeslet, but with a very weak intensity.
5. Conclusions
We determined the hydrodynamic force on a small spherical squirmer in an unsteady, spatially homogeneous flow. We obtained the hydrodynamic force by asymptotic matching of perturbation expansions in the parameter 
$\varepsilon = \sqrt {{Re}_{{p}}\,{Sl}}$, for small particle Reynolds numbers. Our main result (4.1) for the hydrodynamic force generalises the result of Lovalenti & Brady (Reference Lovalenti and Brady1993) from a passive sphere to an active particle, an unsteady spherical squirmer. Equation (4.1) describes how convective inertia changes the kernel of the history force. As explained by Lovalenti & Brady (Reference Lovalenti and Brady1993), convective inertia tends to cause the kernel to decay more rapidly. We believe that the kernels for passive and active particles are identical. We did not demonstrate this analytically, but numerical evaluation for all cases that we considered showed them to be the same. The kernels do not depend upon the particular swimming gait of the squirmer, given by the coefficients 
$B_1(t)$ and 
$B_2(t)$. This information is encoded in 
$\boldsymbol {f}^{(0)}$. For the active particle, this amplitude contains an additional term, compared with the passive sphere, stemming from the active surface-velocity field.
A second difference compared to the result of Lovalenti & Brady (Reference Lovalenti and Brady1993) is that an inhomogeneous part of the inner solution contributes to the hydrodynamic force. Spherical symmetry ensures that this contribution vanishes for a spherical passive particle, but swimming breaks spherical symmetry. This contribution explains why the method used by Redaelli et al. (Reference Redaelli, Candelier, Mehaddi and Mehlig2022) to compute inertial corrections to the hydrodynamic force works in the Saffman limit and for unsteady inertia, but not for the Oseen problem considered here.
Our expression (4.1) for the hydrodynamic force simplifies to known results in two limits. First, when unsteadiness dominates, our result simplifies to that of Wang & Ardekani (Reference Wang and Ardekani2012b) and Redaelli et al. (Reference Redaelli, Candelier, Mehaddi and Mehlig2022), where the history force is determined by the BBO kernel that decays as 
$t^{-1/2}$. Second, when particle inertia is more important than unsteadiness, our expression converges to the steady Oseen approximation obtained by Wang & Ardekani (Reference Wang and Ardekani2012a) and Khair & Chisholm (Reference Khair and Chisholm2014).
In order to illustrate the effects of weak unsteady and convective fluid inertia, we considered a swimmer that suddenly starts its centre-of-mass motion, followed by a sudden stop. We analysed how the disturbance flow created by the swimmer changes as a function of time. During acceleration and deceleration, the disturbance flow is essentially a Stokeslet that decreases more slowly far from the swimmer than the stresslet due to steady swimming. This implies that the swimmer is easier to detect immediately after a sudden start or a sudden stop, because the disturbance can be perceived from further away.
We stress that the theory presented here rests on the assumption that 
$Re_{{p}} \ll 1$. Marine organisms come in many different sizes, and they swim with different speeds and with different swimming gaits. A number of empirical studies have estimated both the particle Reynolds number 
$Re_{{p}}$ and the Strouhal number 
${Sl}$ for different microswimmers. Wadhwa, Andersen & Kiorboe (Reference Wadhwa, Andersen and Kiorboe2014) estimated the Strouhal number for copepods, and concluded that 
${Sl} \sim 1$ for nauplii (with 
$Re_{p} = 5,\ldots, 10$) and for adult copepods that move more vigorously (
$Re_{{p}} \sim 40$). Kiorboe et al. (Reference Kiorboe, Jiang, Goncalves, Nielsen and Wadhwa2014) measured 
$Re_{{p}}$ and 
$Sl$ for copepods in different stages of their evolution, observing 
$Sl \sim 1$, and 
$Re_{{p}}$ values up to 
$70$ (see their table S1). For these values of 
$Re_{{p}}$, our theory most certainly fails. However, there are also organisms that move less vigorously. Mesodinium rubrum, for instance, has particle Reynolds numbers of order of unity (Jiang Reference Jiang2011). Also, cruising copepods tend to swim with smaller particle Reynolds numbers, in the range 
${\sim }0.1\unicode{x2013}0.4$; see tables I and II in Qiu et al. (Reference Qiu, Cui, Climent and Zhao2022). In these cases, small-
$Re_p$ theory may give a qualitatively correct description of thedynamics.
Let us consider in more detail how fluid inertia affects the dynamics of a small neutrally buoyant organism that uses ciliary propulsion to swim. For such organisms, the instantaneous velocity produced by the oscillations of the cilia over the surface of the particle scales as 
$\dot {x} \sim \epsilon _s a \omega$, where 
$\epsilon _s a$ is the amplitude of the oscillation of the cilia (
$\epsilon _s\sim 0.1$ is a small non-dimensional parameter), and 
$\omega$ is the angular frequency of the oscillations. It follows that 
${{Re}}_{{p}} = \epsilon _s a^2 \omega /\nu = \epsilon _s \varepsilon ^2$. So when 
$Re_p$ is unity, unsteady inertia dominates the dynamics, as for instance for Paramecium. In this case, (4.1) shows that the hydrodynamic force is well approximated using the BBO kernel, and adding the instantaneous convective-inertia contribution mentioned above.
When the swimmer is not neutrally buoyant, the particle equation of motion contains an additional term, the buoyancy force. In this case, the amplitude 
$\boldsymbol {f}^{(0)}$ does not tend to zero, even if the swimmer stops swimming, because 
$\boldsymbol {f}^{(0)}$ must balance the buoyancy force. To leading order in 
$Re_p$, the resulting convective-inertia correction is described by our result for the hydrodynamic force.
The disturbance caused by small motile organisms in a marine environment has been described in other ways, by superposing different elementary Stokes solutions. Examples are the impulsive Stokeslet and the impulsive stresslet (Afanasyev Reference Afanasyev2004). Guasto et al. (Reference Guasto, Rusconi and Stocker2012) give examples where this approach fails, because it does not approximate reliably the outer disturbance flow, and they state possible reasons: inertia effect induced by the unsteadiness, buoyancy (giving rise, for example, to Stratlets; Ardekani & Stocker Reference Ardekani and Stocker2010), or a combination of both. At least for small particle Reynolds numbers, our results show how convective and unsteady fluid inertia modify the disturbance flow. At small Reynolds number, the outer flow is universal: shape and swimming gait enter only through an amplitude; the kernels describing the history effect on the outer disturbance flow do not depend on these details. However, for larger particle Reynolds numbers, it remains an open question how history force and disturbance flow depend on the shape and propulsion mechanism of the swimmer.
Funding
The research of B.M. was supported by Vetenskapsrådet, grant no. 2021-4452, and in part also by the Knut and Alice Wallenberg Foundation, grant no. 2019.0079. T.R. and C.E. were supported by funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 834238).
Declaration of interests
The authors report no conflict of interest.
Appendix A
In this appendix, we summarise how the outer solution (3.5) of (3.4) is obtained. These equations are solved by Fourier transform. We use the convention
\begin{equation} \hat{\boldsymbol{w}} = \int \mbox{d}^3 \boldsymbol{r} \,\boldsymbol{w}\exp(- {\rm i} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{r})\quad \mbox{and} \quad \boldsymbol{w} = \frac{1}{8{\rm \pi}^3}\int \mbox{d}^3 \boldsymbol{k} \,\hat{\boldsymbol{w}} \exp( {\rm i} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{r}) . \end{equation}Transforming (3.4) one finds
$$\begin{gather} \boldsymbol{k} \boldsymbol{\cdot} \hat{\boldsymbol{w}}_{{out}} =0, \end{gather}$$
$$\begin{gather}\varepsilon^2\,\frac{\partial \hat{\boldsymbol{w}}_{{out}}}{\partial t} - {\rm i} \varepsilon \left(\frac{\ell_p}{\ell_{O}} \right) (\dot{\boldsymbol{x}}\boldsymbol{\cdot}\boldsymbol{k}) \, \hat{\boldsymbol{w}}_{{out}}={-} {\rm i} \boldsymbol{k} \, \hat{p}_{{out}} - k^2 \hat{\boldsymbol{w}}_{{out}} + \boldsymbol{f}^{(0)} . \end{gather}$$
The pressure 
$\hat {p}_{{out}}$ can be determined by projecting (A2b) along 
$\boldsymbol {k}$, and using (A2a). Substituting the resulting expression for 
$\hat {p}_{{out}}$ into (A2b) yields
\begin{equation} \varepsilon^2\,\frac{\partial \hat{\boldsymbol{w}}_{{out}}}{\partial t} ={-} k^2 \hat{\boldsymbol{w}}_{{out}} + {\rm i} \varepsilon \left(\frac{\ell_p}{\ell_{O}} \right) (\dot{\boldsymbol{x}}\boldsymbol{\cdot}\boldsymbol{k}) \hat{\boldsymbol{w}}_{{out}} + k^2 \,\hat{{\mathbb{G}}} \boldsymbol{\cdot} \boldsymbol{f}^{(0)} , \end{equation}
where 
$\hat {\mathbb{G}}$ is the Fourier transform of the Green tensor of the Stokes equations,
\begin{equation} [\hat{\mathbb{G}}]_{i j} (\boldsymbol{k})= \frac{1}{k^2}\left( \delta_{ij} - \frac{k_ik_j}{k^2} \right),\quad [{\mathbb{G}}]_{i j} (\boldsymbol{r})= \frac{1}{8 {\rm \pi}} \left(\frac{\delta_{i j}}{r} + \frac{r_i r_j}{r^3} \right). \end{equation}The solution of the inhomogeneous differential equation (A3) reads
\begin{equation} \hat{\boldsymbol{w}}_{{out}} = \frac{k^2}{\varepsilon^2} \int_0^t\mbox{d}\tau \, \hat{{\mathbb{G}}} \boldsymbol{\cdot} \boldsymbol{f}^{(0)}(\tau) \exp\left[-\frac{k^2}{\varepsilon^2}\,(t-\tau) + \frac{{\rm i}}{\varepsilon}\left(\frac{\ell_p}{\ell_{O}}\right) \int_{\tau}^t \mbox{d} \tau'\, \boldsymbol{k}\boldsymbol{\cdot} \dot{\boldsymbol{x}} (\tau') \right] . \end{equation}
Here, we assumed that the force 
$\boldsymbol {f}^{(0)}$ vanishes for 
$t\le 0$.
In order to match the outer solution (A5) to the inner solution, we seek an expansion of the outer solution of the form
\begin{equation}
\hat{{\boldsymbol{w}}}_{{out}} =
\boldsymbol{\hat{\mathcal{T}}}_{{\!reg}}^{(0)} +
\varepsilon\,\boldsymbol{\hat{\mathcal{T}}}_{{\!reg}}^{(1)}
+\varepsilon\,
\boldsymbol{\hat{\mathcal{T}}}_{{\!sing}}^{(1)} + O
(\varepsilon^2).
\end{equation}
The terms 
$\boldsymbol{\hat {\mathcal {T}}}_{{reg}}^{(n)}$ correspond to regular parts of the expansion. From (A3), we are led to
$$\begin{gather} \boldsymbol{\hat{\mathcal{T}}}_{{\!reg}}^{(0)} = \hat{{\mathbb{G}}} \boldsymbol{\cdot} \boldsymbol{f}^{(0)}, \end{gather}$$
$$\begin{gather}\boldsymbol{\hat{\mathcal{T}}}_{{\!reg}}^{(1)} = {\rm i} \left(\frac{\ell_p}{\ell_{O}} \right) (\dot{\boldsymbol{x}}\boldsymbol{\cdot}\boldsymbol{k})\,\frac{ \hat{{\mathbb{G}}}\boldsymbol{\cdot}\boldsymbol{f}^{(0)} }{k^2} . \end{gather}$$Transforming back to configuration space gives
$$\begin{gather} \boldsymbol{\mathcal{T}}_{{\!reg}}^{(0)} = {{\mathbb{G}}} \boldsymbol{\cdot} \boldsymbol{f}^{(0)}, \end{gather}$$
$$\begin{gather} \boldsymbol{\mathcal{T}}_{{\!reg}}^{(1)} ={-} \left(\frac{\ell_p}{\ell_{O}}\right)\dot{\boldsymbol{x}}\boldsymbol{\cdot} \boldsymbol{\nabla}\left[\frac{3 r}{32{\rm \pi}} \left( {\unicode{x1D7D9}} - \frac{1}{3}\,\frac{\boldsymbol{r} \otimes \boldsymbol{r}}{r^2}\right) \boldsymbol{\cdot} \boldsymbol{f}^{(0)} \right]. \end{gather}$$
The term 
$\boldsymbol{\hat{\mathcal{T}}}_{{\!sing}}^{(1)}$ is singular in 
$\boldsymbol {k}$-space (Meibohm et al. Reference Meibohm, Candelier, Rosén, Einarsson, Lundell and Mehlig2016), proportional to 
$\delta (\boldsymbol {k})$. The singular term in the expansion (A6) is determined by evaluating the limit (Childress Reference Childress1964; Saffman Reference Saffman1965; Meibohm et al. Reference Meibohm, Candelier, Rosén, Einarsson, Lundell and Mehlig2016)
\begin{equation}
\boldsymbol{\hat{\mathcal{T}}}_{{\!sing}}^{(1)} = \lim_{\varepsilon \to 0}
\frac{\hat{{\boldsymbol{w}}}_{{out}} -
\boldsymbol{\hat{\mathcal{T}}}_{{\!reg}}^{(0)}}{\varepsilon} -
\boldsymbol{\hat{\mathcal{T}}}_{{\!reg}}^{(1)}. \end{equation}
Evaluating the limit, one finds
\begin{equation} \boldsymbol{\hat{\mathcal{T}}}_{{\!sing}}^{(1)} = 8 {\rm \pi}^3\,\boldsymbol{\mathcal{U}}(t)\, \delta(\boldsymbol{k}), \end{equation}with
\begin{align} \boldsymbol{\mathcal{U}}(t) &= \frac{1}{8{\rm \pi}^3}\int \mbox{d}^3\boldsymbol{k} \Bigg\{ - \int_0^t\mbox{d}\tau\,\hat{{\mathbb{G}}} \boldsymbol{\cdot} \frac{\mbox{d} \boldsymbol{f}^{(0)}(\tau)}{\mbox{d} \tau} \exp({-}k^2 (t-\tau) + 2 {\rm i} \sqrt{t-\tau }\, \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{a}(t,\tau) ) \nonumber\\ & \quad {}+ {\rm i} \left(\frac{\ell_p}{\ell_{O}}\right) \int_0^t \mbox{d} \tau\, \hat{{\mathbb{G}}} \boldsymbol{\cdot} \boldsymbol{f}^{(0)}(\tau)\,(\boldsymbol{k}\boldsymbol{\cdot}\dot{\boldsymbol{x}}(\tau)) \exp({-}k^2 (t-\tau) + 2 {\rm i} \sqrt{t-\tau }\,\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{a}(t,\tau)) \Bigg\}. \end{align}
Here, 
$\boldsymbol {a}$ is given in (3.5d). In order to perform the 
$\boldsymbol {k}$ integration in (A13), we follow Lovalenti & Brady (Reference Lovalenti and Brady1993) and express the vectors in spherical coordinates defined around an axis that moves with the displacement vector 
$\boldsymbol {a}$. We therefore write
\begin{equation} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{a} = k\,A(t,\tau)\cos(\theta_a)\quad \mbox{and} \quad \hat{\boldsymbol{q}} = \frac{\boldsymbol{a}(t,\tau)}{A(t,\tau)}, \end{equation}
where 
$\theta _a$ is the angle between 
$\boldsymbol {k}$ and 
$\boldsymbol {a}$, and 
$A(t,\tau )$ is the norm of 
$\boldsymbol {a}$. Performing the 
$\boldsymbol {k}$-integration yields
\begin{align} \boldsymbol{\mathcal{U}}(t) ={-}\frac{1}{6{\rm \pi}} \int_{0}^{t}\mbox{d}\tau \, {{\mathbb{K}}^{(1)}(t,\tau) }\boldsymbol{\cdot} \frac{{\mbox{d}}\boldsymbol{f}^{(0)}(\tau)}{\mbox{d} \tau} - \left(\frac{\ell_p}{\ell_{O}}\right) \frac{1}{6{\rm \pi}} \int_{0}^{t} \mbox{d}\tau\,{{\mathbb{ K}}^{(2)}(t,\tau)} \boldsymbol{\cdot} \boldsymbol{f}^{(0)}(\tau). \end{align}
The expressions for the kernels 
${\mathbb{K}}^{(1)}$ and 
${\mathbb{K}}^{(2)}$ are given in (3.5b) and (3.5c).
