Inertial migration of a neutrally buoyant spheroid in plane Poiseuille flow

Abstract

We study the cross-stream inertial migration of a torque-free neutrally buoyant spheroid, of an arbitrary aspect ratio $\kappa$, in wall-bounded plane Poiseuille flow for small particle Reynolds numbers ($Re_p\ll 1$) and confinement ratios ($\lambda \ll 1$), with the channel Reynolds number, $Re_c = Re_p/\lambda ^2$, assumed to be arbitrary; here $\lambda =L/H$, where $L$ is the semi-major axis of the spheroid and $H$ denotes the separation between the channel walls. In the Stokes limit ($Re_p =0)$, and for $\lambda \ll 1$, a spheroid rotates along any of an infinite number of Jeffery orbits parameterized by an orbit constant $C$, while translating with a time-dependent speed along a given ambient streamline. Weak inertial effects stabilize either the spinning ($C=0$) or tumbling orbit ($C=\infty$), or both, depending on $\kappa$. The asymptotic separation of the Jeffery rotation and orbital drift time scales, from that associated with cross-stream migration, implies that migration occurs due to a Jeffery-averaged lift velocity. Although the magnitude of this averaged lift velocity depends on $\kappa$ and $C$, the shape of the lift profiles are identical to those for a sphere, regardless of $Re_c$. In particular, the equilibrium positions for a spheroid remain identical to the classical Segre–Silberberg ones for a sphere, starting off at a distance of about $0.6(H/2)$ from the channel centreline for small $Re_c$, and migrating wallward with increasing $Re_c$. For spheroids with $\kappa \sim O(1)$, the Jeffery-averaged analysis is valid for $Re_p\ll 1$; for extreme aspect ratio spheroids, the regime of validity becomes more restrictive being given by $Re_p \kappa /\ln \kappa \ll 1$ and $Re_p/\kappa \ll 1$ for $\kappa \rightarrow \infty$ (slender fibres) and $\kappa \rightarrow 0$ (flat disks), respectively.

1. Introduction

In an experimental study of pressure-driven pipe flow of a suspension of neutrally buoyant spheres, Segre & Silberberg (1962a,b) reported the migration of particles to an intermediate annular location. This was termed the tubular pinch effect to indicate that the initially uniform distribution of particles over the pipe cross-section is ‘pinched’ to a narrow annulus with increasing downstream distance. The authors performed experiments for $0.03\leq d_p/D\leq 0.15$ and for $Re$ upto 700, where $d_p$ and $D$ are the particle and pipe diameters, respectively, and $Re$ is the Reynolds number based on $D$ and the mean velocity. The aforementioned pinch occurred at $0.6\times$pipe radius for the smaller $Re$'s, and shifted towards the walls with increasing $Re$. In any unidirectional shearing flow, including pressure-driven flow in the pipe and channel geometries, reversibility of the Stokes equations prohibits cross-stream migration of a neutrally buoyant sphere, and the pinch effect above is a consequence of migration driven by inertial lift forces.

Since the original experiments of Segre and Silberberg, inertial migration in varying flow configurations has been investigated in a number of theoretical (Saffman 1965; Cox & Brenner 1968; Ho & Leal 1974; Vasseur & Cox 1976; Cox & Hsu 1977; Schonberg & Hinch 1989; Asmolov 1999; Matas, Morris & Guazzelli 2009), numerical (Chun & Ladd 2006; Shao, Yu & Sun 2008; Morita, Itano & Sugihara-Seki 2017; Nakayama et al. 2019; Pan, Li & Glowinski 2021) and experimental (Repetti & Leonard 1964; Jeffrey & Pearson 1965; Karnis, Goldsmith & Mason 1966; Tachibana 1973; Aoki, Kurosak & Anzai 1979; Matas, Morris & Guazzelli 2004; Masaeli et al. 2012; Nakayama et al. 2019) studies. From the rheological perspective, the pinch effect in the original pipe-flow experiments is undesirable, since it manifests as an apparent non-Newtonian behaviour for a dilute suspension of spheres. For instance, in a capillary viscometer, depending on the suspension flow rate, the residence time of the particles may or may not be sufficient for inertial migration to be complete; incomplete migration leads to a nonlinear dependence of the inferred viscosity on the shear rate (Segre & Silberberg 1963). On the other hand, inertial migration has recently been used as a tool in microfluidics to separate particles based on size (Di Carlo et al. 2007). Due to the robust fault-tolerant physical effects employed, and high rates of operation (this being a natural consequence of being in the inertial regime), inertial microfluidic systems are expected to have a broad range of applications in continuous bioparticle separation, cell and particle manipulation, and filtration systems (Di Carlo 2009).

A first theoretical explanation of the phenomenon was given by Ho & Leal (1974) who determined the lift force on a sphere in plane Poiseuille flow for $Re_p, Re_c\ll 1$, where $Re_p$ and $Re_c$ are the particle and channel Reynolds numbers, respectively. A more accurate calculation was performed in Vasseur & Cox (1976) using the framework developed earlier in Cox & Brenner (1968), wherein the lift velocity was expressed in terms of a volume integral involving the Green's function for creeping flow in the presence of a pair of plane boundaries (the channel walls). The resulting lift force profile had a pair of zero crossings, symmetrically located on either side of the channel centreline, and that corresponded to stable equilibrium locations. These were interpreted as the analog of the intermediate annular ring observed in the experiments, implicitly pointing to the similarity of the physics governing inertial migration in the channel and pipe geometries. Later, Schonberg & Hinch (1989) calculated the lift velocity of a sphere in plane Poiseuille flow for $Re_p\ll 1$ and for $Re_c$ upto $150$, finding the equilibrium locations to move towards the respective walls with increasing $Re_c$, consistent with the original observations. Asmolov (1999) confirmed these findings, and further extended the calculation to $Re_c=3000$. (This upper bound was regarded as reasonable based on the threshold ($Re_c\sim 11\,544$) for the Tolmein–Schlichting instability. It is now known, however, that the actual transition of plane Poiseuille flow to turbulence has a subcritical character, occurring at much lower $Re_c$'s of $O(2000)$.)

There are numerous instances in microfluidic and other settings where the particles of interest are anisotropic, for example, cancer cells (Suresh 2007), blood cells (Toner & Irimia 2005) or polymeric microstructures (Chung et al. 2008). For anisotropic particles, the inertial lift force, and any equilibria that arise as a result of this force vanishing at specific locations, are expected to depend on particle shape. Motivated by this, and following an earlier study (Hur et al. 2011) on a variety of non-spherical particle geometries, Masaeli et al. (2012) conducted experiments on spheroids of aspect ratio $1\leq \kappa \leq 5$, suspended in pressure-driven flow through a rectangular duct, for $Re_c$ upto $80$ ($Re_c$ defined based on the smaller cross-sectional dimension). The study confirmed the existence of shape-sensitive equilibria, with large aspect ratio spheroids migrating to locations near the channel centreline and those with order unity aspect ratios migrating towards the duct walls, over a range of cross-sectional aspect ratios. The authors also performed numerical simulations of rotating rods in the channel, and the equilibria obtained agreed with the experiments. Recently, Huang, Marson & Larson (2019) performed dissipative particle dynamics simulations to examine inertial migration of relatively large (compared with the channel width) prolate and oblate spheroids in plane Poiseuille flow; for moderate $Re_p$, the equilibrium positions of tumbling prolate and spinning oblate spheroids moved closer to the centreline as the spheroid aspect ratio departed further away from unity. In a later effort, Nizkaya et al. (2020) examined smaller spinning oblate spheroids with aspect ratios in the interval $0.25\leq \kappa \leq 1$ in plane Poiseuille flow, using lattice Boltzmann simulations, and in contrast to the earlier efforts above, found the lift force equilibria to be relatively insensitive to aspect ratio.

Motivated in particular by the experiments of Masaeli et al. (2012), we take a first step towards analysing the inertial migration of a freely rotating neutrally buoyant spheroid in plane Poiseuille flow. Specifically, we calculate the leading-order time-averaged lift velocity for $Re_p\ll 1$ within the framework of a point-particle approximation; $Re_c$, while much larger than $Re_p$, is otherwise arbitrary. Section 2 below presents the governing equations and boundary conditions in the context of the problem formulation. Next, in § 3 we examine the small-$Re_c$ limit where the time-averaged lift velocity is determined semi-analytically using a generalized reciprocal theorem formulation originally used by Ho & Leal (1974), and that is derived in § 3.1. Scaling arguments given in § 3.2 show that the dominant contribution to the lift velocity in this limit comes from scales of the order of the channel width $H$ that is much smaller than the inertial screening length of $O(HRe_c^{-1/2})$. Inertia therefore has a regular character, with the lift velocity being $O(Re_p)$, and its calculation requiring knowledge of only the Stokesian disturbance fields in the confined domain. These Stokesian fields are calculated in § 3.3, and are then used to obtain the time-averaged lift velocity for spheroids in § 3.4. This is followed by a presentation of the results in § 3.5. In § 4 the time-averaged lift velocity is calculated numerically for $Re_c\gtrsim O(1)$, with inertia now acting as a singular perturbation. Following Schonberg & Hinch (1989), this calculation involves a partial Fourier transform of the linearized Navier–Stokes equations, leading to coupled ordinary differential equations (ODEs) – in the transverse coordinate – for the transformed pressure and normal velocity fields, and their numerical solution using a shooting method. A brief outline of the aforementioned calculation procedure is given in § 4.1, with the results obtained discussed in § 4.2. In both §§ 3 and 4, a time averaging is necessary on account of the separation between spheroid rotation and inertial migration time scales, and at leading order in $Re_p$, is based on the Jeffery angular velocity. This, together with the fact that the dominant contributions to the lift velocity come from scales much larger than the particle size for all $Re_c$, lead to the equilibrium locations being independent of the spheroid aspect ratio, and identical to those for spheres; the magnitude of the lift force does depend on aspect ratio. In § 5 we summarize the main results, and briefly discuss their implications for shape-sorting in microfluidic settings.

2. Problem formulation

Figure 1(a) shows a neutrally buoyant spheroid of aspect ratio $\kappa =L/b$ ($L$ and $b$ are the semi-major and minor axes) freely suspended in a wall-bounded plane Poiseuille flow at a distance $d$ from the lower wall; $\kappa <1$ and >1 for oblate and prolate spheroids, respectively. The non-dimensional equations governing the velocity field $\boldsymbol {u}$ and pressure field $p$ are given by

(2.1a)$$\begin{gather} \nabla^2 \boldsymbol{u}-\boldsymbol{\nabla}p =Re_p \left(\frac{\partial\boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}\right), \end{gather}$$

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

where $\boldsymbol {u}$ satisfies the following boundary conditions:

(2.2a)$$\begin{gather} \boldsymbol{u}=\boldsymbol{\varOmega_p}\wedge \boldsymbol{r}\quad \text{for } \boldsymbol{r}\in S_p, \end{gather}$$

(2.2b)$$\begin{gather}\boldsymbol{u}\rightarrow \boldsymbol{u}^\infty\quad \text{for } r_1,r_3\rightarrow \infty\ (r_2\ \text{fixed}), \end{gather}$$

(2.2c)$$\begin{gather}\boldsymbol{u}=-\boldsymbol{U}_p\quad \text{at } r_2=-s\lambda^{-1}, (1-s)\lambda^{-1}. \end{gather}$$

In (2.1a)–(2.2c) all variables are non-dimensionalized using $L$ and the velocity scale $V_c=V_{max}L/H$, this being the order of the velocity change across the ends of the spheroid. The particle Reynolds number $Re_p=V_c L/\nu = V_{max}L^2/H\nu$, $\nu$ being the kinematic viscosity of the fluid. In (2.2a), $\boldsymbol {\varOmega }_p$ is the angular velocity of the spheroid whose surface is denoted by $S_p$. The coordinate system chosen in writing the above equations translates with the spheroid velocity $\boldsymbol {U}_p$, with its origin at the spheroid centre. Thus, $\boldsymbol {u}^\infty$ in (2.2b), the ambient plane Poiseuille flow in this frame, is given by

(2.3)\begin{equation} \boldsymbol{u}^\infty=(\alpha+\beta r_2 +\gamma {r_2}^2)\boldsymbol{1}_1-\boldsymbol{U}_p, \end{equation}

where $\alpha \boldsymbol {1}_1 - \boldsymbol {U}_p$, with $\alpha =4 \lambda ^{-1} s(1-s)$, is the ambient slip velocity at the spheroid centre, $\beta =4(1-2s)$ is the local shear rate that varies linearly across the channel and $\gamma =-4\lambda$ is the constant curvature of the plane Poiseuille profile; $s = d/H$ here being the (non-dimensional) spheroid location. Here $\lambda =L/H$ is the confinement ratio assumed small, so the channel Reynolds number $Re_c=Re_p/\lambda ^2 \gg Re_p$. As will be argued below, at leading order in $Re_p$ and $\lambda$, $\boldsymbol {U}_p$ is along the flow direction.

Figure 1. (a) A neutrally buoyant spheroid with symmetry axis $\boldsymbol {p}$ in plane Poiseuille flow. The position of the spheroid relative to the lab frame ($x_1,x_2,x_3$) is denoted by $\boldsymbol {y}$; $(r_1,r_2,r_3)$ represents the Cartesian frame with origin at the spheroid centre, and translating with it. Schematic (b) shows the body-fixed coordinate system aligned with $\boldsymbol {p}$, along with the polar ($\theta _j$) and azimuthal ($\phi _j$) angles that define the spheroid orientation. The dot-dashed line is the projection of the $\boldsymbol {r}_{b_1}$ axis on the flow-gradient plane.

We now define the disturbance fields $\boldsymbol {u}'=\boldsymbol {u}-\boldsymbol {u}^\infty$, $p'=p-p^\infty$ that satisfy the governing equations

(2.4a)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma}' = \nabla^2 \boldsymbol{u}'-\boldsymbol{\nabla} p' =Re_p\left(\frac{\partial\boldsymbol{u}'}{\partial t}+\boldsymbol{u}' \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}'+\boldsymbol{u}' \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}^\infty+\boldsymbol{u}^\infty \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}'\right), \end{gather}$$

