2.1. Problem set-up

In this paper we develop a unifying theory for the equation of motion of spherical particles of radius $a_p$ and mass $m_p$ immersed in general unsteady incompressible Newtonian flows of fluid density $\rho _f$ (figure 1), placing particular emphasis on fast oscillatory flows, while consistently accounting for particle inertia. The characteristic speed $U^*$ of the unsteady background flow evaluated at the particle position and kinematic viscosity $\nu$ of the fluid define the particle Reynolds number ${Re}_{p}=a_p U^*/\nu$. The particle reaction to the flow importantly also depends on the Stokes number, which we define as $\lambda =a_p^2\omega /(3\nu )$. The unsteady time scale of the flow is written as $1/\omega$, anticipating the oscillatory case, where we can alternatively use the Stokes boundary layer scale $\delta =(2\nu /\omega )^{1/2}$ to write $\lambda =2a_p^2/3\delta ^2$.

Figure 1. (a) Schematic of a spherical particle of radius $a_p$ moving with a velocity $\boldsymbol {U}_p$ as a consequence of the hydrodynamic force $\boldsymbol {F}$ exerted by the surrounding fluid. The undisturbed flow field far away from the particle is denoted by $\boldsymbol {U}$. The hydrodynamic force is generally decomposed into a force due to the undisturbed flow $\boldsymbol {F}^{(0)}$ and the disturbance flow $\boldsymbol {F}^{(1)}$ due to the presence of the particle. (b) The unsteadiness of the flow introduces the Stokes number $\lambda$, which, for oscillatory flows, is a function of the ratio of the particle size to the oscillatory boundary layer thickness $\delta$. The background flow is Taylor-expanded around the particle centre up to the quadratic term.

We follow MR in decomposing the flow around the particle into the given background flow $\boldsymbol {U}$ present without the particle, and the disturbance flow due to the particle's presence. Forces caused by the background flow will carry the superscript ${(0)}$, while those stemming from the disturbance flow will have the superscript ${(1)}$. All forces will be computed for arbitrary $\lambda$ and to first order in $Re_p$, using a regular perturbation expansion. In general flows, such an approach is valid in an inner region, while an outer region (in which inertia reasserts dominance) would have to be treated separately and the complete problem solved by asymptotic matching. However, as shown by Lovalenti & Brady (Reference Lovalenti and Brady1993), an outer region is not present when the oscillatory inertia of the disturbance flow is much greater than its advective inertia. We explain in detail below (see § 5.1) how this condition yields explicit criteria for validity across the range of $\lambda$ values. In oscillatory microfluidics, where we typically have $\lambda \sim 1\unicode{x2013}10$, the resulting criteria are comfortably fulfilled for practically relevant flows, so that it is sufficient to demonstrate the solution by regular expansion. Acoustofluidic devices generally operate at $\lambda \gg 1$, and we will show that particle motion in those cases is also quantitatively covered by our formalism. The case $\lambda \ll 1$ is usually not practically relevant, as the resulting inertial forces on particles become very small.

2.2. Particle motion and fluid flow

Our task is thus to determine explicit expressions of terms in the following equation of motion for the particle velocity $\boldsymbol {U}_p$:

(2.1)\begin{equation} m_p \frac{{\rm d} \boldsymbol{U}_p}{{\rm d}t} = \boldsymbol{F}^{(0)}+\boldsymbol{F}^{(1)} =\boldsymbol{F}^{(0)}_0+{Re}_p\boldsymbol{F}^{(0)}_1+\boldsymbol{F}^{(1)}_0+{Re}_p\boldsymbol{F}^{(1)}_1 + \cdots, \end{equation}

where subscripts denote orders of ${Re}_{p}$. Note that the decomposition of $\boldsymbol {F}^{(0)}$ is exact (there are no terms of higher order in ${Re}_{p}$; cf. MR), while we truncate the expansion of $\boldsymbol {F}^{(1)}$ at first order. Expressions for $\boldsymbol {F}^{(0)}_0$, $\boldsymbol {F}^{(1)}_0$ and part of $\boldsymbol {F}^{(0)}_1$ are contained in the MR equation, and we determine the remaining terms here.

Hydrodynamic force components are computed from the flow field stresses, as $\boldsymbol {F}^{(i)}=(F_S/6{\rm \pi} )\oint _S \boldsymbol {n}\boldsymbol{\cdot } \boldsymbol {\sigma }^{(i)} \,{\rm d}S$ with $i=0,1$, where we use the Stokes drag scale $F_S/6{\rm \pi} =\nu \rho _f a_p U^*$, and the integral is over the particle surface with its outward normal $\boldsymbol {n}$. We use lowercase letters for velocities non-dimensionalized by $U^*$ and it is advantageous in intermediate results to evaluate these velocities in a coordinate system moving with the particle centre, writing $\boldsymbol {w}^{(0)}=\boldsymbol {u}-\boldsymbol {u}_p$ for the undisturbed background flow and $\boldsymbol {w}^{(1)}$ for the disturbance flow. The dimensionless fluid stress tensors are thus written $\boldsymbol {\sigma }^{(i)}=-p^{(i)}\boldsymbol {I} + \boldsymbol {\nabla } \boldsymbol {w}^{(i)} + (\boldsymbol {\nabla } \boldsymbol {w}^{(i)})^{\rm T}$.

The Navier–Stokes equations in the particle frame of reference can be decomposed into background and disturbance components exactly (without approximations),

(2.2a)$$\begin{gather} \nabla^{2} \boldsymbol{w}^{(0)}-\boldsymbol{\nabla} p^{(0)} = 3\lambda \left(\frac{\partial \boldsymbol{w}^{(0)}}{\partial t}+\frac{{\rm d} \boldsymbol{u}_p }{{\rm d}t}\right)+{Re}_{p}\left(\boldsymbol{w}^{(0)} \boldsymbol{\cdot } \boldsymbol{\nabla}\boldsymbol{w}^{(0)}\right)\!, \quad \boldsymbol{\nabla} \boldsymbol{\cdot } \boldsymbol{w}^{(0)} =0, \end{gather}$$

(2.2b)$$\begin{gather}\boldsymbol{w}^{(0)} = \boldsymbol{u}-\boldsymbol{u}_p, \quad \text{as } r \rightarrow \infty, \end{gather}$$

(2.2c)$$\begin{gather} \begin{aligned}\nabla^{2} \boldsymbol{w}^{(1)}-\boldsymbol{\nabla} p^{(1)} &= 3\lambda \frac{\partial \boldsymbol{w}^{(1)}}{\partial t}\nonumber\\ &\quad +{Re}_{p}\left[\boldsymbol{w}^{(0)} \boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{w}^{(1)}+\boldsymbol{w}^{(1)}\boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{w}^{(0)} +\boldsymbol{w}^{(1)}\boldsymbol{\cdot } \boldsymbol{\nabla}\boldsymbol{w}^{(1)}\right]\!, \end{aligned}\end{gather}$$

