1. Introduction
Particle-laden flows are ubiquitous in geophysical and engineering contexts (Balachandar & Eaton Reference Balachandar and Eaton2010). In such flows, the particle dynamics controls how the particles disperse, how they settle, rise or deposit, how they modify the heat/mass transfer efficiency of the carrying flow, at which rate they collide and then possibly fragment or agglomerate, etc. Microscopic descriptions of particle suspensions are challenging, especially because ab initio direct numerical simulations (DNS) resolving the local flow past many particles are extremely demanding, both in terms of numerical schemes and computer time.
An alternative is to use the one-way coupling approximation in which one solves for the fluid motion first, without considering the presence of the particles, and then advances in time an equation of motion for the particles, given a model for the hydrodynamic force and torque acting on them. For very small particles, the force and torque may be obtained by solving the (possibly unsteady) Stokes equation for the disturbance flow caused by the particle in motion, with input data given by the undisturbed fluid-velocity field. In this approximation, the advective terms in the Navier–Stokes equation for the disturbance flow are assumed small compared with the viscous term, so that any finite-Reynolds-number effects are neglected.
This simple approach has two main limitations. First, it considers strictly independent particles, neglecting any interactions within the particle population. This approximation works well for very dilute suspensions, but when two particles come close to each other, their relative motion is affected by hydrodynamic interactions. Second, even for a single particle, it is not known in general how corrections to the above zero-Reynolds-number description affect the particle motion. The larger the particles are, the more important these corrections must become.
The above approximation in which effects of fluid inertia associated with unsteadiness are possibly large while those due to advection are neglected, yields a closed-form expression for the total time-dependent force experienced by a small buoyant spherical particle with radius 
$a$ and mass 
$m_p$ moving at velocity 
${\boldsymbol {v}}_p(t)$ in a uniform flow with velocity 
${\boldsymbol {U}}^\infty (t)$. This is the so-called Basset–Boussinesq–Oseen (BBO) equation (Boussinesq Reference Boussinesq1885; Basset Reference Basset1888; Oseen Reference Oseen1910). Defining the mass of fluid 
$m_f$ corresponding to the particle volume, together with the fluid dynamic and kinematic viscosities 
$\mu$ and 
$\nu$, this equation predicts that the total force acting on the particle at time 
$t$ is
\begin{align} {\boldsymbol{f}}(t) ={\boldsymbol{f}}_b (t) - 6 {\rm \pi}\mu a \left\{{\boldsymbol{u}}_s(t) + \int_0^t \frac{1}{\sqrt{{\rm \pi}\nu (t-\tau)/a^2}} \frac{\mbox{d}{\boldsymbol{u}}_s(\tau)}{{\rm d}\tau}\, \mbox{d}\tau\right\}-\frac{1}{2} m_f \left(\frac{\mbox{d} {\boldsymbol{u}}_s}{\mbox{d}t}\right)(t) , \end{align}
in which 
${\boldsymbol {u}}_s(t)= {\boldsymbol {v}}_p(t) - {\boldsymbol {U}}^{\infty }(t)$ is the slip velocity of the particle with respect to the fluid and 
${\boldsymbol {f}}_b= m_f ({\mbox {d} {\boldsymbol {U}}^\infty }/{\mbox {d} t})+ (m_p - m_f){\boldsymbol {g}}$ may be thought of as the generalized buoyancy force acting on the particle, 
${\boldsymbol {g}}$ denoting gravity. In the expression for 
${\boldsymbol {f}}_b$, the contribution 
$m_f ({\mbox {d} {\boldsymbol {U}}^\infty }/{\mbox {d} t})$ results from the non-uniform pressure distribution at the particle surface induced by the acceleration of the undisturbed flow, and is frequently referred to as the ‘pressure-gradient’ force (Batchelor Reference Batchelor1967). The first term within curly braces is the quasi-steady Stokes drag while the second is the well-known Basset–Boussinesq history force resulting from the unsteady viscous diffusion of the disturbance around the particle arising when the slip velocity changes over time. The last term is the so-called added-mass or virtual-mass force. This force results from the no-penetration condition at the particle surface, which constrains the fluid surrounding the particle to react instantaneously if the particle accelerates with respect to the fluid or vice versa. A counterpart of the BBO equation for a spherical bubble at the surface of which the outer fluid obeys a shear-free condition is also available (Gorodtsov Reference Gorodtsov1975; Morrison & Stewart Reference Morrison and Stewart1976; Yang & Leal Reference Yang and Leal1991). In this case, the 
$6{\rm \pi}$ prefactor weighting the curly braces becomes 
$4{\rm \pi}$ and the 
${\rm \pi} \{\nu (t-\tau )/a^2\}^{-1/2}$ Basset–Boussinesq kernel is changed into 
$2\exp \{{9\nu (t-\tau )/a^2}\}\text {erfc}\{\sqrt {{9\nu (t-\tau )/a^2}}\}$. This kernel yields a finite contribution to the total force at short time, in contrast with the Basset–Boussinesq kernel that diverges as 
$t^{-1/2}$. This difference stresses the fact that history effects are much less severe for a shear-free bubble compared with a rigid sphere, thanks to the slip of the outer fluid along the bubble surface.
The BBO equation was extended to arbitrary non-uniform flows by Gatignol (Reference Gatignol1983) and Maxey & Riley (Reference Maxey and Riley1983), still neglecting any advective effect in the disturbance equation (see table 1 for an overview of the different approximations for the force on a small spherical particle that are discussed in this section). Gatignol (Reference Gatignol1983) also obtained the corresponding equation for the hydrodynamic torque, extending the result previously derived in a uniform flow by Feuillebois & Lasek (Reference Feuillebois and Lasek1978). Once advective effects are assumed negligible in the disturbance dynamics, the non-uniformity of the carrying flow manifests itself in two ways. First, since the background velocity varies with the local position, it is necessary to define the slip velocity at the instantaneous position of the particle centre, 
${\boldsymbol {x}}_p(t)$, which leads to the definition
\begin{equation} {\boldsymbol{u}}_s({\boldsymbol{x}}_p(t))={\boldsymbol{v}}_p(t)- {\boldsymbol{U}}^{\infty}({\boldsymbol{x}}_p(t)).\end{equation}
Computing the disturbance-induced force then reveals that, in addition to 
${\boldsymbol {u}}_s$ and 
${\rm d}{\boldsymbol {u}}_s/{\rm d}t$, the last three terms on the right-hand side of (1.1) involve Faxén corrections proportional to 
$a^2\nabla ^2 U^\infty (\boldsymbol{x}_p(t))$. These corrections are significant if the slip velocity is small (nearly neutrally buoyant particles) and the carrying flow varies significantly and nonlinearly at the particle scale. Second, although advective effects are assumed to have no effect on the disturbance flow, they may be significant in the carrying flow. For this reason, they must be taken into account consistently in the ‘pressure-gradient’ contribution involved in the generalized buoyancy force 
${\boldsymbol {f}}_b$. A simple reasoning shows that this is achieved by replacing the time variation of the fluid velocity 
${\mbox {d} {\boldsymbol {U}}^\infty }/{\mbox {d} t}$ with the Lagrangian acceleration 
$\langle{{\mbox {D} {\boldsymbol {U}}^\infty }/{\text {D} t}}\rangle $, the 
$\langle\ \boldsymbol{\cdot }\ \rangle$ symbol denoting the spatial average over the particle volume. If the fluid acceleration may be considered uniform over this volume or if its spatial variations with respect to the particle centre are odd (such as in a linear flow field), 
$\langle {{\mbox {D} {\boldsymbol {U}}^\infty }/{\mbox {D} t}}\rangle$ reduces to 
${\mbox {D} {\boldsymbol {U}}^\infty }/{\mbox {D} t} |_{{\boldsymbol {x}}_p}$, which is the approximation retained by Maxey & Riley (Reference Maxey and Riley1983). However, to remain consistent with the treatment used for the disturbance, variations of the fluid acceleration within the particle volume must generally be taken into account. The corresponding expansion yields an extra Faxén-type correction, as recognized by Gatignol (Reference Gatignol1983), but also additional quadratic contributions proportional to 
$a^2(\nabla^2 \boldsymbol{U}^\infty) \boldsymbol{\cdot }\boldsymbol {\nabla } \boldsymbol{U}^\infty$ and 
$a^2(\boldsymbol {\nabla } \boldsymbol {U}^\infty \boldsymbol{\colon} \boldsymbol {\nabla }\boldsymbol {\nabla }) \boldsymbol {U}^\infty$ (Rallabandi Reference Rallabandi2021).
Table 1. Review of various approximations beyond the quasi-steady Stokes solution for the hydrodynamic force on a small spherical particle. Values in the 
$O({\cdot })$ symbols compare the magnitude of unsteady and/or advective effects with that of viscous effects. The ratio 
$a/\ell _u=a\|{\boldsymbol {u}}_s\|/\nu$ is the slip Reynolds number based on the relative velocity between the particle and fluid, while 
$\varepsilon =a/\ell _s=(a^2s/\nu )^{1/2}$ is the square root of the shear-based Reynolds number; both Reynolds numbers are assumed small. 
$^{(*)}$This includes uniform shear (e.g. Saffman Reference Saffman1965), solid-body rotation (e.g. Gotoh Reference Gotoh1990), and uniform two-dimensional strain; 
$^{(\text {+})}$Rubinow & Keller (Reference Rubinow and Keller1961) considered the combined effect of particle slip and spin.
As already pointed out, (1.1) is obtained by assuming that (i) unsteady effects are strong enough for the time-derivative term in the Navier–Stokes equation to be comparable with the viscous term; and (ii) advective effects are negligible throughout the flow. However, it has been proved both theoretically and numerically that the latter assumption breaks down when the disturbance reaches distances of the order of the Oseen length scale 
$\ell _u=\nu /\|{\boldsymbol {u}}_s\|$ (Sano Reference Sano1981; Mei, Lawrence & Adrian Reference Mei, Lawrence and Adrian1991; Lovalenti & Brady Reference Lovalenti and Brady1993). Indeed, whatever the ratio 
$a/\ell _u$ (which may be thought of as the particle slip Reynolds number), advective effects cannot be neglected at such a distance, and these effects result in a wake. In this wake, advection being more efficient than viscous diffusion, the disturbance adjusts more quickly to the new slip velocity 
${\boldsymbol {u}}_s(t)$, than in the immediate surroundings of the particle, yielding in most cases a 
$t^{-2}$ long-time decay of the history force instead of the 
$t^{-1/2}$ prediction resulting from the Basset–Boussinesq kernel. Based on these theoretical findings, several authors have proposed approximate extensions of the BBO equation to particles moving at finite Reynolds number in a uniform but possibly time-dependent flow. Mei & Adrian (Reference Mei and Adrian1992), Mei (Reference Mei1994) and Kim, Elghobashi & Sirignano (Reference Kim, Elghobashi and Sirignano1998) designed semi-empirical kernels that recover the correct asymptotic behaviours in the short- and long-time limits. Influence of the Reynolds number was incorporated by introducing empirical functions calibrated against DNS results in the kernel, and replacing Stokes’ expression for the quasi-steady drag by empirical fits based on the standard drag curve. These attempts have proved successful up to slip Reynolds numbers of several tens, even in configurations far from those in which the empirical functions were calibrated, such as the unsteady rise of CO
$_2$ bubbles rapidly dissolving in water (Takemura & Magnaudet Reference Takemura and Magnaudet2004).
Compared with the above picture, extensions of the BBO equation to particles experiencing finite advective effects in non-uniform flows are much less mature, to say the least. This does not have severe consequences regarding the prediction of the drag, which is usually only marginally affected by the local fluid-velocity gradients. However, in many configurations, these gradients are known to induce lift components in the hydrodynamic force. For a sphere, this happens every time the slip velocity is not collinear with one of the eigenvectors of the velocity-gradient tensor. Such lift forces strongly affect the particle dynamics in many situations, such as particle deposition and shear-induced migration in wall-bounded shear flows, or the migration of particles toward high-strain regions or vortex cores in turbulent flows, to mention just a few examples. Most studies devoted to these velocity-gradient-induced lift forces considered a steady framework in which the slip velocity is prescribed and does not vary over time. The best known of these studies is that of Saffman (Reference Saffman1965) who examined the case of a spherical particle immersed in a linear shear flow, with a non-zero slip velocity collinear with the background flow. Assuming the shear-based length scale 
$\ell _s=(\nu /s)^{1/2}$ (with 
$s$ the shear rate) to be much smaller than the Oseen length 
$\ell _u$, Saffman employed the technique of matched asymptotic expansions to calculate the dominant contribution to the shear-induced lift force acting on the particle. This seminal work opened the way to a large number of investigations that considered other canonical linear flows or other orientations of the slip velocity with respect to the flow (see Stone (Reference Stone2000) and Candelier, Mehlig & Magnaudet (Reference Candelier, Mehlig and Magnaudet2019) for reviews). The corresponding predictions for the hydrodynamic forces and torques are widely used. However, they are frequently employed well beyond the original framework in which they were derived and are known to be valid. In particular, Saffman's expression for the shear-induced lift force has routinely been added empirically to the right-hand side of (1.1) to compute the path of particles immersed in non-uniform time-dependent environments barely resembling a stationary linear shear flow; e.g. McLaughlin (Reference McLaughlin1989) and Li & Ahmadi (Reference Li and Ahmadi1992).
Few studies have attempted to derive a rigorous expression for the total hydrodynamic force acting on particles moving in steady linear flows with a time-varying slip velocity. Their common point is the central assumption that the leading-order solution is governed by the quasi-steady Stokes equation. In other words, effects of unsteadiness are assumed to take place over a characteristic time of the order of the inverse shear rate, 
$s^{-1}$, allowing the time rate-of-change term to be treated as a perturbation, similar to the shear-induced advective terms. This is in contrast with the assumption on which (1.1) is grounded, which considers that velocity variations may take place over a characteristic time possibly as short as the viscous time 
$a^2/\nu$. With the former assumption, the first-order corrections to the quasi-steady Stokes force arise at order 
$\varepsilon =a/\ell _s$ and include contributions due to the time variations of the slip velocity as well as to the nonlinear advective interaction of the disturbance with the background velocity field. Miyazaki, Bedeaux & Avalos (Reference Miyazaki, Bedeaux and Avalos1995) employed the so-called induced-force method to derive the 
$O(\varepsilon )$ correction to the force in the case of a linear shear flow, recovering Saffman's prediction in the long-time limit. Asmolov & McLaughlin (Reference Asmolov and McLaughlin1999) made use of the matched asymptotic expansions technique to obtain the time-dependent lift force in the same configuration, in the specific case of a sphere undergoing periodic oscillations. More recently, Candelier et al. (Reference Candelier, Mehlig and Magnaudet2019), hereinafter referred to as CMM, expressed the governing equations for the disturbance in a non-orthogonal coordinate system moving with the undisturbed flow and solved them using matched asymptotic expansions for the three canonical planar linear flows, namely solid-body rotation, uniform straining and uniform shear. At this point, it is important to stress that the problem is nonlinear, owing to the nonlinearity of the advective term. Hence, the solution for an arbitrary linear flow cannot be obtained via a linear superposition of the elementary contributions corresponding to solid-body rotation and uniform straining motion (Candelier & Angilella Reference Candelier and Angilella2006). Candelier et al. (Reference Candelier, Mehlig and Magnaudet2019) expressed the 
$O(\varepsilon )$ correction to the force and torque in the form of a convolution integral involving a tensorial kernel whose components are specific to the linear flow under consideration. This kernel does not depend on the shape of the particle. For this reason, their approach allows the results obtained for a sphere to be straightforwardly extended to arbitrarily shaped particles, simply by performing the dot product of the convolution integral with the appropriate resistance tensor of the particle determined in the creeping-flow limit. Since effects of unsteadiness are assumed to be small compared with viscous effects, the results obtained through this approach are in general not valid at very short times following the introduction of the particle in the flow, i.e. for times 
$t\lesssim a^2/\nu$. In contrast, these results are valid in the intermediate range 
$a^2/\nu \lesssim t\lesssim s^{-1}$ that corresponds to short times with respect to the characteristic time imposed by the velocity gradient. At such ‘short’ times, the kernel is diagonal and its non-zero components behave as 
$t^{-1/2}$, recovering the contribution of the Basset–Boussinesq term in (1.1) to the total force and torque. Corrections to this initial behaviour develop gradually over time, both in the diagonal and off-diagonal components. Their evolution depends on the background flow. For instance, the off-diagonal components, responsible for the lift force, grow as 
$t^{1/2}$ in a uniform shear flow, but only as 
$t^{5/2}$ in a solid-body rotation flow. Each component eventually converges toward a steady-state value for 
$t\gg s^{-1}$ but the duration of the corresponding transient significantly varies from one component to the other. At this point, we need to stress that the validity of the BBO equation in the presence of a linear background flow is limited to times usually significantly shorter than 
$s^{-1}$. For instance, figure 4 in CMM indicates that in a pure shear flow, the lift component that eventually yields the Saffman lift force has already reached 
$20\,\%$ (respectively 
$80\,\%$) of its final value at 
$t=\frac {1}{10}s^{-1}$ (respectively 
$t=s^{-1}$), an effect which is not captured by the BBO approximation. As a consequence, the BBO equation describes, for example, the dynamics of a spherical particle with radius 
$a=100\ \mathrm {\mu }$m moving in a water flow with 
$s=1$ s
$^{-1}$ up to 
$t=a^2/\nu =10^{-2}$ s, but fails to capture the growth of the shear-induced lift that becomes significant for 
$t=\frac {1}{10}s^{-1}=10^{-1}$ s. The situation becomes obviously worse as the shear rate increases, showing that, in this range of particle sizes, the time interval over which the BBO equation is valid does not exceed a few characteristic viscous times.
The approach outlined above yields a consistent prediction of the hydrodynamic force and torque at 
$O(\varepsilon )$ in linear flows, provided the slip velocity does not vary ‘too’ rapidly, and advective effects due to the linear background flow dominate over those due to the slip velocity, i.e. 
$\ell _s\ll \ell _u$. Thus, it may be viewed as the desired rational first-order extension of Stokes’ quasi-steady prediction incorporating finite inertial effects, be they due to unsteadiness or to advection by the linear background flow. However, referring to (1.1), it is clear that this first-order extension does not include any added-mass contribution, owing to the restriction imposed to the time variation of the slip velocity. Such a contribution appears only at next order in the expansion with respect to 
$\varepsilon$.
For reasons of simplicity, many theories for the dynamics of inertial particles in turbulence just use the quasi-steady Stokes approximation without added-mass or history terms (Gustavsson & Mehlig Reference Gustavsson and Mehlig2016). For similar reasons, and also to limit the computational cost, a number of studies devoted specifically to the dynamics of light particles in turbulence retain the added-mass force but ignore any 
$O(\varepsilon )$ contribution, which makes the underlying model formally inconsistent and questions the relevance of its predictions (Babiano et al. Reference Babiano, Cartwright, Piro and Provenzale2000; Bec Reference Bec2003; Calzavarini et al. Reference Calzavarini, Cencini, Lohse and Toschi2008a,Reference Calzavarini, Kerscher, Lohse and Toschib, Reference Calzavarini, Volk, Bourgoin, Leveque, Pinton and Toschi2009; Volk et al. Reference Volk, Calzavarini, Verhille, Lohse, Mordant, Pinton and Toschi2008; Vajedi et al. Reference Vajedi, Gustavsson, Mehlig and Biferale2016). Added-mass effects are known to be central in the dynamics of light particles, not to mention bubbles. Therefore, determining the 
$O(\varepsilon ^2)$ corrections to the above expansion appears as essential to provide a robust model for studying the dynamics of relatively light particles in which all potentially physically relevant building blocks are included. Similarly, the possible rotation of the particle does not have any influence on the 
$O(\varepsilon )$ correction to the force and torque. However, it is well known that, as soon as effects of fluid inertia come into play, a spinning particle translating in a fluid at rest experiences a lift force proportional to the cross-product of the slip velocity and the spinning rate (Rubinow & Keller Reference Rubinow and Keller1961), a clear 
$O(\varepsilon ^2)$ effect. Other cross-effects affecting the particle dynamics may be expected at the same order, due to the possible quadratic terms arising from the various combinations of the slip velocity, spinning rate, strain rate and vorticity of the background flow.
When the added-mass force is mentioned, the known expression for the inviscid hydrodynamic force on a sphere in motion in a general linear flow field comes to mind. In present notation, this expression reads (Auton, Hunt & Prud'homme Reference Auton, Hunt and Prud'homme1988)
\begin{equation} {\boldsymbol{f}}_{AHP} = {\boldsymbol{f}}_b - \frac{1}{2} m_f \left( \left.\frac{\mbox{d} {\boldsymbol{v}}_p}{\mbox{d}t} - \frac{\mbox{D} {\boldsymbol{U}}^\infty}{\mbox{D} t} \right|_{{\boldsymbol{x}}_p} \right) - \frac{1}{2} m_f {\boldsymbol{u}}_s \times ({\boldsymbol{\nabla}} \times {\boldsymbol{U}}^\infty), \end{equation}
with 
${\boldsymbol {f}}_b=m_f({\mbox {D} {\boldsymbol {U}}^\infty }/{\mbox {D} t}) |_{{\boldsymbol {x}}_p}+(m_p-m_f){\boldsymbol {g}}$ in the case of a linear flow. In (1.3) the disturbance-induced force appears to result from the addition of an added-mass force and a shear-induced lift force. The latter was first computed by Auton (Reference Auton1987) for a sphere held fixed in a stationary uniform shear flow, under the condition that the shear-induced velocity at the body scale is weak compared with the slip velocity, i.e. 
$\eta =as / \|{\boldsymbol {u}}_s\|\ll 1$. This condition is needed for the velocity correction induced by the distortion of the background vorticity to remain small compared with the slip-induced velocity at the body surface, which in turn allows the pressure distribution at this surface, hence the force, to be computed at first order with respect to 
$\eta$.
The mathematical form of the added-mass force in (1.3) involves the Lagrangian acceleration of the background flow, in agreement with the expression established by Taylor (Reference Taylor1928) for pure straining (i.e. irrotational) flows. Noting that 
${\mbox {D} {\boldsymbol {U}}^\infty }/{\mbox {D} t} |_{{\boldsymbol {x}}_p}={\mbox {d} {\boldsymbol {U}}^\infty }/{\text {d} t}|_{{\boldsymbol {x}}_p}-{\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {\nabla }}{\boldsymbol {U}}^\infty$, with 
${\mbox {d} {\boldsymbol {U}}^\infty }/{\text {d} t}|_{{\boldsymbol {x}}_p}$ the time derivative of 
${\boldsymbol {U}}_\infty$ following the particle path, allows this force, say 
${\boldsymbol {f}}_{am}$, to be rewritten in the equivalent form
\begin{equation} {\boldsymbol{f}}_{am} ={-} \frac{1}{2} m_f \left( \frac{\mbox{d} {\boldsymbol{v}}_p}{\mbox{d}t} - \left.\frac{\mbox{D} {\boldsymbol{U}}^\infty}{\mbox{D} t} \right|_{{\boldsymbol{x}}_p} \right)\equiv{-} \frac{1}{2} m_f \left( \left.\frac{\mbox{d} {\boldsymbol{u}}_s}{\textrm{d}t}+{\boldsymbol{u}}_s\boldsymbol{\cdot}{\boldsymbol{\nabla}} {\boldsymbol{U}}^\infty\right|_{{\boldsymbol{x}}_p}\right), \end{equation}
with the slip velocity evaluated at the particle centre. The latter form of 
${\boldsymbol {f}}_{am}$ makes it clear that this force may be non-zero even though the slip velocity does not change over time, provided the carrying flow varies in the direction collinear with the slip. According to (1.3), Taylor's expression applies even if the carrying flow has non-zero vorticity. Moreover, Auton's expression for the shear-induced lift force appears to apply even if the carrying flow is unsteady. However, this extension holds only as long as the magnitude of the slip rate of change, 
$\|{\rm d}{\boldsymbol {u}}_s/{\rm d}t\|$, is small compared with 
$\|{\boldsymbol {u}}_s\|^2/a$, in which case the stretching/tilting term in the vorticity equation remains primarily balanced by the advective term. The above conditions limit the validity of (1.3) to weakly unsteady and non-uniform flows, owing to the consequences of the distortion of the upstream vorticity in the three-dimensional flow past a sphere. These limitations do not exist in the two-dimensional case. That is, (1.3) (with the prefactors 
$\frac {1}{2}$ replaced with 
$1$) is an exact solution for the inviscid force acting on a circular cylinder immersed in a general linear time-dependent flow.
The connection between the inviscid force (1.3) and the proper generalization of the visco-inertial force (1.1) to non-uniform flows has been considered by several authors (Maxey & Riley Reference Maxey and Riley1983; Magnaudet, Rivero & Fabre Reference Magnaudet, Rivero and Fabre1995). In particular, Maxey & Riley (Reference Maxey and Riley1983) pointed out that the creeping-flow limit does not allow one to decide whether the form of the added-mass term in (1.1) remains unchanged when the carrying flow is non-uniform, i.e. whether it still involves the time derivative of the slip velocity following the particle, or whether it rather includes the 
$-\frac {1}{2}m_f{\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {\nabla }}{\boldsymbol {U}}^\infty |_{{\boldsymbol {x}}_p}$ contribution resulting from the Lagrangian fluid acceleration in (1.4). Indeed, this term is a quadratic cross-effect in the sense employed above, and consequently appears only at 
$O(\varepsilon ^2)$. Thus, computing the second-order inertial corrections is a necessary step to clarify this issue, as well as to quantify the role of the no-slip condition on the shear-induced lift force by comparing the corresponding prediction with that of (1.3).
Given the open questions identified in this introduction, there is a clear need to calculate all 
$O(\varepsilon ^2)$ contributions to the force and torque acting on a rigid spherical particle immersed in a stationary linear flow. Indeed, as shown above, these contributions contain several physical effects that are not captured at 
$O(\varepsilon )$, and are thus of fundamental importance to better understand how (possibly coupled) effects of slip, spin, unsteadiness and background velocity gradients modify the simpler picture provided by the 
$O(\varepsilon ^0)$ and 
$O(\varepsilon )$ models. Obviously, one intuitively expects second-order contributions to only produce small changes in the particle dynamics since 
$\varepsilon$ is assumed small. However, lift effects on a sphere only arise at 
$O(\varepsilon )$. Therefore, it might be that the sign and magnitude of the 
$O(\varepsilon )$ and 
$O(\varepsilon ^2)$ contributions to the lift force combine in such a way that the force becomes dominated by second-order effects beyond some critical, but still small, 
$\varepsilon$. Indeed, results discussed later will confirm this possibility in some flows, providing an additional confirmation that these effects are worthy of investigation.
In the following, we describe how we computed all 
$O(\varepsilon ^2)$ corrections to the force and torque on a small rigid sphere moving in a steady linear flow with a time-dependent slip velocity. We formulate the problem and the underlying assumptions in § 2. The technique employed to obtain the corresponding contributions makes use of matched asymptotic expansions. The way the outer problem is solved at the required order is a direct extension of the approach developed by CMM to obtain the 
$O(\varepsilon )$ corrections. We summarize this approach in § 3, show how it extends to order 
$O(\varepsilon ^2)$, and insist on the solution of the inner problem that is required to obtain the complete set of second-order contributions. In § 4 we extensively discuss the general results derived in the previous section by considering first the short-term limit (which in the framework of the present theory corresponds to the intermediate range 
$a^2/\nu \ll t \ll s^{-1}$) for the expressions of the force and torque, then the long-term limit 
$t\gg s^{-1}$ in four canonical linear flow configurations. We summarize the main outcomes of the study in § 5, providing in particular the complete expressions for the force and torque up to 
$O(\varepsilon ^2)$ terms in both the short- and long-time limits. Readers mostly interested in applications may directly consider the final results (5.1)–(5.5), together with the specific expressions of the kernel in the different flow configurations, to obtain a complete view of the effects involved in the force and torque balances at the order of approximation considered here.
2. Basic assumptions and disturbance-flow equations
We consider a spherical particle of radius 
$a$ moving freely in a linear flow of the form
\begin{equation} {\boldsymbol{U}}^\infty({\boldsymbol{x}},t) ={\boldsymbol{U}}_0(t)+ {\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{x}} ,\end{equation}
where 
${\boldsymbol {x}}$ is the position vector in the laboratory reference frame whose origin is arbitrary, and 
$\boldsymbol{U}_0$ denotes the (possibly time-dependent) undisturbed velocity 
$\boldsymbol{U}_\infty(\textbf{0},t)$ at 
$\boldsymbol{x}=\textbf{0}$ (see figure 1 in which 
$\boldsymbol{U}_0$ has been set to zero). For reasons discussed in CMM and outlined below in § 3.1, a major simplification is introduced in the calculation of the disturbance induced by the particle in the far field by assuming that 
${\mathbb {A}}$ does not depend upon time, i.e. the undisturbed flow is stationary. However, this puts some restriction on the class of linear flows that can be considered, since 
${\mathbb {A}}$ is also uniform. Indeed, combining these two assumptions implies that the balance for the possibly non-zero vorticity 
${\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty$ reduces to the stretching term 
${\mathbb {A}}\boldsymbol{\cdot }({\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )\equiv {\mathbb {S}}\boldsymbol{\cdot }({\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )$, with 
${\mathbb {S}} = \tfrac {1}{2}({\mathbb {A}} +{\mathbb {A}}^{\sf T})$ the strain-rate tensor, the superscript 
$^{\sf T}$ denoting the transpose. Consequently, this stretching term must be zero, a constraint satisfied by every planar base flow. In contrast, in three-dimensional configurations, this constraint implies that the general solutions obtained below only hold for purely irrotational base flows.
Figure 1. Sketch of a spherical particle moving in a steady linear flow, with some definitions used throughout the paper.
Since the undisturbed fluid velocity is linear in 
$\boldsymbol{x}$, no Faxén corrections can be involved in the loads acting on the particle. Computing the second-order inertial contributions to the force and torque requires solving the Navier–Stokes equation for the disturbance flow 
$\boldsymbol{w}$ caused by the particle motion relative to the undisturbed fluid. To identify the order of magnitude of each term involved in this equation and in the associated boundary conditions, we normalize the fluid and particle translational velocities with a velocity scale 
$u_c$, the angular velocities with 
$\omega _c$, distances with the sphere radius 
$a$, components of 
${\mathbb {A}}$ with a characteristic strain rate 
$s$, pressure with 
$\mu u_c/a$ and time with a characteristic time 
$\tau _c$ over which the relative translational and angular velocities may vary. Using these normalisations, and keeping the previous notations unchanged although all variables are non-dimensional from now on, the equations that govern the disturbance become (CMM)
$$\begin{gather} {\boldsymbol{\nabla}} \boldsymbol{\cdot} {{\boldsymbol{w}}} = {\boldsymbol{0}}, \end{gather}$$
$$\begin{gather}{Re}_s {Sl} \left.\frac{\partial {\boldsymbol{w}}}{\partial t} \right|_{{\boldsymbol{r}}} + {Re}_s ( {\mathbb{A}} \boldsymbol{\cdot} {\boldsymbol{w}} + ({\mathbb{A}} \boldsymbol{\cdot} {\boldsymbol{r}}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}} ) +{Re}_p (-{\boldsymbol{u}}_s \boldsymbol{\cdot}{\boldsymbol{\nabla}} {\boldsymbol{w}} + {\boldsymbol{w}}\boldsymbol{\cdot}{\boldsymbol{\nabla}} {\boldsymbol{w}} ) ={-} {\boldsymbol{\nabla}} p +\boldsymbol{\triangle}{{\boldsymbol{w}}} , \end{gather}$$
$$\begin{gather}{\boldsymbol{w}} \to {\boldsymbol{0}} \quad \mbox{for}\ r \to \infty, \end{gather}$$
with, as shown in figure 1, 
${\boldsymbol {r}}={\boldsymbol {x}}-{\boldsymbol {x}}_p$ the local distance to the particle centre (hence, 
$r=\|{\boldsymbol {r}}\|=1$ at the particle surface), the time derivative in (2.2b) being evaluated at fixed 
${\boldsymbol {r}}$. According to (1.2) and (2.1), the slip velocity in (2.2b) is 
${\boldsymbol {u}}_s(t)= {\boldsymbol {v}}_p(t)-{\mathbb {A}}\boldsymbol{\cdot } {\boldsymbol {x}}_p$. However, the entire problem is left unchanged if, in addition to the linear stationary component 
${\mathbb {A}}\boldsymbol{\cdot }{\boldsymbol {x}}$, the undisturbed flow is assumed to comprise a uniform time-dependent component, say 
${\boldsymbol {U}}_0(t)$, provided the slip velocity is redefined in accordance with (1.2) as 
${\boldsymbol {u}}_s(t)= {\boldsymbol {v}}_p(t)-{\boldsymbol {U}}_0(t)-{\mathbb {A}}\boldsymbol{\cdot } {\boldsymbol {x}}_p$. For a rigid spherical particle, the no-slip condition at the particle surface implies that
\begin{equation} {\boldsymbol{w}} = {\boldsymbol{u}}_s + \frac{{Re}_\omega}{{Re}_p} \boldsymbol{\omega}_p \times {\boldsymbol{r}}- \frac{{Re}_s}{{Re}_p}( \boldsymbol{\varOmega} \times {\boldsymbol{r}} + {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{r}} ) \quad \mbox{for}\ r=1, \end{equation}
with 
$\boldsymbol {\omega }_p$ the particle angular velocity and 
$\boldsymbol {\varOmega } =\tfrac {1}{2}{\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty$ half the undisturbed flow vorticity, the antisymmetric tensor associated with 
${\boldsymbol {\varOmega }}$ being the antisymmetric part of 
${\mathbb {A}}$. The non-dimensional parameters in (2.2) and (2.3) are the rotation, shear and slip Reynolds numbers, respectively, plus a Strouhal number characterizing the magnitude of the time-derivative term in (2.2b), namely
\begin{equation} {Re}_\omega = \frac{a^2 \omega_c}{\nu},\quad {Re}_s= \frac{a^2 s}{\nu},\quad {Re}_p = \frac{a u_c}{\nu} ,\quad {Sl}=\frac{1}{s \tau_c}. \end{equation}
To simplify (2.2) and (2.3), we assume that the particle is small and only weakly positively or negatively buoyant. Therefore, its slip velocity is expected to be small and so is the slip Reynolds number, 
${Re}_p$. For the same reason, if the particle is free to rotate, its angular velocity 
$\boldsymbol {\omega }_p$ is assumed to be close to 
$\boldsymbol {\varOmega }$. We nevertheless keep track of the relative (or slip) angular velocity 
$\boldsymbol {\omega }_s(t)=\boldsymbol {\omega }_p(t)-\boldsymbol {\varOmega }$ in order to allow for a possible transient regime or for a forced particle spin. However, we assume that the magnitude of the particle angular velocity remains of the same order as the shear rate 
$s$ of the undisturbed flow. Having defined
\begin{equation} \varepsilon\equiv\frac{a}{\ell_s}\equiv{Re}_s^{1/2}, \end{equation}we therefore set
$$\begin{gather} u_c\equiv a s,\quad \mbox{so that} \quad {Re}_p=\varepsilon^2 \ll{Re}_s^{1/2}, \end{gather}$$
$$\begin{gather}\omega_c \equiv s \quad \mbox{and} \quad \tau_c \equiv s^{{-}1}, \quad \mbox{so that} \quad {Re}_\omega = {Re}_s=\varepsilon^2 \quad \mbox{and} \quad {Sl}=1. \end{gather}$$
The assumption 
${Re}_p\ll {Re}_s^{1/2}$ corresponds to the limit considered by Saffman in the stationary 
$O(\varepsilon )$ problem (Saffman Reference Saffman1965). It allows effects of slip (but not those due to unsteadiness) to be disregarded in the far field, and we shall show in § 3.1 that this simplification still holds at 
$O(\varepsilon ^2)$. Conditions (2.6b) are less critical and simply allow effects of slip, shear, unsteadiness and spin to all influence the disturbance close to the particle at the retained order of approximation. With these assumptions, the non-dimensional Navier–Stokes equations for the disturbance take the form
$$\begin{gather} {\boldsymbol{\nabla}} \boldsymbol{\cdot} {{\boldsymbol{w}}} = {\boldsymbol{0}} , \end{gather}$$
$$\begin{gather}\varepsilon^2 \frac{\partial {\boldsymbol{w}}}{\partial t} + \varepsilon^2 [ {\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{w}} + ({\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{r}}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}} ] + \varepsilon^2 (-{\boldsymbol{u}}_s \boldsymbol{\cdot}{\boldsymbol{\nabla}} {\boldsymbol{w}} + {\boldsymbol{w}}\boldsymbol{\cdot}{\boldsymbol{\nabla}} {\boldsymbol{w}} ) ={-} {\boldsymbol{\nabla}} {p} + \boldsymbol{\triangle}{{\boldsymbol{w}}}, \end{gather}$$
$$\begin{gather}{\boldsymbol{w}} = {\boldsymbol{u}}_s(t) + [\boldsymbol{\omega}_s(t) \times {\boldsymbol{r}} - {\mathbb{S}}\boldsymbol{\cdot} {\boldsymbol{r}} ]\quad\mbox{for} \ r=1 , \end{gather}$$
$$\begin{gather}{\boldsymbol{w}} \to {\boldsymbol{0}} \quad\mbox{for}\ r \to \infty . \end{gather}$$
The disturbance is assumed to be zero for 
$t<0$. To avoid singularities in the force and torque at 
$t=0$, we assume that the particle is introduced in the flow with zero slip velocity and relative rotation rate (
${\boldsymbol {u}}_s(0)={\boldsymbol {0}},\ \boldsymbol {\omega }_s(0)={\boldsymbol {0}}$), but let both quantities have arbitrary non-zero initial time derivatives (
$({\rm d}{\boldsymbol {u}}_s/{\rm d}t)(0)\neq {\boldsymbol {0}},\ ({{\rm d}\boldsymbol {\omega }_s}/{{\rm d}t})(0)\neq {\boldsymbol {0}}$). In the usual case of freely moving particles, the subsequent evolution of 
${\boldsymbol {u}}_s$ and 
$\boldsymbol {\omega }_s$ is dictated by the overall force and torque balances on the particles. Here in contrast, we are interested in predicting the hydrodynamic loads resulting from arbitrary evolutions of both quantities. The expressions for the force and torque obtained in this way may then be inserted in the overall force and torque balances, which usually involve non-zero external forces (such as the generalized buoyancy force 
${\boldsymbol {f}}_b$ introduced in § 1) and possibly external torques, and the actual evolution of 
${\boldsymbol {u}}_s$ and 
$\boldsymbol \omega _s$ ensues. Starting from (2.7), CMM computed the force and torque to 
$O(\varepsilon )$. They disregarded terms that do not contribute to the loads at this order. This includes the term within square brackets on the right-hand side of (2.7c) because its contribution to the disturbance flow (through a rotlet and a stresslet) vanishes at this order for symmetry reasons, and the advective term due to the slip velocity in (2.7b) which is negligible compared with that due to shear in the far field in the Saffman limit. However, this term contributes to the inner solution at 
$O(\varepsilon ^2)$.
Here our goal is to determine all 
$O(\varepsilon ^2)$ corrections to the force and torque. When commenting on the corresponding results, we shall frequently refer to the original BBO equation (1.1). In the non-dimensional variables defined above, this equation reads
\begin{equation} {\boldsymbol{f}}_{BBO}= {\boldsymbol{f}}_b- 6 {\rm \pi}{\boldsymbol{u}}_s - 6 {\rm \pi}\varepsilon\int_0^t \frac{1}{\sqrt{{\rm \pi} (t-\tau)}} \frac{\mbox{d}{\boldsymbol{u}}_s}{{\rm d}\tau}\, \mbox{d}\tau - \frac{2{\rm \pi}}{3} \varepsilon^2 \frac{\mbox{d} {\boldsymbol{u}}_s}{\mbox{d}t} . \end{equation}
The generalized buoyancy force 
${\boldsymbol {f}}_b$ usually comprises an 
$O(1)$ contribution due to gravity/buoyancy, complemented with contributions due to the Lagrangian acceleration of the undisturbed flow. Noting that 
${\mbox {d} {\boldsymbol {U}}^\infty }/{\text {d} t}|_{{\boldsymbol {x}}_p}\equiv {\boldsymbol {v}}_p\boldsymbol{\cdot }{\boldsymbol {\nabla }}{\boldsymbol {U}}^\infty +\dot {{\boldsymbol {U}}}_0$, with 
$\dot {{\boldsymbol {U}}}_0$ the time derivative of the possible uniform component 
${\boldsymbol {U}}_0(t)$ of 
${\boldsymbol {U}}_\infty$ evaluated in the laboratory frame, one has 
${\mbox {D} {\boldsymbol {U}}^\infty }/{\mbox {D} t} |_{{\boldsymbol {x}}_p}=\dot {{\boldsymbol {U}}}_0+({\boldsymbol {v}}_p-{\boldsymbol {u}}_s)\boldsymbol{\cdot }{\boldsymbol {\nabla }} {\boldsymbol {U}}^\infty$. Since the particle velocity may be much larger than the slip velocity, one concludes that the 
${\mbox {d} {\boldsymbol {U}}^\infty }/{\text {d} t}|_{{\boldsymbol {x}}_p}$ term in the Lagrangian acceleration may contribute up to 
$O(1)$ to the total force, whereas the contribution proportional to 
$-{\boldsymbol {u}}_s\boldsymbol{\cdot }\boldsymbol {\nabla } {\boldsymbol {U}}^\infty$ is of 
$O(\varepsilon ^2)$. Later, we also compare the predictions for the second-order force with their counterpart in the inviscid limit. In non-dimensional form, (1.3) becomes 
${\boldsymbol {f}}_{AHP} = {\boldsymbol {f}}_b+\varepsilon ^2{\boldsymbol {f}}'_{inv}$, where the disturbance-induced force 
${\boldsymbol {f}}'_{inv}$ reads
\begin{equation} {\boldsymbol{f}}'_{inv} ={-} \frac{2{\rm \pi}}{3} \left\{ \left( \left.\frac{\mbox{d} {\boldsymbol{u}}_s}{\mbox{d}t} + {\boldsymbol{u}}_s\boldsymbol{\cdot}\boldsymbol{\nabla} {\boldsymbol{U}}^\infty\right|_{{\boldsymbol{x}}_p} \right)+ {\boldsymbol{u}}_s \times ({\boldsymbol{\nabla}} \times {\boldsymbol{U}}^\infty)\right\} . \end{equation}
In contrast to (2.8), no contribution arises in the disturbance-induced force at 
$O(\varepsilon ^0)$ and 
$O(\varepsilon )$ in the inviscid limit. This is because no vorticity is generated at the surface of the sphere, leading in particular to zero viscous drag (D'Alembert paradox). Keeping in mind that the Saffman lift force is an 
$O(\varepsilon )$ effect in the regime of interest here, the fact that no term exists at this order in (2.9) implies that the inviscid shear lift force (last term on the right-hand side) is one order of magnitude smaller than the dominant lift force present in the low- but finite-Reynolds-number regime.
3. Method and solution
We use matched asymptotic expansions in the spirit of Childress (Reference Childress1964) and Saffman (Reference Saffman1965) to approximate the solution of (2.7) up to 
$O(\varepsilon ^2)$. The 
$O(\varepsilon )$ terms listed below were computed earlier by CMM. Most of the technical steps described in § 3.1 also follow closely the approach developed in that reference.
3.1. Outer problem
Far from the particle, the disturbance-flow equations (2.7) in the outer region simplify thanks to Saffman's assumption, 
${Re}_p \ll \varepsilon \ll 1$. Indeed, the magnitude of the disturbance velocity 
$w$ scales as 
$1/r$ for large 
$r$. In the matching region 
$r\sim \varepsilon ^{-1}$, one then has the following estimates:
\begin{align} -{\boldsymbol{u}}_s \boldsymbol{\cdot}{\boldsymbol{\nabla}} {\boldsymbol{w}} \sim O(\varepsilon^2), \quad {\boldsymbol{w}}\boldsymbol{\cdot}{\boldsymbol{\nabla}} {\boldsymbol{w}} \sim O(\varepsilon^3) ,\quad \frac{ \partial {\boldsymbol{w}}}{\partial t} \sim {\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{w}} + ({\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{r}}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}} \sim O(\varepsilon).\end{align}
Therefore, the first two terms can be dropped in (2.7b) for 
$r \sim \varepsilon ^{-1}$ and beyond, yielding the leading-order momentum equation in the outer region
\begin{align} \varepsilon^2 \left( \frac{ \partial {\boldsymbol{w}}_{{out}}}{\partial t} + {\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{w}}_{{out}} + ({\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{r}}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{out}} \right) ={-} {\boldsymbol{\nabla}} {p}_{{out}} + \boldsymbol{\triangle}{{\boldsymbol{w}}_{{out}}} + ({\boldsymbol{f}}^{(0)} + \varepsilon {\boldsymbol{f}}^{(1)})\delta({\boldsymbol{r}}).\end{align}
In the standard fashion pioneered by Childress (Reference Childress1964), we account for the particle in the outer region through a delta function, 
$\delta (\boldsymbol{r})$. Usually, the strength of the corresponding force, 
${\boldsymbol {f}}^{(0)}$, is taken to be the opposite of the disturbance-induced force experienced by the particle in the creeping-flow limit, 
${\boldsymbol {f}}'^{(0)}$. Here, however, we must allow for corrections of 
$O(\varepsilon )$ since we wish to compute the outer contribution to the force at 
$O(\varepsilon ^2)$. This is why there is a contribution 
$\varepsilon {\boldsymbol {f}}^{(1)}$ on the right-hand side of (3.2). Strictly speaking, the right-hand side of (3.2) should also comprise an additional source term in the form of a dipolar contribution, 
${\mathbb {D}}^{(0)}\boldsymbol{\cdot } \boldsymbol {\nabla } \delta ({\boldsymbol {r}})$, resulting from the strain- and rotation-induced terms in the boundary condition (2.7c). This dipolar term adds a linear correction in the far-field flow, which may be computed after the components of the second-rank tensor 
${\mathbb {D}}^{(0)}$ have been evaluated by matching the leading-order inner and outer flows in the intermediate region 
$r\sim \varepsilon ^{-1}$. This matching procedure yields 
${\mathbb {D}}^{(0)}\boldsymbol{\cdot } \boldsymbol {\nabla } \delta = ({20 {\rm \pi}}/{3}) {\mathbb {S}}\boldsymbol{\cdot } \boldsymbol {\nabla } \delta - 4{\rm \pi} \boldsymbol {\omega }_s\times \boldsymbol {\nabla } \delta$, which corresponds to the stresslet and torque exerted by the particle on the fluid (Batchelor Reference Batchelor1970). Nevertheless, the corresponding far-field correction does not change the force on the particle at any order. Moreover, it may be shown that it only modifies the torque at 
$O(\varepsilon ^3)$ (Meibohm et al. Reference Meibohm, Candelier, Rosen, Einarsson, Lundell and Mehlig2016). For these reasons, we ignore this dipolar term in what follows. We also need to examine whether or not the Oseen-like contribution resulting from the first-order disturbance in the far field, 
$-{\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {\nabla }}{\boldsymbol {w}}_{{out}}^{(1)}$, has to be included in (3.2) to evaluate the second-order correction, 
${\boldsymbol {w}}_{{out}}^{(2)}$. For reasons detailed below, it turns out that this additional forcing term does not add any singular correction to the far-field flow. Therefore, it is sufficient to consider the effect of Oseen-like contributions upon the inner solution to evaluate the modification they introduce on the loads experienced by the particle. This is why no Oseen-like term is included in (3.2). We shall return to this point later in this section.
Equation (3.2) can be solved via Fourier transform (see CMM). In Fourier space, once the pressure has been eliminated with the aid of the continuity equation, the transformed outer momentum equation reads
\begin{align} &\varepsilon^2 \left(
\left.\frac{ \partial
\hat{{\boldsymbol{w}}}_{{out}}}{\partial
t}\right|_{{\boldsymbol{k}}} +
{\mathbb{A}}\boldsymbol{\cdot}
\hat{{\boldsymbol{w}}}_{{out}} - {\boldsymbol{k}}
\boldsymbol{\cdot} ({\mathbb{A}}\boldsymbol{\cdot}
\widehat{{\boldsymbol{\nabla}}})
\hat{{\boldsymbol{w}}}_{{out}} -2 \frac{
({\mathbb{A}}\boldsymbol{\cdot}
\hat{{\boldsymbol{w}}}_{{out}}) \boldsymbol{\cdot}
{\boldsymbol{k}} }{k^2} {\boldsymbol{k}} \right)
\nonumber\\ &\quad ={-} k^2 \hat{{\boldsymbol{w}}}_{{out}}
+ k^2 \hat{{\mathbb{G}}} \boldsymbol{\cdot}
({\boldsymbol{f}}^{(0)} + \varepsilon
{\boldsymbol{f}}^{(1)}), \end{align}
where
\begin{equation} \hat{{\mathbb{G}}} = \frac{1}{k^2}\left({\mathbb{I}} - \frac{{\boldsymbol{k}} \otimes {\boldsymbol{k}}}{k^2}\right), \end{equation}
is the Fourier transform of the Green function 
${\mathbb {G}}({\boldsymbol {r}})$ of the Stokes equation, with 
${\mathbb {I}}$ the unit matrix. Following CMM, (3.3) is expressed in a non-orthogonal coordinate system that moves and deforms with the background flow. This transformation allows the problem to be reduced to a set of ordinary differential equations that are much easier to solve. However, this simplification only holds as long as 
${\mathbb {A}}$ does not depend upon time. Indeed, if 
${\mathbb {A}}$ is time dependent, an extra term arises in the transformed momentum equation, complicating the structure of the general solution and making it much more difficult to obtain. This is why the outer solution described below is only valid when the linear flow is stationary.
Once 
$\hat {{\boldsymbol {w}}}_{{out}}$ is obtained, it may be expanded in terms of generalized functions of 
${\boldsymbol {k}}$ in the form (Meibohm et al. Reference Meibohm, Candelier, Rosen, Einarsson, Lundell and Mehlig2016)
\begin{equation} \hat{{{\boldsymbol{w}}}}_{{out}} = \hat{{\boldsymbol{\mathcal{T}}}}_{{reg}}^{(0)} + \varepsilon( \hat{{\boldsymbol{\mathcal{T}}}}_{{reg}}^{(1)} + \hat{{\boldsymbol{\mathcal{T}}}}_{{sing}}^{(1)}) + \varepsilon^2 (\hat{{\boldsymbol{\mathcal{T}}}}_{{reg}}^{(2)} + \hat{{\boldsymbol{\mathcal{T}}}}_{{sing}}^{(2)})+\cdots . \end{equation}
Inserting this ansatz into (3.3) and expanding in powers of 
$\varepsilon$ provides the regular contributions 
$\hat {{\boldsymbol {\mathcal {T}}}}_{{reg}}^{(n)}$, namely
$$\begin{gather} \hat{{\boldsymbol{\mathcal{T}}}}_{{reg}}^{(0)} = \hat{{\mathbb{G}}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(0)}, \end{gather}$$
$$\begin{gather}\hat{{\boldsymbol{\mathcal{T}}}}_{{reg}}^{(1)} = \hat{{\mathbb{G}}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(1)}, \end{gather}$$
\begin{gather}
\hat{{\boldsymbol{\mathcal{T}}}}_{{reg}}^{(2)}
={-}\frac{1}{k^2}\hat{{\mathbb{G}}}\boldsymbol{\cdot}
\frac{\mbox{d}{\boldsymbol{f}}^{(0)}}{\mbox{d}t} -
\frac{1}{k^4} {\mathbb{A}}\boldsymbol{\cdot}
{\boldsymbol{f}}^{(0)} \nonumber\\ \quad -
\frac{1}{k^6}\left\{2[({\mathbb{A}}\boldsymbol{\cdot}
{\boldsymbol{k}}) \boldsymbol{\cdot} {\boldsymbol{k}}]
{\mathbb{I}} - ({\mathbb{A}}-{\mathbb{A}}^{\sf T})
\boldsymbol{\cdot} {\boldsymbol{k}}\otimes{\boldsymbol{k}}
- {\boldsymbol{k}}\otimes{\boldsymbol{k}}
\boldsymbol{\cdot} {\mathbb{A}} \right\}\boldsymbol{\cdot}
{\boldsymbol{f}}^{(0)} \nonumber\\ \quad + \frac{2}{k^8}
[({\mathbb{A}} \boldsymbol{\cdot} {\boldsymbol{k}})
\boldsymbol{\cdot} {\boldsymbol{k}} ]
{\boldsymbol{k}}\otimes{\boldsymbol{k}} \boldsymbol{\cdot}
{\boldsymbol{f}}^{(0)} .
\end{gather}
Transforming these contributions from Fourier space back to the configuration space yields
$$\begin{gather} {{\boldsymbol{\mathcal{T}}}}_{{reg}}^{(0)} = {{\mathbb{G}}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(0)}, \end{gather}$$
$$\begin{gather}{{\boldsymbol{\mathcal{T}}}}_{{reg}}^{(1)} = {{\mathbb{G}}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(1)}, \end{gather}$$
\begin{gather}
{{\boldsymbol{\mathcal{T}}}}_{{reg}}^{(2)} = \frac{3
r}{32{\rm \pi}} \left( {\mathbb{I}} - \frac{1}{3}
\frac{{\boldsymbol{r}} \otimes
{\boldsymbol{r}}}{r^2}\right) \boldsymbol{\cdot}
\frac{\mbox{d} {\boldsymbol{f}}^{(0)}}{\mbox{d}t}
\nonumber\\ \quad + \frac{r }{32{\rm \pi}} \left( 3 \left(
{\mathbb{I}} - \frac{1}{3} \frac{{\boldsymbol{r}} \otimes
{\boldsymbol{r}}}{r^2}\right) \boldsymbol{\cdot}
{\mathbb{A}} - ({\mathbb{A}}-{\mathbb{A}}^{\sf T})
\boldsymbol{\cdot}\left( {\mathbb{I}} +
\frac{{\boldsymbol{r}} \otimes {\boldsymbol{r}}}{r^2}
\right) + 2\frac{{\mathbb{A}} \boldsymbol{:} {\boldsymbol{r}}
\otimes {\boldsymbol{r}}}{r^2} {\mathbb{I}}\right)
\boldsymbol{\cdot} {\boldsymbol{f}}^{(0)} \nonumber\\ \quad
- \frac{r}{96{\rm \pi} } \left( 2{\mathbb{S}} + 2
\frac{{\boldsymbol{r}} \otimes {\boldsymbol{r}}}{r^2}
\boldsymbol{\cdot} {\mathbb{S}} + 2 {\mathbb{S}}
\boldsymbol{\cdot} \frac{{\boldsymbol{r}} \otimes
{\boldsymbol{r}}}{r^2} + \frac{{\mathbb{A}} \boldsymbol{:}
{\boldsymbol{r}} \otimes {\boldsymbol{r}}}{r^2}
\left({\mathbb{I}} - \frac{{\boldsymbol{r}} \otimes
{\boldsymbol{r}}}{r^2} \right) \right)
\boldsymbol{\cdot}{\boldsymbol{f}}^{(0)},
\end{gather}
with
\begin{equation} {\mathbb{G}}({\boldsymbol{r}})= \frac{1}{8{\rm \pi}} \left(\frac{{\mathbb{I}}}{r} + \frac{{\boldsymbol{r}}\otimes{\boldsymbol{r}}}{r^3}\right). \end{equation}
The expansion (3.5) suggests that there are singular terms in 
$\boldsymbol{k}$ space, 
$\hat {{\boldsymbol {\mathcal {T}}}}_{{sing}}^{(n)}$, that cannot be obtained in this way. These terms, which are concentrated at 
${\boldsymbol {k}}={\boldsymbol {0}}$, correspond to uniform flow corrections in the far field, resulting from the presence of the particle. As is well known, these corrections are directly related to the 
$1/r$ divergence of the advective contribution 
${\mathbb {A}}\boldsymbol{\cdot } {\boldsymbol {w}} + ({\mathbb {A}}\boldsymbol{\cdot } {\boldsymbol {r}}) \boldsymbol{\cdot } {\boldsymbol {\nabla }} {\boldsymbol {w}}$ based on the leading-order solution. These singular terms may be computed by evaluating in the sense of generalized functions the successive limits (Meibohm et al. Reference Meibohm, Candelier, Rosen, Einarsson, Lundell and Mehlig2016)
\begin{equation}
\hat{\boldsymbol{\mathcal{T}}}_{{sing}}^{(1)} =
\lim_{\varepsilon \to 0}
\frac{\hat{{{\boldsymbol{w}}}}_{{out}} -
\hat{\boldsymbol{\mathcal{T}}}_{{reg}}^{(0)}}{\varepsilon} -
\hat{\boldsymbol{\mathcal{T}}}_{{reg}}^{(1)},
\end{equation}
and
\begin{equation}
\hat{\boldsymbol{\mathcal{T}}}_{{sing}}^{(2)} = \lim_{\varepsilon \to 0}
\frac{\hat{{{\boldsymbol{w}}}}_{{out}} -
\hat{\boldsymbol{\mathcal{T}}}_{{reg}}^{(0)} - \varepsilon
\hat{\boldsymbol{\mathcal{T}}}_{{reg}}^{(1)} -\varepsilon
\hat{\boldsymbol{\mathcal{T}}}_{{sing}}^{(1)}}{\varepsilon^2} -
\hat{\boldsymbol{\mathcal{T}}}_{{reg}}^{(2)} .
\end{equation}
The outer solution is obtained from (3.3), allowing the successive singular contributions to be evaluated. The term 
$\hat {{\boldsymbol {\mathcal {T}}}}_{{sing}}^{(1)}$ was calculated in this way by CMM who showed that it takes the form of a convolution product between a tensorial kernel, 
${\mathbb {K}}(t)$, and the time derivative of the leading-order forcing term in (3.3), 
${\boldsymbol {f}}^{(0)}$. That is,
\begin{equation} \hat{{\boldsymbol{\mathcal{T}}}}_{{sing}}^{(1)} ={-} 8 {\rm \pi}^3\left[ \int_0^{t} {\mathbb{K}}(t-\tau) \boldsymbol{\cdot} \frac{\mbox{d} {\boldsymbol{f}}^{(0)}}{\mbox{d}\tau}\, \mbox{d} \tau \right] \delta({\boldsymbol{k}}) \equiv 8 {\rm \pi}^3 {\boldsymbol{\mathcal{U}}}_1(t) \delta({\boldsymbol{k}}) . \end{equation}
The 
$[\mathbb{K]}_i^j(t)$ component of the kernel expresses how an instantaneous change in the 
$i$th component of the slip velocity (hence, in that of the force exerted by the particle on the fluid) influences the 
$j$th component of the uniform velocity correction in the far field at later time. The characteristic time scale over which this influence manifests itself is that imposed by the shear, which is larger than the viscous time scale by a factor of 
$O(\varepsilon ^{-2})$. Since 
$\hat {{\boldsymbol {\mathcal {T}}}}_{{sing}}^{(1)}$ contributes to 
$\hat {{\boldsymbol {w}}}_{{out}}$ at 
$O(\varepsilon )$ according to (3.5), it may be concluded that time variations in the slip velocity (hence, in 
${\boldsymbol {f}}^{(0)}$) already affect the force on the particle at 
$O(\varepsilon )$ through the far-field velocity correction 
${\boldsymbol {\mathcal {U}}}_1(t)$. Conversely, as (3.7c) shows, these time variations only affect the inner solution at 
$O(\varepsilon ^2)$. Candelier et al. (Reference Candelier, Mehlig and Magnaudet2019) computed the kernel 
$\mathbb{K}(t)$ and discussed its properties for the three canonical planar linear flows. At very short time, comparable with the viscous time scale (
$t\sim \varepsilon ^2$), or in the absence of fluid-velocity gradients, 
$\mathbb{K}$ reduces to the Basset–Boussinesq kernel, i.e. 
$\mathbb{K}(t)\rightarrow ({\rm \pi} t)^{-1/2}{\mathbb {I}}$. The second-order singular term takes the form
\begin{equation} \hat{{\boldsymbol{\mathcal{T}}}}_{{sing}}^{(2)} ={-}8 {\rm \pi}^3 \left[ \int_0^{t} {\mathbb{K}}(t-\tau) \boldsymbol{\cdot} \frac{\mbox{d} {\boldsymbol{f}}^{(1)}}{\mbox{d}\tau}\, \mbox{d} \tau \right] \delta({\boldsymbol{k}}) \equiv 8 {\rm \pi}^3 {\boldsymbol{\mathcal{U}}}_2(t) \delta({\boldsymbol{k}}) . \end{equation}
Remarkably, this term involves the same kernel, 
$\mathbb{K}(t)$, as 
$\hat {{\boldsymbol {\mathcal {T}}}}_{{sing}}^{(1)}$. Since 
${\boldsymbol {f}}^{(1)}(t)=-6{\rm \pi} {\boldsymbol {\mathcal {U}}}_1(t)$ results from the convolution product between the kernel 
${\mathbb {K}}(t)$ and the time derivative of the force 
${\boldsymbol {f}}^{(0)}$, 
$\hat {{\boldsymbol {\mathcal {T}}}}_{{sing}}^{(2)}$ has the form of a double convolution product. In the configuration space, the above two singular terms give rise to two spatially uniform but time-dependent velocity corrections, namely
\begin{equation} {{\boldsymbol{\mathcal{T}}}}_{{sing}}^{(1)} = {\boldsymbol{\mathcal{U}}}_1(t) \quad \mbox{and} \quad {{\boldsymbol{\mathcal{T}}}}_{{sing}}^{(2)} = {\boldsymbol{\mathcal{U}}}_2(t) . \end{equation}
Equations (3.7) and (3.13a,b) are the desired solutions of the outer problem. The velocity field resulting from the superposition of these solutions serves as the outer boundary condition for the inner problem, since the inner and outer solutions must match at 
$r \sim \varepsilon ^{-1}$. Note that we also evaluated 
${{\boldsymbol {\mathcal {T}}}}_{{sing}}^{(2)}$ with the Oseen-like term included in (3.2). The far-field disturbance resulting from this term comprises odd and even functions of 
${\boldsymbol {k}}$, but the former is one order lower with respect to 
$\varepsilon$ than the latter, owing to the assumption that the Saffman length is shorter than the Oseen length by a factor of 
$O(\varepsilon )$ (see appendix A in CMM for details on the evaluation of the contributions to 
${\boldsymbol {w}}_{{out}}$ in 
${\boldsymbol {k}}$ space). Conversely, and for the same reason, the even part is one order lower with respect to 
$\varepsilon$ than the odd one in the case of the disturbance associated with the linear flow. Since odd terms eventually integrate to zero, it turns out that, at the order of approximation we need to consider, the Oseen-like term does not contribute to 
${{\boldsymbol {\mathcal {T}}}}_{{sing}}^{(2)}$ whereas terms associated with the linear flow do. In contrast, the Oseen-like term produces contributions in the regular term 
${{\boldsymbol {\mathcal {T}}}}_{{reg}}^{(2)}$. However, the evaluation of 
${{\boldsymbol {\mathcal {T}}}}_{{reg}}^{(2)}$ is not needed, since this term merely matches the second-order inner solution for 
$r\sim \varepsilon ^{-1}$. Therefore, only the contributions induced in this inner solution by the Oseen-like term 
$-{\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {\nabla }}{\boldsymbol {w}}_{{in}}^{(1)}$ and the companion term 
${\boldsymbol {w}}_{{in}}^{(1)}\boldsymbol{\cdot }{\boldsymbol {\nabla }}{\boldsymbol {w}}_{{in}}^{(1)}$ need to be considered (see below) to evaluate the loads on the body at 
$O(\varepsilon ^2)$.
3.2. Inner problem
At order 
$\varepsilon ^0$, the problem to solve is
$$\begin{gather} {\boldsymbol{\nabla}} \boldsymbol{\cdot} {{\boldsymbol{w}}}_{{in}}^{(0)}= {\boldsymbol{0}}, \quad - {\boldsymbol{\nabla}} {p}_{{in}}^{(0)} + \boldsymbol{\triangle} {{\boldsymbol{w}}}_{{in}}^{(0)} = {\boldsymbol{0}}, \end{gather}$$
$$\begin{gather}{{\boldsymbol{w}}}_{{in}}^{(0)} = {\boldsymbol{u}}_s +\boldsymbol{\omega}_s\times {\boldsymbol{r}} - {\mathbb{S}}\boldsymbol{\cdot} {\boldsymbol{r}}\ \mbox{at}\ r=1, \quad {{\boldsymbol{w}}}_{{in}}^{(0)} \sim {\boldsymbol{\mathcal{T}}}^{(0)} = {\mathbb{G}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(0)} \quad\mbox{for}\ r \to \infty. \end{gather}$$In a linear flow, the solution of this standard Stokes problem is known to be (Kim & Karrila Reference Kim and Karrila1991)
\begin{equation} {{\boldsymbol{w}}}_{{in}}^{(0)} = {\mathbb{G}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(0)} + \frac{1}{6}{\boldsymbol{\triangle}} ( {\mathbb{G}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(0)} ) + \frac{\boldsymbol{\omega}_s\times {\boldsymbol{r}}}{r^3} -\frac{{\mathbb{S}}\boldsymbol{\cdot} {\boldsymbol{r}}}{r^5}-\frac{5}{2 r^5}\left(1-\frac{1}{r^2}\right) {\boldsymbol{r}}\otimes{\boldsymbol{r}}\boldsymbol{\cdot} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{r}} . \end{equation}
Only the first term on the right-hand side contributes to the disturbance-induced force 
${\boldsymbol {f}}'^{(0)}=-{\boldsymbol {f}}^{(0)}$ on the particle, which is of course nothing but the Stokes drag corresponding to the slip velocity 
${\boldsymbol {u}}_s$, namely
\begin{equation} {\boldsymbol{f}}'^{(0)}(t) ={-}6{\rm \pi} {\boldsymbol{u}}_s(t) .\end{equation}Similarly, only the rotlet (third term on the right-hand side of (3.16)) contributes to the torque for a sphere, leading to
\begin{equation} {\boldsymbol{\tau}}'^{(0)}(t) ={-}8{\rm \pi} \boldsymbol{\omega}_s(t).\end{equation}
At order 
$\varepsilon$, the problem to solve is
\begin{equation} {\boldsymbol{\nabla}} \boldsymbol{\cdot} {{\boldsymbol{w}}}_{{in}}^{(1)}= {\boldsymbol{0}} , \quad - {\boldsymbol{\nabla}} {p}_{{in}}^{(1)} + \boldsymbol{\triangle} {{\boldsymbol{w}}}_{{in}}^{(1)} = {\boldsymbol{0}}, \end{equation}subject to the boundary conditions
\begin{equation} {{\boldsymbol{w}}}_{{in}}^{(1)} = {\boldsymbol{0}} \quad \mbox{at}\ r=1, \quad {{\boldsymbol{w}}}_{{in}}^{(1)} \sim {\boldsymbol{\mathcal{U}}}^{(1)}(t)+{\boldsymbol{\mathcal{T}}}^{(1)}_{{reg}} \quad \mbox{for}\ r \to \infty . \end{equation}This is a homogeneous Stokes problem. Its solution may be sought in the form
\begin{equation} {{\boldsymbol{w}}}_{{in}}^{(1)} = {\boldsymbol{\mathcal{U}}}^{(1)}(t) + {{\boldsymbol{v}}}_{{in}}^{(1)} . \end{equation}
Since 
${\boldsymbol {\mathcal {U}}}^{(1)}(t)$ is uniform (see (3.13a,b)), the problem for 
${{\boldsymbol {v}}}_{{in}}^{(1)}$ becomes
$$\begin{gather} {\boldsymbol{\nabla}} \boldsymbol{\cdot} {{\boldsymbol{v}}}_{{in}}^{(1)}= {\boldsymbol{0}} , \quad - {\boldsymbol{\nabla}} {p}_{{in}}^{(1)} + \boldsymbol{\triangle} {{\boldsymbol{v}}}_{{in}}^{(1)} = {\boldsymbol{0}}, \end{gather}$$
$$\begin{gather} {{\boldsymbol{v}}}_{{in}}^{(1)} ={-}{\boldsymbol{\mathcal{U}}}^{(1)}(t) \quad \mbox{at}\ r=1, \quad {{\boldsymbol{v}}}_{{in}}^{(1)} \sim {\mathbb{G}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(1)} \quad \mbox{for}\ r\to \infty . \end{gather}$$
This problem is formally identical to the Stokes problem for a sphere moving in a fluid at rest. Therefore, the solution for 
${{\boldsymbol {w}}}_{{in}}^{(1)}$ is readily obtained as
\begin{equation} {{\boldsymbol{w}}}_{{in}}^{(1)} ={\boldsymbol{\mathcal{U}}}^{(1)}(t) + {\mathbb{G}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(1)} + \tfrac{1}{6}{\boldsymbol{\triangle}} ({\mathbb{G}} \boldsymbol{\cdot} {\boldsymbol{f}}^{(1)}) . \end{equation}
This implies that the contribution 
${\boldsymbol {f}}'^{(1)}=-{\boldsymbol {f}}^{(1)}$ to the force on the particle is the Stokes drag corresponding to the slip velocity 
$-{\boldsymbol {\mathcal {U}}}^{(1)}$, namely
\begin{equation} {\boldsymbol{f}}'^{(1)}(t) = 6{\rm \pi} {\boldsymbol{\mathcal{U}}}^{(1)}(t) ={-}6 {\rm \pi}\int_0^t {\mathbb{K}}(t-\tau) \boldsymbol{\cdot} \frac{\mbox{d} {\boldsymbol{f}}^{(0)}}{\mbox{d}\tau}\, \mbox{d}\tau, \end{equation}
where the second equality stems from (3.11). Since 
${\boldsymbol {f}}'^{(1)}$ is generally not collinear to 
${\boldsymbol {f}}'^{(0)}$, the possibly non-zero component 
${\boldsymbol {f}}'^{(1)}-({{\boldsymbol {f}}'^{(1)}\boldsymbol{\cdot }{\boldsymbol {f}}'^{(0)}}/{{\boldsymbol {f}}'^{(0)} \boldsymbol{\cdot }{\boldsymbol {f}}'^{(0)}}){\boldsymbol {f}}'^{(0)}$ represents a lift force. In the case of a sphere, and more generally of a body with a fore-aft symmetry, this is even the leading-order lift contribution, since symmetry considerations indicate that no lift force can exist at 
$O(\varepsilon ^0)$ for such bodies (Bretherton Reference Bretherton1962).
For an arbitrarily shaped particle, the 
$O(\varepsilon )$ torque comprises a contribution proportional to the right-hand side of (3.25). However, 
${\boldsymbol {\mathcal {U}}}^{(1)}$ is a uniform velocity field. Therefore, just as 
${\boldsymbol {u}}_s$ cannot induce any torque on a spherical particle in the creeping-flow limit, no 
$O(\varepsilon )$ torque may result from 
${\boldsymbol {\mathcal {U}}}^{(1)}$, implying that
\begin{equation} {\boldsymbol{\tau}}'^{(1)}(t) = \boldsymbol{0}.\end{equation}Obviously, this conclusion would not hold if the shape of the particle were such that its translational and rotational dynamics are coupled in the creeping-flow limit.
The 
$O(\varepsilon ^2)$ problem is more complicated, owing to the 
$O(\varepsilon ^2)$ terms on the left-hand side of (2.7b). More specifically, one has to solve
$$\begin{gather} {\boldsymbol{\nabla}} \boldsymbol{\cdot} {{\boldsymbol{w}}}_{{in}}^{(2)}= {\boldsymbol{0}} , \end{gather}$$
$$\begin{gather}\frac{ \partial {\boldsymbol{w}}_{{ in}}^{(0)}}{\partial t} + {\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{w}}_{{in}}^{(0)} + ({\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{r}}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{in}}^{(0)} - {\boldsymbol{u}}_s \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{in}}^{(0)} + {\boldsymbol{w}}_{{in}}^{(0)} \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{in}}^{(0)} ={-} {\boldsymbol{\nabla}} {p}_{{in}}^{(2)} + \boldsymbol{\triangle} {{\boldsymbol{w}}_{{in}}^{(2)}}, \end{gather}$$together with the boundary conditions
\begin{equation} {{\boldsymbol{w}}}_{{in}}^{(2)} = {\boldsymbol{0}} \quad \mbox{at}\ r=1, \quad {{\boldsymbol{w}}}_{{in}}^{(2)} \sim {\boldsymbol{\mathcal{U}}}^{(2)}(t)+{\boldsymbol{\mathcal{T}}}^{(2)}_{{reg}} \quad \mbox{for} \ r \to \infty . \end{equation}
The above problem is inhomogeneous. We seek its solution in the form of a particular (forced) solution, to which we add the uniform contribution 
${\boldsymbol {\mathcal {U}}}_2(t)$, plus a solution of the homogeneous problem. In other words, we set
\begin{equation} {p}_{{in}}^{(2)} = p_p + p_h , \quad {{\boldsymbol{w}}_{{in}}^{(2)}} = {\boldsymbol{w}}_{p}+ {\boldsymbol{\mathcal{U}}}_2(t) + {\boldsymbol{w}}_{h} . \end{equation}
We obtain the formal solutions corresponding to 
${\boldsymbol {w}}_{p}$ and 
${\boldsymbol {w}}_{h}$ using Maple
$^{\circledR}$. Here we just outline the main steps of the procedure but do not provide the final expressions since they are extremely lengthy. Nevertheless, the corresponding routines may be obtained on request from the authors. Moreover, a pivotal technical step in the building of these inner solutions turned out to be the generic determination of the Fourier transforms of functions involving negative or positive powers of 
$r$. Since this aspect may be of interest to some readers, we provide the corresponding results in supplementary material available at https://doi.org/10.1017/jfm.2022.1015.
To determine the particular solution 
${\boldsymbol {w}}_{p}$, we first compute the Fourier transform (
${\boldsymbol {\mathcal {F}}}$) of the associated pressure in the form
\begin{align} \hat{p}_p = \frac{\texttt{i}}{k^2} \left({\boldsymbol{k}}\boldsymbol{\cdot} {\boldsymbol{\mathcal{F}}}\left(\frac{\partial {\boldsymbol{w}}_{{in}}^{(0)}}{\partial t} + {\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{w}}_{{in}}^{(0)} + ({\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{r}}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{in}}^{(0)} - {\boldsymbol{u}}_s \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{in}}^{(0)} + {\boldsymbol{w}}_{{in}}^{(0)} \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{in}}^{(0)} \right)\right) , \end{align}
with 
$\texttt {i}^2=-1$. This allows us to obtain the particular solution for the velocity in Fourier space as
\begin{align} \hat{{\boldsymbol{w}}}_p ={-} \hat{{\mathbb{G}}} \boldsymbol{\cdot} {\boldsymbol{\mathcal{F}}}\left( \frac{\partial {\boldsymbol{w}}{{in}}^{(0)}}{\partial t} + {\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{w}}_{{in}}^{(0)} + ({\mathbb{A}}\boldsymbol{\cdot} {\boldsymbol{r}}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{in}}^{(0)} - {\boldsymbol{u}}_s \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{in}}^{(0)} + {\boldsymbol{w}}_{{in}}^{(0)} \boldsymbol{\cdot} {\boldsymbol{\nabla}} {\boldsymbol{w}}_{{in}}^{(0)} \right). \end{align}The inverse Fourier transform then yields the particular solution in the configuration space as
\begin{equation} {\boldsymbol{w}}_p = {\boldsymbol{\mathcal{F}}}^{{-}1} ( \hat{{\boldsymbol{w}}}_p) \quad \mbox{and} \quad p_p = {\boldsymbol{\mathcal{F}}}^{{-}1} ( \hat{p}_p) .\end{equation}
Note that, in line with the comments made in § 3.1, the contribution of the Oseen-like terms resulting from the leading-order solution 
${\boldsymbol {w}}_{{in}}^{(0)}$ is accounted for in (3.28)–(3.33a,b). Next, we consider the homogeneous solution 
$\boldsymbol{w}_h$. It satisfies
\begin{equation} {\boldsymbol{\nabla}} \boldsymbol{\cdot} {\boldsymbol{w}}_{h} = 0 ,\quad - {\boldsymbol{\nabla}} {p}_h + \boldsymbol{\triangle} {\boldsymbol{w}}_{h} = {\boldsymbol{0}} , \end{equation}together with the boundary conditions
\begin{equation} {\boldsymbol{w}}_{h} ={-} {\boldsymbol{w}}_p - {\boldsymbol{\mathcal{U}}}_2(t) \quad \mbox{at}\ r=1,\quad {\boldsymbol{w}}_{h} \sim {\boldsymbol{0}} \quad\mbox{for} \ r \to \infty. \end{equation}
As can be seen, the particular solution computed previously now appears as a boundary condition on the particle surface, together with the uniform velocity correction 
$\boldsymbol{\mathcal {U}}_2(t)$ resulting from the outer solution. We finally obtain the solution of (3.35) using Lamb's expansion (Lamb Reference Lamb1932, art. 336). Note that, unlike the calculation of 
${\boldsymbol {w}}_{{out}}$, the evaluation of the inner solution is left unchanged if 
${\mathbb {A}}$ is time dependent. Therefore, the second-order inner contributions to the force and torque discussed below are valid even if the underlying linear flow is not stationary.
3.3. Second-order force and torque
Using the successive solutions described above, the second-order contributions to the force and torque acting on the particle may be computed from the standard definitions
\begin{equation}
{{\boldsymbol{f}}'}^{(2)} = \int_\mathcal{S}
\boldsymbol{\mathbb{\sigma}}^{(2)}_{{in}} \boldsymbol{\cdot}
{\boldsymbol{n}} \mbox{d} \mathcal{S}, \quad
{\boldsymbol{\tau}'}^{(2)} = \int_\mathcal{S}
{\boldsymbol{r}} \times \boldsymbol{\mathbb{\sigma}}^{(2)}_{{in}}
\boldsymbol{\cdot} {\boldsymbol{n}} \mbox{d} \mathcal{S}.
\end{equation}
Here 
${\boldsymbol {n}}$ is the outward unit normal to the particle, and the second-order stress tensor 
$\mathbb{\sigma} ^{(2)}_{{in}}$ is defined as 
$\boldsymbol{\mathbb{\sigma}}^{(2)}_{{in}} = - p_{{in}}^{(2)} \one + 2 {\mathbb {S}}_{{in}}^{(2)}$, where 
${\mathbb {S}}_{{in}}^{(2)}$ is the symmetric part of the velocity-gradient tensor based on the velocity field 
${\boldsymbol {w}}_{{in}}^{(2)}$. The final result for the second-order force and torque resulting from the inner solution is found to be
$$\begin{gather} {\boldsymbol{f}}'^{(2)}(t) = \frac{16{\rm \pi}}{3} \frac{\mbox{d}{\boldsymbol{u}}_s}{\mbox{d}t} + \frac{37{\rm \pi}}{12} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s + \frac{\rm \pi}{3} (4\boldsymbol{\varOmega}+3\boldsymbol{\omega_s} )\times {\boldsymbol{u}}_s + 6{\rm \pi} {\boldsymbol{\mathcal{U}}}_2(t) , \end{gather}$$
$$\begin{gather}\boldsymbol{\tau}'^{(2)}(t) ={-}\frac{8{\rm \pi}}{3}\frac{\mbox{d}\boldsymbol{\omega}_s }{\mbox{d}t} -\frac{8{\rm \pi}}{3} {\mathbb{S}}\boldsymbol{\cdot} \boldsymbol{\varOmega} + {\rm \pi}{\mathbb{S}}\boldsymbol{\cdot} \boldsymbol{\omega}_s . \end{gather}$$
Equations (3.37) and (3.38) are the main results of the paper. For the reason mentioned above, the expression for the torque and the inner contribution to the force are valid for an arbitrary linear background flow, i.e. an arbitrary combination of a uniform straining motion and a solid-body rotation, both of which may possibly be time dependent. In contrast, for the reason explained in § 3.1, we were only able to obtain the far-field uniform corrections 
${\boldsymbol {\mathcal {U}}}_1(t)$ and 
${\boldsymbol {\mathcal {U}}}_2(t)$ in the presence of a steady linear component of the carrying flow, which restricts the general result for the force on the particle to this subclass of flows. It must also be stressed that, because of the nonlinearity of the outer problem, the kernel 
${\mathbb {K}}$ must be computed separately for each given linear flow. This kernel is already known explicitly for pure shear, solid-body rotation and planar elongation (see CMM).
The quadratic contributions in (3.37) and (3.38) involve all possible combinations of 
${\mathbb {S}}$, 
$\boldsymbol\varOmega$, 
$\boldsymbol{\omega}_s$ and 
${\boldsymbol {u}}_s$ allowed by symmetry constraints, with the exception of 
$\boldsymbol \varOmega \times \boldsymbol {\omega }_s=\boldsymbol \varOmega \times \boldsymbol {\omega }_p$ in (3.38). Note that 
${\boldsymbol {u}}_s$ does not appear in (3.38). This is because the torque is an axial (or pseudo-) vector while 
${\boldsymbol {u}}_s$ is a polar (or true) vector, and no axial vector can be formed by combining quadratically 
${\boldsymbol {u}}_s$ with 
$\boldsymbol \varOmega$, 
$\boldsymbol {\omega }_s$ or 
${\mathbb {S}}$. Note also that 
${\boldsymbol {\mathcal {U}}}_2$ does not contribute to the second-order torque, for reasons identical to those already discussed in connection with 
$\boldsymbol {\tau }'^{(1)}$. The presence of the term 
${\rm \pi} \boldsymbol {\omega }_s\times {\boldsymbol {u}}_s$ in (3.37) is noticeable. This contribution, which represents a lift force, is the extension to linear flows of that obtained by Rubinow & Keller (Reference Rubinow and Keller1961) for a sphere rotating and translating in a fluid at rest. This force may be thought of as the visco-inertial analogue of the inviscid Magnus lift force. It results from the coupling of the translational and rotational velocity dynamics operated in the hydrodynamic force by advective effects, a feature that does not exist in the creeping-flow regime, owing to the geometrical symmetries of the particle.
In § 3 we showed that the 
$O(\varepsilon ^2)$ problem is inhomogeneous, i.e. non-zero terms arise on the left-hand side of (3.28). The solution of the problem is linear with respect to these terms. In order to determine this full solution, one can therefore solve a succession of ‘elementary’ problems, considering first, for instance, only the unsteady term 
$\partial _t {\boldsymbol {w}}^{(0)}_{{in}}$, then only terms 
${\mathbb {A}}\boldsymbol{\cdot }{\boldsymbol {w}}_{{in}}^{(0)} + ({\mathbb {A}}\boldsymbol{\cdot } {\boldsymbol {r}}) \boldsymbol{\cdot } {\boldsymbol {\nabla }} {\boldsymbol {w}}_{{in}}^{(0)}$ and so on. The full solution is obtained by summing up these partial solutions. With this procedure, one can trace back which contribution to the force is due to the pressure gradient, which is of viscous origin, etc. Only viscous stresses contribute to the torque since 
${\boldsymbol {r}}$ and 
${\boldsymbol {n}}$ are collinear for a sphere. The origin of the various contributions to the second-order force and torque is summarized in table 2.
Table 2. Summary of all contributions to the second-order force and torque, the former expressed with respect to 
${\boldsymbol {f}}^{(0)} = 6{\rm \pi} {\boldsymbol {u}}_s$. Terms arising from the regular and singular parts of the solution are given in the upper and lower tables, respectively. We also distinguish between the contributions to the force induced by the viscous part of the stress tensor (e.g. 
${\boldsymbol {f}}_\nu$) and those induced by the pressure (e.g. 
${\boldsymbol {f}}_p$). The total second-order hydrodynamic force is the sum of these two contributions.
4. Discussion
4.1. The 
$O(\varepsilon ^2)$ torque
At leading order, the angular velocity of a torque-free spherical particle immersed in a linear flow is dictated by the vorticity of the undisturbed carrying flow. Indeed, for such a particle, (3.18) and (3.26) imply that 
$\boldsymbol {\omega }_p = \boldsymbol {\varOmega } + O(\varepsilon ^2)$. Now, considering the long-time limit of (3.38), the 
$O(\varepsilon ^2)$ disturbance-induced torque on such a particle is
\begin{equation} \boldsymbol{\tau}'^{(2)} ={-}\frac{8{\rm \pi}}{3} {\mathbb{S}}\boldsymbol{\cdot} \boldsymbol{\varOmega} . \end{equation}
This second-order torque is non-zero only if the base flow is three dimensional and has a non-zero vorticity (which implies that it is unsteady for the reason discussed at the beginning of § 2), since 
${\mathbb {S}}\boldsymbol{\cdot } \boldsymbol {\varOmega }$ is a vortex-stretching term vanishing in a two-dimensional flow. Note that computing the spatial average of 
${\boldsymbol {r}}\times ({\mbox {D} {\boldsymbol {U}}^\infty }/{\text {D} t})$ over the particle volume reveals that the base flow generally brings a complementary contribution to the second-order torque, namely
\begin{equation} \boldsymbol{\tau}_b={-}\frac{8{\rm \pi}}{15}\varepsilon^2 {\mathbb{S}}\boldsymbol{\cdot} \boldsymbol{\varOmega} . \end{equation}
Adding (4.1) and (4.2), one concludes that the total second-order torque resulting from the velocity gradients of a general three-dimensional linear flow is 
$-({16{\rm \pi} }/{5}) {\mathbb {S}}\boldsymbol{\cdot } \boldsymbol {\varOmega }$, as derived by Candelier, Einarsson & Mehlig (Reference Candelier, Einarsson and Mehlig2016) for a neutrally buoyant particle. Since the translational dynamics does not influence the torque for the symmetry reasons discussed above, this result is unchanged if the particle is not neutrally buoyant, as (3.38) indicates. Finally, setting the total torque 
$\boldsymbol {\tau }'^{(0)}+\varepsilon \boldsymbol {\tau }'^{(1)}+\varepsilon ^2 \boldsymbol {\tau }'^{(2)}+\boldsymbol {\tau }_b$ to zero yields the angular velocity of a torque-free particle as
\begin{equation} \boldsymbol{\omega}_p \approx \boldsymbol{\varOmega} - \tfrac{2}{5}\varepsilon^2 {\mathbb{S}}\boldsymbol{\cdot} \boldsymbol{\varOmega} . \end{equation}No history (or memory) contribution appears in the expression of the second-order torque provided by (3.38). This is surprising at first glance, since the problem reduces to the unsteady Stokes equation when unsteady effects are large enough for advective terms to become negligible near the particle. In this limit, Feuillebois & Lasek (Reference Feuillebois and Lasek1978) (see also Gatignol Reference Gatignol1983) computed the torque acting on a spinning sphere, showing that
\begin{align} \boldsymbol{\tau}'(t)={-}8 {\rm \pi}\boldsymbol{\omega}_s(t) - \frac{8{\rm \pi}}{3} \int_0^t \frac{\mbox{d} \boldsymbol{\omega}_s}{\mbox{d}\tau} \left[ \frac{1}{\sqrt{\dfrac{{\rm \pi} (t-\tau)}{\varepsilon^2}}} - \exp\left(\frac{t-\tau}{\varepsilon^2}\right) \mbox{erfc}\left(\sqrt{\frac{t-\tau}{\varepsilon^2}}\right) \right]\,\mbox{d}\tau . \end{align}
The first term on the right-hand side corresponds to the leading-order torque (3.18). To compare the second term in (4.4) with the prediction (3.38), it must be borne in mind that all terms on the left-hand side of (2.2b) are assumed small compared with unity, so that the present theory is valid provided 
${Sl} \ll {Re}_s^{-1} = \varepsilon ^{-2}$. We show in Appendix A that in the limit 
$\varepsilon \rightarrow 0$, the term in square brackets in (4.4) tends toward 
$\varepsilon ^2 \delta (t)$, with 
$\delta (t)$ the delta function (see (A7)). In the limit of small 
$\varepsilon$, the memory term in (4.4) therefore tends towards the first term on the right-hand side of (3.38), namely
\begin{align} - \frac{8{\rm \pi}}{3} \int_0^t \frac{\mbox{d} \boldsymbol{\omega}_s}{\mbox{d}\tau} \left[ \frac{1}{\sqrt{\dfrac{{\rm \pi} (t-\tau)}{\varepsilon^2}}} - \exp\left(\frac{t-\tau}{\varepsilon^2}\right) \mbox{erfc}\left(\sqrt{\frac{t-\tau}{\varepsilon^2}}\right) \right]\,\mbox{d}\tau \to - \frac{8{\rm \pi}}{3} \varepsilon^2\frac{\mbox{d} \boldsymbol{\omega}_s}{\mbox{d}t} . \end{align}
This shows that (3.38) and (4.4) are consistent, and the first contribution on the right-hand side of the former is what is left of the history torque at 
$O(\varepsilon ^2)$ when 
$\varepsilon$ is small. In other words, the memory term derived by Feuillebois & Lasek (Reference Feuillebois and Lasek1978) converges – after a short transient of the order of the viscous time scale – toward the expression obtained here, in which only the instantaneous angular acceleration 
${\mbox {d} \boldsymbol {\omega }_s}/{\mbox {d}t}$ appears. This observation also explains why the theory describing the 
$O({Re}_s)$ effects on the angular dynamics of particles immersed in a stationary shear flow (Einarsson et al. Reference Einarsson, Candelier, Lundell, Angilella and Mehlig2015) does not contain a memory term.
4.2. Second-order force at short times
We now examine the force acting on the particle at short time with respect to the time scale of the background flow, i.e. non-dimensional times such that 
$\varepsilon ^2\ll t\ll 1$. Within this time interval, the vorticity generated at the particle surface at 
$t=0$ has already diffused several radii away from the particle but has not yet entered the wake located at distances 
$r\gtrsim \varepsilon ^{-1}$. At order 
$\varepsilon$, it is known that unsteady inertial effects are responsible for the existence of a history force, the expression of which takes the form of a convolution product between a time-dependent kernel 
${\mathbb {K}}(t)$ and the particle relative acceleration (or more precisely the force 
${\boldsymbol {f}}'^{(0)}$ acting on the particle in the creeping-flow limit). As discussed in § 3.1, a term with a similar structure exists at 
$O(\varepsilon ^2)$, namely 
${\boldsymbol {\mathcal {U}}}_2$, the kernel involved (twice) in the expression of this second-order term being the same as that involved in the 
$O(\varepsilon )$-memory term. Both kernels depend on the specific undisturbed flow under consideration. However, at short times (
$\varepsilon ^2\ll t\ll 1$) and for any linear flow, 
${\mathbb {K}}$ is closely approximated by the Basset–Boussinesq kernel, namely
\begin{equation} 6 {\rm \pi}{\mathbb{K}}(t) \sim \frac{{\mathbb{I}}}{\sqrt{{\rm \pi} t}} . \end{equation}
In this case, 
${\boldsymbol {\mathcal {U}}}_2(t)$ simplifies to
\begin{equation} {\boldsymbol{\mathcal{U}}}_2(t) = \int_0^t \frac{1}{6 {\rm \pi}} \frac{1}{\sqrt{{\rm \pi} (t-\tau)}} \frac{\mbox{d}}{\mbox{d}\tau} \left( \int_0^\tau \frac{1}{\sqrt{{\rm \pi} (\tau-\tau')}} \frac{\mbox{d} {\boldsymbol{f}}^{(0)}}{\mbox{d}\tau'} \right) \,\mbox{d}\tau= \frac{1}{6 {\rm \pi}} \frac{\mbox{d} {\boldsymbol{f}}^{(0)}}{\mbox{d} t } = \frac{\mbox{d} {\boldsymbol{u}}_s}{\mbox{d} t } . \end{equation}
This can be shown using the definition of the fractional derivative 
$\mbox {d}^{1/2}/\mbox {d} t^{1/2}$ or, equivalently, Laplace transform. Remarkably, the double convolution product reduces to a simple ‘local’ term (with respect to time) expressible in the form of a time derivative of the relative velocity. Using (4.7), (3.37) simplifies to
\begin{equation} {\boldsymbol{f}}'^{(2)} ={-} \frac{2{\rm \pi}}{3} \frac{\mbox{d}{\boldsymbol{u}}_s}{\mbox{d}t} + \frac{37{\rm \pi}}{12} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s + \frac{\rm \pi}{3} (4\boldsymbol{\varOmega}+3\boldsymbol{\omega_s} )\times {\boldsymbol{u}}_s . \end{equation}
The first term on the right-hand side corresponds to the added-mass force if the undisturbed flow is uniform; see (2.8). All quadratic terms in (4.8) come from the inner solution. In contrast, the added-mass term results mostly from 
${\boldsymbol {\mathcal {U}}}_2(t)$, i.e. from the second-order outer solution, as its sign differs from that of the corresponding term provided by the inner solution in (3.37). This situation contrasts with that found in (3.38) for the torque, for which the contribution proportional to the time derivative 
${\mbox {{d}}\boldsymbol {\omega }_s}/{\mbox {{d}}t}$ results entirely from the inner solution. Therefore, one has to conclude that the inner and outer regions of the disturbance may both contribute to the ‘remains’ of history effects, the region providing the dominant contribution to the corresponding time-derivative term differing according to the load component under consideration.
Equation (4.8) may be recast in the form
$$\begin{gather} {\boldsymbol{f}}'^{(2)} ={-} \frac{2{\rm \pi}}{3} \left( \left.\frac{{\rm d} {\boldsymbol{v}}_p}{{\rm d}t} - \frac{\mbox{D} {\boldsymbol{U}}^\infty}{\mbox{D} t} \right|_{{\boldsymbol{x}}_p}\right) + \frac{15{\rm \pi}}{4} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s + {\rm \pi}\left( \boldsymbol{\varOmega} + \boldsymbol{\omega}_p \right)\times{\boldsymbol{u}}_s , \end{gather}$$
using again the identity 
$\text {d}{\boldsymbol {U}}^\infty /\text {d}t|_{{\boldsymbol {x}}_p}=\text {D}{\boldsymbol {U}}^\infty /\text {D}t|_{{\boldsymbol {x}}_p}+{\mathbb {A}}\boldsymbol{\cdot } {\boldsymbol {u}}_s$ implying 
$\mbox {d}{\boldsymbol {u}}_s/\mbox {d}t = {\rm d}{\boldsymbol {v}}_p/{\rm d}t - \mbox {D}{\boldsymbol {U}}^\infty /\mbox {D}t|_{{\boldsymbol {x}}_p} - {\mathbb {A}}\boldsymbol{\cdot } {\boldsymbol {u}}_s$. Equation (4.9) for 
$\boldsymbol{f}'^{(2)}$ allows a direct comparison with the inviscid form (2.9) of the disturbance-induced force. In (4.9) the added-mass term (first term on the right-hand side) is identical to that known in the inviscid limit. From (4.9) and (2.9), the difference 
${\boldsymbol {f}}'^{(2)}-{\boldsymbol {f}}'_{inv}$ evaluates to
\begin{equation} {\boldsymbol{f}}'^{(2)}-{\boldsymbol{f}}'_{inv}= \frac{15{\rm \pi}}{4} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s + {\rm \pi}(\boldsymbol{\varOmega}-{\boldsymbol{\nabla}}\times {\boldsymbol{U}}^\infty )\times {\boldsymbol{u}}_s + {\rm \pi}\boldsymbol{\omega}_p\times{\boldsymbol{u}}_s . \end{equation}
Note that in (4.10) the prefactor of the inviscid lift force is 
${\rm \pi}$ instead of 
${2{\rm \pi} }/{3}$ in (2.9). This is because we have to employ the short-time expression for this force to remain consistent with that considered for 
${\boldsymbol {f}}'^{(2)}$. It has been shown (Legendre & Magnaudet Reference Legendre and Magnaudet1998) that in this limit, more precisely for 
$t\ll u_c/as\equiv {Re}_s/{Re}_p\sim 1$, the lift coefficient (usually defined as 
${3}/{4{\rm \pi} }$ times the above prefactor) is 
$\frac {3}{4}$ instead of 
$\frac {1}{2}$ in the steady state, which yields the above result. With this prefactor for the inviscid lift force, the second term on the right-hand side of (4.10) may be simplified as 
$({{\rm \pi} }/{2}) {\boldsymbol {u}}_s \times ({\boldsymbol {\nabla }}\times {\boldsymbol {U}}^\infty )$ using the identity 
$\boldsymbol {\varOmega }=\frac {1}{2}( {\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )$, but the original form appears more suitable for the discussion that follows. Not surprisingly, 
${\boldsymbol {f}}'^{(2)}$ does not coincide with 
${\boldsymbol {f}}'_{inv}$, although the two added-mass terms do. The last contribution on the right-hand side of (4.10) gives an immediate clue to understand where the differences come from. Indeed, the particle rotation does not matter in the inviscid limit, since (i) rotation does not displace any fluid in the specific case of a sphere, and (ii) the fluid is free to slip at the particle surface. Therefore, no force proportional to 
$\boldsymbol {\omega }_p$ can exist in this limit, which is in stark contrast with the visco-inertial regime considered here, in which the no-slip condition forces the surrounding fluid to follow the particle rotation, yielding the non-zero Rubinow–Keller lift force.
The no-slip condition is also responsible for the contribution 
$({15{\rm \pi} }/{4}) {\mathbb {S}} \boldsymbol{\cdot } {\boldsymbol {u}}_s$ that has no counterpart in the inviscid limit, an indication that it results from the influence of the ambient strain on the disturbance generated in the sphere vicinity by the no-slip condition. Note that in general this term may contribute to both the drag and lift components of the total force. The previous argument regarding the role of the no-slip condition holds for the contribution 
$-{\rm \pi} {\boldsymbol {u}}_s \times \boldsymbol {\varOmega }$ but the presence of the inviscid force 
${\rm \pi} {\boldsymbol {u}}_s \times ({\boldsymbol {\nabla }}\times {\boldsymbol {U}}^\infty )$ makes the point a little bit more subtle. As is well known, the inviscid shear lift force results from the distortion of the vorticity 
${\boldsymbol {\nabla }}\times {\boldsymbol {U}}^\infty$ contained in the background flow as it is transported along the curved streamlines around the (inviscid) sphere. Lighthill (Reference Lighthill1956) gave an illuminating quantitative description of how this inviscid stretching/tilting mechanism results in the generation of a pair of counter-rotating streamwise vortices in the wake of a sphere immersed in a linear shear flow. This mechanism subsists qualitatively at a finite Reynolds number. However, when the no-slip condition applies at the sphere surface, the local kinematic structure of the undisturbed flow, characterized by 
${\mathbb {S}}$ and 
$\boldsymbol {\varOmega }$, influences the velocity disturbance required to satisfy this condition; hence, the near-surface vorticity. The term 
$-{\rm \pi} {\boldsymbol {u}}_s \times \boldsymbol {\varOmega }$ in (4.10), or at least a part of it, results from this effect. It may be that, owing to the identity 
$\boldsymbol {\varOmega }=\frac {1}{2}( {\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )$, a non-zero part of this term rather results from what remains at a low but finite Reynolds number of the aforementioned stretching/tilting mechanism of the upstream vorticity, and that the two combine to yield the prefactor 
$-{\rm \pi}$. In any case, there is no reason for the entire finite-Reynolds-number contribution to be identical to the inviscid lift force. Noting that the magnitude of the former is smaller than that of the latter (
$-{{\rm \pi} }/{2}$ instead of 
$-{\rm \pi}$ once expressed with respect to 
${\boldsymbol {u}}_s \times ({\boldsymbol {\nabla }}\times {\boldsymbol {U}}^\infty )$), we may even conclude that the alteration of the near-surface disturbance by the rotational component 
$\boldsymbol {\varOmega }$ of the undisturbed velocity gradient yields a contribution to the force whose sign is opposite to that of the inviscid lift force.
4.3. Stationary limit: pure straining flows
In the stationary limit, the slip velocity between the particle and the fluid no longer varies. Similarly, 
${{\rm d}{\boldsymbol {f}}^{(0)} }/{{\rm d}t}$ tends to zero and the singular term in (3.37) becomes
\begin{equation} 6{\rm \pi} {\boldsymbol{\mathcal{U}}}_2 \to -(6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} {\boldsymbol{f}}^{(0)} , \end{equation}
where 
$6{\rm \pi} \bar {{\mathbb {K}}}$ may be thought of as the steady resistance tensor induced by fluid inertia effects. The second-order force resulting from the flow disturbance then becomes
\begin{equation} {\boldsymbol{f}}'^{(2)} = \frac{37 {\rm \pi}}{12} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s - \frac{\rm \pi}{6} {\boldsymbol{u}}_s \times ( {\boldsymbol{\nabla}} \times {\boldsymbol{U}}^\infty) - {\rm \pi}{\boldsymbol{u}}_s\times\boldsymbol{\omega}_p - 6{\rm \pi} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} {\boldsymbol{u}}_s.\end{equation}Again, this expression is general in that it is valid in any linear flow. The actual difficulty to use it in practice is that the explicit expression for the long-time kernel is generally not known, except in some canonical configurations.
We first specialize the above result to the case of a planar extensional flow, defined as
\begin{equation} {\mathbb{A}} = \left(\begin{array}{@{}ccc@{}} 1 & 0 & 0\\ 0 & -1 & 0 \\ 0 & 0 & 0\\ \end{array}\right) .\end{equation}
The corresponding steady-state kernel was computed by CMM. However, a technical difficulty (a non-convergence of a three-dimensional integral to be evaluated in Fourier space) prevented the computation of the component of the kernel corresponding to the compressional direction 
${\boldsymbol {e}}_2$ beyond a certain time (
$t\approx 30$). The final result obtained in this reference reads
\begin{equation} 6{\rm \pi} \bar{{\mathbb{K}}} \approx \left(\begin{array}{@{}ccc@{}} 0.901 & 0 & 0\\ 0 & [\bar{{\mathbb{K}}}]_2^2 & 0 \\ 0 & 0 & 0.420\\ \end{array} \right) ,\end{equation}
with 
$[\bar {{\mathbb {K}}}]_2^2$ presumably negative and larger than 
$1.48$ in absolute value. Fortunately, the diagonal nature of 
$\bar {{\mathbb {K}}}$ leaves the 
${\boldsymbol {e}}_1$ component of the resulting force unaffected by this uncertainty, even at 
$O(\varepsilon ^2)$. Let us assume that the particle is at rest at the location 
${\boldsymbol {x}}_p=(1,x_{p2},x_{p3})$. This implies that 
${\boldsymbol {u}}_s=-{\boldsymbol {e}}_1+x_{p2}{\boldsymbol {e}}_2$, so that the first-order correction to the 
${\boldsymbol {e}}_1$ component of the disturbance-induced force becomes 
$\boldsymbol{f}'^{(1)}\boldsymbol{\cdot }{\boldsymbol {e}}_1=6{\rm \pi} {\boldsymbol {e}}_1\boldsymbol{\cdot }(6{\rm \pi} \bar {{\mathbb {K}}}) \boldsymbol{\cdot }{\boldsymbol {e}}_1\approx 16.98$. As (3.38) shows, the undisturbed flow does not provide any contribution to the second-order torque because 
$\boldsymbol {\varOmega }\equiv \boldsymbol {0}$. Therefore, a torque-free sphere does not rotate in the present flow whatever its transverse position 
$x_{p2}$, unless it has been given an initial non-zero rotation. If one sets 
$\boldsymbol {\omega }_p=\boldsymbol {0}$ in (3.37), the second-order correction to the 
${\boldsymbol {e}}_1$ component of the disturbance-induced force is found to be
\begin{equation} \boldsymbol{f}'^{(2)}\boldsymbol{\cdot}{\boldsymbol{e}}_1= \frac{37{\rm \pi}}{12} {\boldsymbol{e}}_1\boldsymbol{\cdot}({\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s) + 6{\rm \pi} {\boldsymbol{e}}_1\boldsymbol{\cdot}(6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} (6{\rm \pi} \bar{{\mathbb{K}}})\boldsymbol{\cdot}{\boldsymbol{u}}_s\approx5.615.\end{equation}
Another straining motion of interest is the uniaxial axisymmetric flow characterized by a constant stretching rate along the symmetry axis and a uniform compression in the plane perpendicular to it. This flow was not considered by CMM. The corresponding kernel is computed in Appendix B, for 
${\mathbb {S}}=-({\boldsymbol {e}}_1{\boldsymbol {e}}_1+{\boldsymbol {e}}_2{\boldsymbol {e}}_2)+2{\boldsymbol {e}}_3{\boldsymbol {e}}_3$. The steady-state result (B3) indicates that the axial component of the kernel is 
$6{\rm \pi} [\bar {{\mathbb {K}}}]_3^3=6{\rm \pi} {\boldsymbol {e}}_3\boldsymbol{\cdot }\bar {{\mathbb {K}}} \boldsymbol{\cdot }{\boldsymbol {e}}_3\approx 1.26$. Therefore, considering a particle held fixed on the flow axis at the axial location 
$x_{p3}=\frac {1}{2}$ (which implies 
${u}_s=-\boldsymbol{e}_3$), the first-order correction to the force reads 
$\boldsymbol{f}'^{(1)}= -6{\rm \pi} (6{\rm \pi} \bar {{\mathbb {K}}}) \boldsymbol{\cdot } {\boldsymbol {u}}_s\approx 23.75\, {\boldsymbol {e}}_3$, while the second-order correction is
\begin{equation} \boldsymbol{f}'^{(2)}= \frac{37{\rm \pi}}{12}{\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s -6{\rm \pi} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} {\boldsymbol{u}}_s\approx 10.55 {\boldsymbol{e}}_3 . \end{equation}In Appendix A we establish the counterpart of (4.16) in the case of a spherical bubble at the surface of which the surrounding flow obeys a shear-free condition instead of the no-slip condition in (2.3). The corresponding result reads
\begin{equation} \boldsymbol{f}'^{(2)}= \frac{4{\rm \pi}}{3} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s -4{\rm \pi} (4{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} (4{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} {\boldsymbol{u}}_s\approx0.49 {\boldsymbol{e}}_3.\end{equation}
Comparing (4.16) and (4.17) reveals that the inner and outer contributions to 
$\boldsymbol{f}'^{(2)}$ both differ. In particular, the former is more than twice as small in the case of a bubble. This is a clear confirmation of the prominent role played by the vorticity produced at the particle surface in all contributions to the hydrodynamic force.
The above comments call for a comparison between the predictions for 
$\boldsymbol{f}'^{(2)}$ and those for the inviscid force 
${\boldsymbol {f}}'_{inv}$ given by (2.9). In a pure straining flow, 
${\boldsymbol {f}}'_{inv}$ reduces to the added-mass force, i.e. 
${\boldsymbol {f}}'_{inv}= ({2{\rm \pi} }/{3}) ({\mbox {D} {\boldsymbol {U}}^\infty }/{\mbox {D} t} |_{{\boldsymbol {x}}_p}-({\mbox {d} {\boldsymbol {v}}_p}/{\mbox {d}t}) )$. If the particle is held fixed, 
${\boldsymbol {U}}^\infty |_{{\boldsymbol {x}}_p}=-{\boldsymbol {u}}_s$. Hence, in the case of the planar flow (4.13), setting again 
${\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {e}}_1=-1$, the 
${\boldsymbol {e}}_1$ component of the inviscid force is 
${\boldsymbol {f}}'_{inv}\boldsymbol{\cdot }{\boldsymbol {e}}_1=({2{\rm \pi} }/{3}) {\boldsymbol {e}}_1\boldsymbol{\cdot }{\mathbb {S}}\boldsymbol{\cdot }{\boldsymbol {e}}_1={2{\rm \pi} }/{3}$. Similarly, in the uniaxial straining flow, 
${\boldsymbol {f}}'_{inv}\boldsymbol{\cdot }{\boldsymbol {e}}_3={4{\rm \pi} }/{3}$. Remarkably, although the inner contribution to 
$\boldsymbol {f}'^{(2)}$ and 
${\boldsymbol {f}}'_{inv}$ are both directly proportional to 
${\mathbb {S}}\boldsymbol{\cdot }{u}_s$, the former is negative for both types of straining flow and whatever the nature of the particle, while the latter is always positive. This observation sheds light on the debate summarized in § 1 regarding the proper form of the added-mass force in the low- but finite-Reynolds-number limit. Indeed, from a theoretical point of view, this turns out to be a ‘non-question’ since the dynamic boundary condition at the particle surface, by nature a viscous effect, is found to produce contributions to the force proportional to 
${\mathbb {S}}\boldsymbol{\cdot }{u}_s$, just as the added mass does if Taylor's expression involving the Lagrangian derivative of 
${\boldsymbol {U}}^\infty$ holds. Having noticed this, our only hope to disentangle both effects stood in the splitting of the inner contributions to 
$\boldsymbol{f}'^{(2)}$ into viscous and pressure drag components in (4.16) and (4.17), from which a comparison of all contributions for a solid sphere and a spherical bubble immersed in the same flow is possible. However, nothing emerged from this comparison, as the ratio of the viscous and pressure components (once the possible 
${4{\rm \pi} }/{3}$ added-mass contribution had been removed from the latter) differs for the two types of particle. Based on this unsuccessful attempt, we have to conclude that whether the added-mass force involves the fluid acceleration 
${\mbox {D} {\boldsymbol {U}}^\infty }/{\mbox {D} t}|_{{\boldsymbol {x}}_p}$ or the derivative of 
${\boldsymbol {U}}^\infty$ following the particle motion, i.e. 
${\mbox {d} {\boldsymbol {U}}^\infty }/{\mbox {d} t}|_{{\boldsymbol {x}}_p}$, is undecidable at a low but finite Reynolds number.
Things are different at moderate-to-large Reynolds numbers, as the direct simulations of Magnaudet et al. (Reference Magnaudet, Rivero and Fabre1995) showed. Their study considered the flow experienced by a rigid sphere or a spherical bubble held fixed in an axisymmetric uniaxial (or biaxial) straining flow over a wide range of Reynolds numbers. Whatever the Reynolds number 
${Re}_p$ (which was varied from 
$0.05$ to 
$150$), the magnitude of the fluid acceleration and the nature of the particle, the results revealed the existence of two series of contributions dependent on both strain and slip effects. One, found in the pressure drag, has a magnitude independent of 
${Re}_p$ and equal, to numerical accuracy, to that of the inviscid force predicted by (2.9). In contrast, the other, found in both the pressure and viscous drag components, depends on 
${Re}_p$ in the same way as the corresponding drag component in the case the body is immersed in a uniform flow. Its origin was readily identified by considering the vorticity and pressure distributions at the surface of the particle. In particular, compared with the same particle in a uniform stream, the fluid acceleration was observed to decrease the surface vorticity on the front half and to increase it on the rear half. These findings establish that these contributions are a direct consequence of the changes in the vorticity distribution at the particle surface, which result from the strain component of the base flow. To summarize, these simulations made it clear that (i) added-mass effects due to the advective acceleration 
${\boldsymbol {U}}^\infty \boldsymbol{\cdot }\boldsymbol {\nabla } {\boldsymbol {U}}^\infty \equiv -{\mathbb {S}}\boldsymbol{\cdot }{\boldsymbol {u}}_s$ exist in the considered configuration whatever 
${Re}_p$ and are unaffected by the Reynolds number; and (ii) the local strain characterizing the undisturbed flow modifies the near-surface vorticity distribution resulting from the dynamic boundary condition (no slip or shear free), which in turn yields additional contributions to the drag that depend linearly on 
${\mathbb {S}}$ but experience large variations with 
${Re}_p$. Apart from their different behaviours with respect to 
${Re}_p$, there is no way to separate the various 
${\mathbb {S}}$-dependent contributions in the overall drag. The theoretical low- but finite-Reynolds-number results obtained above in the planar and uniaxial straining flows are just another illustration of this coexistence of two categories of strain-dependent contributions to the drag resulting from two totally different physical mechanisms. Unlike the aforementioned simulations in which 
${Re}_p$ was varied by several orders of magnitude, these contributions appear at the same 
$O(\varepsilon ^2)$ order in the present asymptotic theory and, for this reason, cannot be disentangled.
4.4. Stationary limit: vortical flows
We now turn to the case of a linear shear flow in which the velocity-gradient tensor 
${\mathbb {A}}$ takes the form
\begin{equation} {\mathbb{A}} = \left(\begin{array}{@{}ccc@{}} 0 & 0 & 1\\ 0 & 0 & 0 \\ 0 & 0 & 0\\ \end{array}\right). \end{equation}
Saffman (Reference Saffman1965) computed the leading-order 
$O(\varepsilon )$-lift force in this flow for a sphere translating steadily along the direction of the undisturbed stream and possibly rotating perpendicular to it. Saffman also computed the 
$O(\varepsilon ^2)$ contributions to the lift force arising from the inner region of the disturbance. In contrast, Saffman did not consider the 
$O(\varepsilon ^2)$ force arising from the singular term 
${\boldsymbol {\mathcal {U}}}_2(t)$. Setting the slip velocity 
${\boldsymbol {u}}_s$ to 
$-{\boldsymbol {e}}_1$ and the angular velocity of the particle to 
$\boldsymbol {\omega }_p= \omega _p {\boldsymbol {e}}_2$, Saffman's second-order result (his (2.17)) reads
\begin{equation} {{\boldsymbol{f}}'}_S^{(2)} = \left( {\rm \pi}\omega_p -\frac{11{\rm \pi}}{8} \right){\boldsymbol{e}}_3 . \end{equation}
The present theory allows us to compute all 
$O(\varepsilon ^2)$ contributions, including those resulting from 
${\boldsymbol {\mathcal {U}}}_2(t)$. Moreover, the orientation of the slip velocity may be kept arbitrary. Evaluating all terms in (3.37) leads to
\begin{equation} {\boldsymbol{f}}'^{(2)} = \left(\begin{array}{@{}ccc@{}} 0 & 0 & \dfrac{41{\rm \pi}}{24} + {\rm \pi}\omega_p\\ 0 & 0 & 0 \\ \dfrac{11 {\rm \pi}}{8} - {\rm \pi}\omega_p & 0 & 0\\ \end{array}\right) \boldsymbol{\cdot} {\boldsymbol{u}}_s - 6{\rm \pi} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} {\boldsymbol{u}}_s .\end{equation}
The stationary kernel 
$\bar {{\mathbb {K}}}$ is known to be (Miyazaki et al. Reference Miyazaki, Bedeaux and Avalos1995, CMM)
\begin{equation} 6{\rm \pi} \bar{{\mathbb{K}}} \approx \left(\begin{array}{@{}ccc@{}} 0.0737 & 0 & 0.9436\\ 0 & 0.5766 & 0 \\ 0.3425 & 0 & 0.3269\\ \end{array} \right) ,\end{equation}
so that the leading-order correction to the Stokes drag is 
${\boldsymbol {f}}'^{(1)}=-6{\rm \pi} (6{\rm \pi} \bar {{\mathbb {K}}})\boldsymbol{\cdot }{\boldsymbol {u}}_s$. In particular, the leading-order lift component 
${\boldsymbol {f}}'^{(1)}\boldsymbol{\cdot }{\boldsymbol {e}}_3$ in the configuration considered by Saffman is 
${\boldsymbol {f}}'^{(1)}\boldsymbol{\cdot }{\boldsymbol {e}}_3=6{\rm \pi} (6{\rm \pi} {\boldsymbol {e}}_3\boldsymbol{\cdot }\bar {{\mathbb {K}}} \boldsymbol{\cdot }{\boldsymbol {e}}_1)\approx 6.456$.
Using (4.21), (4.20) provides the next-order force as
\begin{equation} {\boldsymbol{f}}'^{(2)} = \left(\begin{array}{@{}ccc@{}} 0 & 0 & \dfrac{41{\rm \pi}}{24} + {\rm \pi}\omega_p\\ 0 & 0 & 0 \\ \dfrac{11 {\rm \pi}}{8} - {\rm \pi}\omega_p & 0 & 0\\ \end{array}\right) \boldsymbol{\cdot} {\boldsymbol{u}}_s - \left(\begin{array}{@{}ccc@{}} 6.1942 & 0 & 7.1252\\ 0 & 6.2669 & 0 \\ 2.5862 & 0 & 8.1062\\ \end{array}\right) \boldsymbol{\cdot} {\boldsymbol{u}}_s .\end{equation}If the particle does not rotate, this reduces to
\begin{equation} {\boldsymbol{f}}'^{(2)} ={-} \left(\begin{array}{@{}ccc@{}} 6.1942 & 0 & 1.7583\\ 0 & 6.2669 & 0 \\ -1.7335 & 0 & 8.1062\\ \end{array}\right) \boldsymbol{\cdot} {\boldsymbol{u}}_s .\end{equation}
The off-diagonal components of 
$\boldsymbol{f}'^{(2)}$ yield the second-order lift force 
${\boldsymbol {f}}'^{(2)}_{L}\approx -1.7583({\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {e}}_3) {\boldsymbol {e}}_1+1.7335({\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {e}}_1){\boldsymbol {e}}_3$. In the specific case 
${\boldsymbol {u}}_s=- {\boldsymbol {e}}_1$, (4.23) implies that 
${\boldsymbol {f}}'^{(2)}\boldsymbol{\cdot }{\boldsymbol {e}}_3=-1.7335$, to be compared with 
${\boldsymbol {f}}_S^{(2)}\boldsymbol{\cdot }{\boldsymbol {e}}_3=-\frac {11}{8}{\rm \pi} \approx -4.32$ according to (4.19). Hence, the actual magnitude of the second-order lift force is approximately 
$2.5$ times smaller than Saffman's incomplete prediction. It may be noted that the two prefactors involved in 
${\boldsymbol {f}}'^{(2)}_{L}$ have very close values, albeit with opposite signs, so that on a purely empirical basis the second-order lift force may be approximated as 
${\boldsymbol {f}}'^{(2)}_{L}\approx 1.746\,{\boldsymbol {u}}_s\times ({\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )$. This expression may be compared with the prediction for the inviscid disturbance-induced force 
${\boldsymbol {f}}'_{inv}$ in (2.9), which in the present case reduces to the shear lift force 
${\boldsymbol {f}}'_{inv}=-({2{\rm \pi} }/{3}){\boldsymbol {u}}_s\times ({\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )$. Although the two expressions have the same structure and the two prefactors have a comparable magnitude (since 
${2{\rm \pi} }/{3}\approx 2.094$), they have opposite signs, a clear indication that the mechanism responsible for 
${\boldsymbol {f}}'^{(2)}_{L}$ is totally different from the Lighthill mechanism discussed in § 4.2.
Similarly, if the particle rotates freely, the torque-free condition implies that 
$\|\boldsymbol {\omega }_s\| =O(\varepsilon ^2)$, i.e. 
$\boldsymbol {\omega }_p\approx \boldsymbol {\varOmega }=\frac {1}{2}{\boldsymbol {e}}_2$, so that (4.22) simplifies to
\begin{equation} \boldsymbol{f}'^{(2)} = \left(\begin{array}{@{}ccc@{}} -6.1942 & 0 & -0.1876\\ 0 & -6.2669 & 0 \\ 0.1626 & 0 & -8.1061\\ \end{array} \right)\boldsymbol{\cdot} {\boldsymbol{u}}_s .\end{equation}
Now the off-diagonal components of 
$\boldsymbol{f}'^{(2)}$ yield the second-order lift force 
${\boldsymbol {f}}'^{(2)}_{L}\approx -0.1876({\boldsymbol {u}}_s\boldsymbol{\cdot } {\boldsymbol {e}}_3){\boldsymbol {e}}_1+0.1626({\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {e}}_1){\boldsymbol {e}}_3$. Again, the two prefactors have very close values, so that this contribution may be empirically approximated in the form 
${\boldsymbol {f}}'^{(2)}_{L}\approx 0.175{\boldsymbol {u}}_s\times ({\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )$. The prefactors involved in 
${\boldsymbol {f}}'^{(2)}_{L}$ are fairly small compared with those found in the non-rotating case because the contribution brought by the singular term 
$- 6{\rm \pi} (6{\rm \pi} \bar {{\mathbb {K}}}) \boldsymbol{\cdot } (6{\rm \pi} \bar {{\mathbb {K}}})\boldsymbol{\cdot }{\boldsymbol {u}}_s$ nearly balances the sum of the regular terms. In particular, with 
${\boldsymbol {u}}_s=- {\boldsymbol {e}}_1$, (4.24) implies that 
${\boldsymbol {f}}'^{(2)}_L=-0.1626{\boldsymbol {e}}_3$. In contrast, with the relevant angular velocity 
${\omega }_p=\frac {1}{2}$, Saffman's partial result (4.19) yields 
${\boldsymbol {f}}_S^{(2)}=-\frac {7}{8}{\rm \pi} {\boldsymbol {e}}_3\approx -2.75{\boldsymbol {e}}_3$, which overpredicts the actual second-order lift force by more than one order of magnitude. The above two examples indicate that Saffman's incomplete result for the second-order correction to the lift force is grossly inaccurate. Therefore, practitioners using point-particle models incorporating the steady-state form of the Saffman lift force should either consider only the corresponding leading-order prediction or make use of (4.24) or (4.23) to estimate the second-order correction, depending on whether the particle rotates or not.
Last, we turn to the solid-body rotation flow characterized by the velocity-gradient tensor
\begin{equation} {\mathbb{A}} = \left(\begin{array}{@{}ccc@{}} 0 & -1 & 0\\ 1 & 0 & 0 \\ 0 & 0 & 0\\ \end{array}\right) .\end{equation}
Since 
${\mathbb {S}}\equiv \boldsymbol {0}$, there is no source term in (3.38), so that the torque-free condition implies that 
$\boldsymbol {\omega }_s=\boldsymbol {0}$ if the particle is not given an initial differential rotation. Therefore, once the slip velocity no longer varies, (3.37) reduces to
\begin{equation} {\boldsymbol{f}}'^{(2)} = \frac{4{\rm \pi}}{3} \boldsymbol{\varOmega}\times {\boldsymbol{u}}_s- 6{\rm \pi} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} {\boldsymbol{u}}_s .\end{equation}For this flow, the steady-state kernel may be obtained in closed form and reads (Gotoh Reference Gotoh1990; Miyazaki Reference Miyazaki1995, CMM)
\begin{equation} 6{\rm \pi} \bar{{\mathbb{K}}} = \left(\begin{array}{@{}ccc@{}} \dfrac{3\sqrt{2}(19+9\sqrt{3})}{280} & -\dfrac{3\sqrt{2}(19-9\sqrt{3})}{280} & 0\\ \dfrac{3\sqrt{2}(19-9\sqrt{3})}{280} & \dfrac{3\sqrt{2}(19+9\sqrt{3})}{280} & 0 \\ 0 & 0 & \dfrac{4}{7} \end{array} \right).\end{equation}
The off-diagonal terms of (4.27) yield the leading-order rotational lift force 
${\boldsymbol {f}}'^{(1)}_{L}=({9{\rm \pi} \sqrt {2}(19-9\sqrt {3})}/{140})[({\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {e}}_2){\boldsymbol {e}}_1- ({\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {e}}_1){\boldsymbol {e}}_2]$ 
$\approx 0.980[({\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {e}}_2){\boldsymbol {e}}_1-({\boldsymbol {u}}_s \boldsymbol{\cdot }{\boldsymbol {e}}_1){\boldsymbol {e}}_2]$. Making use of (4.26), the second-order disturbance-induced force is found to be
\begin{align} {\boldsymbol{f}}'^{(2)} &=\frac{4{\rm \pi}}{3} {\boldsymbol{e}}_3\times {\boldsymbol{u}}_s-6{\rm \pi} \left(\begin{array}{ccc} \dfrac{1539\sqrt{3}}{9800} & -\dfrac{531}{9800} & 0\\ \dfrac{531}{9800} & \dfrac{1539\sqrt{3}}{9800} & 0 \\ 0 & 0 & \dfrac{16}{49}\\ \end{array} \right)\boldsymbol{\cdot}{\boldsymbol{u}}_s\nonumber\\ &\approx \left(\begin{array}{@{}ccc@{}} -5.127 & -3.167 & 0\\ 3.167 & -5.127 & 0\\ 0 & 0 & -6.155\\ \end{array} \right) \boldsymbol{\cdot}{\boldsymbol{u}}_s. \end{align}
In particular, the second-order rotational lift force resulting from the off-diagonal terms of (4.28) is 
${\boldsymbol {f}}'^{(2)}_{L}\approx -3.167[({\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {e}}_2) {\boldsymbol {e}}_1-({\boldsymbol {u}}_s\boldsymbol{\cdot }{\boldsymbol {e}}_1){\boldsymbol {e}}_2]\approx -1.584 {\boldsymbol {u}}_s\times ({\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )$, since 
${\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty =2{\boldsymbol {e}}_3$ according to the definition (4.25) of the velocity-gradient tensor.
The prefactor involved in 
${\boldsymbol {f}}'^{(2)}_{L}$ is significantly larger than that found in 
${\boldsymbol {f}}'^{(1)}_{L}$ and of opposite sign. The total lift force at the present order of approximation being 
${\boldsymbol {f}}'_{L} = \varepsilon {\boldsymbol {f}}'^{(1)}_{L} +\varepsilon ^2{\boldsymbol {f}}'^{(2)}_{L}$, this force changes sign and becomes dominated by the second-order contribution beyond 
$\varepsilon \approx 0.31$, i.e. beyond a critical shear Reynolds number 
${Re}_s\approx 0.096$. Hence, under many practical conditions, one may expect a small particle held fixed in a solid-body rotation flow to experience a lift force of opposite sign to the one predicted by the leading-order estimate 
${\boldsymbol {f}}'^{(1)}_{L}$. This surprising feature contrasts with what happens in the pure shear configuration. Indeed, (4.24) indicates that the total lift force experienced by a torque-free particle with a slip velocity 
${\boldsymbol {u}}_s=-{\boldsymbol {e}}_1$ is in that case 
${\boldsymbol {f}}'_{L} =(6.456\varepsilon -0.1626\varepsilon ^2){\boldsymbol {e}}_3$. Again, the first- and second-order contributions have opposite signs but the prefactor of the 
$O(\varepsilon ^2)$ term is much smaller than that of the 
$O(\varepsilon )$ term, so that 
${\boldsymbol {f}}'_{L}$ never changes sign within the domain of validity of the present theory.
We may again compare 
${\boldsymbol {f}}'^{(2)}_{L}$ with the inviscid force predicted by (2.9), noting that in the present case 
${\mbox {D} {\boldsymbol {U}}^\infty }/{\mbox {D} t}|_{{\boldsymbol {x}}_p}\equiv \frac {1}{2}{\boldsymbol {u}}_s \times ({\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )$ if the particle is at rest. For this reason, the added-mass and shear-induced lift forces in (2.9) combine in the form of an ‘extended’ lift force 
${\boldsymbol {f}}'_{inv}=-({{\rm \pi} }/{3}) {\boldsymbol {u}}_s \times ({\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty )$. This inviscid force and the second-order lift force 
${\boldsymbol {f}}'^{(2)}_{L}$ are seen to have the same sign, which contrasts with what was noted above in the pure shear configuration. However, the magnitude of 
${\boldsymbol {f}}'^{(2)}_{L}$ (
$3.167$) is larger than that of 
${\boldsymbol {f}}'_{inv}$ (
$\approx 2.094$), which once again underlines the role of the no-slip condition in the structure of the disturbance flow. That 
${\boldsymbol {f}}'^{(2)}_{L}$ and 
${\boldsymbol {f}}'_{inv}$ keep the same sign in a solid-body rotation flow but have opposite signs in a linear shear flow emphasizes the crucial importance of the alteration introduced in the near-surface disturbance by the presence of a non-zero strain in the latter case. More specifically, this change of sign suggests that the surface vorticity associated with this strain-induced disturbance contributes more to the lift force than any other mechanism affecting the vorticity distribution around the particle, especially those resulting from the stretching and tilting of the upstream vorticity 
${\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty$.
5. Summary and concluding remarks
In this work, we computed the second-order inertial corrections to the creeping-flow approximation of the time-dependent hydrodynamic force and torque acting on a small rigid sphere immersed in a general steady linear flow. Our primary motivation was to obtain a consistent approximation for the force and torque in which all fundamental physical effects resulting from fluid inertia, such as shear- and spin-induced lift and added mass, are accounted for to 
$O(\varepsilon ^2)$. To achieve this goal, we used and extended the methodology developed by CMM that employs asymptotic matching expansions and formulates the problem in a coordinate system co-moving with the carrying flow. Some of the second-order corrections derived here were already known, but others were not. Computing these corrections in a systematic fashion is much more difficult than obtaining the first-order corrections. This is on the one hand because the inner problem at 
$O(\varepsilon ^2)$ is inhomogeneous, and on the other hand because the second-order singular contribution brought by the outer solution results from a double convolution product.
Our results for the second-order force and torque are summarized in (3.37) and (3.38). However, it is only after the singular outer velocity correction 
${\boldsymbol {\mathcal {U}}}_2(t)$ has been explicitly computed that one can fully appreciate how the different contributions combine in the second-order force, and of which physical effect they account for. This is especially true regarding terms involving the relative acceleration between the particle and fluid. In particular, the classical added-mass force in a uniform flow is found to result from the sum of an inner and an outer contribution evaluated in the short-time limit 
$t\ll 1$, and it is dominated by the latter. That the outer contribution dominates indicates that the majority of the fluid instantaneously displaced when the particle or the fluid accelerates is located far from the body, i.e. at distances larger than the Saffman length. The expression (3.38) for the torque also comprises a contribution involving the instantaneous relative angular acceleration, as if there were a non-zero rotational added-mass effect. This is of course not the case given the point-symmetric geometry of the considered particle, and this contribution is just due to the ‘remains’ of the exact history torque beyond the initial stage that extends over a few viscous time units. Indeed, it is important to bear in mind that the present theory is not designed to deal with very large levels of unsteadiness. For this reason, it does not in general provide the exact form of the kernel at very short times, i.e. 
$t\lesssim \varepsilon ^2$. Another example of this limitation is observed in the case of the force acting on a spherical bubble, for which the total contribution proportional to the instantaneous relative acceleration is found to be the sum of the added-mass force and of a (dominant) term left by the early decay of the exact history kernel. In the present theory, history effects appear at first and second orders through the far-field velocity corrections 
${\boldsymbol {\mathcal {U}}}_1(t)$ and 
${\boldsymbol {\mathcal {U}}}_2(t)$, respectively. While these outer contributions are crucial in the inertial corrections to the force, they do not affect the torque in the case of a spherical (or spheroidal) particle, for symmetry reasons.
In (3.38) two quadratic contributions to the second-order inertial torque are identified. We showed that the one proportional to 
${\mathbb {S}}\boldsymbol{\cdot }\boldsymbol \varOmega$ is identical to that computed by Candelier et al. (Reference Candelier, Einarsson and Mehlig2016) for a neutrally buoyant particle. It is non-zero only in three-dimensional linear flows whose structure generates a uniform vortex stretching, which makes it relevant to the motion of small particles in turbulence. The contribution proportional to 
${\mathbb {S}}\boldsymbol{\cdot }\boldsymbol \omega _s$ was apparently not identified so far. It reveals that a spinning particle immersed in a linear flow with a non-zero strain component aligned with the spin axis experiences an inertial torque. This correction makes the total torque 
$\boldsymbol \tau '=\boldsymbol \tau '^{(0)}+\varepsilon ^2\boldsymbol \tau '^{(2)}$ weaker (respectively stronger) than the primary viscous torque if the particle spins about the elongational (respectively compressional) axis of the straining flow. It induces a torque component orthogonal to the primary spinning direction otherwise. For instance, the total torque on a particle spinning with the angular velocity 
$\omega _p{\boldsymbol {e}}_1$ in the linear shear flow 
${\boldsymbol {U}}^\infty =x_3{\boldsymbol {e}}_1$ is 
$\boldsymbol \tau '=-8{\rm \pi} \omega _p({\boldsymbol {e}}_1-({\varepsilon ^2}/{16}){\boldsymbol {e}}_3)$.
The second-order force (3.37) involves three quadratic contributions. The one corresponding to the spin-induced lift force originally computed by Rubinow & Keller (Reference Rubinow and Keller1961) is well known. Our results show that this force subsists in a linear flow, provided the slip rotation rate 
$\boldsymbol {\omega }_s=\boldsymbol {\omega }_p-\boldsymbol \varOmega$ is substituted for the particle spin rate 
$\boldsymbol {\omega }_p$. The other two contributions involve the kinematic structure of the background flow, through 
${\mathbb {S}}$ and 
$\boldsymbol \varOmega$, and the particle slip velocity, 
${\boldsymbol {u}}_s$. The contribution proportional to 
$\boldsymbol \varOmega \times {\boldsymbol {u}}_s$ is a pure lift force while that proportional to 
${\mathbb {S}}\boldsymbol{\cdot }{\boldsymbol {u}}_s$ may comprise drag and lift components. We compared the prefactors of the corresponding terms with those found in inviscid flow in the weak shear limit (
$as/u_c\ll 1$). We also compared the prefactors of the 
${\mathbb {S}}\boldsymbol{\cdot }{\boldsymbol {u}}_s$ contribution determined in an axisymmetric uniaxial flow for a rigid particle and a spherical bubble, respectively. These comparisons shed light on the central role played by the dynamic boundary condition at the particle surface. To a large extent, the relative magnitude of these contributions is determined by the changes imposed to the vorticity distribution at the particle surface by the strain and/or rotation component of the base flow. These changes modify the spatial distribution of the fluid momentum in the region close to the particle, i.e. at distances less than 
$\ell _s$, which in turn results in a net force on it. The velocity gradients 
${\mathbb {S}}$ and 
$\boldsymbol \varOmega$ enter the kernel 
${\mathbb {K}}$ in a subtle manner (this is what makes it flow dependent), and hence also the far-field velocity correction 
${\boldsymbol {\mathcal {U}}}_2$. For a given flow, the second-order force on a rigid particle tends to 
$({37{\rm \pi} }/{12}) {\mathbb {S}} \boldsymbol{\cdot } {\boldsymbol {u}}_s + ({4{\rm \pi} }/{3}) \boldsymbol {\varOmega }\times {\boldsymbol {u}}_s+ 6{\rm \pi} {\boldsymbol {\mathcal {U}}}_2$ in the long-term limit. We found that, for a given slip velocity, the inner and outer contributions have opposite signs, and the sign of their sum depends on the considered flow, as the examination of several canonical configurations in § 4 revealed. Moreover, the total lift force is obtained by summing the first- and second-order lift components. The two may have opposite signs, especially in vortical flows. Although the first-order term is necessarily dominant from an asymptotic point of view, we observed that in a solid-body rotation flow the prefactor of the second-order component is significantly larger than that of the first-order one, making the total lift force change sign beyond a modest value of the shear Reynolds number.
Having isolated the 
$({37{\rm \pi} }/{12}) {\mathbb {S}} \boldsymbol{\cdot } {\boldsymbol {u}}_s$ contribution to the second-order force in straining flows led us to re-examine the old yet still open question whether one should use the Lagrangian derivative, 
${\mbox {D} {\boldsymbol {U}}^\infty }/{\mbox {D} t}|_{{\boldsymbol {x}}_p}$, or the time derivative following the particle, 
${\mbox {d} {\boldsymbol {U}}^\infty }/{\mbox {d} t}|_{{\boldsymbol {x}}_p}$, in the expression for the added-mass force in non-uniform flows. Indeed, selecting the former, which is known to be the proper choice in the inviscid limit, introduces a contribution 
$-({2{\rm \pi} }/{3}){\mathbb {S}} \boldsymbol{\cdot } {\boldsymbol {u}}_s$ in the second-order force in straining flows, while no such contribution appears if the latter is chosen. However, the above question turned out to be undecidable at the present order of approximation, owing to the dominant influence of the dynamic boundary condition on the strength of the 
${\mathbb {S}} \boldsymbol{\cdot } {\boldsymbol {u}}_s$ term, as underlined by the opposite signs of the above two prefactors. Indeed, the flow configurations examined in § 4 did not bring any convincing clue that might allow us to conclude that the 
${37{\rm \pi} }/{12}$ prefactor should rather be interpreted as 
$-{2{\rm \pi} }/{3}+{15{\rm \pi} }/{4}$, with the first and second contributions corresponding respectively to the added-mass effect induced by the no-penetration condition and the ‘vortical’ effect resulting from the dynamic boundary condition. Higher-order expansions, or more likely DNS results covering a wide range of Reynolds number, such as those of Magnaudet et al. (Reference Magnaudet, Rivero and Fabre1995), may help to disentangle the two types of contributions. At the moment, the latter reference provides the clearest heuristic evidence that the inviscid form of the added-mass term remains valid at finite Reynolds number.
In view of applications, it appears useful to summarize the complete expressions for the force and torque resulting from the present theory in dimensional variables. We remind the reader that this theory assumes that the background flow has the form 
${\boldsymbol {U}}_\infty ({\boldsymbol {x}},t)={\boldsymbol {U}}_0(t)+{\mathbb {S}}\boldsymbol{\cdot }{\boldsymbol {x}}+ \boldsymbol \varOmega \times {\boldsymbol {x}}$ (hence, 
${\boldsymbol {\nabla }} \times {\boldsymbol {U}}^\infty =2\boldsymbol \varOmega$), with 
${\mathbb {S}}$ and 
$\boldsymbol \varOmega$ uniform and time independent. Actually, the latter restriction applies to the results for the force, while those for the torque remain valid even through 
${\mathbb {S}}$ and 
$\boldsymbol \varOmega$ depend upon time. This restriction also implies that, in three-dimensional flows, the results for the force are only valid in pure straining flows (
$\boldsymbol \varOmega ={\boldsymbol {0}}$). With these restrictions in mind, using the definitions introduced in §§ 1 and 2, the total force and torque, obtained by summing the creeping-flow result with the first- and second-order inertial corrections and the force/torque resulting from the undisturbed flow, read
\begin{align} {\boldsymbol{f}}(t)&={\boldsymbol{f}}_b(t)-6{\rm \pi}\mu a\left\{{\boldsymbol{u}}_s+\left(\frac{a^2s}{\nu}\right)^{1/2} {\int_0^t} 6{\rm \pi} {\mathbb{K}}(s(t-\tau))\boldsymbol{\cdot} \frac{\textrm{d}{\boldsymbol{u}}_s}{\textrm{d}\tau}\,{\rm d}\tau\right\}\nonumber\\ &\quad +m_f\left\{4 \frac{\mbox{d}{\boldsymbol{u}}_s}{\mbox{d}t} + \frac{37}{16} {\mathbb{S}}\boldsymbol{\cdot} {\boldsymbol{u}}_s- \frac{1}{2} {\boldsymbol{u}}_s \times \left( {\boldsymbol{\nabla}} \times{\boldsymbol{U}}^\infty\right)- \frac{3}{4} {\boldsymbol{u}}_s\times\boldsymbol{\omega}_s\right. \nonumber\\ &\quad -\left. \frac{9}{2} {\int_0^t} 6{\rm \pi} {\mathbb{K}}(s(t-\tau))\boldsymbol{\cdot} \frac{\textrm{d}}{\textrm{d}\tau}\left( \int_0^\tau 6{\rm \pi} {\mathbb{K}}(s(t-\tau))\boldsymbol{\cdot} \frac{\textrm{d}{\boldsymbol{u}}_s}{\textrm{d}\tau'}\,{\rm d}\tau'\right) \, \textrm{d}\tau\right\}, \end{align}
\begin{align} \boldsymbol{\tau}(t) &={-}8{\rm \pi}\mu a^3\boldsymbol{\omega}_s -m_fa^2\left(2\frac{\mbox{d}\boldsymbol{\omega}_s }{\mbox{d}t} + \frac{6}{5}{\mathbb{S}}\boldsymbol{\cdot} ( {\boldsymbol{\nabla}} \times {\boldsymbol{U}}^\infty) -\frac{3}{4}{\mathbb{S}}\boldsymbol{\cdot} \boldsymbol{\omega}_s\right), \end{align}with
\begin{align} {\boldsymbol{f}}_b(t)&=(m_p-m_f){\boldsymbol{g}}+m_f\left.\frac{\mbox{D} {\boldsymbol{U}}^\infty}{\mbox{D} t} \right|_{{\boldsymbol{x}}_p(t)}\nonumber\\ &=(m_p-m_f){\boldsymbol{g}}+m_f\left\{\dot{{\boldsymbol{U}}}_0+{\mathbb{S}} \boldsymbol{\cdot}{\boldsymbol{U}}^\infty|_{{\boldsymbol{x}}_p(t)}+\frac{1}{2} ( {\boldsymbol{\nabla}} \times {\boldsymbol{U}}^\infty) \times{\boldsymbol{U}}^\infty|_{{\boldsymbol{x}}_p(t)}\right\} . \end{align}
Since 
${\mathbb {K}}(st)\rightarrow ({1}/{6{\rm \pi} })({{\mathbb {I}}}/{\sqrt {{\rm \pi} st}})$ for 
$st\rightarrow 0$, the dimensional force in the short-term limit is
\begin{align} {\boldsymbol{f}}({t\ll s^{{-}1}})&={\boldsymbol{f}}_b(t)-6{\rm \pi}\mu a\left\{{\boldsymbol{u}}_s+\int_0^t \frac{1}{\sqrt{{\rm \pi}\nu(t-\tau)/a^2}} \frac{\textrm{d}{\boldsymbol{u}}_s}{\textrm{d}\tau}\,{\rm d}\tau\right\}\nonumber\\ &\quad +m_f\left\{-\frac{1}{2} \frac{\mbox{d}{\boldsymbol{u}}_s}{\mbox{d}t} +\frac{37}{16} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s - \frac{1}{2} {\boldsymbol{u}}_s \times ( {\boldsymbol{\nabla}} \times {\boldsymbol{U}}^\infty) - \frac{3}{4}{\boldsymbol{u}}_s\times\boldsymbol{\omega}_s\right\} . \end{align}
This is the original BBO equation (1.1) enriched with the last three terms that account for the interaction of the slip velocity with the background velocity gradients and the particle rotation. Note that (5.4) remains valid for times of 
$O(a^2/\nu )$ or even shorter, since the short-time limit of the history kernel and that of the added-mass force coincide with those found in the BBO approximation. This is a bonus specific to the case of the force on a rigid spherical particle. The same does not hold for a bubble, nor for the torque on a spherical particle, because in these cases the BBO kernel involves a term proportional to 
$\exp [{9\nu (t-\tau )/a^2}]\text {erfc}[{9\nu (t-\tau )/a^2}]^{1/2}$ that is not captured by the present theory. A similar term is known to exist for drops of arbitrary viscosity (Gorodtsov Reference Gorodtsov1975; Yang & Leal Reference Yang and Leal1991) and non-spherical particles (Lawrence & Weinbaum Reference Lawrence and Weinbaum1986). Obtaining expressions for the loads that incorporate finite-Reynolds-number effects and are uniformly valid in the limit 
$t\rightarrow 0$ even when the exact kernel does not reduce to the Basset–Boussinesq form 
$[\nu (t-\tau )/a^2]^{-1/2}$ requires a theory fundamentally different from the one described here, the starting point of which is the unsteady Stokes equation. Since, instead of the simple form (3.16), the solution to this equation is non-local in time, a good part of the methodology outlined in § 3 no longer applies and it is doubtful that results may be obtained in closed form in the time domain. Nevertheless, the problem is worthy of consideration since in emerging fields, such as acoustofluidics, the particle motion is driven by a high-frequency acoustic field, which corresponds to high levels of unsteadiness with 
${Re}_s\text {Sl}=O(1)$ in (2.2b), while small but finite advective effects may have a significant influence (Agarwal et al. Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021).
For 
$st\rightarrow \infty$, 
${\mathbb {K}}(st)\rightarrow \bar {{\mathbb {K}}}$, so that in the long-term limit one has
\begin{align} {\boldsymbol{f}}({t\gg s^{{-}1}})&={\boldsymbol{f}}_b-6{\rm \pi}\mu a\left\{{\boldsymbol{u}}_s+\left(\frac{a^2s}{\nu}\right)^{1/2} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} {\boldsymbol{u}}_s\right\} \nonumber\\ &\quad +m_f\left\{\frac{37}{16} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s - \frac{1}{2} {\boldsymbol{u}}_s \times ( {\boldsymbol{\nabla}} \times {\boldsymbol{U}}^\infty) - \frac{3}{4}{\boldsymbol{u}}_s\times\boldsymbol{\omega}_s\right.\notag\\ &\quad - \left. \frac{9}{2}s(6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} (6{\rm \pi} \bar{{\mathbb{K}}}) \boldsymbol{\cdot} {\boldsymbol{u}}_s\right\} . \end{align}
In (5.1) and (5.5) the velocity-gradient scale 
$s$ is defined as 
$s=(({1}/{\mathcal {D}!}){\mathbb {A}}\colon {\mathbb {A}})^{1/2}$, with 
$\mathcal {D}=1,2$ or 
$3$ depending on whether the base flow is one, two or three dimensional.
Besides the fact that the background flow is assumed stationary and all Reynolds numbers involved have to be small, the present theory suffers from two main limitations. First, the condition 
${Re}_p\ll {Re}_s^{1/2}$ requires the slip velocity to be very small, making the theory essentially suitable for weakly positively or negatively neutrally buoyant particles. In particular, our second-order results do not comprise the classical Oseen correction to the drag, which results from advective terms that are negligible in the outer region under the above condition. Beyond the drag increase they induce (Kaplun & Lagerstrom Reference Kaplun and Lagerstrom1957; Proudman & Pearson Reference Proudman and Pearson1957), these terms are known to decrease the strength of the Saffman lift force (Asmolov Reference Asmolov1990; McLaughlin Reference McLaughlin1991). For instance, this force is reduced by 
$25\,\%$ when 
${Re}_p\approx {Re}_s^{1/2}$, and the larger the ratio 
${Re}_p/{Re}_s^{1/2}$ the smaller the lift force. Consequently, finite slip effects affect the kernel components involved in the first- and second-order inertial corrections, and extending our theory to incorporate these effects is crucial to make it applicable to situations involving large fluid-to-particle density ratios.
Second, the kernel is only known for the simplest linear flows. Since the problem is nonlinear, one cannot just determine the general kernel by linearly superposing the ‘elementary’ kernels for the canonical flows that were considered in § 4. To make the present theory useful in applications involving arbitrary linear flows, it is necessary to render the dependence of the kernel with respect to the components of the velocity-gradient tensor 
${\mathbb {A}}$ explicit. Closed-form expressions of 
${\mathbb {K}}$ in a general linear flow may presumably be obtained only in the short- and long-term limits, while semi-empirical fits will certainly be required to obtain approximate expressions valid for arbitrary times. This, we think, is the price to pay for making the results (5.1)–(5.5) applicable in practice to a wide range of configurations of interest in sheared suspensions and possibly in turbulent flows.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2022.1015.
Funding
B.M. was supported by Vetenskapsrådet (grant no. 2021-4452) and by the Knut-and-Alice Wallenberg Foundation (grant no. 2019.0079).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Second-order force on a spherical bubble
We complemented the results for a rigid sphere with those for a spherical bubble. Our initial hope was that this might allow us to disentangle inertial effects resulting from the particle shape (i.e. added-mass effects) from those resulting from the vorticity generated at the particle surface. Indeed, comparing a spherical bubble to a rigid sphere leaves the no-penetration condition unchanged, whereas the no-slip condition is changed into a shear-free condition for the bubble, assuming that the viscosity of the gas that fills it is negligibly small and the gas–liquid interface is free of any contamination. For technical reasons related to the determination of the second-order inner solution, we did not succeed to solve this problem for a general linear flow. Nevertheless, we managed to solve it when the bubble moves, with a possibly time-dependent slip velocity, in the axisymmetric straining flow considered in Appendix B.
The steps involved in the calculation of the first- and second-order force disturbance on a spherical bubble are similar to those described in § 3. Only the boundary condition at the particle surface differs. In the 
$O(\varepsilon )$ problem, the first boundary condition in (3.23) is replaced by
\begin{equation} {{\boldsymbol{v}}}_{{in}}^{(1)} \boldsymbol{\cdot} {\boldsymbol{n}} ={-}{\boldsymbol{\mathcal{U}}}^{(1)}(t) \boldsymbol{\cdot} {\boldsymbol{n}},\quad {\boldsymbol{n}}\times(\boldsymbol\sigma^{(1)}\boldsymbol{\cdot} {\boldsymbol{n}})= {\boldsymbol{0}} \ \mbox{at}\ r=1, \end{equation}
while at order 
$\varepsilon ^2$, instead of the first boundary condition in (3.29) one has
\begin{equation} {\boldsymbol{w}}_{h} \boldsymbol{\cdot} {\boldsymbol{n}} ={-} ({\boldsymbol{w}}_p + {\boldsymbol{\mathcal{U}}}_2(t) ) \boldsymbol{\cdot} {\boldsymbol{n}} , \quad {\boldsymbol{n}}\times(\boldsymbol\sigma_{h}\boldsymbol{\cdot} {\boldsymbol{n}})={-} {\boldsymbol{n}}\times(\boldsymbol\sigma_{p}\boldsymbol{\cdot} {\boldsymbol{n}})\ \mbox{at}\ r=1 . \end{equation}
In (A1), 
$\boldsymbol \sigma ^{(1)}\boldsymbol{\cdot }{\boldsymbol {n}}$ denotes the first-order viscous traction resulting from the velocity field 
${{\boldsymbol {v}}}_{{in}}^{(1)}$, while in (A2a,b), 
$\boldsymbol{\sigma}_{h}\boldsymbol{\cdot } {\boldsymbol {n}}$ and 
$\boldsymbol{\sigma}_{p}\boldsymbol{\cdot } {\boldsymbol {n}}$ denote the second-order viscous tractions resulting from the velocity components 
${\boldsymbol {w}}_{h}$ and 
${\boldsymbol {w}}_{p}$ defined in (3.30a,b), respectively.
Carrying out the steps described in § 3, the total disturbance force acting on the bubble is found to be
\begin{align} {\boldsymbol{f}}'&= {\boldsymbol{f}}'^{(0)} +\varepsilon{\boldsymbol{f}}'^{(1)}+\varepsilon^2{\boldsymbol{f}}'^{(2)}={-}4{\rm \pi} {\boldsymbol{u}}_s - 4{\rm \pi} \varepsilon \int_0^t {\mathbb{K}}(t-\tau) \frac{\mbox{d} {\boldsymbol{f}}^{(0)}}{\mbox{d}\tau}\, \mbox{d} \tau \nonumber\\ &\quad + \varepsilon^2 \left(2{\rm \pi} \frac{\mbox{d} {\boldsymbol{u}}_s}{\mbox{d}t}+\frac{4{\rm \pi}}{3} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s - 4{\rm \pi} \int_0^t {\mathbb{K}}(t-\tau) \frac{\mbox{d}}{\mbox{d}\tau}\left[ \int_0^\tau 4{\rm \pi} {\mathbb{K}}(\tau-\tau') \frac{\mbox{d} {\boldsymbol{f}}^{(0)}}{\mbox{d}\tau'}\, \mbox{d}\tau'\right]\, \mbox{d}\tau \right) . \end{align} As for a rigid sphere, the different terms in (A3) may be split into pressure and viscous contributions. The result is detailed in table 3. Note that the 
$({4{\rm \pi} }/{3}) {\mathbb {S}} \boldsymbol{\cdot } {\boldsymbol {u}}_s$ contribution was only determined for a bubble located on the symmetry axis of the axisymmetric straining flow defined by (B1). However, whatever the bubble position and the velocity-gradient tensor 
${\mathbb {A}}$ characterizing the linear straining flow, the stationary second-order force resulting from the inner solution must have the vector form 
${\mathbb {S}} \boldsymbol{\cdot } {\boldsymbol {u}}_s$, similar to (3.37) for a rigid sphere. This is why the particular result derived above remains valid for a bubble standing at an arbitrary position in an arbitrary linear straining flow.
Table 3. Summary of all contributions to the second-order force on a bubble moving along the axis of an axisymmetric elongational flow, expressed with respect to 
${\boldsymbol {f}}^{(0)} = 4{\rm \pi} {\boldsymbol {u}}_s$. Terms arising from the regular and singular parts of the solution are given in the upper and lower tables, respectively. In the last row of the former, note that 
$\frac {3029}{6200} -\frac {281}{9300}-\frac {1}{8}=\frac {1}{3}$, which yields the 
$({4{\rm \pi} }/{3}){\mathbb {S}}\boldsymbol{\cdot }{\boldsymbol {u}}_s$ contribution in (A3). We also distinguish between the contributions to the force induced by the viscous part of the stress tensor (e.g. 
${\boldsymbol {f}}'_\nu$) and those induced by the pressure (e.g. 
${\boldsymbol {f}}'_p$). The total second-order hydrodynamic force is the sum of these two contributions.
Apart from the fact that the prefactors 
$6{\rm \pi}$ are replaced by 
$4{\rm \pi}$ for a bubble, (A3) differs from (3.37) through the prefactors of the second-order contributions 
${\mbox {d} {\boldsymbol {u}}_s}/{\mbox {d}t}$ (
$2{\rm \pi}$ instead of 
${16{\rm \pi} }/{3}$) and 
${\mathbb {S}} \boldsymbol{\cdot } {\boldsymbol {u}}_s$ (
${4{\rm \pi} }/{3}$ instead of 
${37{\rm \pi} }/{12}$). The significance of the latter difference is discussed in § 4. In what follows we examine the origin of the former. Indeed, once the inner solution and the short-time contribution of the kernel are combined, the 
${\mbox {d} {\boldsymbol {u}}_s}/{\mbox {d}t}$ term should provide the same added-mass force 
$- ({2{\rm \pi} }/{3}) \varepsilon ^2({\mbox {d} {\boldsymbol {u}}_s}/{\mbox {d}t})$ as that found in (4.8) for a rigid sphere, the added-mass force being independent of the dynamic boundary condition.
Since the kernel depends on the boundary condition at the particle surface only through the prefactors 
$6{\rm \pi}$ or 
$4{\rm \pi}$, one still has 
$6 {\rm \pi}{\mathbb {K}}(t) = {{\mathbb {I}}}/{\sqrt {{\rm \pi} t}}$ at short times. Inserting this asymptotic expression in (A3) yields
\begin{align} {\boldsymbol{f}}'(\varepsilon^2\ll t\ll1) ={-}4{\rm \pi} {\boldsymbol{u}}_s - \frac{8{\rm \pi}}{3} \varepsilon \int_0^t \frac{1}{\sqrt{{\rm \pi} (t-\tau)}} \frac{\mbox{d} {\boldsymbol{u}}_s}{\mbox{d}\tau}\, \mbox{d} \tau + \varepsilon^2\left( \frac{2 {\rm \pi}}{9} \frac{\mbox{d} {\boldsymbol{u}}_s}{\mbox{d}t}+\frac{4{\rm \pi}}{3} {\mathbb{S}} \boldsymbol{\cdot} {\boldsymbol{u}}_s\right) . \end{align} Surprisingly, the prefactor of the 
${\mbox {d} {\boldsymbol {u}}_s}/{\mbox {d}t}$ term is 
${2 {\rm \pi}}/{9}$ instead of 
$-{2{\rm \pi} }/{3}$ as expected. To understand this difference, one has to compare (A4) (considering from now on 
${\mathbb {S}}={\boldsymbol {0}}$) with the exact solution obtained by Gorodtsov (Reference Gorodtsov1975) and Yang & Leal (Reference Yang and Leal1991) for a bubble translating unsteadily in a fluid at rest in the limit of negligibly small advective effects, namely
\begin{equation} {\boldsymbol{f}}_{G} ={-}4{\rm \pi} \left\{{\boldsymbol{u}}_s +2\int_0^t \exp \left( \frac{9 (t-\tau)}{\varepsilon^2}\right) \mbox{erfc}\left( \sqrt{\frac{9 (t-\tau)}{\varepsilon^2}}\right) \frac{\mbox{d} {\boldsymbol{u}}_s}{\mbox{d}\tau}\, \mbox{d} \tau \right\}- \frac{2{\rm \pi}}{3} \varepsilon^2\frac{\mbox{d} {\boldsymbol{u}}_s}{\mbox{d}t} . \end{equation}
Similar to the BBO equation (2.8), the three terms on the right-hand side respectively represent the quasi-steady drag, the visco-inertial history effect and the added-mass force, the latter resulting exclusively from the pressure contribution to the total force. At first glance, (A5) confirms that the second-order prediction (A4) differs from the exact result. However, the present theory is designed to work up to large but ‘not too large’ levels of unsteadiness. In particular, it is only supposed to work for times much larger than 
$\varepsilon ^2$. In this limit, the kernel involved in (A5) is such that, at leading order,
\begin{equation} \exp \left( \frac{9 t}{\varepsilon^2}\right) \mbox{erfc}\left( \sqrt{\frac{9 t}{\varepsilon^2}}\right) \sim \frac{1}{3} \frac{\varepsilon }{\sqrt{{\rm \pi} t}}.\end{equation} To find out the next term in the expansion of this kernel for large 
$t/\varepsilon ^2$, one may start from the mathematical result 
$\int _0^\infty \{({1}/{\sqrt {{\rm \pi} x}})-\exp ( x) \mbox {erfc}( \sqrt { x}) \}{{\rm d}\kern0.06em x}=1$. Moreover, in the sense of generalized functions, one knows that 
$1=\int _{-\infty }^{+\infty }\delta (t)\, {\rm d} t$, with 
$\delta (t)$ the Dirac delta function. Then, setting 
$x=t/\zeta$ and defining the function 
$F_\zeta (t)$ such that 
$F_\zeta (t)\equiv 0$ for 
$t<0$ and 
$F_\zeta (t)=({1}/{\zeta })\{({1}/{\sqrt {{\rm \pi} t/\zeta }})-\exp ( t/\zeta ) \mbox {erfc}( \sqrt { t/\zeta })\}$ for 
$t>0$, one finds that (Boccara Reference Boccara1997, pp. 53–54)
\begin{equation} \lim_{\zeta\rightarrow0}F_\zeta (t)=\delta(t).\end{equation}
Finally, setting 
$\zeta =\varepsilon ^2/9$, one obtains
\begin{equation} \lim_{t/\varepsilon^2 \gg 1} \exp \left( \frac{9 t}{\varepsilon^2}\right) \mbox{erfc}\left( \sqrt{\frac{9 t}{\varepsilon^2}}\right) =\frac{1}{3} \varepsilon \frac{1}{\sqrt{{\rm \pi} t}} -\frac{1}{9}\varepsilon^2 \delta(t)+O(\varepsilon^3) . \end{equation}
Inserting this result into (A5), one recovers (A4). Therefore, in contrast to the case of a rigid sphere, the history force acting on a bubble contains a term of 
$O(\varepsilon ^2)$. For this reason, after a very short transient, the 
$O(\varepsilon ^2)$ correction to the force acting on the bubble results from a combination of added-mass and visco-inertial history effects. The latter is even dominant since the resulting prefactor, 
${2{\rm \pi} }/{9}$, has an opposite sign to that of the added-mass force, 
$-{2{\rm \pi} }/{3}$. Similar to the problem encountered with the terms 
${\mathbb {S}}\boldsymbol{\cdot }{\boldsymbol {u}}_s$ and 
$\boldsymbol \varOmega \times {\boldsymbol {u}}_s$ in § 4, one faces a situation in which, under a certain asymptotic condition, a term with a given mathematical structure encapsulates two effects produced by two distinct physical mechanisms. A similar situation was identified by Magnaudet (Reference Magnaudet2003) for a rigid particle moving unsteadily close to a wall: also in this case, the contribution proportional to the instantaneous relative acceleration between the particle and fluid results from a combination of added-mass and history effects (see (14) of this reference and the discussion that follows this equation).
Appendix B. Kernel in an axisymmetric straining flow
In CMM only planar linear flows were considered. However, another configuration of fundamental interest is the so-called uniaxial straining flow corresponding to the velocity-gradient tensor
\begin{equation} {\mathbb{A}} = \left(\begin{array}{@{}ccc@{}} -1 & 0 & 0\\ 0 & -1 & 0 \\ 0 & 0 & 2\\ \end{array}\right) .\end{equation} The forces experienced by a spherical rigid sphere or a bubble held fixed in this flow were considered numerically by Magnaudet et al. (Reference Magnaudet, Rivero and Fabre1995). We computed the kernel corresponding to this flow using the techniques described by CMM. Figure 2 shows how the non-zero components of 
${\mathbb {K}}$ reach their stationary value. At short times, these components behave as
\begin{align} 6{\rm \pi} {\mathbb{K}}(t) \sim \frac{1}{\sqrt{{\rm \pi} t}}\left(\begin{array}{@{}ccc@{}} 1- \dfrac{7}{10 } t +\dfrac{2}{21 }t^{2}\cdots & 0 & 0 \\ 0 & 1- \dfrac{7}{10 } t +\dfrac{2}{21 }t^{2}\cdots & 0\\ 0 & 0 & 1 + \dfrac{7}{5 } t - \dfrac{19}{210 }t^{2}\cdots\\ \end{array}\right) .\end{align}
That is, deviations from the Basset–Boussinesq limit grow initially as 
$t^{1/2}$. Up to terms of 
$O(t^{1/2})$, the short-time evolution of the kernel may be written in the simple form 
$6{\rm \pi} {\mathbb {K}}(t) \sim ({1}/{\sqrt {{\rm \pi} t}})({\mathbb {I}}+\frac {7}{10}{\mathbb {S}}t+\cdots )$. This form was shown to be universal, i.e. independent of the specific linear flow considered, by Bedeaux & Rubi (Reference Bedeaux and Rubi1987) and Miyazaki et al. (Reference Miyazaki, Bedeaux and Avalos1995). According to figure 2, the axial and radial components reach their long-term asymptote at 
$t\approx 1$ and 
$t\approx 4$, respectively. Interestingly, the radial component changes sign beyond 
$t\approx 2.2$. This means that the long-time 
$O(\varepsilon )$ correction to the drag is negative, so that the drag of a particle translating in the radial direction (along which fluid elements are compressed) is decreased compared with its creeping-flow value. In the stationary regime one has
\begin{align} 6{\rm \pi}\bar{{\mathbb{K}}} \approx\left(\begin{array}{ccc} -0.097 & 0 & 0 \\ 0 & -0.097 & 0\\ 0 & 0 & 1.260 \\ \end{array}\right).\end{align}
Since the 
$O(\varepsilon ^2)$ correction to the drag is proportional to 
$\bar {{\mathbb {K}}}\boldsymbol{\cdot }\bar {{\mathbb {K}}}$, this correction is found to be significant along the flow axis but negligibly small in the radial direction.
Figure 2. Evolution of the kernel 
${\mathbb {K}}(t)$ in a uniaxial straining flow. Solid line: axial component 
${\mathbb {K}}_3^3$; dotted line: radial components 
${\mathbb {K}}_1^1={\mathbb {K}}_2^2$.
We also computed the kernel associated with the so-called biaxial straining flow associated with the velocity-gradient tensor
\begin{equation} {\mathbb{A}} = \left(\begin{array}{@{}ccc@{}} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & -2\\ \end{array}\right) .\end{equation}At short times, the components of this kernel behave as
\begin{equation} 6{\rm \pi} {\mathbb{K}}(t) \sim \frac{1}{\sqrt{{\rm \pi} t}}\left(\begin{array}{@{}ccc@{}} 1+ \dfrac{7}{10 } t +\dfrac{2}{21 }t^{2}\cdots & 0 & 0 \\ 0 & 1+ \dfrac{7}{10 } t +\dfrac{2}{21 }t^{2}\cdots & 0\\ 0 & 0 & 1 - \dfrac{7}{5 } t - \dfrac{19}{210 }t^{2}\cdots\\ \end{array}\right) .\end{equation}
Comparing (B2) and (B5) indicates that even-order terms are identical while odd-order terms have the same magnitude but opposite signs, as expected from the universal short-time form of 
${\mathbb {K}}(t)$ mentioned above.
The corresponding steady-state kernel is found to be
\begin{equation} 6{\rm \pi}\bar{{\mathbb{K}}} \approx\left(\begin{array}{@{}ccc@{}} 0.9544 & 0 & 0 \\ 0 & 0.9544 & 0\\ 0 & 0 & -0.8622 \\ \end{array}\right) .\end{equation}
Here, according to figure 3, the axial component changes sign for 
$t\approx 0.7$, so that inertial effects decrease the axial drag at longer times. Comparing (B6) and (B3) makes it clear that the steady-state value of each component is totally different in the two flows, which underlines the nonlinear nature of advective effects.
Figure 3. Evolution of the kernel 
${\mathbb {K}}(t)$ in a biaxial straining flow. Solid line: axial component 
${\mathbb {K}}_3^3$; dotted line: radial components 
${\mathbb {K}}_1^1={\mathbb {K}}_2^2$.
Open access