(2.4b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'=0, \end{gather}$$

where $\boldsymbol {u}'$ satisfies

(2.5a)$$\begin{gather} \boldsymbol{u}'=\boldsymbol{U}_p+\boldsymbol{\varOmega}_p\wedge \boldsymbol{r}-(\alpha+\beta r_2 + \gamma r_2^2)\boldsymbol{1}_1\quad \text{for } \boldsymbol{r}\in S_p, \end{gather}$$

(2.5b)$$\begin{gather}\boldsymbol{u}'\rightarrow 0\quad \text{for } r_1,r_3\rightarrow \infty\ (r_2\ \text{fixed}), \end{gather}$$

(2.5c)$$\begin{gather}\boldsymbol{u}'=0\quad \text{at } r_2=-s\lambda^{-1}, (1-s)\lambda^{-1}. \end{gather}$$

Unlike a sphere, one needs to include the unsteady terms in the inertial acceleration on account of spheroid rotation even in the Stokes limit ($Re_p = 0$). We analyse the inertial migration problem defined above in the limit $\lambda, Re_p\ll 1$ with $Re_c=Re_p/\lambda ^2$ arbitrary. It is also assumed that $s$, $(1-s)\gg \lambda$, implying that the analysis is restricted to the spheroid being at distances from either wall that are much larger than $O(L)$.

Before embarking on a detailed analysis for weak inertia, it is worth summarizing the nature of spheroid motion in plane Poiseuille flow in the Stokes limit. Since the local linear flow approximation for the plane Poiseuille profile is simple shear flow, the neutrally buoyant spheroid must rotate along Jeffery orbits (Jeffery 1922), this rotation being characterized by the polar ($\theta _j$) and azimuthal ($\phi _j$) angles of the spheroid symmetry axis (see figure 1b) being functions of time; these equations are given later in § 3 (see (3.19a,b)). Jeffery rotation is best described in $(C,\tau$) coordinates, where $C\in [0,\infty )$ is the orbit constant and $\tau$ is the orbit phase that changes at a constant albeit $\kappa$-dependent rate from $0$ to $2{\rm \pi}$ over a single period (Leal & Hinch 1971; Dabade, Marath & Subramanian 2016); the two limiting orbits correspond to the spinning ($C=0$) and tumbling ($C=\infty$) modes. The Jeffery period is independent of $C$, being given by $T_{jeff}=2{\rm \pi} (\beta V_{max}/H)^{-1} (\kappa +\kappa ^{-1})$, where $\beta$ accounts for the linearly varying shear rate. The curvature of the plane Poiseuille profile only affects the spheroid translation velocity that is now a function of its orientation, and thence, of time. A neutrally buoyant spheroid in an unbounded plane Poiseuille flow continues to move along a given ambient streamline, with a speed that ranges from a maximum, corresponding to the flow-aligned orientation ($\phi _j=0$, ${\rm \pi}$ and ${{\rm \pi} }/{2}$, ${3{\rm \pi} }/{2}$ for $\kappa >1$ and <1, respectively), to a minimum when the spheroid is aligned orthogonal to the flow direction ($\phi _j={{\rm \pi} }/{2}$, ${3{\rm \pi} }/{2}$ and $0$, ${\rm \pi}$ for $\kappa >1$ and ${<}1$, respectively) (Chwang 1975). For very large (small) $\kappa$, the spheroid spends an increasing fraction of a Jeffery period in the flow-aligned orientation, and the translational motion thereby acquires an increasingly jerky character owing to the spheroid abruptly slowing down during the brief periods of misalignment.

In presence of walls, the leading-order correction to the motion above arises from interaction of the spheroid with time-dependent image stresslets induced by each wall, resulting in an $O(\lambda ^2)$ lateral velocity component even in the Stokes limit, with an additional $O(\lambda ^3)$ correction to the Jeffery angular velocity. Note that this image-stresslet interaction does not lead to transverse motion for a sphere, as may be seen from the stresslet orientation (along the local extensional axis) and the associated purely radial velocity field; as already mentioned in § 1, this must be so, independent of $\lambda$, due to reversibility constraints. Thus, for $\lambda \ll 1$, the centre of mass of a neutrally buoyant spheroid in wall-bounded plane Poiseuille flow exhibits a small-amplitude oscillatory motion about an ambient streamline, the amplitude being $O(\lambda ^2)$. In other words, unlike a sphere, Stokesian reversibility does not preclude an instantaneous lift force for an anisotropic particle. Evidence for such oscillating-cum-tumbling spheroid trajectories, for finite $\lambda$, is available from earlier computations, with there being a transition from tumbling (rotation) to angular oscillations beyond a threshold $\lambda$ close to unity, corresponding to sufficiently narrow channels (Sugihara-Seki 1993, 1996; Staben, Zinchenko & Davis 2003, 2006). However, Stokesian reversibility still forbids a net migration of the spheroid in the transverse direction, and to allow for such a motion, one needs inertia. It is this net cross-stream migration, for small $Re_p$ and $\lambda$, that is analysed in the following sections. We obtain the time-averaged motion of a spheroid in this limit, and since the wall and inertia-induced modification of the primary translational and rotational motion are weak in the aforementioned limit, the time average corresponds to an average over a Jeffery period.

It is well known that the Stokes equations do not provide a uniformly valid leading-order approximation for small but finite $Re_p$, and that in general one requires a matched asymptotic expansion approach (Proudman & Pearson 1957) involving an inner expansion in the neighbourhood of the particle, and an outer expansion at distances of the order of an inertial screening length, to calculate inertial corrections. The screening length for the present problem is $L Re_p^{-1/2}$ (Saffman 1965), or equivalently, $HRe_c^{-1/2}$. For small $Re_c$, the screening length is larger than $H$, so the channel walls lie in the inner Stokesian region where fluid inertia may be treated as a regular perturbation. The inertial migration problem for a sphere in this limit was investigated by Ho & Leal (1974), and in a series of papers by Brenner, Cox and collaborators (Cox & Brenner 1968; Vasseur & Cox 1976; Cox & Hsu 1977). We examine the analogous problem for a spheroid in this limit in § 3 below. When $Re_c\gtrsim O(1)$, the inertial screening length is of $O(H)$ or smaller, and the solution procedure involves the outer expansion, the leading-order term of which satisfies the linearized Navier–Stokes equations. This limit was first examined for a sphere in Schonberg & Hinch (1989), and we examine the same for a spheroid in § 4.

3. The inertial lift velocity for $Re_c\ll 1$

Although inertial effects are a regular perturbation for small $Re_c$, following Ho & Leal (1974), we do not calculate this perturbation explicitly, and instead use the generalized reciprocal theorem that relates the velocity and stress fields of the problem of interest with those of a simpler test problem whose solution is known. The problem of interest ($\boldsymbol {u}',\boldsymbol {\sigma }'$) corresponds to the motion of a torque-free neutrally buoyant spheroid in a wall-bounded plane Poiseuille flow, taking into account the inertial acceleration of the suspending fluid. Since the quantity of interest is the inertial lift velocity, the test problem ($\boldsymbol {u}^t,\boldsymbol {\sigma }^t$) is taken to be a torque-free spheroid, in a quiescent ambient between parallel walls, acted on by an arbitrarily oriented unit force; the test spheroid has the same instantaneous orientation as the one in the actual problem. The governing equations and boundary conditions for the actual problem have already been given in (2.4a,b) and (2.5a–c). Those for the test problem are

(3.1a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}^t=\nabla^2 \boldsymbol{u}^t-\boldsymbol{\nabla} p^t = 0 , \end{gather}$$

(3.1b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^t=0, \end{gather}$$

with the boundary conditions

(3.2a)$$\begin{gather} \boldsymbol{u}^t=\boldsymbol{U}_p^t+\boldsymbol{\varOmega}_p^t\wedge\boldsymbol{r}\quad \text{for } \boldsymbol{r}\in S_p, \end{gather}$$

(3.2b)$$\begin{gather}\boldsymbol{u}^t\rightarrow 0\quad \text{for } r_1,r_3 \rightarrow \infty\ (r_2\ \text{fixed}), \end{gather}$$

(3.2c)$$\begin{gather}\boldsymbol{u}^t=0\quad \text{at } r_2=-s\lambda^{-1}, (1-s)\lambda^{-1}, \end{gather}$$

where $\boldsymbol {U}_p^t$ and $\boldsymbol {\varOmega }_p^t$ are the spheroid translational and angular velocities.

3.1. The generalized reciprocal theorem

To derive the generalized reciprocal theorem identity, we contract (2.4a) with $\boldsymbol {u}^t$ and (3.1a) with $\boldsymbol {u}'$, and subtract the resulting expressions (Ho & Leal 1974) to obtain

(3.3)\begin{equation} (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}')\boldsymbol{\cdot} \boldsymbol{u}^t-(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}^t) \boldsymbol{\cdot} \boldsymbol{u}'=Re_p\,\boldsymbol{u}^t\boldsymbol{\cdot}\boldsymbol{f}, \end{equation}

where $\boldsymbol {f}=({\partial \boldsymbol {u}'}/{\partial t}+\boldsymbol {u}'\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}'+\boldsymbol {u}'\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}^\infty +\boldsymbol {u}^\infty \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}')$. Integrating (3.3) over the fluid volume ($V^F$) between the channel walls gives

(3.4)\begin{align} \int_{V^F} \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\sigma}'\boldsymbol{\cdot} \boldsymbol{u}^t-\boldsymbol{\sigma}^t\boldsymbol{\cdot}\boldsymbol{u}')\,{\rm d}V + \int_{V^F} (\boldsymbol{\sigma}^t\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}'-\boldsymbol{\sigma}' \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^t)\,{\rm d}V = Re_p\int_{V^F}\boldsymbol{u}^t\boldsymbol{\cdot}\boldsymbol{f}\,{\rm d}V. \end{align}

Using (2.4b) and (3.1b), along with the definitions of $\boldsymbol {\sigma }'$ and $\boldsymbol {\sigma }^t$, the second integral on the left-hand side of (3.4) can be shown to be identically zero. Applying the divergence theorem to the first integral yields

(3.5)\begin{equation} -\int_{S_p+S_w+S_\infty} (\boldsymbol{\sigma}'\boldsymbol{\cdot}\boldsymbol{u}^t-\boldsymbol{\sigma}^t\boldsymbol{\cdot} \boldsymbol{u}')\boldsymbol{\cdot}\boldsymbol{n}\,{\rm d}S = Re_p\int_{V^F}\boldsymbol{u}^t\boldsymbol{\cdot}\boldsymbol{f} \,{\rm d}V, \end{equation}

where $\boldsymbol {n}$ is the unit normal to all bounding surfaces pointing into the fluid domain. The set of bounding surfaces include the particle surface ($S_p$), the channel walls ($S_w$) and the surface at infinity ($S_\infty$); the latter can be thought of as the curved surface of a cylinder of radius $R$, with its axis along the gradient direction, in the limit $R\to \infty$. For $R\gg H$, the disturbance fields in the actual and test problems are exponentially small (see (A25) and (3.32)), and therefore, the integral over $S_\infty$ in (3.5) is vanishingly small. The integral $\int _{S_w} (\boldsymbol {\sigma }'\boldsymbol {\cdot }\boldsymbol {u}^t-\boldsymbol {\sigma }^t\boldsymbol {\cdot }\boldsymbol {u}')\boldsymbol {\cdot } \boldsymbol {n}\,{\rm d}S$ can be shown to be identically zero on using the no-slip conditions (2.5c) and (3.2c) on the channel walls. Thus, (3.5) reduces to

(3.6)\begin{equation} -\int_{S_p} (\boldsymbol{\sigma}'\boldsymbol{\cdot}\boldsymbol{u}^t-\boldsymbol{\sigma}^t\boldsymbol{\cdot} \boldsymbol{u}')\boldsymbol{\cdot}\boldsymbol{n}\,{\rm d}S =Re_p\int_{V^F}\boldsymbol{u}^t \boldsymbol{\cdot}\boldsymbol{f}\,{\rm d}V. \end{equation}

Applying the no-slip conditions (2.5a) and (3.2a) to (3.6),

(3.7)\begin{align} &-\boldsymbol{U}_p^t\boldsymbol{\cdot}\int_{S_p} \boldsymbol{\sigma}'\boldsymbol{\cdot}\boldsymbol{n}\,{\rm d}S - \boldsymbol{\varOmega}_p^t\boldsymbol{\cdot}\int_{S_p}\boldsymbol{r}\wedge(\boldsymbol{\sigma}' \boldsymbol{\cdot}\boldsymbol{n})\,{\rm d}S + \boldsymbol{\varOmega}_p \boldsymbol{\cdot}\int_{S_p} \boldsymbol{r}\wedge(\boldsymbol{\sigma}^t\boldsymbol{\cdot} \boldsymbol{n})\,{\rm d}S\nonumber\\ &\quad +\boldsymbol{U}_p\boldsymbol{\cdot}\int_{S_p} \boldsymbol{\sigma}^t\boldsymbol{\cdot} \boldsymbol{n}\,{\rm d}S-\int_{S_p} (\boldsymbol{\sigma}^t\boldsymbol{\cdot}\boldsymbol{n})\boldsymbol{\cdot}(\alpha+\beta r_2+\gamma r_2^2)\boldsymbol{1}_1\,{\rm d}S=Re_p\int_{V^F}\boldsymbol{u}^t\boldsymbol{\cdot}\boldsymbol{f} \,{\rm d}V. \end{align}

The hydrodynamic force and torque experienced by the particle may be written as $\int _{S_p} \boldsymbol {\sigma }'\boldsymbol {\cdot }\boldsymbol {n}\,{\rm d}S=Re_p\,{\rm d}\boldsymbol {U}_p/{\rm d}t$ and $\int _{S_p} \boldsymbol {r}\wedge (\boldsymbol {\sigma }'\boldsymbol {\cdot }\boldsymbol {n})\,{\rm d}S=Re_p\, {\rm d}(\boldsymbol {I}_p\boldsymbol {\cdot }\boldsymbol {\varOmega }_p)/{\rm d}t$, $\boldsymbol {I}_p$ being the spheroid moment of inertia tensor. Next, to $O(Re_p)$, one may replace $\boldsymbol {u}'$ by $\boldsymbol {u}_s$ in the volume integral in (3.7), $\boldsymbol {u}_s$ being the corresponding Stokesian approximation. The validity of such a replacement requires the resulting volume integral to be convergent, this being related to inertia acting as a regular perturbation. Detailed arguments in this regard are given in the next subsection. Furthermore, noting that $\int _{S_p} \boldsymbol {r}\wedge (\boldsymbol {\sigma }^t\boldsymbol {\cdot }\boldsymbol {n})\,{\rm d}S = 0$ owing to the spheroid being torque free in the test problem, one obtains

(3.8)\begin{align} &-Re_p\,\boldsymbol{U}_p^t\boldsymbol{\cdot}\frac{{\rm d}\boldsymbol{U}_p}{{\rm d}t} -Re_p\, \boldsymbol{\varOmega}_p^t\boldsymbol{\cdot}\frac{{\rm d}(\boldsymbol{I}_p\boldsymbol{\cdot} \boldsymbol{\varOmega}_p)}{{\rm d}t}+\boldsymbol{U}_p\boldsymbol{\cdot}\int_{S_p} \boldsymbol{\sigma}^t\boldsymbol{\cdot}\boldsymbol{n}\,{\rm d}S\nonumber\\ &\quad-\int_{S_p} (\boldsymbol{\sigma}^t\boldsymbol{\cdot}\boldsymbol{n})\boldsymbol{\cdot}(\alpha+\beta r_2+\gamma r_2^2)\boldsymbol{1}_1\,{\rm d}S=Re_p\int_{V^F}\boldsymbol{u}^t \boldsymbol{\cdot}\boldsymbol{f}_s\,{\rm d}V, \end{align}

where $\boldsymbol {f}_s$ denotes the approximation of $\boldsymbol {f}$ based on replacing $\boldsymbol {u}'$ by $\boldsymbol {u}_s$, being given by

(3.9)\begin{equation} \boldsymbol{f}_s=\frac{\partial\boldsymbol{u}_s}{\partial t}+\boldsymbol{u}_s\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_s+\boldsymbol{u}_s\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}^\infty+\boldsymbol{u}^\infty\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_s. \end{equation}

While $\boldsymbol {u}_s$ above is the Stokesian disturbance field induced by an arbitrarily oriented neutrally buoyant spheroid in a bounded plane Poiseuille flow, as will be seen in § 3.2, at leading order in $\lambda$, the spheroid may be replaced by the corresponding stresslet singularity, with the quadrupolar component of the disturbance field (that arises in response to the quadratic part of the ambient flow) being neglected.

One may now use the small-$Re_p$ expansions, $\boldsymbol {U}_p=\boldsymbol {U}_{p0}+Re_p \boldsymbol {U}_{p1}+O(Re_p^2)$ and $\boldsymbol {\varOmega }_p=\boldsymbol {\varOmega }_{p0}+Re_p \boldsymbol {\varOmega }_{p1}+O(Re_p^2)$, for the spheroid translational and angular velocities, in (3.8), and obtain the following relations at successive orders in $Re_p$.

1. (i) $O(1)$:

   (3.10)\begin{equation} \boldsymbol{U}_{p0}\boldsymbol{\cdot}\int_{S_p} \boldsymbol{\sigma}^t\boldsymbol{\cdot} \boldsymbol{n}\,{\rm d}S=\int_{S_p} (\boldsymbol{\sigma}^t\boldsymbol{\cdot} \boldsymbol{n})\boldsymbol{\cdot}(\alpha+\beta r_2+\gamma r_2^2) \boldsymbol{1}_1\,{\rm d}S; \end{equation}

2. (ii) $O(Re_p)$:

   (3.11)\begin{equation} -\boldsymbol{U}_p^t\boldsymbol{\cdot}\frac{{\rm d}\boldsymbol{U}_{p0}}{{\rm d}t}- \boldsymbol{\varOmega}_p^t\boldsymbol{\cdot}\frac{{\rm d}(\boldsymbol{I}_p\boldsymbol{\cdot} \boldsymbol{\varOmega}_{p0})}{{\rm d}t}+\boldsymbol{U}_{p1}\boldsymbol{\cdot} \int_{S_p} \boldsymbol{\sigma}^t\boldsymbol{\cdot}\boldsymbol{n}\,{\rm d}S= \int_{V^F}\boldsymbol{u}^t\boldsymbol{\cdot}\boldsymbol{f}_s\,{\rm d}V .\end{equation}

3.1.1. The $O(1)$ problem