(2.2d)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot } \boldsymbol{w}^{(1)} = 0, \end{gather}$$

(2.2e)$$\begin{gather}\boldsymbol{w}^{(1)} = \boldsymbol{u}_p-\boldsymbol{U}, \quad \text{on } r=1 \quad \text{and} \quad \boldsymbol{w}^{(1)} = 0, \quad \text{as } r \rightarrow \infty. \end{gather}$$

2.3. Forces from background flow

Both the $O(1)$ and $O({Re}_{p})$ components of $\boldsymbol {F}^{(0)}$ in (2.1) can be evaluated directly using the divergence theorem and the above Navier–Stokes equations valid for the background flow $\boldsymbol {w}^{(0)}$. In laboratory coordinates (see MR; Rallabandi Reference Rallabandi2021) they read

(2.3a,b)\begin{equation} \boldsymbol{F}^{(0)}_0= \frac{F_S}{6{\rm \pi}} \int_{V} \left( 3\lambda \partial_t \boldsymbol{u}\right) \,{\rm d}V,\quad \boldsymbol{F}^{(0)}_1=\frac{F_S}{6{\rm \pi}} \int_{V} \left(\boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{u}\right) \,{\rm d}V. \end{equation}

To make further progress, we need to evaluate forces due to successive orders of $\boldsymbol {w}^{(1)}$, which are ultimately also derived from the given background field $\boldsymbol {u}$. To render our solution strategy analytically tractable, we expand $\boldsymbol {u}$ around the leading-order particle position $\boldsymbol {r}_{p_0}$ into spatial moments of alternating symmetry,

(2.4)\begin{equation} \boldsymbol{u}=\boldsymbol{u}|_{\boldsymbol{r}_{p_0}} + \boldsymbol{r}\boldsymbol{\cdot } \boldsymbol{E} + \boldsymbol{r}\boldsymbol{r}:\boldsymbol{G}+\cdots,\end{equation}

where $\boldsymbol {E}=(a_p/L_\varGamma )\boldsymbol {\nabla } \boldsymbol {u}|_{\boldsymbol {r}_{p_0}}$ and $\boldsymbol {G}=\tfrac {1}{2}(a_p^2/L_\kappa ^2)\boldsymbol {\nabla }\boldsymbol {\nabla } \boldsymbol {u}|_{\boldsymbol {r}_{p_0}}$ are time-dependent, with gradient $L_\varGamma$ and curvature $L_\kappa$ length scales. Such an expansion is valid for $a_p/L_\varGamma \ll 1$ and $a_p/L_\kappa \ll 1$, conditions readily satisfied in microfluidic scenarios. Based on this, (2.3a,b) was recently evaluated analytically by Agarwal et al. (Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021) and Rallabandi (Reference Rallabandi2021), showing that an $O({Re}_p)$ contribution from $\boldsymbol {F}^{(0)}_1$ had been missed in MR, while being, in fact, of the same order as other terms in the original MR equation.

2.4. Disturbance flow: zeroth order

The Navier–Stokes equations for the disturbance flow at $O({Re}_p^0)$ read

(2.5a)$$\begin{gather} \nabla^{2} \boldsymbol{w}_0^{(1)}-\boldsymbol{\nabla} p_0^{(1)} =3\lambda \frac{\partial \boldsymbol{w}_0^{(1)}}{\partial t}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot } \boldsymbol{w}_0^{(1)} =0, \end{gather}$$

(2.5b)$$\begin{gather}\boldsymbol{w}_0^{(1)} = \boldsymbol{u}_{p_0}-\boldsymbol{u}, \quad \text{on } r=1 \quad \text{and} \quad \boldsymbol{w}_0^{(1)} = 0 ,\quad \text{as } r \rightarrow \infty. \end{gather}$$

Unlike MR, where the solution to this unsteady Stokes equation (2.5) was not explicitly needed to compute the force resulting from it, our present approach does require expressions for $\boldsymbol {w}^{(1)}_0$ to compute the full $O({Re}_p)$ force. This is accomplished by substituting the expansion (2.4) into (2.5). Each spatial moment gives rise to a linear equation with known solutions, the sum of which yields the general expression (see Landau & Lifshitz Reference Landau and Lifshitz1959; Pozrikidis et al. Reference Pozrikidis1992)

(2.6)\begin{equation} \boldsymbol{w}^{(1)}_0 = \boldsymbol{\mathcal{M}}_D \boldsymbol{\cdot } \boldsymbol{u}_s - \boldsymbol{\mathcal{M}}_Q \boldsymbol{\cdot }\left(\boldsymbol{r}\boldsymbol{\cdot } \boldsymbol{E}\right) -\boldsymbol{\mathcal{M}}_O \boldsymbol{\cdot }\left(\boldsymbol{r}\boldsymbol{r}:\boldsymbol{G}\right) + \cdots, \end{equation}

where $\boldsymbol {u}_s = \boldsymbol {u}_{p_0}-\boldsymbol {u}\vert _{\boldsymbol {r}_{p_0}}$ is the slip velocity and $\boldsymbol {\mathcal {M}}_{D,Q,O}(r,\lambda )$ are mobility tensors with known spatial dependence. For oscillatory flows, their dependence on the Stokes number $\lambda$ is known analytically. Explicit expressions for these tensors are given in Appendix A.

2.5. Disturbance flow: first order using a reciprocal theorem

Fast oscillatory particle motion can give rise to large disturbance flow gradients, so that terms involving $\boldsymbol {\nabla }\boldsymbol {w}^{(1)}$ on the right-hand side of (2.2d) are not necessarily negligible compared with the viscous diffusion term on the left-hand side, and $O({Re}_p)$ force terms in $\boldsymbol {F}^{(1)}$ become important.

With $\boldsymbol {w}^{(1)}_0$ explicitly known, the equations at $O({Re}_p)$ read

(2.7a)$$\begin{gather} \nabla^{2} \boldsymbol{w}^{(1)}_1-\boldsymbol{\nabla} p^{(1)}_1 =\boldsymbol{\nabla} \boldsymbol{\cdot } \boldsymbol{\sigma}^{(1)}_1=3\lambda \frac{\partial \boldsymbol{w}^{(1)}_1}{\partial t}+\boldsymbol{f}_0, \quad \boldsymbol{\nabla} \boldsymbol{\cdot } \boldsymbol{w}^{(1)}_1 =0, \end{gather}$$

(2.7b)$$\begin{gather}\boldsymbol{w}^{(1)}_1 =-\boldsymbol{u}_{p_1}, \quad \text{on } r=1 \quad \text{and} \quad \boldsymbol{w}^{(1)}_1 = 0, \quad \text{as }\quad r \rightarrow \infty , \end{gather}$$

with $\boldsymbol {f}_0=\boldsymbol {w}^{(0)}_0\boldsymbol{\cdot } \boldsymbol {\nabla } \boldsymbol {w}^{(1)}_0+\boldsymbol {w}^{(1)}_0\boldsymbol{\cdot } \boldsymbol {\nabla } \boldsymbol {w}^{(0)}_0 +\boldsymbol {w}^{(1)}_0\boldsymbol{\cdot } \boldsymbol {\nabla }\boldsymbol {w}^{(1)}_0$ as the leading-order nonlinear forcing.

In order to compute the force $\boldsymbol {F}^{(1)}_1$, we do not solve for the flow field $\boldsymbol {w}^{(1)}_1$ in (2.7) but instead employ a reciprocal relation in the Laplace domain. The reciprocal theorem infers the force from the known stress of a test flow $\boldsymbol {u}'$, which is here chosen to be a dipolar unsteady Stokes flow around the particle with arbitrary directionality $\boldsymbol {e}$; see Appendix B for a detailed derivation. We obtain for the magnitude of the force along ${\boldsymbol {e}}$,

(2.8)\begin{equation} \boldsymbol{e}\boldsymbol{\cdot } \boldsymbol{F}^{(1)}_1 = \frac{F_S}{6{\rm \pi}}\mathcal{L}^{-1}\left\{\int_{S_p}\frac{\hat{\boldsymbol{u}}_{p_1}}{\hat{u}'}\boldsymbol{\cdot } \left( \hat{\boldsymbol{\sigma}}'\boldsymbol{\cdot } \boldsymbol{n}\right)\,{\rm d}S - {Re}_p\int_V \frac{\hat{\boldsymbol{u}}' \boldsymbol{\cdot } \hat{\boldsymbol{f}}_0}{\hat{u}'}\,{\rm d}V \right\}\!,\end{equation}

where the hat denotes the Laplace transform and $\mathcal {L}^{-1}$ is the inverse Laplace transform. When applied at $O(1)$, this reciprocal-theorem strategy similarly yields

(2.9)\begin{equation} \boldsymbol{e}\boldsymbol{\cdot } \boldsymbol{F}^{(1)}_0 =F^{(1)}_0= \frac{F_S}{6{\rm \pi}}\mathcal{L}^{-1}\left\{\int_{S_p}\frac{\hat{\boldsymbol{w}}_0^{(0)}}{\hat{u}'}\boldsymbol{\cdot } \left( \hat{\boldsymbol{\sigma}}'\boldsymbol{\cdot } \boldsymbol{n}\right)\,{\rm d}S\right\} , \end{equation}

which is precisely the force expression obtained by MR. Since the variable in the overall equation of motion (2.1) is the unexpanded particle velocity $\boldsymbol {u}_{p}$, we make the substitution $\boldsymbol {w}^{(0)}_0=\boldsymbol {w}^{(0)} - {Re}_p \boldsymbol {u}_{p_1}+O({Re}_p^2)$. Adding (2.9) and (2.8), the $O({Re}_p)$ term in (2.9) exactly cancels the first term in (2.8) and produces a correction term that is $O({Re}_p^2)$. The net force on the particle due to its disturbance flow then reads

(2.10a)$$\begin{gather} \boldsymbol{e}\boldsymbol{\cdot } \boldsymbol{F}^{(1)} = \boldsymbol{e}\boldsymbol{\cdot } \left(\boldsymbol{F}^{(1)}_0+{Re}_p\boldsymbol{F}^{(1)}_1\right) +O({Re}_p^2), \end{gather}$$

(2.10b)$$\begin{gather}\boldsymbol{e}\boldsymbol{\cdot } \boldsymbol{F}^{(1)}_0= \frac{F_S}{6{\rm \pi}}\mathcal{L}^{-1}\left\{\int_{S_p}\frac{\hat{\boldsymbol{w}}^{(0)}}{\hat{u}'}\boldsymbol{\cdot } \left( \hat{\boldsymbol{\sigma}}'\boldsymbol{\cdot } \boldsymbol{n}\right)\,{\rm d}S\right\}\!, \end{gather}$$

(2.10c)$$\begin{gather}\boldsymbol{e}\boldsymbol{\cdot } \boldsymbol{F}^{(1)}_1= \frac{F_S}{6{\rm \pi}}\mathcal{L}^{-1}\left\{- {Re}_p\int_V \frac{\hat{\boldsymbol{u}}' \boldsymbol{\cdot } \hat{\boldsymbol{f}}_0}{\hat{u}'}\,{\rm d}V \right\}\!, \end{gather}$$

where we have also replaced $\boldsymbol {w}^{(0)}_0$ by $\boldsymbol {w}^{(0)}$ in $\boldsymbol {f}_0$, resulting in an error that is again $O({Re}_p^2)$. Only certain products in $\boldsymbol {f}_0$ are non-vanishing when the angular integration is performed due to alternating symmetry of terms in the background flow field multipole expansion (2.4). These non-zero terms are conveniently labelled by the multipole orders involved in the product

(2.11)\begin{equation} {Re}_p\frac{F_S}{6{\rm \pi}}\mathcal{L}^{-1}\left\{-\int_V \frac{\hat{\boldsymbol{u}}' \boldsymbol{\cdot } \hat{\boldsymbol{f}}_0}{\hat{u}'}\,{\rm d}V\right\}= F_{\sigma \varGamma}^{(1)}+ F_{\varGamma\kappa}^{(1)} + \cdots.\end{equation}

Here, $F_{\sigma \varGamma }^{(1)}$, $F_{\varGamma \kappa }^{(1)}$ are the inertial force contributions obtained by successive contractions of adjacent tensors involving $\boldsymbol {u}_s$ (index $\sigma$), $\boldsymbol {E}$ (index $\varGamma$), $\boldsymbol {G}$ (index $\kappa$) and so on. The volume integral is tedious but straightforward to compute since all the integrations resulting from the leading-order velocity fields (2.4) and (2.6) are convergent. The evaluation of the Laplace transforms can be performed analytically if the flow has harmonic time dependence. This is not a severe restriction as arbitrary time dependences can be decomposed into harmonic contributions. Rectified forces resulting from averaging time-periodic processes can then be computed separately by frequency. Such cases of oscillatory flow are the most relevant in practical applications and are focused on in § 3 below. To simplify notation, we therefore assume a single oscillatory frequency $\omega$ in the following, without loss of generality.