For the $O(1)$ problem defined by (3.10), we choose $\int _{S_p} \boldsymbol {\sigma }^t\boldsymbol {\cdot } \boldsymbol {n}\,{\rm d}S=\boldsymbol {1}_1$, corresponding to the test spheroid translating due to a unit flow-aligned force. In the absence of boundaries, the induced velocity field is known in terms of vector spheroidal harmonics (Dabade, Marath & Subramanian 2015), and the surface force density to be used on the right-hand side is given by (Kushch & Sangani 2003)

(3.12)\begin{equation} \boldsymbol{\sigma}^t\boldsymbol{\cdot} \boldsymbol{n} =-P^t \boldsymbol{1}_\xi+2\left(\frac{(\xi^2-1)^{1/2}}{\xi_0 (\xi^2-\eta^2)^{1/2}} \frac{\partial \boldsymbol{u}^t}{\partial \xi}+\frac{1}{2}\boldsymbol{1}_\xi\wedge\boldsymbol{\nabla}\wedge\boldsymbol{u}^t\right), \end{equation}

where $\xi \geq \xi _0>1$ is the coordinate characterizing a family of confocal spheroids with $\xi _0$ representing the spheroid surface and $\boldsymbol {1}_\xi$ is the unit normal to the spheroid surface; note that $\kappa =\xi _0/(\xi _0^2-1)^{1/2}$ and $\kappa =(\xi _0^2-1)^{1/2}/\xi _0$ for prolate and oblate spheroids, respectively. Evaluating the surface integral in (3.10), using (3.12), one obtains

(3.13)\begin{equation} \boldsymbol{U}_{p0}= \left[ \alpha+\frac{\gamma}{3\kappa^2}[\cos^2\phi_j+(\cos^2\theta_j+\kappa^2\sin^2\theta_j)\sin^2\phi_j] \right]\boldsymbol{1}_1, \end{equation}

which is Faxen's law for the translation of a force-free spheroid truncated to second order, the term proportional to $\gamma$ arising from the curvature of the ambient flow; higher-order terms in the expansion are zero for an ambient quadratic flow (Happel & Brenner 2012). (3.13) pertains to a spheroid of given orientation, and a time trajectory requires knowing $\theta _j$ and $\phi _j$ as functions of time. For this purpose, one may use Faxen's law relating the torque and angular velocity of a spheroid. This is again an infinite expansion in the general case, but with only the leading-order term, proportional to the ambient velocity gradient, surviving for a quadratic flow. Thus, spheroid rotation remains identical to that in an ambient linear flow, and application of the torque-free constraint yields

(3.14)\begin{align} \boldsymbol{\varOmega}_{p0}&=-\frac{\beta(\kappa^2-1)}{4(\kappa^2+1)}\cos\phi_j\sin 2\theta_j\boldsymbol{1}_1+ \frac{\beta(\kappa^2-1)}{4(\kappa^2+1)}\sin\phi_j\sin 2\theta_j\boldsymbol{1}_2\nonumber\\ &\quad +\frac{\beta}{2}\left[-1+\frac{(\kappa^2-1)}{(\kappa^2+1)}\cos 2\phi_j\sin^2 \theta_j\right] \boldsymbol{1}_3. \end{align}

The equations governing Jeffery orbits arise as components of (3.14) in the body-aligned coordinate system, and are given below in (3.18a,b). With $\theta _j$ and $\phi _j$ defined in this manner, $\boldsymbol {U}_{p0}$ as defined in (3.13) exhibits the time dependence described earlier in § 2. Substituting $\boldsymbol {U}_{p0}$ in (2.3), the ambient flow in the reference frame chosen for the reciprocal theorem is given by

(3.15)\begin{equation} \boldsymbol{u}^\infty=\left[ (\beta r_2 +\gamma r_2^2)-\frac{\gamma}{3\kappa^2}[\cos^2\phi_j+(\cos^2\theta_j+\kappa^2\sin^2\theta_j)\sin^2\phi_j] \right] \boldsymbol{1}_1. \end{equation}

Note that the $O(\lambda ^2)$ centre-of-mass oscillations mentioned in § 2 can only arise from $\boldsymbol {U}_{p0}$ having a component along the gradient direction, and this requires an expression for $\boldsymbol {\sigma }^t \boldsymbol {\cdot } \boldsymbol {n}$ that includes the influence of the plane boundaries.

3.1.2. The $O(Re_p)$ problem

At $O(Re_p)$, one chooses the test force as $\int _{S_p} \boldsymbol {\sigma }^t\boldsymbol {\cdot }\boldsymbol {n}\, {\rm d}S=-\boldsymbol {1}_2$, which leads to

(3.16)\begin{align} V_p=Re_p\,\boldsymbol{U}_{p1}\boldsymbol{\cdot}\boldsymbol{1}_2 &=-Re_p\int_{V^F}\boldsymbol{u}^{t2}\boldsymbol{\cdot}\left(\frac{\partial\boldsymbol{u}_s}{\partial t}+\boldsymbol{u}_s\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_s+\boldsymbol{u}_s \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^\infty+\boldsymbol{u}^\infty \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_s\right) {\rm d}V\nonumber\\ &\quad -Re_p\,\boldsymbol{U}_p^{t2}\boldsymbol{\cdot} \frac{{\rm d}\boldsymbol{U}_{p0}}{{\rm d}t}-Re_p\,\boldsymbol{\varOmega}_p^{t2} \boldsymbol{\cdot}\frac{{\rm d}(\boldsymbol{I}_p\boldsymbol{\cdot}\boldsymbol{\varOmega}_{p0})}{{\rm d}t} \end{align}

for the instantaneous lift velocity of a neutrally buoyant spheroid in plane Poiseuille flow, where the volume integral is a function of the changing spheroid orientation via $\boldsymbol {u}^{t2}$ and $\boldsymbol {u}_s$. The additional superscript ‘$2$’ for the test problem quantities indicate the choice of the test force orientation above.

As already mentioned, for small confinement ratios, a spheroid in wall-bounded plane Poiseuille flow rotates along Jeffery orbits (Jeffery 1922) in the Stokes limit. For small but finite $Re_p$, there is an additional $O(Re_p)$ drift across orbits that stabilizes the tumbling mode ($C=\infty$) for prolate spheroids, the spinning mode ($C=0$) for oblate spheroids with $0.14<\kappa <1$ and, depending on the initial orientation, either the spinning or tumbling mode for oblate spheroids with $\kappa <0.14$ (Einarsson et al. 2015; Dabade et al. 2016). The time scale for rotation along a Jeffery orbit is $T_{jeff}\sim O(H/V_{max})$ for $\kappa \sim O(1)$. The inertia-driven orbital drift above occurs on longer time scales $t_{drift} \sim O(Re_p^{-1}H/V_{max})$. Scaling arguments in § 3.2 below show that the volume integral in (3.16) is the dominant contribution to the cross-stream migration, being $O(1)$ for $\lambda \ll 1$, so the time scale for inertial migration is $t_{lift} \sim O(Re_p^{-1}\lambda ^{-1}H/V_{max})$. Since $t_{lift}/t_{drift}\sim \lambda ^{-1}\gg 1$, for purposes of the migration calculation, one may assume that the spheroid has settled into its stable Jeffery orbit. Moreover, since $t_{lift}/T_{jeff}\sim Re_p^{-1}\lambda ^{-1} \gg 1$, the leading-order migration is the result of an orientation-averaged lift velocity, the average being over the orientations sampled in the stabilized Jeffery orbit. As indicated below, this average may be approximated based on the Jeffery angular velocity for small $Re_p$. Superposed on this orientation-averaged migration trajectory would be small-amplitude oscillations of $O(Re_p\lambda )$, corresponding to the fluctuations in the instantaneous lift velocity arising from the rapidly changing spheroid orientation, as may be established formally using the method of multiple scales.

Owing to the aforementioned time scale separation in the limit $\lambda, Re_p \ll 1$, it is of interest to determine the Jeffery-averaged rather than the instantaneous lift velocity. Averaging both sides of (3.16) over a Jeffery period, one obtains

(3.17)\begin{equation} \langle V_p\rangle=-Re_p\int_{V^F}\left\langle\boldsymbol{u}^{t2}\boldsymbol{\cdot} \left(\frac{\partial\boldsymbol{u}_s}{\partial t}+\boldsymbol{u}_s\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_s+\boldsymbol{u}_s\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^\infty+\boldsymbol{u}^\infty \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}_s\right)\right\rangle {\rm d}V, \end{equation}

where the averaging operation is defined as $\langle.\rangle = ({1}/{T_{jeff}})\int _0^{T_{jeff}} (.)\,{\rm d}t$, and pertains to the inertially stabilized Jeffery orbit (that is, $C$ is either $0$ or $\infty$) for times longer than $O(Re_p^{-1}H/V_{max})$.

In writing down (3.17), we have neglected the terms on the right-hand side of (3.16) involving the translational and angular accelerations in the actual problem in the Stokesian limit. These depend on $\lambda$ at leading order, and are therefore asymptotically smaller than the Jeffery-averaged volume integral that, as mentioned above (see the paragraph following (3.16)), is independent of $\lambda$ for $\lambda \ll 1$. For a spinning spheroid, both acceleration terms are trivially zero because, independent of $\lambda$, a spinning spheroid (like a sphere) translates with a constant speed along an ambient streamline, while rotating at a uniform rate about the ambient vorticity direction. For a spheroid rotating in any other Jeffery orbit, both terms do lead to non-zero Jeffery-averaged contributions. The first of these terms is non-zero because of a correlation between the time-periodic variations of ${{\rm d}\boldsymbol {U}_{p0}}/{{\rm d}t}$ along the flow direction, and the analogous variation of the flow-directed component of $\boldsymbol {U}_p^{t2}$, the latter arising due to the periodically varying test spheroid orientation for a fixed force (along $-\boldsymbol {1}_2$). The acceleration in the actual problem is only $O(\lambda ^2)$, however, as may be seen from the Faxen's translation law given earlier. The smallness of the second term is because $\boldsymbol {\varOmega }_p^{t2}$, at leading order, is driven by the velocity gradient associated with the equal and oppositely directed image Stokeslet, and is again $O(\lambda ^2)$; the angular acceleration in the actual problem, ${{\rm d}(\boldsymbol {I}_p\boldsymbol {\cdot }\boldsymbol {\varOmega }_{p0})}/{{\rm d}t}$, is $O(1)$. Thus, both acceleration terms in (3.16) are $O(\lambda ^2)$ smaller than the volume integral.

The averaging operation in (3.17) corresponds to a fixed $C$, and is naturally evaluated in $(C,\tau )$ coordinates (Leal & Hinch 1971; Dabade et al. 2016). Here, $C={\tan \theta _j(\kappa ^2\sin ^2\phi _j+\cos ^2\phi _j)^{1/2}}/{\kappa }$ and $\tau =\tan ^{-1}[{1}/{\kappa \tan \phi _j}]$, where $\theta _j$ and $\phi _j$ are the polar and azimuthal angles characterizing the spheroid orientation. The latter is denoted by the unit vector $\boldsymbol {p}$ in figure 1(b) with $\boldsymbol {p}=\sin \theta _j\cos \phi _j\boldsymbol {1}_1+\sin \theta _j\sin \phi _j\boldsymbol {1}_2+\cos \theta _j\boldsymbol {1}_3$. The equations

(3.18a)$$\begin{gather} \frac{{\rm d}\phi_j}{{\rm d}t} =\beta \left(-\frac{1}{2}+\frac{\kappa^2-1}{2 (\kappa^2+1)}\cos 2\phi_j\right), \end{gather}$$

(3.18b)$$\begin{gather}\frac{{\rm d}\theta_j}{{\rm d}t} =\beta\frac{(\kappa^2-1)}{4(\kappa^2+1)}\sin 2\theta_j\sin 2\phi_j \end{gather}$$

describe rotation along Jeffery orbits, and are obtained as the individual Cartesian components of the Faxen angular velocity relation, (3.14), given above ($\boldsymbol {\varOmega }_{p0} \boldsymbol {\cdot } \boldsymbol {1}_{r_{b_2}} = {{\rm d}\theta _j}/{{\rm d}t}$, $\boldsymbol {\varOmega }_{p0} \boldsymbol {\cdot } \boldsymbol {1}_{r_{b_1}} = -\sin \theta _j({{\rm d}\phi _j}/{{\rm d}t})$). In terms of $C$ and $\tau$, (3.18a) and (3.18b) take the form ${{\rm d}C}/{{\rm d}t}=0$ and ${{\rm d}\tau }/{{\rm d}t}={\beta }/{\kappa +\kappa ^{-1}}$. The former must be the case by definition, while the latter may be used to transform the time-averaged integral in (3.17) into a $\tau$-averaged one. Inverting the definitions above, one may write $\theta _j$ and $\phi _j$ as

(3.19a)$$\begin{gather} \phi_j=\sin^{-1}\left[\frac{\cos\tau}{(\kappa^2-(\kappa^2-1) \cos^2\tau)^{1/2}}\right], \end{gather}$$

(3.19b)$$\begin{gather}\theta_j=\cos^{-1}\left[\frac{1}{(1+C^2\kappa^2-C^2(\kappa^2-1)\cos^2\tau)^{1/2}}\right], \end{gather}$$

which are used below in evaluating $\langle V_p \rangle$.

For $Re_p$ fixed, the Jeffery-averaged description of inertial migration, given by (3.17), breaks down for sufficiently small or large $\kappa$. This is because the leading-order Jeffery angular velocity becomes small in these limits, being $O(\kappa ^{-2})$ for fibres ($\kappa \gg 1$) close to flow alignment, and $O(\kappa ^2)$ for flat disks ($\kappa \ll 1$) close to alignment with the gradient–vorticity plane. As a result, the $O(Re_p)$ inertial correction becomes comparable in magnitude to the Jeffery contribution, leading to a slow down and eventual arrest of rotation (Marath & Subramanian 2017), first shown by Subramanian & Koch (2005) for the case of a slender fibre. Herein, we confine ourselves to analysing the Jeffery-averaged approximation, only noting that it remains valid provided $Re_p\,\kappa /\ln \kappa \ll 1\ (Re_p/\kappa \ll 1)$ for $\kappa \gg 1 (\kappa \ll 1)$, an increasingly restrictive assumption for extreme aspect ratio particles (Subramanian & Koch 2005; Marath & Subramanian 2017). The consequence of an inertia-induced slow down at leading order, for the said particles, will be analysed in a later communication.

A final point worth mentioning is that, for thin oblate spheroids with $\kappa < 0.14$, the aforementioned inertial drift time scale of $O(Re_p^{-1}H/V_{max})$ only applies in the absence of stochastic orientation fluctuations. In the presence of such fluctuations, either of a thermal origin or otherwise (for instance, due to pair-hydrodynamic interactions or weak turbulence; see Subramanian & Marath 2022), there is a barrier-hopping time associated with the eventual equilibration between the numbers of spinning and tumbling spheroids in a manner independent of the initial orientation distribution (Dabade et al. 2016; Marath, Dwivedi & Subramanian 2017; Marath & Subramanian 2018). This time scale increases exponentially with decreasing amplitude of the fluctuations, becoming much longer than the nominal drift time scale, and likely comparable to the time scale for inertial migration. Under these conditions, the migration dynamics will have a probabilistic rather than deterministic character, being described by a kinetic equation for the probability density that is a function of both $C$ and the transverse channel coordinate ($s$) – we briefly revisit this issue when calculating the lift velocity for arbitrary $C$ later in this section.

3.2. Scaling analysis and the point-particle formulation

The dominant contribution to the volume integral in (3.17), for $\lambda \ll 1$, can arise from either scales of $O(L)$ or those of $O(H)$, the inertial screening length being irrelevant in the small-$Re_c$ limit. In order to assess the relative importance of these contributions, we consider the intermediate asymptotic interval $1\ll r\ll \lambda ^{-1}$ ($r$ is measured in units of $L$), in which case $\boldsymbol {u}^{t2}\sim 1/r$, $\boldsymbol {u}_s\sim \beta /r^2+\gamma /r^3$, corresponding to the Stokeslet scaling for the test velocity field, and the stresslet-cum-force-quadrupole scaling for the velocity field in the actual problem. Using these $r$ scalings along with $\boldsymbol {u}^\infty \sim \beta r + \gamma r^2$ for the ambient flow, and separating the estimates for the linear and nonlinear parts of the integrand in (3.17), one obtains