When the particle is neutrally buoyant, the first term $F_{\sigma \varGamma }^{(1)}$ vanishes so that the leading term is $F_{\varGamma \kappa }^{(1)}$, which was derived in Agarwal et al. (Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021) as an unexpected inertial force for density-matched particles. This term (involving the product $\boldsymbol {E} : \boldsymbol {G}$) has no analogue in the previous literature and for completeness, we reproduce it as follows for harmonic oscillatory flows $\boldsymbol {U}$:

(2.12)\begin{equation} F_{\varGamma \kappa}^{(1)}=m_f a_p^2 \left[\boldsymbol{\nabla}\boldsymbol{U}:\boldsymbol{\nabla}\left(\boldsymbol{\nabla} \boldsymbol{U}\right)\right]\boldsymbol{\cdot } \boldsymbol{e} \mathcal{F}^{(1)}_1. \end{equation}

The $\lambda$-dependent dimensionless function $\mathcal {F}^{(1)}_1$ results from the volume integration, which also yields $m_f$, the mass of fluid displaced by the particle, via $(4{\rm \pi} a_p^3/3){Re}_p F_S/(6{\rm \pi} ) = m_f a_p^2 (U^*)^2$. In the next section, we follow a similar strategy for non-neutrally buoyant particles.

2.6. Disturbance flow: evaluation of $F_{\sigma \varGamma }^{(1)}$

Non-neutrally buoyant particles have a slip velocity and thus a non-zero $F_{\sigma \varGamma }^{(1)}$, involving the product $\boldsymbol {u_s} \boldsymbol{\cdot } \boldsymbol {E}$. Appropriate to fast harmonic oscillatory flows, we approximate the background flow as a potential flow with a given single frequency. The slip velocity as a linear response can then be generally decomposed into an in-phase and an out-of-phase component with respect to the background flow, i.e. $\boldsymbol {u}_s(\boldsymbol {r}_p,t)=\boldsymbol {u}_{s}^{I}(\boldsymbol {r}_p,t)+\boldsymbol {u}_{s}^{O}(\boldsymbol {r}_p,t)$. The corresponding force is written as

(2.13)\begin{equation} \frac{F_{\sigma \varGamma}^{(1)}}{{Re}_p F_S/(6{\rm \pi})}=\frac{4{\rm \pi}}{3}(\boldsymbol{u}_{s}^{I}\boldsymbol{\cdot }\boldsymbol{E}) \boldsymbol{\cdot } \boldsymbol{e}\, \mathcal{{G}}_1(\lambda) + \frac{4{\rm \pi}}{3}(\boldsymbol{u}_{s}^{O}\boldsymbol{\cdot }\boldsymbol{E}) \boldsymbol{\cdot } \boldsymbol{e} \mathcal{{G}}_2(\lambda) ,\end{equation}

where the $\mathcal {{G}}_1$ and $\mathcal {{G}}_2$ terms are explicit outcomes of the volume integration in (2.11) and capture the $\lambda$-dependence of the in- and out-of-phase contributions, respectively. For fast oscillatory background flows, we can replace the in-phase component with $\boldsymbol {u}_{s}\boldsymbol{\cdot }\boldsymbol {E}$ and the out-of-phase component with $\partial _t \boldsymbol {u}_{s}\boldsymbol{\cdot }\boldsymbol {E}$ (see Appendix C for details), resulting in

(2.14)\begin{align} F_{\sigma \varGamma}^{(1)} &= \frac{4{\rm \pi}}{3}\rho_f a_p^2 U^{*^2} \left( \boldsymbol{u}_s\boldsymbol{\cdot }\boldsymbol{E} \boldsymbol{\cdot } \boldsymbol{e} \mathcal{G}_1(\lambda) + \partial_t{\boldsymbol{u}_s}\boldsymbol{\cdot }\boldsymbol{E} \boldsymbol{\cdot } \boldsymbol{e} \mathcal{G}_2(\lambda) \right)\nonumber\\ &= m_f \left[ (\boldsymbol{U}_p -\boldsymbol{U})\boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{U}\right]\boldsymbol{\cdot } \boldsymbol{e} \mathcal{G}_1(\lambda) + m_f \left[ \partial_t(\boldsymbol{U}_p -\boldsymbol{U}) \boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{U}\right] \boldsymbol{\cdot } \boldsymbol{e} \frac{\mathcal{G}_2(\lambda)}{\omega}. \end{align}

While the exact, lengthy expressions for the universal functions $\mathcal {G}_{1,2}$ are given in Appendix C, an excellent uniformly valid solution can be constructed by simply adding the leading orders of the small and large $\lambda$ expansions of $\mathcal {G}_1$ (analogous to the function $\mathcal {F}$ in Agarwal et al. (Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021)). Taylor expansion in both the viscously dominated limit ($\lambda \to 0$) and the inviscid limit ($\lambda \to \infty$) obtains

(2.15a,b)\begin{equation} \mathcal{G}_1^{v} =-\frac{63}{80}\sqrt{\frac{3}{2\lambda}}+O(1),\quad \mathcal{G}_1^{i}=-\frac{1}{2} + O(1/\sqrt{\lambda}), \end{equation}

from which the following uniformly valid result is constructed:

(2.16)\begin{equation} \mathcal{G}_1^{uv}(\lambda) \approx-\left(\frac{1}{2} +\frac{63}{80}\sqrt{\frac{3}{2\lambda}}\right)\!. \end{equation}

Figure 2(a,b) illustrates that this simple two-term expression agrees very well with the full result (C3) over the entire range of the parameter $\lambda$, with a maximum error of ${\sim }6\,\%$.

Figure 2. (a) Plot of the in-phase force function $\mathcal {G}_1$. The uniformly valid expression (purple dashed) closely tracks the full solution (red). Also displayed are the viscous (green) and inviscid (blue) limit asymptotes. (b) The magnitude of the percentage error between the uniformly valid and full solutions is small throughout the entire range of $\lambda$, with a maximum error of ${\sim }6\,\%$. (c) Plot of the out-of-phase force function $\mathcal {G}_2$ (red) together with its viscous (green) and inviscid (blue) limit expressions.

A Taylor expansion of the out-of-phase term $\mathcal {G}_2$, in the viscous and inviscid limits, respectively, results in

(2.17a,b)\begin{equation} \mathcal{G}_2^v = \frac{3}{16}\sqrt{\frac{3}{2\lambda}} +O(1),\quad \mathcal{G}_2^i =-\frac{57}{40}\sqrt{\frac{3}{2\lambda}} + O(1/\lambda). \end{equation}

Both of the above expansions have a $O(1/\sqrt {\lambda })$ leading-order term and a simple, two-term approximation fails. In the following, we use the full expression (C4), noting that the contribution from this term is small in most practical situations, i.e. when $\lambda \gtrsim 1$.