1. (i) $\boldsymbol {u}^{t2}\boldsymbol {\cdot }({\partial \boldsymbol {u}_s}/{\partial t}+ \boldsymbol {u}_s\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u}^\infty + \boldsymbol {u}^\infty \boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}_s)\sim {\beta ^2}/{r^3}+{\beta \gamma }/{r^2}+{\gamma \beta }/{r^4}+{\gamma ^2}/{r^3}$ (linear in $\boldsymbol {u}_s$),

2. (ii) $\boldsymbol {u}^{t2}\boldsymbol {\cdot }(\boldsymbol {u}_s\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}_s)\sim {\beta ^2}/{r^6}+{\beta \gamma }/{r^7}+{\gamma ^2}/{r^8}$ (nonlinear in $\boldsymbol {u}_s$).

Note that, over the range of scales under consideration, wall effects only contribute at a smaller order in $\lambda$, and hence, use of the unbounded domain estimates above. Next, using ${\rm d}V\sim O(r^2\,{\rm d}r)$, one obtains the following estimates for the contributions of the different terms to the lift velocity integral:

(3.20a)$$\begin{gather} V_p^{linear}\sim Re_p\left(\beta^2\ln r+\beta\gamma r+\dfrac{\gamma\beta}{r}+\gamma^2\ln r\right), \end{gather}$$

(3.20b)$$\begin{gather}V_p^{nonlinear}\sim Re_p\left(\frac{\beta^2}{r^3}+\dfrac{\beta\gamma}{r^4}+\dfrac{\gamma^2}{r^5}\right). \end{gather}$$

The algebraically growing contribution in (3.20a) will be dominated by scales of $O(H)$ ($r\sim \lambda ^{-1}$, the outer region), while contributions in (3.20a,b) that decay with $r$ will be dominated by length scales of $O(L)$ ($r\sim O(1)$, the inner region). The algebraic growth with $r$ will be cutoff for $r\gtrsim \lambda ^{-1}$ from the more rapid decay of the disturbance velocity fields due to wall-induced screening (recall that this more rapid decay was used in neglecting the integral over $S_\infty$ in (3.5)). Likewise, the apparent divergence for $r\to 0$, for the algebraically decaying terms, will be cutoff at $r \sim O(1)$ by the finite size of the spheroid. The terms proportional to $\ln r$ in (3.20a) imply the dominance of the matching interval $1 \ll r \ll \lambda ^{-1}$, resulting in a leading-order contribution proportional to $\ln \lambda ^{-1}$, with logarithmically smaller contributions arising from the inner and outer regions. Use of these cutoffs leads to the estimates

(3.21a)$$\begin{gather} V_p^{linear}\sim Re_p\left[\beta^2(O(1)+\ln\lambda^{-1})+\beta\gamma\lambda^{-1}+\gamma\beta+ \gamma^2(O(1)+\ln\lambda^{-1})\right], \end{gather}$$

(3.21b)$$\begin{gather}V_p^{nonlinear}\sim Re_p(\beta^2+\beta\gamma+\gamma^2), \end{gather}$$

for the contributions from terms linear and nonlinear in $\boldsymbol {u}_s$. Using the definitions of $\beta$ and $\gamma$ given below (2.3), and reverting to the dimensional form, imply a leading-order lift velocity of the form $V_{max}^2 L^3/(\nu H^2)[\ln (H/L)+O(1)]$. Here, the contribution to the $O(1)$ term within brackets arises from the $\beta \gamma \lambda ^{-1}$ term in (3.21a) pertaining to the outer region, and from the $\beta ^2$ terms in both (3.21a) and (3.21b) pertaining, respectively, to the outer and inner regions; the logarithmically larger lift contribution arises from the $\beta ^2\ln \lambda ^{-1}$ term in (3.21a) that pertains to the matching region. However, the matching-region $\beta ^2\ln \lambda ^{-1}$ contribution and the inner-region $\beta ^2$ contribution must both vanish by symmetry because they cannot involve boundaries, and therefore, relate to the time-averaged lift on a neutrally buoyant spheroid in an unbounded simple shear flow. Since such a spheroid cannot exhibit a net lateral drift, only the outer-region $\beta ^2$ contribution survives. This and the outer-region $\beta \gamma$ contribution both yield a dimensional lift velocity of $V_{max}^2L^3/(\nu H^2)$; in dimensionless terms, the Jeffery-averaged volume integral in (3.17) is $O(1)$, with the scaled lift velocity being $O(Re_p)$. More detailed arguments along these lines given in Anand & Subramanian (2023) show that the next order contribution to the volume integral for a neutrally buoyant sphere arises from the inner region, involves an additional factor of $\lambda$, but leads to a qualitative alteration of the lift velocity profiles for large $Re_c$. Owing to the acceleration terms in (3.16) only being $O(\lambda ^2)$, the aforementioned inner-region contribution will also be relevant to a neutrally buoyant spheroid.

Owing to the dominant contribution arising from length scales of $O(H)$ implied by the above arguments, the spheroids in the actual and test problems can be replaced by the corresponding point singularities. Thus, $\boldsymbol {u}_s$ and $\boldsymbol {u}^{t2}$ may be approximated as being induced by a time-dependent stresslet ($\boldsymbol {u}_{str}$) and a Stokeslet due to a point force directed along the negative gradient direction ($\boldsymbol {u}_{St}$), respectively. The Jeffery-averaged lift velocity given by (3.17), at leading order in $\lambda$, may therefore be written as

(3.22)\begin{equation} \langle V_p\rangle =-Re_p\int_{V^F+V^P}\left\langle\boldsymbol{u}_{St}\boldsymbol{\cdot} \left(\dfrac{\partial\boldsymbol{u}_{str}}{\partial t}+ \boldsymbol{u}_{str} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^\infty+\boldsymbol{u}^\infty\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_{str} \right)\right\rangle {\rm d}\boldsymbol{r}, \end{equation}

where, on account of the subdominant nature of the length scales of $O(L)$, the domain of integration has now been extended to include the particle volume $V^P$. Thus, (3.22) is the leading-order point-particle approximation for the lift velocity for $Re_c\ll 1$.

In (3.22), $\boldsymbol {u}_{St}$ is time independent since the unit force in the test problem points in a fixed direction, despite the changing orientation of the test spheroid. Furthermore, $\langle \partial \boldsymbol {u}_{str}/\partial t\rangle =0$ for rotation along Jeffery orbits. Also noting that the ambient flow is steady in the chosen non-rotating reference frame, (3.22) reduces to

(3.23)\begin{equation} \langle V_p\rangle =-Re_p\int_{V^F+V^P}\boldsymbol{u}_{St} \boldsymbol{\cdot}(\langle\boldsymbol{u}_{str}\rangle\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}^\infty+ \boldsymbol{u}^\infty\boldsymbol{\cdot}\boldsymbol{\nabla}\langle \boldsymbol{u}_{str}\rangle)\, {\rm d}\boldsymbol{r}. \end{equation}

Using $r_2 \sim O(\lambda ^{-1})$ in (3.15) on account of the outer-region dominance, the Faxen's correction to spheroid translation turns out to be $O(\lambda ^2)$ smaller than the terms linear and quadratic in $r_2$, and one may therefore approximate the ambient flow in (3.23) as $\boldsymbol {u}^\infty \approx (\beta r_2 +\gamma r_2^2)\boldsymbol {1}_1$.

With the spheroid volume neglected, the volume integral in (3.23) is most easily evaluated by Fourier transforming the flow ($r_1$) and vorticity ($r_3$) coordinates, the partial Fourier transform being defined as

(3.24)\begin{equation} \hat{f}(k_1,r_2,k_3) =\int_{-\infty}^\infty \int_{-\infty}^\infty {\rm d}r_1 \,{\rm d}r_3\exp({\iota(k_1 r_1+k_3 r_3)})f(r_1,r_2,r_3). \end{equation}

Applying the convolution theorem along these coordinates (Arfken, Weber & Harris 2011), and substituting the above approximate form of $\boldsymbol {u}^\infty$ in (3.23), one obtains

(3.25)\begin{align} \langle V_p\rangle &=-\frac{Re_p}{4{\rm \pi}^2}\int_{-s\lambda^{-1}}^{(1-s)\lambda^{-1}} \,{\rm d}r_2\int {\rm d}\boldsymbol{k}_\perp\,\hat{\boldsymbol{u}}_{St}(-\boldsymbol{k}_\perp,r_2;y_2) \boldsymbol{\cdot} \left[\langle\hat{\boldsymbol{u}}_{str}\rangle(\boldsymbol{k}_\perp,r_2;y_2)\boldsymbol{\cdot} \boldsymbol{1}_2(\beta+2\gamma r_2)\boldsymbol{1}_1\right.\nonumber\\ &\quad \left.-\iota k_1(\beta r_2 +\gamma r_2^2) \langle\hat{\boldsymbol{u}}_{str}\rangle(\boldsymbol{k}_\perp,r_2;y_2)\right], \end{align}

in terms of the partial Fourier transforms of the Stokeslet and stresslet velocity fields; here, $\boldsymbol {k}_\perp \equiv (k_1,k_3)$ and $y_2 =s\lambda ^{-1}$ is the transverse distance of the stresslet from the lower wall. The expression for $\hat {\boldsymbol {u}}_{St}$ is derived in Appendix A, while that for $\langle \hat {\boldsymbol {u}}_{str}\rangle$ is developed in the next subsection.

3.3. Solution for $(\langle \boldsymbol{u}_{str}\rangle, \langle \,p_{str}\rangle )$

The Jeffery-averaged disturbance fields ($\langle \boldsymbol {u}_{str}\rangle, \langle \,p_{str}\rangle$), that appear in the point-particle approximation of the inertial lift velocity given by (3.23), satisfy

(3.26a)$$\begin{gather} \nabla^2 \langle\boldsymbol{u}_{str}\rangle-\boldsymbol{\nabla}\langle \,p_{str}\rangle= \beta\langle\boldsymbol{S}(\boldsymbol{p})\rangle\boldsymbol{\cdot} \boldsymbol{\nabla}\delta(\boldsymbol{r}), \end{gather}$$

(3.26b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\boldsymbol{u}_{str}\rangle=0, \end{gather}$$

with the boundary conditions

(3.27a)$$\begin{gather} \langle\boldsymbol{u}_{str}\rangle=0 \quad\text{at } r_2=-s\lambda^{-1}, (1-s)\lambda^{-1}, \end{gather}$$

(3.27b)$$\begin{gather}\langle\boldsymbol{u}_{str}\rangle\rightarrow 0\quad \text{for } r_1,r_3\rightarrow\infty (r_2\ \text{fixed}), \end{gather}$$

where the stresslet singularity on the right-hand side of (3.26a) approximates the torque-free neutrally buoyant spheroid, in an unbounded simple shear flow, on scales much larger than $O(L)$. The stresslet coefficient $\boldsymbol {S}$ is a function of the spheroid orientation $\boldsymbol {p}$, being given by (Marath & Subramanian 2017)

(3.28)\begin{align} \boldsymbol{S}(\boldsymbol{p})&= A_1 \frac{3}{2}(\boldsymbol{E}:\boldsymbol{pp})\left(\boldsymbol{pp}-\frac{\boldsymbol{I}}{3}\right) + A_2 ((\boldsymbol{I}-\boldsymbol{pp})\boldsymbol{\cdot}\boldsymbol{E}\boldsymbol{\cdot} \boldsymbol{pp}+\boldsymbol{pp}\boldsymbol{\cdot}\boldsymbol{E}\boldsymbol{\cdot}(\boldsymbol{I}-\boldsymbol{pp}))\nonumber\\ &\quad +A_3\left((\boldsymbol{I}-\boldsymbol{pp})\boldsymbol{\cdot}\boldsymbol{E} \boldsymbol{\cdot}(\boldsymbol{I}-\boldsymbol{pp})+(\boldsymbol{I}-\boldsymbol{pp}) \frac{\boldsymbol{E}:\boldsymbol{pp}}{2}\right), \end{align}

where $E_{ij}={\beta }/{2}(\delta _{i1}\delta _{j2}+\delta _{i2}\delta _{j1})$ is the rate of strain tensor associated with the local simple shear. The simple shear flow may be resolved into an axisymmetric extension aligned with $\boldsymbol {p}$, longitudinal planar extensions in a pair of orthogonal planes containing $\boldsymbol {p}$, and a pair of transverse planar extensions in the plane perpendicular to $\boldsymbol {p}$ (Subramanian & Koch 2006; Dabade et al. 2016). The $A_i$'s in (3.28) are the $\kappa$-dependent stresslet amplitudes corresponding to the aforesaid component flows, there being only three of these (in contrast to the five component flows) owing to the axisymmetry of the spheroidal geometry (Kim & Karrila 1991; Marath & Subramanian 2017). For $\kappa >1$, these amplitudes are given by

(3.29)$$\begin{gather} A_1=-\frac{16{\rm \pi}(\kappa^2-1)^{5/2}}{9\kappa^3[-3(\kappa^2-1)^{1/2}\kappa+2\kappa^2 \cosh^{-1}(\kappa)+\cosh ^{-1}(\kappa)]}, \end{gather}$$

(3.30)$$\begin{gather}A_2=-\frac{16{\rm \pi}(\kappa^2-1)^3}{3\kappa^2(\kappa^2+1)(\kappa^4+\kappa^2-3(\kappa^2-1)^{1/2}\kappa\cosh^{-1}(\kappa)-2)}, \end{gather}$$

(3.31)$$\begin{gather}A_3=-\frac{32{\rm \pi}(\kappa^2-1)^3}{3\kappa^3(2\kappa^5-7\kappa^3+3(\kappa^2-1)^{1/2}\cosh^{-1}(\kappa)+5\kappa)}. \end{gather}$$

The corresponding expressions for an oblate spheroid ($\kappa <1$) may be obtained by first substituting $\kappa =\xi _0/(\xi _0^2-1)^{1/2}$ in terms of the coordinate $\xi _0\ ({>}1)$ labelling the surface of the spheroid, and then using the transformation $d\to -\iota d, \xi _0\to \iota (\xi _0^2-1)^{1/2}$ in the dimensional form of the stresslet obtained from multiplying $\boldsymbol {S}$ above by $\mu L^3 V_{max}/H$ with $L= d\xi _0$ (Marath & Subramanian 2017). Note that $\boldsymbol {p}$, and thence $\boldsymbol {S}$, is a function of time owing to rotation along Jeffery orbits. In the limit of a sphere, $A_1 = A_2 = A_3 = -{20{\rm \pi} }/{3}$, and $\boldsymbol {S}$ reduces to $-20{\rm \pi} \boldsymbol {E}/3$, independent of time, corresponding to the stresslet induced by a freely rotating sphere that leads to the well-known Einstein coefficient in the suspension viscosity.

While one may, in principle, solve for ($\langle \boldsymbol {u}_{str}\rangle, \langle \,p_{str}\rangle$) in a manner similar to that for the bounded domain Stokeslet fields (see Appendix A) by defining the disturbance fields as the sum of an unbounded domain contribution, and a contribution that accounts for the confinement induced by channel walls, one may also obtain them directly as a gradient of the bounded domain Stokeslet fields, derived in Appendix A; the gradient is with respect to the location $\boldsymbol {y}$ of the Stokeslet (Swan & Brady 2010). Therefore,

(3.32)\begin{equation} \langle\boldsymbol{u}_{str}\rangle=\beta\langle\boldsymbol{S}\rangle:\frac{\partial\boldsymbol{J}}{ \partial\boldsymbol{y}}, \end{equation}

where the bounded domain Stokeslet field corresponding to a point force $\boldsymbol {F}$ is given by $\boldsymbol {J}\boldsymbol {\cdot } \boldsymbol {F}$; the Fourier transform ($\hat {\boldsymbol {J}}$) of the tensor $\boldsymbol {J}$ is defined in Appendix A (see text below (A25)). Since it is $\langle \hat {\boldsymbol {u}}_{str} \rangle$ that appears in the final expression for the point-particle approximation of the lift velocity, (3.25), we Fourier transform (3.32) to obtain

(3.33)\begin{equation} \langle\hat{u}_{{str},i}\rangle=\beta \langle S_{jm}\rangle\hat{N}_{ijm}, \end{equation}

where the third-order tensor $\hat {N}_{ijm}$ is defined as

(3.34)\begin{equation} \hat{N}_{ijm}=\iota\hat{J}_{im}(k_1 \delta_{j1}+k_3 \delta_{j3})+\dfrac{\partial\hat{J}_{im}}{\partial y_2}\delta_{j2}. \end{equation}

In (3.33), $\langle \boldsymbol {S}\rangle ={(1}/{2{\rm \pi} })\int _0^{2{\rm \pi} } \boldsymbol {S}\,{\rm d}\tau$ denotes the time-averaged stresslet corresponding to a neutrally buoyant spheroid rotating along a given Jeffery orbit in an unbounded simple shear flow (with unit velocity gradient). The expression for $\langle \boldsymbol {S}\rangle$ is given below in the next subsection.

3.4. The time-averaged lift velocity $\langle V_p \rangle$

Using (3.33) for the stresslet velocity field in (3.25), the Jeffery-averaged lift velocity takes the form

(3.35)\begin{align} \langle V_p\rangle &=-\frac{Re_p}{4{\rm \pi}^2}\int_{-s\lambda^{-1}}^{(1-s)\lambda^{-1}} \,{\rm d}r_2\int {\rm d}\boldsymbol{k}_\perp\,\hat{\boldsymbol{u}}_{St}(-\boldsymbol{k}_\perp,r_2;y_2)\nonumber\\ &\quad \boldsymbol{\cdot} \left[\beta\boldsymbol{1}_2\boldsymbol{\cdot}\hat{\boldsymbol{N}}: \langle\boldsymbol{S}\rangle(\boldsymbol{k}_\perp,r_2;y_2)(\beta+2\gamma r_2)\boldsymbol{1}_1\right.\nonumber\\ &\quad \left.-\iota k_1(\beta r_2 +\gamma r_2^2) \beta\hat{\boldsymbol{N}}:\langle\boldsymbol{S}\rangle(\boldsymbol{k}_\perp,r_2;y_2)\right]. \end{align}

Note that the averaging operation ($\langle.\rangle$) in (3.35) is independent of the shear rate since it involves the scaled Jeffery angular velocity that is only a function of $\kappa$ and $\phi _j$ for the inertially stabilized orbits; the shear-rate dependence only comes in at $O(Re_p^{{3}/{2}}$) (Marath & Subramanian 2017). Thus, $\langle \boldsymbol {S}\rangle$ must be linear in $\boldsymbol {E}$, with a prefactor that is a function of $C$ and $\kappa$. For simple shear flow, this implies that $\langle S_{12}\rangle =\langle S_{21}\rangle$ are the only non-zero components of the Jeffery-averaged stresslet. The following expressions for these components may be obtained by starting from (3.28), with $\boldsymbol {p}$ expressed in terms of $C$ and $\tau$ using (3.19a,b), followed by an integration over $\tau$:

(3.36)\begin{align} \langle S_{12}\rangle&=\frac{1}{4(\kappa^2-1)^2[(C^2+1)(C^2\kappa^2+1)]^{1/2}}\left\{3 A_1 \kappa^2\left(2+C^2(\kappa^2+1)\right.\right.\nonumber\\ &\quad \left.-2 [(C^2+1) (C^2 \kappa^2+1)]^{1/2}\right)\nonumber\\ &\quad +2 A_2(\kappa^2+1)\left(\kappa^2\left([(C^2+1) (C^2 \kappa^2+1)]^{1/2}-2 C^2-1\right)\right.\nonumber\\ &\quad\left.+[(C^2+1) (C^2 \kappa^2+1)]^{1/2}-1\right)\nonumber\\ &\quad +A_3 \left(\kappa^4(C^2+2)+\kappa^2(-2 [(C^2+1) (C^2 \kappa^2+1)]^{1/2}\right.\nonumber\\ &\quad \left.\left.+C^2-2)+2\right)\right\} \end{align}

for prolate spheroids and

(3.37)\begin{align} \langle S_{12}\rangle &= \frac{1}{4(\kappa^2-1)^2(C^2+1)(C^2\kappa^2+1)}\left\{3 A_1 \kappa^2\left(-2C^4\kappa^2 +C^2(\kappa^2+1)\right.\right.(-2\nonumber\\ &\quad\left.+[(C^2+1)(C^2\kappa^2+1)]^{1/2})+2 (-1+[(C^2+1)(C^2\kappa^2+1)]^{1/2})\right)\nonumber\\ &\quad +2 A_2(\kappa^2+1)\left(C^4 (\kappa^2+\kappa^4)-(1+\kappa^2)(-1+[(C^2+1)(C^2\kappa^2+1)]^{1/2})\right.\nonumber\\ &\quad\left.+C^2(1+\kappa^2(2+\kappa^2-2[(C^2+1)(C^2\kappa^2+1)]^{1/2}))\right)\nonumber\\ &\quad +A_3 \left(2[(C^2+1)(C^2\kappa^2+1)]^{1/2}+\kappa^2(-2-2C^4\kappa^2\right.\nonumber\\ &\quad+2(\kappa^2-1)[(C^2+1)(C^2\kappa^2+1)]^{1/2})\nonumber\\ &\quad\left.\left.+C^2(1+\kappa^2)(-2+[(C^2+1)(C^2\kappa^2+1)]^{1/2})\right) \right\} \end{align}

for oblate spheroids. For both (3.36) and (3.37), $\lim _{\kappa \rightarrow 1} \langle S_{12}\rangle =-10{\rm \pi} /3$ independent of $C$, which yields the Einstein coefficient ($5/2$) in the $O(\phi )$ contribution to the suspension viscosity, $\phi$ being the sphere volume fraction; the $C$ independence arises owing to the sphere orientation being a degenerate degree of freedom. At the other extreme, one has $\lim _{\kappa \to \infty } \langle S_{12}\rangle |_{C=\infty }=-2{\rm \pi} /(3\kappa \ln \kappa )$ in (3.36). Here, the factor $1/\ln \kappa$ arises from viscous slender body theory (Subramanian & Koch 2005), while the additional $\kappa ^{-1}$ factor reflects the probability of the occurrence of non-aligned orientations that contribute dominantly to the viscosity of a dilute non-interacting suspension of slender fibres (Leal & Hinch 1971; Dabade et al. 2016).

For times much longer than $O(Re_p^{-1}H/V_{max})$, only the inertially stabilized orbits are relevant to cross-stream migration. As mentioned earlier, for prolate spheroids, this is the tumbling orbit ($C=\infty$), in which case, $\langle S_{12}\rangle$ in (3.36) reduces to

(3.38)\begin{equation} \langle S_{12}\rangle|_{C=\infty} = \frac{ (3 A_1+A_3)\kappa+2 A_2 (\kappa ^2+1)}{4 (\kappa +1)^2}. \end{equation}

For oblate spheroids, one can have either the spinning ($C=0$) or tumbling ($C=\infty$) orbits depending on $\kappa$, and accordingly, $\langle S_{12}\rangle$ in (3.37) reduces to

(3.39)$$\begin{gather} \langle S_{12}\rangle|_{C=\infty} = \frac{ (3 A_1+A_3)\kappa+2 A_2 (\kappa ^2+1)}{4 (\kappa +1)^2}, \end{gather}$$

(3.40)$$\begin{gather}\langle S_{12}\rangle|_{C=0} = \frac{A_3}{2}. \end{gather}$$

Equations (3.38)–(3.40) will be used to determine $\langle V_p \rangle$ below. In light of the above, one may write (3.33) in the form

(3.41)\begin{equation} \langle\hat{u}_{{str},i}\rangle=\beta \langle S_{12}\rangle\,(\hat{N}_{i12}+\hat{N}_{i21}), \end{equation}

with the Jeffery-averaged lift velocity given by

(3.42)\begin{align} \langle V_p\rangle &=-\frac{Re_p\langle S_{12}\rangle} {4{\rm \pi}^2} \int_{-s\lambda^{-1}}^{(1-s)\lambda^{-1}} {\rm d}r_2\int {\rm d}\boldsymbol{k}_\perp\,\hat{u}_{{St},i} (-\boldsymbol{k}_\perp,r_2;y_2) \nonumber\\ &\quad \times\big[\beta(\beta+2\gamma r_2)(\hat{N}_{212}+\hat{N}_{221})(\boldsymbol{k}_\perp,r_2;y_2)\delta_{i1}\nonumber\\ &\quad -\iota \beta k_1(\beta r_2 +\gamma r_2^2) (\hat{N}_{i12}+\hat{N}_{i21}) (\boldsymbol{k}_\perp,r_2;y_2)\big]. \end{align}

Since the dominant contributions to the integral in (3.42) come from scales of $O(H)$ (or, equivalently, $k_\perp \sim O(H^{-1})$), it is natural to transform to rescaled variables $\boldsymbol {k}_\perp =\boldsymbol {k}_\perp ''\lambda, r_2=r_2''/\lambda$. The original disturbance fields may be written in terms of the rescaled ones as $\hat {\boldsymbol {u}}_{St} =\hat {\boldsymbol {u}}_{St}^{\prime \prime } /\lambda$, $\hat {\boldsymbol {N}}=\hat {\boldsymbol {N}}''$, and furthermore, noting that $y_2''= s$ and ${\rm d}\boldsymbol {k}_\perp \,{\rm d}r_2\equiv \lambda \, {\rm d}\boldsymbol {k}_\perp ''\,{\rm d}r_2''$, one obtains

(3.43)\begin{align} \langle V_p\rangle &=-\frac{Re_p\langle S_{12}\rangle} {4{\rm \pi}^2} \int_{-s}^{1-s} {\rm d}r_2''\int {\rm d}\boldsymbol{k}''_\perp\,\hat{u}_{{St},i}^{\prime\prime} (-\boldsymbol{k}''_\perp,r_2'';s) \nonumber\\ &\quad\times \big[\beta(\beta+2\gamma'' r_2'')(\hat{N}_{212}''+\hat{N}_{221}'') (\boldsymbol{k}''_\perp,r_2'';s)\delta_{i1}\nonumber\\ &\quad -\iota \beta k_1''(\beta r_2'' +\gamma'' r_2''^2) (\hat{N}_{i12}''+\hat{N}_{i21}'') (\boldsymbol{k}''_\perp,r_2'';s)\big], \end{align}

where $\gamma ''=\gamma \lambda ^{-1}=-4$, and $\langle S_{12} \rangle$ is given by (3.38), (3.39) or (3.40) depending on $\kappa$. The essential consequence of the rescaling above is to show that the volume integral in (3.43) is independent of $\lambda$. In fact, the integral is only a function of $s$ (the particle location in the channel; see figure 1), with the additional dependence on $\kappa$ entirely contained in $\langle S_{12} \rangle$. Accounting for the scaling $V_{max}\lambda$ used for $V_p$, (3.43) shows that the point-particle framework leads to a lift velocity of $O(V_{max}\lambda Re_p)$, which was used earlier in § 3.1 to estimate the inertial migration time scale ($t_{lift}$).

A consequence of the integral in (3.43) being independent of $\kappa$ is that, within the Jeffery-averaged framework used, the inertial lift velocity of a sphere, and the time-averaged lift velocity of a spheroid, only differ by a multiplicative function of $C$ and $\kappa$ in general, and only by a function of $\kappa$ if one assumes the spheroid to have settled onto the inertially stabilized Jeffery orbit – this multiplicative function is the ratio of the Jeffery-averaged spheroid stresslet to the sphere stresslet, being equal to $-{3\langle S_{12} \rangle }/{10{\rm \pi} }$. An immediate outcome of this proportionality relation is that the shapes of the lift velocity profiles for the two cases must be identical, with the zero crossings in particular (that correspond to the Segre–Silberberg equilibria for a sphere) being identical. In other words, a change in $\kappa$ only affects the magnitude of the inertial lift, the equilibrium positions being unaffected. In light of the obvious importance of the above conclusion for shape-sorting applications, it is worth documenting the reasons that lead to the simple proportionality relation.

1. (i) The dominant contribution from the linearized inertial terms on scales of $O(H)$. This outer-region dominance led to the point-particle approximation, and thence, to the time independence of the test velocity field in the reciprocal theorem volume integral.

2. (ii) The asymptotic smallness, for small confinement ratios, of the additional contributions arising from the translational and angular accelerations associated with the motion of the neutrally buoyant spheroid in the Stokes limit,

3. (iii) The approximation of the time averaging based on the $Re_p-$independent Jeffery angular velocity, which leads to the time-averaged spheroid stresslet being proportional to $\boldsymbol {E}$ (just as the sphere stresslet).

The Jeffery-averaged lift velocity in (3.43) may be written in the compact form

(3.44)\begin{equation} \langle V_p\rangle=Re_p\langle S_{12}\rangle(\kappa)(\beta^2 F(s)+\beta\gamma''G(s)), \end{equation}

where the functions $F(s)$ and $G(s)$ are given by

(3.45)\begin{align} F(s)&=-\frac{1} {4{\rm \pi}^2} \int_{-s}^{1-s} {\rm d}r_2''\int {\rm d}\boldsymbol{k}''_\perp\, \hat{u}_{{St},i}^{\prime\prime} (-\boldsymbol{k}''_\perp,r_2'';s) \nonumber\\ &\quad \times\big[(\hat{N}_{212}''+\hat{N}_{221}'') (\boldsymbol{k}''_\perp,r_2'';s)\delta_{i1}-\iota k_1'' r_2'' (\hat{N}_{i12}''+\hat{N}_{i21}'') (\boldsymbol{k}''_\perp,r_2'';s)\big], \end{align}

(3.46)\begin{align} G(s)&=-\frac{1} {4{\rm \pi}^2} \int_{-s}^{1-s} {\rm d}r_2''\int {\rm d}\boldsymbol{k}''_\perp\, \hat{u}_{{St},i}^{\prime\prime} (-\boldsymbol{k}''_\perp,r_2'';s) \nonumber\\ &\quad\times\big[2 r_2''(\hat{N}_{212}''+\hat{N}_{221}'') (\boldsymbol{k}''_\perp,r_2'';s)\delta_{i1}-\iota k_1''r_2''^2 (\hat{N}_{i12}''+\hat{N}_{i21}'') (\boldsymbol{k}''_\perp,r_2'';s)\big], \end{align}

and satisfy $F(s)=-F(1-s)$ and $G(s)=G(1-s)$, consistent with the antisymmetry of the lift velocity profile about the channel centreline.

In (3.44) the $\beta ^2$ contribution characterizes the effect of the asymmetrically located walls on the disturbance stresslet, and leads to migration away from the walls. In contrast, the $\beta \gamma ''$ contribution, which characterizes the interaction of the stresslet with the ambient profile curvature, causes migration away from the channel centreline; note that this contribution is still influenced by the walls, given that the dominant scales are $O(H)$. The Segre–Silberberg equilibria emerge from a balance between these two opposing effects. Substituting $\beta =4(1-2s)$ and $\gamma ''=-4$ in (3.44), the final expression for the time-averaged lift velocity of a neutrally buoyant spheroid of an arbitrary aspect ratio $\kappa$, for $Re_c\ll 1$, is given by

(3.47)\begin{equation} \langle V_p\rangle=Re_p\langle S_{12}\rangle(\kappa)\left[16(1-2s)^2 F(s)-16(1-2s)G(s)\right]. \end{equation}

The Fourier integrals in (3.45) and (3.46) can be calculated using plane polar coordinates, $k_1''=k_\perp ''\cos \phi$, $k_3''=k_\perp ''\sin \phi$, with ${\rm d}\boldsymbol {k}_\perp ''\equiv k_\perp ''\,{\rm d}k_\perp ''\,{\rm d}\phi$; here $k_\perp ''\in [0,\infty )$ and $\phi \in [0,2{\rm \pi} ]$. The integrals over the transverse coordinate $r_2''$ and the azimuthal angle $\phi$ may be done analytically, reducing (3.45) and (3.46) to the one-dimensional integrals

(3.48)$$\begin{gather} F(s)=\int_0^\infty {\rm d}k_\perp'' \dfrac{k_\perp''{\rm e}^{-k_\perp'' (27 s+16)} I(k_\perp'',s)}{48 {\rm \pi}\big({\rm e}^{2 k_\perp''}-1\big) \left[-2 {\rm e}^{2 k_\perp''} \big(2 k_\perp''^2+1\big)+ {\rm e}^{4 k_\perp''}+1\right]^2}, \end{gather}$$

(3.49)$$\begin{gather}G(s)=\int_0^\infty {\rm d}k_\perp'' \dfrac{{\rm e}^{-k_\perp'' (27 s+16)} J(k_\perp'',s)}{192 {\rm \pi}k_\perp''^2 \big({\rm e}^{2 k_\perp''}-1\big) \left[-2 {\rm e}^{2 k_\perp''} \big(2 k_\perp''^2+1\big)+{\rm e}^{4 k_\perp''}+1\right]^2}, \end{gather}$$

with $I(k_\perp '',s)$ and $J(k_\perp '',s)$ defined in Appendix B. The $k_\perp ''$ integrals above are evaluated numerically for various $s$, using Gauss–Legendre quadrature, after replacing the infinite interval with a finite one – $(0,K_{max})$. Numerical convergence, especially near the channel walls, depends sensitively on $K_{max}$. As the spheroid approaches either wall, the contribution to the $k_\perp ''$ integrals comes from progressively smaller scales (compared with $H$) or, equivalently, from larger and larger $k_\perp ''$. In fact, the relevant scale changes from $H$ ($k_\perp ''\sim O(1)$) to the distance from the wall – either $s$ (lower wall) or $1-s$ (upper wall). As a result, the choice of any finite $K_{max}$, however large, will only lead to converged lift velocities down to a certain non-zero $s$.

To analyse the lift close to the lower wall, for example, one defines a rescaled wavenumber $k_\perp ''=k_w/s$. Next, expanding the integrands in (3.48) and (3.49) for $s\to 0$ with $k_w$ fixed, the near-wall lift velocity takes the form

(3.50)\begin{equation} \lim_{s\to 0}\,\langle V_p\rangle=\langle V_p\rangle^{wall}=-\frac{Re_p\langle S_{12}\rangle(\kappa)}{3{\rm \pi}} \int_0^\infty {\rm d}k_w\,{\rm e}^{-2 k_w} k_w (3 k_w^2-2 k_w+3). \end{equation}

The limiting value near the upper wall is equal in magnitude, but of an opposite sign, as must be the case by symmetry. The integral in (3.50) can be evaluated analytically, and yields

(3.51)\begin{equation} \langle V_p\rangle^{wall}={\pm} \dfrac{11 Re_p\langle S_{12}\rangle(\kappa)}{24{\rm \pi}}, \end{equation}

with the upper and lower signs pertaining to the corresponding channel wall. For a sphere, (3.51) reduces to $55 Re_p/36$, a value originally given by Vasseur & Cox (1976), albeit without an accompanying explanation. In the results shown in the next subsection, we have chosen $K_{max}=10^4$ for the numerical evaluation of (3.48) and (3.49), to ensure accuracy down to $s=0.001$ without the aid of a large-$k_\perp ''$ asymptote. This choice also enables a close approach to the near-wall limiting value above, and the sufficiency of $200$ quadrature points in obtaining converged results for arbitrary $\kappa$.

The sign of the limiting value in (3.51) points to a wall-induced repulsion. Furthermore, (3.51) arises entirely from the $\beta ^2$ contribution, with the $\beta \gamma ''$ contribution being asymptotically small, implying that the leading-order repulsion may be determined independently by considering a sphere or spheroid moving parallel to a single plane wall subject to a linear shearing flow. The near-wall estimates of the said contributions may be obtained by considering a spheroid, regarded as the relevant point singularity, close to either wall. For the lower wall, this would mean $s\ll 1$, with an inner region characterized by distances of the order of the spheroid size ($\lambda < r''\ll s$), and an outer region characterized by distances of the order of the spheroid-wall separation ($r'' \gtrsim s$); recall that $r''$ was measured in units of $H$. In accordance with (3.20), the lift velocity integral grows on scales pertaining to the matching interval, $\lambda \ll r'' \ll s$, with the dominant contributions being of the form $O(\beta ^2)[\ln r'' + O(1)]$ and $O(\beta \gamma '' r'')$. The logarithmic term is again zero due to symmetry, and the linear growth is cutoff for distances larger than $O(s)$ owing to the more rapid decay arising from wall-induced cancellation of the leading-order singularity. For $s \ll r'' \ll 1$, the actual velocity field exhibits an $O(1/r''^3)$ quadrupolar decay owing to wall-induced stresslet cancellation; the test velocity field is also $O(1/r''^3)$ over this range of distances owing to the wall-induced cancellation of both the point force and force dipole for a Stokeslet of a perpendicular orientation. This more rapid decay ensures convergence of the volume integral, and using $r'' \sim s$ in the growing terms above yields the limiting forms of $O(\beta ^2)$ and $O(\beta \gamma '' s)$ close to the lower wall.

A final important point concerns the approach of the near-wall lift velocity to a finite value, rather than zero, as one approaches the channel walls. The latter must be the case on account of the eventual diverging lubrication resistance associated with the thin intervening fluid layer between the particle and the wall, and the discrepancy is on account of the regime of validity of the present calculation. As stated in the problem definition in § 2, the near-wall limit above pertains to spheroid-wall separations that, although much smaller than $O(H)$, are nevertheless much larger than $L$. We comment further on the nature of the lift profile for small separations in § 4.

3.5. Results and discussion

In figure 2 we compare the lift velocity profile for a sphere, obtained by taking $\langle S_{12} \rangle = -{10{\rm \pi} }/{3}$ in (3.47), with profiles plotted using the data from table 4 of Ho & Leal (1974) and digitizing figure 8 of Vasseur & Cox (1976). There is a clear mismatch between our profile and the Ho–Leal one, an aspect that we comment on in the conclusion. The result of Vasseur & Cox (1976) shows good agreement with the present calculation throughout the channel. The Segre–Silberberg equilibria, corresponding to the zeros of $\langle V_p\rangle$, are located at $s\approx 0.182$ and $0.818$; these are at a (dimensional) distance of $0.636\times H/2$ from the centreline, which agrees well with the intermediate annulus at $\sim 0.6\times$pipe radius in the original experiments (Segre & Silberberg 1962a). Note that the lift velocity profiles in both our and the Vasseur & Cox (1976) analysis asymptote to the wall values $\pm 55/36$ (the horizontal dashed magenta lines in figure 2) given in the previous subsection. As already argued therein, $H$ is no longer the relevant scale (at leading order) close to the wall. This implies that $Re_c$ should also not be relevant, and the near-wall lift values should therefore be independent of $Re_c$. This is validated in the next section where the lift velocity profiles are seen to approach (3.51) for $s \to 0,1$ even for $Re_c\gtrsim O(1)$, although this approach occurs in a shrinking neighbourhood of the wall with increasing $Re_c$. The inset in figure 2 plots the $\beta ^2$ and $\beta \gamma ''$ contributions to the inertial lift, confirming that the latter curvature-induced contribution becomes vanishingly small at the walls, consistent with the scaling arguments given at the end of § 3.4.

Figure 2. Comparison of the small-$Re_c$ lift profile for a sphere obtained from (3.46) with profiles extracted from Ho & Leal (1974) and Vasseur & Cox (1976). The vertical lines mark the Segre–Silberberg equilibria ($s\approx 0.182$ and $0.818$) on either side of the centreline; the horizontal dotted lines denote the near-wall limiting values given by (3.51). The inset shows the profiles of the component $\beta ^2$ and $\beta \gamma ''$ contributions.

As mentioned after (3.43), while the magnitude of the Jeffery-averaged lift velocity is sensitive to $\kappa$, the equilibrium locations remain the same as those for a sphere. This aspect ratio insensitivity of the equilibria, for spinning oblate spheroids, has been observed in Nizkaya et al. (2020). The lack of a $\kappa$ dependence is seen from figures 3 and 4, which show that the inertial lift profiles for both tumbling prolate spheroids and spinning/tumbling oblate spheroids of different $\kappa$ have the same zero crossings. In figure 3(a) the inertial lift for a tumbling prolate spheroid at a fixed $s$ is seen to decrease with increasing $\kappa$, consistent with the decreasing magnitude of the disturbance velocity field. To better examine the large-$\kappa$ limit, in figure 3(b) we plot the lift profiles normalized by the large-$\kappa$ scaling of $\langle S_{12} \rangle$ found earlier. The scaled profiles approach a $\kappa -$independent limiting form for $\kappa \to \infty$, although the approach is non-monotonic – the scaled profile with the maximum amplitude corresponds to $\kappa = 5$ – reflecting the non-monotonic approach of the scaled stresslet, $\kappa \ln \kappa \langle S_{12} \rangle$, to its infinite-$\kappa$ limit of $-{2{\rm \pi} }/{3}$.

Figure 3. (a) The small-$Re_c$ lift profiles for tumbling prolate spheroids of various aspect ratios. (b) Lift profiles in (a) rescaled using the infinite $\kappa$ $\langle S_{12}\rangle$ scaling. The thick black curve denotes the limiting $\kappa$-independent profile for $\kappa \gtrsim 1000$.

Figure 4. The small-$Re_c$ lift profiles for (a) spinning oblate spheroids of $0.14\leq \kappa <1$, and (b) spinning/ tumbling oblate spheroids of $\kappa <0.14$, for $Re_c\ll 1$. The inset in (b) shows the collapsed profiles on rescaling with the $\kappa \to 0$ limit of $\langle S_{12}\rangle$.

In figure 4(a) the inertial lift velocity of a spinning oblate spheroid is seen to decrease in magnitude as $\kappa$ decreases from $1$ to $0.14$, again consistent with the expected decrease in the magnitude of the disturbance velocity field. At the latter $\kappa$, a reversal of the inertia-induced orbital drift leads to both tumbling and spinning modes being stabilized, with the respective basins of attraction being demarcated by a pair of unstable limit cycles on the unit sphere (Einarsson et al. 2015; Dabade et al. 2016). In figure 4(b) we therefore plot the inertial lift profiles of both spinning and tumbling oblate spheroids for $0 < \kappa < 0.14$. The lift profiles for spinning spheroids approach a finite limiting form for $\kappa \rightarrow 0$, while those for tumbling spheroids are much smaller in magnitude. The inertial lift in the latter case, in fact, goes to zero on account of $\langle S_{12} \rangle$ being $O(\kappa )$ for $\kappa \to 0$; the inset in figure 4(b) shows the tumbling oblate spheroid lift profiles, when scaled by $\kappa$, approaching a finite limiting form in this limit. The $O(\kappa )$ scaling above reflects both the small $O(\kappa )$ probability of occurrence of thin oblate spheroids with symmetry axes not closely aligned with the gradient–vorticity plane, with $O(1)$ stresslet coefficients, and the $O(\kappa )$ stresslet coefficient of such spheroids oriented in the neighbourhood of the gradient–vorticity plane. Unlike the prolate case, both the aligned and non-aligned phases of a tumbling thin oblate spheroid make comparable $O(\kappa )$ contributions to $\langle S_{12} \rangle$.

The invariance of the equilibrium locations associated with the lift profiles in figures 3 and 4, with changing $\kappa$, implies that inertia-induced shape sorting of spheroidal particles is precluded in sufficiently long channels. All spheroids regardless of $\kappa$ will end up migrating to the same pair of Segre–Silberberg locations in such channels. However, one may achieve separation in channels short enough for the residence time of $O(L_c/V_{max})$ to be comparable or smaller than $t_{lift}$; this requires $L_c \lesssim HRe_p^{-1}\lambda ^{-1}$, $L_c$ being the channel length. Shape sorting of spheroids would occur in these channels owing to differential rates of migration: for instance, with $L$ (and, therefore, $Re_p$) remaining the same, prolate spheroids with $\kappa \sim O(1)$ will migrate to the equilibrium locations relatively rapidly, while those with larger $\kappa$ remain close to their initial positions due to weaker lift forces. A similar scenario would prevail for oblate spheroids provided $\kappa \gtrsim 0.14$. For $\kappa \lesssim 0.14$, the disparity in magnitudes of the lift velocities for spinning and tumbling oblate spheroids would lead instead to sorting of spheroids with the same $\kappa$, but in different (tumbling vis-a-vis spinning) orientation modes.

For even shorter channels with $L_c \lesssim O(H Re_p^{-1})$, the residence time is no longer enough for a spheroid to be able to migrate to its stable Jeffery orbit. For residence times of $O(L_c/V_{max}) \lesssim O(H/V_{max}Re_p^{-1}) \gg H/V_{max}$, the cross-stream migration is still determined, at leading order, by a Jeffery-averaged lift velocity, but one pertaining to orbits other than the tumbling or spinning mode. Thus, $\langle V_p \rangle$ is now also a function of $C$, being given by (3.44) with $\langle S_{12} \rangle$ given by its full form involving both $C$ and $\kappa$; see (3.36) and (3.37). The $C$-dependent lift velocity profiles for $\kappa =2$ and $\kappa =0.5$ are shown in figures 5(a) and 5(b), respectively. Since, for a fixed $\kappa$, the amplitude of the disturbance field is the largest for a tumbling prolate spheroid and a spinning oblate spheroid, the largest amplitude profiles in figures 5(a) and 5(b) correspond to $C = \infty$ and $C = 0$, respectively. Determining single-spheroid trajectories now requires solving the following system of coupled ODEs in $s$ and $C$:

(3.52)$$\begin{gather} \frac{{\rm d}s}{{\rm d}t_s} = \lambda\langle S_{12}\rangle(C,\kappa)\left[16(1-2s)^2 F(s)-16(1-2s)G(s)\right], \end{gather}$$

(3.53)$$\begin{gather}\frac{{\rm d}C}{{\rm d}t_s} = \frac{\beta C}{2{\rm \pi}} \left(\sum_{n=1}^{6} I_n(C,\kappa) F_n^f (\kappa)+\sum_{n=1}^{4} J_n(C,\kappa) G_n^f(\kappa)\right). \end{gather}$$

Here $t_s = Re_pt$ is a slow time variable, and the functions $F_n^f$, $G_n^f$, $I_n$ and $J_n$ have been defined in Dabade et al. (2016) and Marath & Subramanian (2018). Note that the $\lambda$ prefactor in (3.52) indicates that spheroids only migrate by a distance of the order of their own size in such short channels.

Figure 5. Lift velocity profiles for prolate and oblate spheroids, of the indicated aspect ratios, rotating in different Jeffery orbits; $C=0$ and $\infty$ correspond to the spinning and tumbling modes. Results are shown for (a) $\kappa =2$ and (b) $\kappa =0.5$.

At higher volume fractions, the trajectory of a given inertial spheroid will acquire a stochastic character on account of occasional pair interactions that may be modelled via a scattering kernel involving pre- and post-interaction orbit constants (Marath et al. 2017). Such interactions are particularly important for oblate spheroids with $\kappa \lesssim 0.14$ where fluctuations in $C$, and thence, changes in $\langle V_p \rangle (C,\kappa )$, occur on time scales that can be much longer than the nominal drift time scale of $O(Re^{-1}H/V_{max})$, and that characterize barrier-hopping events equilibrating the tumbling and spinning orientations. In the presence of such interactions, the transient evolution of an initial distribution of spheroids, on time scales long compared with $T_{jeff}$, is formally described by a probability density $P(s,C,t_s)$ that satisfies a kinetic equation of the form

(3.54)\begin{align} &\frac{\partial P}{\partial t_s} + \frac{\partial}{\partial s}(\dot{s}P) + \frac{\partial}{\partial C}(\dot{C}P) \nonumber\\ &\quad = nL^3\displaystyle\int {\rm d}C' \displaystyle\int {\rm d}\boldsymbol{r}_\perp (s-s') \displaystyle\int {\rm d}\hat{C}\,{\rm d}\hat{C}'[\hat{P}\hat{P}'{\mathcal{K}}(\hat{C},\hat{C}'|C,C';s) - P P'], \end{align}