3.1. Generalized Auton correction

We first comment on the limiting case of the well-known correction to MR due to Auton (Auton et al. Reference Auton, Hunt and Prud'Homme1988). The equation of motion derived by MR deviated from previous versions in the form of the convective term in (3.1b), using $m_f (\boldsymbol {U}\boldsymbol{\cdot }\boldsymbol {\nabla } \boldsymbol {U})$ instead of $m_f (\boldsymbol {U}_p\boldsymbol{\cdot }\boldsymbol {\nabla } \boldsymbol {U})$ – the values of these two derivatives can differ substantially when the Reynolds number is not small. Similarly, Auton et al. (Reference Auton, Hunt and Prud'Homme1988) showed that in the limit of potential flows, the added mass term should read $\tfrac {1}{2}m_f ({{\rm d}\boldsymbol {U}_p}/{{\rm d}t}-{{\rm D}\boldsymbol {U}}/{{\rm D}t})$ instead of $\tfrac {1}{2}m_f ({{\rm d}\boldsymbol {U}_p}/{{\rm d}t}-{{\rm d}\boldsymbol {U}}/{{\rm d}t})$. Again, these two expressions are identical in the zero Reynolds number limit employed by MR, but in flows with substantial inertial effects, they can differ significantly.

Our formalism naturally addresses these concerns through the rigorous treatment of the disturbance flow around the particle. The first term on the right-hand side of (3.1d), involving $\mathcal {G}_1$, modifies the added mass term in (3.1c) and reproduces the Auton correction (Auton et al. Reference Auton, Hunt and Prud'Homme1988) in the inviscid, potential flow limit ($\lambda \to \infty$, ${Re}_p\ll 1$), modifying ${{\rm d} \boldsymbol {U}}/{{\rm d}t}$ to ${{\rm D} \boldsymbol {U}}/{{\rm D}t}$, or explicitly,

(3.2)$$\begin{gather} -\frac{1}{2}m_f \frac{{\rm d}}{{\rm d}t}\left[\boldsymbol{U}_p-\boldsymbol{U}\right] + m_f \left[ (\boldsymbol{U}_p -\boldsymbol{U})\boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{U}\right] \mathcal{G}_1(\lambda)\nonumber\\ \approx-\frac{1}{2}m_f \left[\frac{{\rm d}}{{\rm d}t}\boldsymbol{U}_p-\frac{{\rm D}}{{\rm D}t}\boldsymbol{U}\right] - \frac{63}{80}\sqrt{\frac{3}{2\lambda}} m_f \left[ (\boldsymbol{U}_p -\boldsymbol{U})\boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{U}\right]\!, \end{gather}$$

where we use the simple two-term approximation (2.16) for $\mathcal {G}_1$. Thus, instead of heuristically modifying the added mass term, our approach rigorously derives its dependence on $\lambda$. Note that in most practically relevant oscillatory microfluidic flows, the value of $\lambda$ is $O(1\unicode{x2013}10)$, so that the contribution from the second term of (3.2) – capturing the effect of viscous streaming around the particle – results in the inertial force being quite large due to the $1/\sqrt {\lambda }$ scaling.

We note that the second term in (3.1d), involving $\mathcal {G}_2$, arises due to the out-of-phase component of the slip velocity and thus characterizes diffusion of vorticity from the particle. This term is analogous to the Basset–Boussinesq history force and contributes most prominently when $\lambda \sim O(1)$, while it is subdominant for both small and large $\lambda$.

3.2. Time scale separation and connection to acoustofluidics

Equation (3.1) describes unsteady particle dynamics as an integral equation containing a history integral, which can be explicitly evaluated in special cases, particularly for particles executing purely oscillatory motion. In more general settings, where there is a superposition of slower rectified or transport fluid flows – with a clear separation of scales from the fast oscillatory motion – we can still find an explicit, analytical evaluation of the memory integral by employing the method of multiple scales. This approach results in a simple overdamped equation of motion for the particle that captures the slow dynamics accurately, as outlined in the following (see Appendix D for details).

For flows induced by a localized oscillating source with curvature scale $a_b$, amplitude $\epsilon a_b$ and angular frequency $\omega$, we non-dimensionalize our equation with $a_b$, $\epsilon a_b \omega$ and $1/\omega$ as characteristic length, velocity and time scales, respectively. Equation (3.1) then reads

(3.3)\begin{align} \lambda \left( \hat{\kappa} + 1\right)\frac{{\rm d}^2 \boldsymbol{r}_p}{{\rm d} t^2} &= \epsilon\lambda \frac{\partial \boldsymbol{u}}{\partial t} + \frac{2\lambda}{3} \epsilon^2 \boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{u} - \frac{\lambda}{3}\epsilon^2 \lambda \boldsymbol{U}_p \boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{u} -\left(\frac{{\rm d} \boldsymbol{r}_p}{{\rm d} t}-\epsilon\boldsymbol{u}\right)\nonumber\\ &\quad + \sqrt{\frac{3\lambda}{\rm \pi}}\int_{-\infty}^t \frac{{\rm d}/{\rm d}\tau \left[{\rm d}\boldsymbol{r}_p(\tau)/{\rm d}\tau - \epsilon\boldsymbol{u}(\boldsymbol{r}_p(\tau),\tau)\right]}{\sqrt{t-\tau}}\,{\rm d}\tau\nonumber\\ &\quad + \frac{2\lambda}{3}\epsilon \mathcal{G}_1 \left(\frac{d \boldsymbol{r}_p}{d t}-\epsilon\boldsymbol{u}\right)\boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{u} + \frac{2\lambda}{3}\epsilon\, \mathcal{G}_2 \partial_t \left(\frac{{\rm d} \boldsymbol{r}_p}{{\rm d} t}-\epsilon\boldsymbol{u}\right)\boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{u} \nonumber\\ &\quad + \frac{2\lambda}{3}\epsilon^2 \alpha^2 \mathcal{F} \boldsymbol{\nabla}\boldsymbol{u}:\boldsymbol{\nabla} \boldsymbol{\nabla} \boldsymbol{u}, \end{align}

where $\hat {\kappa }=2/3({\rho _p}/{\rho _f}-1)$ is a dimensionless measure of density difference, $\alpha =a_p/a_b$ is the relative particle size and ${{\rm d} \boldsymbol {r}_p}/{{\rm d}t} =\epsilon \boldsymbol {u}_p$. As in Agarwal et al. (Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021), we write $\mathcal {F}=\mathcal {F}^{(0)}_1+\mathcal {F}^{(1)}_1$.

We employ standard techniques of time scale separation (see Appendix D) to obtain the leading-order overdamped equation of particle motion. Briefly, the fast oscillatory dynamics in (3.3) are time-averaged over the oscillation period and the resulting equation describes the dynamics of the leading-order mean particle position $\boldsymbol {r}_{p_0}$ on the slow time scale $T=\epsilon ^2 t$,

(3.4)\begin{equation} \frac{{\rm d} \boldsymbol{r}_{p_0}}{{\rm d} T}=\frac{\hat{\kappa}\lambda}{(\hat{\kappa}+1)} \mathcal{G}(\lambda) \left\langle \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{u}\right\rangle + \frac{2\lambda}{3}\alpha^2 \mathcal{F}(\lambda) \left\langle \boldsymbol{\nabla} \boldsymbol{u}:\boldsymbol{\nabla} \boldsymbol{\nabla} \boldsymbol{u} \right\rangle ,\end{equation}

with $\mathcal {F}(\lambda )\approx \tfrac {1}{3} + \tfrac {9}{16}\sqrt {{3}/{2\lambda }}$ derived in Agarwal et al. (Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021) and

(3.5)\begin{equation} \mathcal{G}(\lambda)=\frac{(\hat{\kappa}+1)(2(1-\mathcal{G}_1)(d+\hat{\kappa})\lambda^2 +c\left(2\lambda\mathcal{G}_2-3\right))}{3(c^2+(d+\hat{\kappa})^2\lambda^2)},\end{equation}

where $c = 1+\sqrt {3\lambda /2}$ and $d=1+\sqrt {3/(2\lambda )}$ are expressions resulting from the integration of the history force term (cf. Appendix D for details).

The first term on the right-hand side of (3.4) can be rewritten as $F_{R}\mathcal {G}(\lambda )$, where $F_{R} = ({\hat {\kappa }\lambda }/{(\hat {\kappa }+1)}) \langle \boldsymbol {u} \boldsymbol{\cdot } \boldsymbol {\nabla } \boldsymbol {u}\rangle$ is a time-averaged force formally identical to the acoustic radiation force induced by an incident sound field with velocity field $\boldsymbol {u}$ (Bruus Reference Bruus2012). In the acoustofluidic context, this velocity field may be caused by an oscillating object (bubble) excited by a primary acoustic wave. The resulting force from the bubble on a distant particle is then often denoted as the secondary radiation force (Doinikov & Zavtrak Reference Doinikov and Zavtrak1996). As the acoustic formalism is based on the assumption of inviscid flow, $\mathcal {G}(\lambda )$ generalizes the far-field inviscid $F_{R}$ to include viscous effects that, as shown below, can change the resulting particle motion quantitatively and qualitatively. Note that $\mathcal {G}(\lambda \to \infty )=1$, recovering the inviscid case, while the viscous limit depends on the density contrast, $\mathcal {G}(\lambda \to 0)=-(1+\hat {\kappa })$.

In the next section, we specialize (3.4) to the simplest case of a background flow induced by a volumetrically oscillating object – a situation commonly encountered in many practical microfluidic set-ups involving acoustically excited microbubbles – and compare our results with DNS.

4.1. Simulation approach and results

To computationally simulate the relevant flow scenarios, we employ an axisymmetric formulation of the incompressible Navier–Stokes equations (see Appendix E). Figure 3(a) presents the simulation set-up. A spherical particle of radius $a_p$ is initially released with zero velocity at a distance $r_{p0}$ from the oscillating monopole. It is thus exposed to the model flow of frequency $\omega$ and velocity amplitude $\epsilon \omega$ (the nominal source size $a_b$ is normalized to 1). We choose $\epsilon =0.01$, $a_p=0.05$, $\omega = 16{\rm \pi}$ throughout, and $r_{p0}=2$ unless otherwise stated. The fluid viscosity is determined from the corresponding values of $\lambda$ in each simulation. The top half of the panel in figure 3(a) shows representative streamlines of the instantaneous, near-radial flow, while the bottom half shows time-averaged streamlines, highlighting the ensuing steady, rectified flow pattern. Varying the ratio of particle density to fluid density in figure 3(c–e) while keeping all other parameters constant shows that the direction of this rectified flow reverses, but not for matching densities – rather, the flow pattern loses directionality around $\rho _p/\rho _f\approx 0.95$.

Figure 3. Direct numerical simulation of the prototypical problem: (a) a spherical particle of radius $a_p$ is exposed to an oscillating monopole placed at a distance $r_p$ from the particle centre with primary flow velocity $U^*$. Top half of the panel are instantaneous streamlines (the colourbar is flow speed in units of $U^*$); bottom half of the panel are the time-averaged streamlines (the colourbar is steady flow speed in units of $\epsilon U^*$). (b) Particle coordinate as a function of time. For three different density contrasts, we show the full oscillatory dynamics (see inset for a close-up) as well as the steady particle motion (averaged once per oscillation cycle). (c–e) Time-averaged flow fields around the particle for the three cases of (b).

Accordingly, the particle motion in the simulation (figure 3b) reverses direction: particles lighter than $\approx 0.95\rho _f$ are repelled over time, while those of greater density are attracted towards the monopole (which includes the density-matched case).

4.2. Comparison of particle trajectories

A comparison between unsteady DNS dynamics and predictions from the unsteady theory equation (3.3) is possible, although it entails evaluating the non-local Basset memory integral, which is computationally expensive and typically not of relevance in applications. For a clearer and more practical validation, in figure 4 we focus on comparing time-averaged DNS dynamics and predictions from the analytically derived equation (4.1) for the rectified steady dynamics, which is easily integrated in time.

Figure 4. Comparison of theoretical particle motion with DNS. (a–d) Time-averaged dynamics from the theory using (4.1) with the full analytical expressions for $\mathcal {G}$ and $\mathcal {F}$ agree with DNS (magenta) for the entire range of $\lambda$ and density contrast values (all results are for $r_{p0}=2$). Two exemplary density contrast and $\lambda$ combinations are displayed, $\rho _p/\rho _f=1.1$ ($\hat {\kappa }=0.067)$ and $\rho _p/\rho _f=0.9$ ($\hat {\kappa }=-0.067)$. The classical MR equation solutions (green) fail to even qualitatively capture the particle repulsion in (a), and (b–d) otherwise strongly underestimate the force. The inviscid formalism of Agarwal, Rallabandi & Hilgenfeldt (Reference Agarwal, Rallabandi and Hilgenfeldt2018) (light blue) has similar, though quantitatively less severe, shortcomings. (e) Best-fits of $\mathcal {G}(\lambda )$ to (4.1) are extracted from DNS and show excellent agreement with the full theory (3.5), for both heavier ($\rho _p/\rho _f=1.1$, red) and lighter ($\rho _p/\rho _f=0.9$, teal) particles.

Figure 4(a–d) depict examples of such averaged radial dynamics for different density ratios and different Stokes numbers $\lambda$, all with $r_{p0}=2$. Across a wide range of parameters, DNS dynamics (magenta) and analytical results (red) from the uniformly valid asymptotic expressions of $\mathcal {G}(\lambda )$ are found to be in very good agreement. Predictions from the classical MR equation (green) instead deviate significantly in all cases and, for some parameter combinations (see figure 4a), even misidentify the direction of the particle motion. The theory of Agarwal et al. (Reference Agarwal, Rallabandi and Hilgenfeldt2018) (light blue), which relies on inviscid flow throughout, also misses important force contributions and shows deviations similar in nature to those of MR, though quantitatively smaller. Only properly accounting for particle inertia successfully reproduces the range of numerically observed behaviours.

To illustrate the success of (4.1) over the entire range of practically relevant $\lambda$ values, figure 4(e) condenses all results by extracting a $\mathcal {G}(\lambda )$ value from best-fitting (4.1) to the numerically simulated particle trajectories (see Appendix F for details), given the previously established accuracy of the $\mathcal {F}(\lambda )$ function (Agarwal et al. Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021). Both for heavier ($\rho _p/\rho _f=1.1$, red) and lighter particles ($\rho _p/\rho _f=0.9$, teal), the analytical equation yields excellent agreement with the simulated rectified drift of the particle, indicating that it captures the key physical mechanisms at play. We note here that even for $\lambda =20$, there are significant deviations of $\mathcal {G}(\lambda )$ from its inviscid asymptotic value of 1, showing that viscous effects remain important in quantitative device design even at large Stokes numbers.

This validation demonstrates the utility of our theoretical framework in predicting the dynamics of solid particles in oscillatory flows, as each individual DNS simulation incurs a large computational cost up to $\sim$24–48 core hours on a single node on the Expanse supercomputer (see Appendix E), while the theory ODE is trivial to solve.

4.3. Particles at large distances: connection to acoustofluidics

Acoustofluidics has been a fruitful field of study aiming to manipulate fluid and particles using acoustic waves (Bruus Reference Bruus2012; Friend & Yeo Reference Friend and Yeo2011). As mentioned above, our framework specializes to the far-field acoustofluidic secondary radiation force when the distance between the particle and the oscillating source is large, $r_{p_0}\gg 1$. In this case, the force on the particle is the first term of (3.4), i.e. the nominal inviscid acoustic radiation force $F_{R}$ multiplied by the Stokes-number-dependent factor $\mathcal {G}(\lambda )$. That such a $\lambda$-dependence exists has been known in acoustofluidics, and several approaches have been used to quantify it. We compile these predictions in figure 5(a) for reference.

Figure 5. (a) Stokes number dependence of the overall dimensionless inertial force magnitude $\mathcal {G}$, representing the ratio between acoustofluidic forces (limit of large distance between source and particle) to the radiation force $F_R$. Lines are results from different theories, symbols from DNS, all for $\rho _p/\rho _f = 1.1$ ($\hat {\kappa }=0.067$), $\epsilon =0.01$. The DNS values are best fits of $\mathcal {G}$ given the full expression for $\mathcal {F}$ in (4.1). The present work (red line) is in excellent agreement with all DNS data, while both the Agarwal et al. (Reference Agarwal, Rallabandi and Hilgenfeldt2018) (light blue) and MR formalisms (green) significantly underestimate the forces. (b) Contour plots for steady particle velocity at $r_p=2$ with varying $\lambda$ and $\rho _p/\rho _f$. The solid red line marks the transition from attraction to repulsion. Solid circles indicate simulation outcomes with blue and red circles representing attraction and repulsion, respectively.

Predictions using the MR equation fail to correctly reproduce the inviscid limit ($\lambda \to \infty$) due to the incorrect form of the fluid acceleration in the added mass term (see the discussion of the Auton correction in § 3.1). The formalism of Settnes & Bruus (Reference Settnes and Bruus2012) instead misses the opposite viscous limit ($\lambda \to 0$), as it ignores viscosity completely. In previous work (Agarwal et al. Reference Agarwal, Rallabandi and Hilgenfeldt2018), the present authors heuristically combined the leading-order inviscid and viscous effects. This simplified formalism agrees exactly with the much more elaborate theory of Doinikov (Reference Doinikov1994) in both the viscously dominated ($\lambda \to 0$) and the inviscid limits ($\lambda \to \infty$), while quantitative discrepancies remain in the intermediate $\lambda$ regime, where the $\mathcal {G}(\lambda )$ of Doinikov (Reference Doinikov1994) is larger than that of Agarwal et al. (Reference Agarwal, Rallabandi and Hilgenfeldt2018).

The theory of the present work agrees with the previously established viscous and inviscid limits, and makes new predictions in the intermediate $\lambda$ range, with values between those of Doinikov (Reference Doinikov1994) and Agarwal et al. (Reference Agarwal, Rallabandi and Hilgenfeldt2018). The DNS data in figure 5(a) demonstrate that our theory is in excellent agreement with the forces observed in a full Navier–Stokes simulation, significantly improving on all previous approaches. The relative error between our analytical predictions and the DNS is ${\approx } 5\,\%\unicode{x2013}10\,\%$ across the simulation range $1\leq \lambda \leq 20$. This range is restricted by computational cost on both ends, as the boundary layer scale becomes large for $\lambda \ll 1$ and the size of the computational domain must be increased, while for $\lambda \gg 1$ the extremely thin boundary layer requires ever greater mesh refinement. The diverging computational cost underscores again the advantage of using analytical theory to predict particle behaviour at high Stokes numbers and into the regime of acoustofluidic applications.

Our results reaffirm that viscous effects can significantly affect the behaviour of particles in acoustofluidic systems, and have important implications for the design and optimization of microfluidic devices that utilize acoustic waves for particle manipulation.

4.4. Transition from attraction to repulsion

Equation (4.1) predicts that particles in monopolar oscillatory flows can exhibit equilibrium positions (at finite $r$) where the net force acting on the particle is zero. Setting ${{\rm d} r_p}/{{\rm d} T}=0$ in (4.1) obtains the critical radial position (in units of the particle radius $a_p$) as

(4.2)\begin{equation} r_{p_c}= \sqrt{-\frac{(\hat{\kappa}+1)\mathcal{F}(\lambda)}{\hat{\kappa}\mathcal{G}(\lambda)}}. \end{equation}

In most practically relevant situations, $\lambda \gtrsim O(1)$, and thus $\mathcal {G}>0$ (cf. figure 5a). A real $r_{p_c}$ then exists if the particle is lighter than the surrounding medium ($\hat {\kappa }<0$). Such an equilibrium position is necessarily unstable, as the repulsive term in (4.1) decays more slowly. Thus, for light particles and $\lambda \gtrsim O(1)$ this model predicts a critical radial distance below which the particle is always attracted towards the oscillating source. In a practical set-up, a particle can be transported into this attractive range by streaming flows or other appropriately designed flow fields. Thus, $r_{p_c}$ is an important quantity to consider in the design of microfluidic devices that make use of acoustically excited microbubbles to selectively trap particles (cf. Chen et al. Reference Chen, Fang, Merritt, Strack, Xu and Lee2016; Zhang et al. Reference Zhang, Song, Bai, Guo, Feng and Arai2021a,Reference Zhang, Song, Bai, Jia, Song, Guo and Fengb).

Conversely, given a certain distance from the oscillating object, attraction or repulsion of a particle can be designed by adjusting density contrast or Stokes number (oscillation frequency). Figure 5(b) plots the isolines of the right-hand side of (4.1) as a function of the parameters $\lambda$ and $\hat {\kappa }$, for a fixed $r_{p}=2$. The red line is the zero contour separating attractive from repulsive dynamics. Particles of density equal or higher than fluid density are always attracted towards the source, while light particles ($\rho _p/\rho _f<1$) are repelled above a threshold Stokes number. Comparison with DNS data (circles in figure 5b) confirms these predictions. The sign change of $\mathcal {G}(\lambda )$ at $\lambda \ll 1$ complicates this picture (in principle, repulsion can be achieved even for heavier particles), although force magnitudes in this regime are typically too small to be practically relevant.

5.1. Avoiding effects of outer-flow inertia

The results obtained in this study show that particle motion can be described quantitatively by inertial forcing terms. Often, such computations are complicated by a transition between a viscous-dominated inner flow volume (near the particle) and an inertia-dominated outer volume, necessitating an asymptotic matching of the two limits, such as for the Oseen (Oseen Reference Oseen1910) and Saffman (Saffman Reference Saffman1965) problems. Our formalism, however, only employs an inner-solution expansion and still obtains accurate predictions. This can be rationalized by invoking the analysis of Lovalenti & Brady (Reference Lovalenti and Brady1993), who showed that an outer region is not present when the magnitude of oscillatory inertia in the disturbance flow $\partial \boldsymbol {w}^{(1)}/\partial t$ is much larger than that of the advective term ${\boldsymbol {f}}$, i.e. the characteristic unsteady time scale $\omega ^{-1}$ is shorter than the convective inertial time scale $\nu /(U^*w^{(0)})^2$, where $w^{(0)}$ is the dimensionless velocity scale of the fluid in the particle reference frame. This ensures that vorticity cannot diffuse to the distance of the Oseen wake during the time scale of unsteadiness. For neutrally buoyant particles, $w^{(0)} = O(\alpha )$, while for non-neutrally buoyant particles, $w^{(0)} = O(\hat {\kappa })$. With our definition of ${Re}_p$, this translates into $\lambda \gg \max (\alpha ^2,\kappa ^2) {Re}_p^2$, interpreting $\lambda \propto \omega$ as a measure of unsteadiness. However, in the scenarios of oscillatory microfluidics most relevant to our analysis, ${Re}_p$ itself scales with $\omega$ as well as with the oscillation amplitude $\epsilon$, so that the criterion becomes

(5.1)\begin{equation} \epsilon^2 \lambda \ll \min(\alpha^2/\hat{\kappa}^2,1). \end{equation}

As long as the density contrast between the particle and fluid is small, or $|\hat {\kappa }|\ll 1$, (5.1) is easily satisfied in most experimental situations, and it reverts to the criterion $\epsilon ^2\lambda \ll 1$ established for neutrally buoyant particles (Agarwal et al. Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021). An interesting point to note is that the density-dependent condition $\epsilon ^2 \lambda \ll \alpha ^2/\hat {\kappa }^2$ can be rewritten in the $a_p$-independent form $\epsilon \ll \delta _S/(a_b\hat {\kappa })$. This is because the leading term of the background flow field expansion at the particle position contains no information about the particle length scale.

5.2. Magnitude and practical relevance of inertial effects

In figure 5(a), we illustrated how our formalism, in agreement with DNS, predicts much stronger inertial forces than either MR (which emphasizes viscous effects) or Agarwal et al. (Reference Agarwal, Rallabandi and Hilgenfeldt2018), which treats the background flow as inviscid. For particles typically encountered in microfluidic applications involving biological cells, the density difference tends to be around $5\,\%$, while the size parameter is $\alpha \lesssim 0.2$ and $\lambda \gtrsim 1$. A practically useful metric to quantify the effect of the inertial force acting on the particle by a localized oscillating source is the time needed for radial displacement of a particle diameter. In most particle manipulation strategies, $r_p\gtrsim a_b$ and, upon solving (4.1) with these nominal parameter values, we find that our formalism predicts a time scale of ${\sim }10\ {\rm ms}$ compared with ${\sim }50\ {\rm ms}$ predicted by the inviscid formalism. This translates to much more efficient design strategies for sorting particles based on size or density. The MR formalism predicts a time scale of ${\sim }500\ {\rm ms}$, which is off by more than one order of magnitude and severely underestimates the performance of oscillatory microfluidic set-ups.

For these prototypical cases where particles are close to the interface of the oscillating object ($r_p\gtrsim a_b$), the major contribution to inertial forces is due to $F_{\varGamma \kappa }$, as discussed in Agarwal et al. (Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021), However, since $F_{\varGamma \kappa }$ decays more strongly with the distance from the source than $F_{\sigma \varGamma }$, the density contrast dependent force can easily become comparable in magnitude, resulting in the rich behaviour of attraction and repulsion separated by the critical (and tunable) distance $r_{p_c}$ as described in § 4.4. Thus, the present work suggests new avenues for particle trapping/sorting relying on density contrast; some of these ideas will be explored in future publications.

In a microfluidic set-up, the oscillatory flow is induced around an obstacle, e.g. a cylinder or bubble of radius $a_b$, and as mentioned above, particles in typical applications will approach quite closely to the interface of this obstacle. We have not accounted for effects due to such a nearby boundary in this analysis. In Agarwal et al. (Reference Agarwal, Rallabandi and Hilgenfeldt2018), we demonstrated the existence of a stable fixed-point position when the particle is in very close proximity to the boundary. This stable equilibrium is a consequence of the repulsive lubrication force near the interface balancing the attractive force discussed here. As long as $|\hat {\kappa }|\ll 1$, which is the case in an overwhelming majority of practical applications, the conclusions of Agarwal et al. (Reference Agarwal, Rallabandi and Hilgenfeldt2018) are not affected by the present findings, i.e. a particle attracted to an oscillating obstacle is expected to come to rest at a stable equilibrium distance that is extremely small compared with the interface scale, and typically even compared with the particle scale.