where $P' \equiv P(s,C')$, $\hat {P} \equiv P(s,\hat {C})$, $\hat {P}' \equiv P(s,\hat {C}')$; $\dot {s}$ and $\dot {C}$ being given by (3.52) and (3.53), respectively. The right-hand side of (3.54) denotes a Boltzmann-type kernel involving the pre-($[\hat {C},\hat {C}']$) and post-($[C,C']$) interaction orbit constant pairs.

At the end of § 3.4, we highlighted the finite value of the inertial lift velocity attained even as the particle approaches either wall. This runs counter to one's expectation of the lift velocity vanishing in this limit due to the diverging resistance associated with the thin lubricating layer of fluid between the particle and the wall. As mentioned therein, the discrepancy arises due to the present analysis only being valid for $s,1-s\gg \lambda$, a restriction that comes from treating the spheroid as a point (stresslet) singularity. Our analysis has to be supplemented by one that accounts for the finite size of the spheroid and is therefore valid for $s,1-s \sim O(\lambda )$, which in turn would connect to the lubrication regime corresponding to $s-\lambda,(1-s)+\lambda \ll \lambda$. Such a connection is possible for the case of a sphere based on results available in the literature. The inertial lift on a sphere, in the presence of a single plane boundary subject to a linear shearing flow, has been evaluated numerically for $s \sim O(\lambda )$ by Cherukat & Mclaughlin (1994), with the aid of an integral expression obtained using bispherical coordinates. The authors found the lift force to always have a repulsive character (that is, to be directed away from the wall). The numerical value asymptoted to the near-wall limit (${\approx }55{\rm \pi} Re_p/6$) of the two-wall (channel) problem mentioned above for $s\gg \lambda$, while asymptoting to a different finite value ($9.22 Re_p$) when the sphere touches the wall ($s=\lambda$). The latter value was shown to agree with the inertial lift force acting on a stationary non-rotating sphere, in contact with a plane boundary, evaluated by Leighton & Acrivos (1985). These authors examined a non-rotating sphere in light of results obtained earlier (Goldman, Cox & Brenner 1967) that showed that the angular velocity of a torque-free sphere must decrease to zero logarithmically in the limit of a vanishing sphere-wall separation on account of the lubrication resistance associated with the relative tangential motion in the narrow gap; the finiteness of the force implies that the dominant contributions to the lift arise from the fluid domain outside the thin gap. When combined with the known $O(s-\lambda )^{-1}$ divergence of the translational resistance for normal approach towards the wall, one concludes that the inertial lift velocity for a sphere must start from $55 Re_p/36$ for $\lambda \ll s \ll 1$, and eventually approach zero linearly in the limit $s-\lambda \ll \lambda$. In the reciprocal theorem formulation used here, the approach to zero would appear via the divergence of the test problem resistance coefficient. To our knowledge, analogous results for an arbitrary aspect ratio spheroid when $s, 1-s \sim O(\lambda )$ are not available, and will involve more effort. On one hand, there is no spheroidal analog of the bispherical coordinate system (Cherukat & Mclaughlin 1994) that would enable the derivation of a formal expression for the inertial lift for arbitrary spheroid-wall separations. On the other hand, the additional orientation degrees of freedom make the spheroid problem substantially more complicated even in the Stokes limit. Apart from an increased period, the nature of the orientation dynamics itself changes from Jeffery rotation to eventually pole vaulting (Mody & King 2005; Stover & Cohen 1990) for a sufficiently close approach towards either wall, with the changes in orientation intimately coupling to corresponding changes in position (Pozrikidis 2005).

4. The inertial lift velocity for $Re_c\gtrsim O(1)$

Herein, we calculate the inertial lift profiles for $Re_c \gtrsim O(1)$ using a numerical shooting method employed originally by Schonberg & Hinch (1989), and later for higher $Re_c$'s by Asmolov (1999), both for the case of a sphere. For $Re_c\gtrsim O(1)$, the inertial screening length ($H Re_c^{-1/2}$) is of the order of the channel width or smaller, and one must solve the governing equations (2.4a,b) with the boundary conditions (2.5a–c) using a matched asymptotic expansions approach. The inner region is characterized by scales of $O(L)$, and the outer region by scales of $O(HRe_c^{-1/2})$, with the channel width $H$ also entering via the wall boundary conditions. Although we use a shooting method below to evaluate the lift velocity, it is worth noting that the reciprocal theorem formulation in § 3.1 remains valid for any $Re_c$, provided one uses the finite-$Re_c$ disturbance velocity field ($\boldsymbol {u}'$) in the volume integral. Thus, the scaling arguments in § 3.2 may be extended to finite $Re_c$, and in doing so, one finds that the dominant contribution to the inertial lift continues to come from scales asymptotically larger than $O(L)$. For $Re_c\ll 1$ in the previous section, this outer-region dominance meant solving the Stokes equations in a domain confined by plane channel walls with the spheroid approximated as a stresslet singularity. In the present section this means solving the linearized Navier–Stokes equations at leading order, driven by the same stresslet singularity (since $Re_p$ is assumed small). Note that $Re_c$, interpreted as the square of the ratio of the two outer-region length scales ($H$ and $HRe_c^{-{1}/{2}}$), appears as a parameter in the governing equations below, and the inertial lift velocity is therefore now a function of $Re_c$.

4.1. The finite-$Re_c$ formulation

The Stokesian disturbance field due to a freely suspended particle in an ambient linear flow only decays algebraically, as $1/r^2$, at large distances, a fact already used in the scaling arguments in § 3.2. The algebraic decay implies that distinct outer-region expansions for both the velocity and pressure fields become necessary on scales of $O(HRe_c^{-1/2})$, with the leading terms satisfying the linearized Navier–Stokes and continuity equations. To write down these equations, one transforms to outer coordinates using $\boldsymbol {r}=Re_p^{-1/2}\boldsymbol {R}$, which corresponds to using $HRe_c^{-{1}/{2}}$ (rather than $L$ as in § 2) as the relevant length scale. The Stokesian rates of decay in the inner region suggest the scalings $\boldsymbol {u}'=Re_p\boldsymbol {U}$ and $p'=Re_p^{3/2} P$ for the leading-order terms in the outer expansions, where $\boldsymbol {U}$ and $P$ satisfy

(4.1a)\begin{gather} \frac{\partial^2 U_i}{\partial R_m^2}-\frac{\partial P}{\partial R_i}-\frac{\partial U_i}{\partial t}-U_2(\beta+2\gamma''R_2 Re_c^{-1/2})\delta_{i1}-(\beta R_2+\gamma''R_2^2 Re_c^{-1/2})\frac{\partial U_i}{\partial R_1}\nonumber\\ =\beta S_{im}\frac{\partial\delta(\boldsymbol{R})}{\partial R_m}, \end{gather}

(4.1b)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}=0. \end{gather}

The Faxen correction contributes at a higher order in $\lambda$ and, therefore, $u^\infty _i\approx Re_p^{-1/2}(\beta R_2+\gamma ''\,Re_c^{-1/2} R_2^2) \delta _{i1}$ has been used in (4.1a). Equations (4.1a,b) must be supplemented by the conditions

(4.2a)$$\begin{gather} U_i \sim \frac{3\beta R_i R_j S_{jm}R_m}{4{\rm \pi} R^5} \quad \text{for } \boldsymbol{R}\rightarrow 0, \end{gather}$$

(4.2b)$$\begin{gather}\boldsymbol{U}=0\quad \text{at } R_2=-s\,Re_c^{1/2}, (1-s)\,Re_c^{1/2}, \end{gather}$$

where (4.2b) denotes the no-slip conditions on the channel walls, while (4.2a) is the requirement of matching to the stresslet velocity field that arises as the far-field form of the inner-region Stokesian field. In (4.1a) and (4.2a), $\boldsymbol {S}$ is the tensorial amplitude defined in (3.28).

As explained in § 3.2, the separation between the migration and drift time scales implies one needs only solve for the Jeffery-averaged lift velocity, and towards this end, we average (4.1a,b) and (4.2a,b) over a single period of the stable Jeffery orbit. One obtains

(4.3a)\begin{gather} \frac{\partial^2\langle U_i\rangle}{\partial R_m^2} -\frac{\partial\langle P\rangle}{\partial R_i}-\langle U_2\rangle(\beta+2\gamma''R_2 Re_c^{-1/2})\delta_{i1}-(\beta R_2+\gamma''R_2^2 Re_c^{-1/2})\frac{\partial\langle U_i\rangle}{\partial R_1}\nonumber\\ =\beta\langle S_{12}\rangle\left[\delta_{i1}\frac{\partial\delta(\boldsymbol{R})}{\partial R_2}+\delta_{i2} \frac{\partial\delta(\boldsymbol{R})}{\partial R_1}\right], \end{gather}

(4.3b)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \langle\boldsymbol{U}\rangle=0, \end{gather}

where $\langle \boldsymbol {U}\rangle$ satisfies

(4.4a)$$\begin{gather} \langle U_i\rangle \sim \frac{3 \beta \langle S_{12}\rangle R_1 R_2 R_i}{4{\rm \pi} R^5}\quad\text{for } \boldsymbol{R}\rightarrow 0, \end{gather}$$

(4.4b)$$\begin{gather}\langle\boldsymbol{U}\rangle=0\quad \text{at } R_2=-s\,Re_c^{1/2}, (1-s)\,Re_c^{1/2}. \end{gather}$$

Here, $\langle S_{12}\rangle (\kappa )$ corresponds to the stresslet averaged over the relevant stable Jeffery orbit, and has been given in (3.38)–(3.40). Due to the linearity of (4.3a,b) and (4.4a,b), and the fact that the Jeffery-averaged spheroid stresslet tensor differs from that for a sphere only by a scalar multiplicative factor, $\langle \boldsymbol {U} \rangle$ and $\langle P \rangle$ differ from their spherical analogs only by $-{3\langle S_{12}\rangle (\kappa )}/{10{\rm \pi} }$, corresponding to the ratio of the aforementioned stresslets. This proportionality relation must hold for any linear functional of the disturbance fields, and in particular, for the lift velocity that is a linear functional of $\langle U_2 \rangle$; see (4.8) below. Thus, similar to the case of $Re_c\ll 1$, the Jeffery-averaged lift profiles at a given finite $Re_c$, for an arbitrary aspect ratio spheroid, have the same shape as those for a sphere at the same $Re_c$. It follows that the associated pair of equilibria are identical to those for a sphere regardless of $Re_c$, and as for a sphere (Schonberg & Hinch 1989), must migrate wallward with increasing $Re_c$. It is worth reiterating that the requirement for a Jeffery-averaged analysis to remain valid becomes restrictive for extreme aspect ratio spheroids. The regime of validity, originally stated after (3.19) and expressed in terms of $Re_c$, is given by $Re_c\lambda ^2 \kappa /\ln \kappa \ll 1$ and $Re_c\lambda ^2/\kappa \ll 1$ for $\kappa \gg 1$ and $\kappa \ll 1$, respectively; these point to the restriction becoming more severe with increasing $Re_c$. The implications of the finite-$Re_c$ Jeffery-averaged analysis above, for shape sorting, remain the same as those discussed in § 3.5.

For purposes of completeness, we now follow along the lines of Schonberg & Hinch (1989), and briefly present the manner in which inertial lift is determined for $Re_c \gtrsim O(1)$. After implementing the partial Fourier transform defined in (3.24), followed by some algebraic manipulation, one obtains the coupled ODEs for $\widehat {\langle P \rangle}$ and $\widehat {\langle U_2 \rangle}$,

(4.5a)$$\begin{gather} \frac{{\rm d}^2 \widehat{\langle P\rangle}}{{\rm d}R_2^2}-k_\perp^2\widehat{\langle P\rangle}= 2\iota k_1 \widehat{\langle U_2\rangle}(\beta+2\gamma''R_2 Re_c^{-1/2}), \end{gather}$$

(4.5b)$$\begin{gather}\frac{{\rm d}^2\widehat{\langle U_2\rangle}}{{\rm d}R_2^2}-k_\perp^2\widehat{\langle U_2\rangle}= \frac{{\rm d}\widehat{\langle P\rangle}}{{\rm d}R_2}-\iota k_1 \widehat{\langle U_2\rangle}(\beta R_2+\gamma''R_2^2 Re_c^{-1/2}), \end{gather}$$

with the conditions,

(4.6a)$$\begin{gather} \widehat{\langle U_2\rangle} \sim\frac{\iota \beta\langle S_{12}\rangle k_1 |R_2|\,{\rm e}^{-k_\perp |R_2|}}{2}\quad \text{for } k_1,k_3\to\infty \text{ and } R_2\rightarrow 0, \end{gather}$$

(4.6b)$$\begin{gather}\widehat{\langle U_2\rangle}=\frac{{\rm d}\widehat{\langle U_2\rangle}}{{\rm d}R_2}=0\quad \text{at } R_2=-s\,Re_c^{1/2},R_2=(1-s)\,Re_c^{1/2}, \end{gather}$$

where $k_\perp ^2=k_1^2+k_3^2$ as before. The delta-function forcing in (4.3a) leads to the following jump conditions across the particle location ($R_2=0$):

(4.7a)$$\begin{gather} \widehat{\langle P\rangle}^+(k_1,0^+,k_3)-\widehat{\langle P\rangle}^-(k_1,0^-,k_3) =2\iota k_1 \beta\langle S_{12}\rangle, \end{gather}$$

(4.7b)$$\begin{gather}\frac{{\rm d}\widehat{\langle P\rangle}^+}{{\rm d}R_2}(k_1,0^+,k_3)- \frac{{\rm d}\widehat{\langle P\rangle}^-}{{\rm d}R_2}(k_1,0^-,k_3) =0, \end{gather}$$

(4.7c)$$\begin{gather}\widehat{\langle U_2\rangle}^+(k_1,0^+,k_3)-\widehat{\langle U_2\rangle}^-(k_1,0^-,k_3) =0, \end{gather}$$

(4.7d)$$\begin{gather}\frac{{\rm d}\widehat{\langle U_2\rangle}^+}{{\rm d}R_2}(k_1,0^+,k_3)- \frac{{\rm d}\widehat{\langle U_2\rangle}^-}{{\rm d}R_2}(k_1,0^-,k_3) =\iota k_1 \beta \langle S_{12}\rangle. \end{gather}$$

Here the superscripts ‘$+$’ and ‘$-$’ denote the limiting values attained on approaching $R_2 = 0$ from the regions $0< R_2\leq (1-s)\,Re_c^{1/2}$ and $-s\,Re_c^{1/2}\leq R_2<0$, respectively; these conditions are derived in Appendix C. The limiting form of $\langle \boldsymbol {U} \rangle$ in the matching region ($\boldsymbol {R}\ll 1$) is the sum of the singular stresslet contribution given in (4.4a), and a uniform flow along the gradient (cross-stream) direction that is a consequence of fluid inertia. The neutrally buoyant spheroid being force free is convected by this uniform flow that therefore equals the inertial lift velocity, and may be determined from the limit of the inverse transform for $\boldsymbol {R}\to 0$,

(4.8)\begin{equation} \langle V_p\rangle=\frac{Re_p}{4{\rm \pi}^2}{\rm Re}\left\{\int_{-\infty}^\infty\int_{-\infty}^\infty \,\widehat{\langle U_2\rangle}^\pm(k_1,0^\pm,k_3)\,{\rm d}k_1\, {\rm d}k_3\right\}. \end{equation}

Here, ${\rm Re}\{.\}$ denotes taking the real part that eliminates the purely imaginary stresslet contribution. As indicated, one may use either $\widehat {\langle U_2\rangle }^-$ or $\widehat {\langle U_2\rangle }^+$ on account of continuity; see (4.7c). The governing ODEs (4.5a,b), the boundary conditions (4.6a,b) and the jump conditions (4.7a–d) can be solved using the shooting technique described in Appendix A of Schmid, Henningson & Jankowski (2002), a brief description of which is given in Appendix D.

The integral in (4.8) is evaluated numerically using a two-dimensional Gauss-Legendre quadrature over a circle of a sufficiently large radius $K_m$ in the $k_1\unicode{x2013}k_3$ plane. An analytical large-$k_\perp$ asymptote was calculated using the steps outlined in Hogg (1994), and added to the numerical integral, to improve convergence. The analysis involves expanding $\widehat {\langle U_2\rangle }$ and $\widehat {\langle P\rangle }$ in inverse powers of $k_\perp$, an ansatz valid only in an $O(k_\perp ^{-1})$ neighbourhood of $R_2 = 0$. As a result, satisfaction of (4.6a,b) is replaced by a far-field decay requirement for ${R_2k_\perp \gg 1}$. Including only the exponentially decaying solutions on either side of $R_2 = 0$, satisfying the jump conditions (4.7a–d), and then performing the inverse Fourier transform, one obtains

(4.9)\begin{equation} \langle V_p\rangle^{far\ field}\approx-\frac{9\beta\gamma''\langle S_{21}\rangle Re_p}{64{\rm \pi} K_m Re_c^{1/2}}, \end{equation}

at leading order. As in § 3, the choice of $K_m$ is dictated by the distance of the particle from the walls. The relevant scale sufficiently close to either wall is still the spheroid-wall separation ($s$ or $1-s$) – the near-wall lift remains the same as that for $Re_c\ll 1$, being the far-field lift experienced by a particle in the presence of a single plane boundary. While one needs to keep increasing $K_m$ with approach to either wall (that is, for sufficiently small $s$ or $1-s$), to obtain a converged result, this increase is only necessary once $s$ or $(1-s)$ becomes less than $O(Re_c^{-1/2})$. Thus, to capture the wall-induced repulsion for finite $Re_c$, one needs to ensure $K_m Re_c^{1/2} s_{min}, K_m Re_c^{1/2} (1-s)_{min}\gg 1$. In our calculations we chose $K_m=200$ for $0.1\leq Re_c\leq 10$ to ensure accurate lift velocities down to $s, 1-s\approx 0.03$. For $Re_c>10$, accurate lift profiles were obtained down to $s, 1-s\approx 0.05$ for $K_m=40$.

4.2. Results and discussion

We begin with figure 6 that shows the inertial lift profiles for a sphere for $0.5 \leq Re_c \leq 75$, along with the limiting small-$Re_c$ profile given by (3.47); only profiles in the half-channel have been plotted owing to their antisymmetry about the centreline. The profiles for $Re_c=1$ and $75$ exhibit good agreement with the data extracted from Schonberg & Hinch (1989). Interestingly, the profiles for $Re_c = 0.5, 1$ and $10$ compare closely with the small-$Re_c$ limiting form. The inset shows the approach of the finite-$Re_c$ profiles to the near-wall value given by (3.51). Figure 7 shows that the sphere lift profiles for higher $Re_c$ agree well with those extracted from Asmolov (1999) all the way upto $Re_c=3000$; note that our profiles have been continued to smaller $s$ to emphasize the approach to the common near-wall limiting value mentioned above. The profiles highlight the rapid decrease in the lift magnitude for $Re_c\gtrsim O(10)$, reflective of weakening particle-wall interactions. Inertia-induced faster decay of the velocity field, on scales larger than $O(HRe_c^{-{1}/{2}})$, is responsible for the reduced influence of the walls with increasing $Re_c$. For $Re_c \gtrsim 300$, the lift profiles begin to exhibit an intermediate region of oppositely signed curvature. It has recently been shown that by including finite-size effects this intermediate region is pushed towards the zero-lift line, eventually leading to the emergence of new equilibria closer to the centreline for sufficiently large $Re_c$ (Anand & Subramanian 2023). The inset in figure 7 shows the collapsed lift profiles on a log-log scale, with the rescaled abscissa highlighting the $O(Re_c^{-{1}/{2}})$ neighbourhood of the wall where the lift profiles begin to rise towards to near-wall limit.

Figure 6. Comparison of finite-$Re_c$ lift velocity profiles for a sphere with the semi-analytical small-$Re_c$ profile, and the data from Schonberg & Hinch (1989). The inset shows the approach of the finite-$Re_c$ profiles towards the wall asymptote given by (3.51).

Figure 7. Comparison of lift velocity profiles for a sphere at higher $Re_c$ with Asmolov (1999) (dashed curves). The inset shows the collapse and eventual approach of the lift profiles, for large $Re_c$'s, towards the wall asymptote (horizontal dotted line) near the wall (zero crossings appear as dips to negative infinity).

Figure 8 shows the magnitude of the sphere lift velocity, as a function of $Re_c$, at different locations on either side of the Segre–Silberberg equilibrium ($s_{eq} = 0.182$). In accordance with the above discussion, for all $s$ values considered, the lift velocity starts off on a small-$Re_c$ plateau that extends until $Re_c \approx 10$. Thus, the small-$Re_c$ approximation remains a good approximation well beyond $Re_c$'s of order unity. In fact, recent results of Hood, Lee & Roper (2015) show that, for a square duct, the small-$Re_c$ analysis appears to remain valid until a duct Reynolds number of $O(80)$. For $s=0.3$ and $0.4$, with increasing $Re_c$, the lift velocity in figure 8 directly transitions from the plateau to an eventual algebraic decrease, this being typical for all $s>s_{eq}$. On the other hand, the lift velocity magnitude for $s=0.1$ (and for all $s< s_{eq}$) exhibits a non-monotonic variation with an intermediate zero crossing that, for $s=0.1$, is at $Re_c \approx 300$. This is due to the Segre–Silberberg equilibrium crossing the given $s$ in course of its wallward movement (with increasing $Re_c$). The lift velocity increases again at larger $Re_c$, but to a value smaller than the small-$Re_c$ plateau, finally transitioning to a steeper algebraic decrease.

Figure 8. Sphere lift velocity (magnitude), as a function of $Re_c$, for different $s$. In all cases, $|V_p|/Re_p$ exhibits a plateau (dashed lines) until $Re_c \approx 10$, transitioning to an algebraic decrease for sufficiently large $Re_c$. For $s < s_{eq}\ (=0.182)$, the transition is preceded by a zero crossing that appears as a sharp dip to negative infinity. The dotted lines are empirical fits to the large-$Re_c$ behaviour, and highlight the $s$-dependent decay exponent.

As mentioned earlier, the Jeffery-averaged lift profiles for a spheroid of an arbitrary aspect ratio may simply be obtained by multiplying the sphere lift profile at the same $Re_c$ by the ratio of the stresslets. Consequently, features pertaining to the $Re_c$ dependence mentioned above, including the range of validity of the small-$Re_c$ approximation, remain true for spheroids. Figures 9(a) and 9(b) show the inertial lift profiles for tumbling prolate spheroids over a range of $\kappa$, for $Re_c = 300$, with and without the large-$\kappa$ stresslet scaling. Figures 10(a) and 10(b) shows the inertial lift profiles for both spinning ($0.14\leq \kappa <1$) and spinning/tumbling ($0<\kappa <0.14)$ oblate spheroids, again for $Re_c = 300$, with and without the small-$\kappa$ stresslet scaling. Note that the latter scaling is only used for tumbling oblate spheroids in the inset of figure 10(b), since $\langle S_{12}\rangle$ remains of order unity for spinning spheroids even as $\kappa \to 0$. As in § 3, the scaled lift profiles in figure 9(b), and in the inset of 10(b), approach $\kappa$-independent limiting forms for $\kappa \to \infty$ and $0$, respectively. The magnitudes of the lift velocity profiles reflect that of the disturbance field (via the Jeffery-averaged stresslet), and therefore, decrease with increasing (decreasing) $\kappa$ for tumbling prolate (oblate) spheroids. Although not shown, the finite-$Re_c$ lift profiles may be evaluated for arbitrary $C$, and the $C$ dependence again correlates to the magnitude of the disturbance velocity field, being the largest for tumbling prolate and spinning oblate spheroids.

Figure 9. (a) Lift velocity profiles for tumbling prolate spheroids of various aspect ratios for $Re_c=300$. (b) Lift profiles in (a) rescaled using the $\kappa \to \infty$ limit of $\langle S_{12}\rangle$. The thick black curve denotes the limiting $\kappa -$independent profile for $\kappa \gtrsim 1000$.

Figure 10. (a) Lift velocity profiles for spinning oblate spheroids in the interval $0.14\leq \kappa <1$ for $Re_c=300$. (b) Lift profiles for spinning (solid lines) and tumbling (dashed lines) spheroids in the interval $0<\kappa <0.14$. The inset shows the profiles rescaled using the $\kappa \to 0$ limit of $\langle S_{12}\rangle$.

5. Conclusions

The primary result of this manuscript is to show that the inertial lift velocity profiles for neutrally buoyant spheroids in plane Poiseuille flow, for sufficiently small $Re_p$, can be calculated within a Jeffery-averaged approximation. Within this framework, spheroid lift profiles differ from those for a sphere only by a multiplicative factor, regardless of $Re_c$, which leads to the original Segre–Silberberg equilibria being dependent on $Re_c$ but not on the aspect ratio $\kappa$; this $\kappa$ independence is broadly consistent with the findings of Nizkaya et al. (2020). The multiplicative factor above is the ratio of the Jeffery-averaged spheroid stresslet to the sphere stresslet, and is a function of $C$ and $\kappa$ in the general case; for times relevant to cross-stream migration, and that are much longer than those characterizing the inertia-induced orientation drift, it is only a function of $\kappa$ since $C$ now corresponds to the inertially stabilized orbit. While the exact $\kappa$ dependence of this stresslet ratio is evidently non-trivial, for the moderate aspect ratio spinning oblate spheroids considered in Nizkaya et al. (2020), this dependence is approximately linear, consistent with the simpler $O(L^3b)$ lift force scaling proposed therein. The Jeffery-averaged approximation only requires $Re_p \ll 1$ for spheroids with $\kappa \sim O(1)$; the inertial correction to the rotation period occurs at $O(Re_p^{3/2})$ for these cases (Marath & Subramanian 2017). However, this approximation becomes increasingly restrictive for spheroids with asymptotically large and small $\kappa$, and with increasing $Re_c$. Deviations from this approximation owing to inertia-induced slow down of rotation and eventual arrest, that lead to $\kappa$-dependent equilibria for extreme aspect ratio spheroids, will be reported in a future investigation. Herein, we only note that accounting for the rotational slow down leads to the equilibria migrating back towards the centreline beyond a certain $\kappa$-dependent $Re_c$ threshold, this being consistent with both experiments (Masaeli et al. 2012) and computations (Chen, Pan & Chang 2012).

En-route to deriving the Jeffery-averaged lift on a spheroid, we presented, in some detail, the inertial lift profiles for a sphere in plane Poiseuille flow. While most of these results are known from earlier literature, and are spread across multiple efforts (Ho & Leal 1974; Vasseur & Cox 1976; Schonberg & Hinch 1989; Asmolov 1999), each of these only pertains to specific ranges of $Re_c$. Moreover, there are important features of the lift profile not reported in any of these efforts. We have expressed the small-$Re_c$ sphere lift velocity, obtained within the framework of a point-particle approximation, in terms of a one-dimensional Fourier integral, with detailed expressions for the integrands given in Appendix B. These expressions allow one to analytically determine the near-wall lift value. The arguments leading to this limiting value also show that it must be independent of $Re_c$, numerical evidence for which is provided by the finite-$Re_c$ lift profiles obtained using a shooting technique. Furthermore, the near-wall limiting lift was also identified with the far-field limit of the lift acting on a finite-sized sphere moving parallel to a single plane wall subject to a linear shearing flow. This connection allows one, in principle, to construct a uniformly valid lift profile across the entire channel. We also established the surprisingly large range of validity (up to $Re_c \approx 10$) of the small-$Re_c$ approximation derived first by Ho & Leal (1974) and Vasseur & Cox (1976).

The scaling arguments in § 3.2 are a crucial element of the overall analysis, and clearly show that the dominant scales contributing to the inertial lift remain much greater than the particle size, regardless of $Re_c$. This leads to use of a point-particle framework to obtain the leading-order approximation for the inertial lift, with scaling arguments identifying both competing contributions due to profile curvature and wall-shear-induced repulsion. The first lift force calculation for a sphere, by Ho & Leal (1974) for $Re_c\ll 1$, used the same reciprocal theorem formulation as that given here. Although in a less transparent form, the said authors, via scaling arguments, did recognize the dominant contribution of the linearized inertial terms, on scales of $O(H)$, to the reciprocal theorem volume integral. The effect of confining plane boundaries was modelled differently, and in a manner not readily generalizable to a spheroid. Rather than directly consider the relevant point singularity (Stokeslet or stresslet) between plane parallel walls, the authors used a partial Fourier representation of the full form of the unbounded domain velocity field, including both the finite-size terms and the quadrupolar disturbance induced by the quadratic component of the ambient flow, in order to derive the wall-induced contribution (the first reflection). For evaluation of the final volume integral, the physical space integration appears to have been done analytically, with the partial Fourier integral done numerically. Although more circuitous, carrying out the physical space integration should have led to a residual Fourier integral identical to the one obtained directly here from use of the convolution theorem. While the inaccuracy of the resulting profile, and the errors in the equilibrium locations, have been pointed out in the context of figure 2, the detailed calculational procedure in Ho & Leal (1974) is correct, and the origin of the errors remains uncertain. In contrast to the intuitive scaling arguments here, the later small-$Re_c$ lift calculation of Vasseur & Cox (1976) made use of a velocity field that emerged from a formal matched asymptotics expansions approach, developed in earlier papers by Cox & Brenner (1967, 1968).

The scaling arguments in § 3.2 can be generalized in several directions. First, the arguments apply virtually unchanged to pipe Poiseuille flow. Thus, for small pipe Reynolds numbers ($Re$), the dominant contributions to the inertial lift must arise from scales of the order of the pipe radius, and a calculation of this lift would involve approximating the particle as a stresslet within a cylindrical domain. The required Stokesian velocity fields should either be obtainable from that known for a Stokeslet in this domain(Liron & Shahar 1978), or be derivable in a manner similar to that in Appendix A (adapted to cylindrical coordinates). Such an analysis would yield the radius of the Segre–Silberberg annulus for $Re\to 0$. Although finite-$Re$ lift force profiles for pipe Poiseuille flow have been computed earlier (Matas et al. 2009), the profiles differ significantly for the two lowest $Re$'s examined ($1$ and $30$), preventing one from inferring the range of validity of a small-$Re$ approximation. The small-$Re$ analysis in § 3 can also be generalized to ducts of non-circular cross-sections, provided one has the confined Stokeslet field for the relevant cross-sectional geometry. The latter may be derivable for rectangular or elliptical cross-sections, using the procedure in Appendix A, using the Greens function of the (two-dimensional) Laplacian. The broken symmetry for the duct case should then allow for the prediction of discrete inertial equilibria in the transverse plane, and their variation with cross-sectional aspect ratio.

Importantly, the scaling arguments in § 3.2 allow for the systematic incorporation of finite-size effects that modify the leading-order point-particle estimate of the inertial lift velocity for both spheres and anisotropic particles in plane Poiseuille flow. As mentioned in § 3.2, and discussed in more detail in Anand & Subramanian (2023), the finite-size contributions involve an additional factor of $\lambda$, and pertain to the inner region (scales of the order of the particle size), implying that they arise independently of the outer-region point-particle contribution. It is shown in Anand & Subramanian (2023) that, for spheres, these contributions invariably become important for sufficiently large $Re_c$, leading to the emergence of a pair of new equilibria closer to the channel centreline, consistent with the results of recent experiments. In contrast to the point-particle analysis presented here, calculation of the finite-size contributions involves consideration of the nonlinear terms, and also requires knowing the disturbance velocity field induced by a sphere in an ambient quadratic flow. The analogous calculation for a spheroid will proceed along similar lines. Although more complicated, the disturbance velocity field in an ambient quadratic flow, for a spheroid of an arbitrary aspect ratio and orientation, may be constructed using a superposition of the appropriate vector spheroidal harmonics (Dabade et al. 2015, 2016). Crucially, the nonlinear inertial terms as well as the time dependence of the test velocity field, that arise in the inner region, imply that there is no longer a simple proportionality relationship between the finite-size contribution for a sphere and the Jeffery-averaged version of the same for a spheroid. As a result, one expects the incorporation of finite-size contributions to lead to $\kappa$-dependent equilibria, for neutrally buoyant spheroids, even within a Jeffery-averaged framework.

The above role of finite-size effects will be examined in a separate communication. Herein, we only mention a recent calculation of shape-selective lift forces for spheroids in an unbounded ambient quadratic flow by Bagge et al. (2021). The lift force in this case arises due to particle inertia alone, that is, for a finite Stokes number ($St$) and $Re = 0$; the force is directed towards increasing shear rate. For small $St$, the resulting lift velocity is of $O(St\,\gamma L^2)$, as also expected from dimensional considerations, with $St = \rho _p \gamma L^3/\mu$ in terms of the ambient profile curvature. Using $\rho _p = \rho _f$ for neutrally buoyant particles, and $\gamma = V_{max}/H^2$, one finds the lift velocity to be $O(Re_p V_{max}\lambda ^3)$, or $O(\lambda ^2)$ smaller than the point-particle contribution calculated here, and also $O(\lambda )$ smaller than the finite-size inner-region contribution described above. The contribution in question stems from the term in (3.16) that couples the spheroid velocity in the test problem to the translational acceleration of the spheroid in the actual problem. For this contribution to therefore be relevant, one requires sufficiently massive spheroids with $St \sim O(Re_p/\lambda ^2) \gg Re_p$; although, even in this limit, one needs to account for the angular acceleration contribution in (3.16) that explicitly invokes the presence of boundaries.