Density-contrast induced inertial forces on particles in oscillatory flows

Abstract Oscillatory flows have become an indispensable tool in microfluidics, inducing inertial effects for displacing and manipulating fluid-borne objects in a reliable, controllable and label-free fashion. However, the quantitative description of such effects has been confined to limit cases and specialized scenarios. Here we develop an analytical formalism yielding the equation of motion of density-mismatched spherical particles in oscillatory background flows, generalizing previous work. Inertial force terms are systematically derived from the geometry of the flow field together with analytically known Stokes number dependences. Supported by independent, first-principles direct numerical simulations, we find that these forces are important even for nearly density-matched objects such as cells or bacteria, enabling their fast displacement and separation. Our formalism thus consistently incorporates particle inertia into the Maxey–Riley equation, and in doing so provides a generalization of Auton's modification to added mass, as well as recovering the description of acoustic radiation forces on particles as a limiting case.


Introduction
One of the most fundamental problems in fluid dynamics that has evaded a general solution is describing the unsteady motion of particles immersed in a prescribed background flow.Most analytical attempts work under the severe assumption of unsteady Stokes flows, for which symmetry-breaking inertial effects are neglected (see Michaelides (1997) for a brief overview).Following the pioneering contributions of Basset (1888), Boussinesq (1885) and Oseen (1927) for uniform background flows, resulting in what is now termed the BBO equation, effects of flow non-uniformity were taken into account by Gatignol (1983) and Maxey & Riley (1983) (MR).Because of its comprehensive, rigorous and systematic character, MR has been widely used to describe hydrodynamic forces on particles over the past 40 years, although it still operates strictly in the limit of vanishing inertial effects.
The Maxey-Riley equation for spherical particles (MR equation) assumes the validity of the unsteady Stokes assumption, which implies that (i) the particle Reynolds number based on a typical difference velocity between particle speed and background flow must be small, and (ii) the background flow gradients must be small compared with viscous momentum diffusion.These assumptions do constrain the applicability of MR in a number of situations.One of the most glaring shortcomings was pointed out by Leal (1992), and concerns the incompatibility of MR with the experimentally observed phenomenon of lateral migration of particles due to lift forces caused by inertial effects.Subsequent work aimed at the development of equations valid at finite particle Reynolds numbers has yielded specialized results, for example for steady flow (Ho & Leal 1974;Martel & Toner 2014;Hood, Lee & Roper 2015) or for forces occurring in acoustic fields (Baudoin & Thomas 2020;Rufo et al. 2022).
The advent of oscillatory microfluidics (Lutz, Chen & Schwartz 2003;Marmottant & Hilgenfeldt 2003;Thameem, Rallabandi & Hilgenfeldt 2017;Mutlu, Edd & Toner 2018;Zhang et al. 2020Zhang et al. , 2021aZhang et al. ,b, 2023) ) has since introduced the use of much stronger particle inertia effects, enabling fast and high-throughput particle manipulation.Yet again, quantitative modelling and prediction of such effects has been largely lacking, with experimental results often explained qualitatively, and/or by appealing to analogies with theories that are not obviously applicable, such as particle forces in acoustofluidics (Gor'kov 1962;Friend & Yeo 2011;Devendran, Gralinski & Neild 2014;Chen et al. 2016;Nadal & Lauga 2016;Collins et al. 2019;Wu et al. 2019).Given the versatility and richness of microfluidic flows, what is needed is a fundamental understanding of inertial hydrodynamic forces acting on particles immersed in a general unsteady background flow, that is, a true generalization of MR.
In a first step towards such a generalization, Agarwal et al. (2021) rigorously described inertial forces on density-matched particles.Whereas MR does not predict any net force on neutrally buoyant particles immersed in unsteady fluid flows, Agarwal et al. (2021) showed that such a force can be very significant and is often dominant in oscillatory microfluidics.In the present paper, we augment that formalism to include finite density contrast between particle and fluid (a relevant scenario in microfluidics), thus completing the consistent generalization of MR.In our theory, density-contrast dependent contributions to inertial forces specialize to the well-known Auton correction (Auton, Hunt & Prud'Homme 1988) in the potential flow limit, but continue to play a significant role in the presence of unsteady viscous effects.In a different specific limit our framework establishes a quantitative connection with acoustofluidic formulae for radiation forces on particles, which, as has been remarked above, have up to now been used in an unsystematic way in oscillatory microfluidics.
The organization of this paper is as follows.In § 2 we describe the general theoretical formalism for inertial forces and their evaluation for oscillatory flows.In § 3, we develop an explicit time-averaged equation of motion for spherical particles, and in § 4 we rigorously compare its predictions with direct numerical simulations (DNS) as well as with existing theories in specialized limits.Section 5 discusses the validity and importance of the present approach in practical situations, while § 6 draws conclusions.985 A33-2

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 ρ 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 ν of the fluid define the particle Reynolds number Re p = a p U * /ν.The particle reaction to the flow importantly also depends on the Stokes number, which we define as λ = a 2 p ω/(3ν).The unsteady time scale of the flow is written as 1/ω, anticipating the oscillatory case, where we can alternatively use the Stokes boundary layer scale δ = (2ν/ω) 1/2 to write λ = 2a 2 p /3δ 2 .We follow MR in decomposing the flow around the particle into the given background flow 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 λ 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 (1993), 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 λ values.In oscillatory microfluidics, where we typically have λ ∼ 1-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 λ 1, and we will show that particle motion in those cases is also quantitatively covered by our formalism.The case λ 1 is usually not practically relevant, as the resulting inertial forces on particles become very small.

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 U p : where subscripts denote orders of Re p .Note that the decomposition of F (0) is exact (there are no terms of higher order in Re p ; cf.MR), while we truncate the expansion of F (1) at first order.Expressions for 0 and part of 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 (i) dS with i = 0, 1, where we use the Stokes drag scale F S /6π = νρ f a p U * , and the integral is over the particle surface with its outward normal 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 w (0) = u − u p for the undisturbed background flow and w (1)  for the disturbance flow.The dimensionless fluid stress tensors are thus written σ (i) The Navier-Stokes equations in the particle frame of reference can be decomposed into background and disturbance components exactly (without approximations), 0) • ∇w (0) , ∇ • w (0) = 0, (2.2a) ∇ 2 w (1) − ∇p (1) = 3λ ∂w (1)  ∂t + Re p w (0) • ∇w (1) + w (1) • ∇w (0) + w (1) • ∇w (1) , (2.2c) w (1) = u p − U, on r = 1 and w (1) = 0, as r → ∞.
(2.2e) 2.3.Forces from background flow Both the O(1) and O(Re p ) components of F (0) in (2.1) can be evaluated directly using the divergence theorem and the above Navier-Stokes equations valid for the background flow w (0) .In laboratory coordinates (see MR; Rallabandi 2021) they read To make further progress, we need to evaluate forces due to successive orders of w (1) , which are ultimately also derived from the given background field u.To render our solution strategy analytically tractable, we expand u around the leading-order particle position r p 0  Agarwal et al. (2021) and Rallabandi (2021), showing that an O(Re p ) contribution from F (0) 1 had been missed in MR, while being, in fact, of the same order as other terms in the original MR equation.

Disturbance flow: zeroth order
The Navier-Stokes equations for the disturbance flow at O(Re 0 p ) read w (1) 0 = u p 0 − u, on r = 1 and w (1) 0 = 0, as r → ∞. (2.5b) 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 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 1959;Pozrikidis et al. 1992) where u s = u p 0 − u| r p 0 is the slip velocity and M D,Q,O (r, λ) are mobility tensors with known spatial dependence.For oscillatory flows, their dependence on the Stokes number λ 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 ∇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 F (1) become important.
With w (1) 0 explicitly known, the equations at O(Re p ) read as the leading-order nonlinear forcing.
In order to compute the force F (1) 1 , we do not solve for the flow field 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 u , which is here chosen to be a dipolar unsteady Stokes flow around the particle with arbitrary directionality e; see Appendix B for a detailed derivation.We obtain for the magnitude of the force along e, e • F (1)  1 = where the hat denotes the Laplace transform and L −1 is the inverse Laplace transform.When applied at O(1), this reciprocal-theorem strategy similarly yields e • F (1) 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 u p , we make the substitution w . 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 2 p ).The net force on the particle due to its disturbance flow then reads where we have also replaced w (0) 0 by w (0) in f 0 , resulting in an error that is again O(Re 2 p ).Only certain products in 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 Here, F (1) σ Γ , F (1) Γ κ are the inertial force contributions obtained by successive contractions of adjacent tensors involving u s (index σ ), E (index Γ ), G (index κ) 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 ω in the following, without loss of generality.
When the particle is neutrally buoyant, the first term F (1) σ Γ vanishes so that the leading term is F for density-matched particles.This term (involving the product E : G) has no analogue in the previous literature and for completeness, we reproduce it as follows for harmonic oscillatory flows U: (2.12) The λ-dependent dimensionless function F (1)  1 results from the volume integration, which also yields m f , the mass of fluid displaced by the particle, via (4πa In the next section, we follow a similar strategy for non-neutrally buoyant particles. 2.6.Disturbance flow: evaluation of F (1)   σ Γ Non-neutrally buoyant particles have a slip velocity and thus a non-zero F (1) σ Γ , involving the product u s • 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. u s (r The corresponding force is written as where the G 1 and G 2 terms are explicit outcomes of the volume integration in (2.11) and capture the λ-dependence of the in-and out-of-phase contributions, respectively.For fast oscillatory background flows, we can replace the in-phase component with u s • E and the out-of-phase component with ∂ t u s • E (see Appendix C for details), resulting in (2.14) While the exact, lengthy expressions for the universal functions 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 λ expansions of G 1 (analogous to the function F in Agarwal et al. (2021)).Taylor expansion in both the viscously dominated limit (λ → 0) and the inviscid limit (λ → ∞) obtains from which the following uniformly valid result is constructed: 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 λ, with a maximum error of ∼6 %.(c) Plot of the out-of-phase force function G 2 (red) together with its viscous (green) and inviscid (blue) limit expressions.
A Taylor expansion of the out-of-phase term G 2 , in the viscous and inviscid limits, respectively, results in Both of the above expansions have a O(1/ √ λ) 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 λ 1.

Equation of motion for a particle immersed in an oscillatory flow
We now collect all the force contributions from (2.3a,b), (2.10c), (2.12), (2.14), and combine them with the results of MR and Agarwal et al. (2021).We use dimensional variables for easier physical interpretation.The following is the equation of motion for the velocity U p of a rigid spherical particle immersed in an oscillatory background flow field U, taking into account all force terms up to O(Re p ): 1 .
(3.1d) Here, we have dropped the contraction with e in (2.12) and (2.14) to derive F (1) 0 , since the direction e is arbitrary (cf. the equivalent argument in MR).Equation (3.1b) includes the background flow force term missing from MR mentioned in § 2.3, proportional to F (0) 1 = 1/5.As we are modelling the oscillatory background flow as potential (cf.§ 2.6), Faxén terms proportional to ∇ 2 U are absent.For the same reason of irrotational background flow, particle rotation is neglected.Note that the scales of all the inviscid and inertial force terms use m f , while the viscous force terms contain ν explicitly.In the following, we point out that (3.1), while containing new physics, encompasses a number of earlier results as special cases, clarifying connections between them.

Generalized Auton correction
We first comment on the limiting case of the well-known correction to MR due to Auton (Auton et al. 1988).The equation of motion derived by MR deviated from previous versions in the form of the convective term in (3.1b), using m f (U • ∇U) instead of m f (U p • ∇U) -the values of these two derivatives can differ substantially when the Reynolds number is not small.Similarly, Auton et al. (1988) showed that in the limit of potential flows, the added mass term should read 1 2 m f (dU p /dt − DU/Dt) instead of 1 2 m f (dU p /dt − dU/dt).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 G 1 , modifies the added mass term in (3.1c) and reproduces the Auton correction (Auton et al. 1988) in the inviscid, potential flow limit (λ → ∞, Re p 1), modifying dU/dt to DU/Dt, or explicitly, where we use the simple two-term approximation (2.16) for G 1 .Thus, instead of heuristically modifying the added mass term, our approach rigorously derives its dependence on λ.Note that in most practically relevant oscillatory microfluidic flows, the value of λ is O(1-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/ √ λ scaling.We note that the second term in (3.1d), involving 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 λ ∼ O(1), while it is subdominant for both small and large λ.

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 985 A33-9 https://doi.org/10.1017/jfm.2024.303Published online by Cambridge University Press 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 a b and angular frequency ω, we non-dimensionalize our equation with a b , a b ω and 1/ω as characteristic length, velocity and time scales, respectively.Equation (3.1) then reads where κ = 2/3(ρ p /ρ f − 1) is a dimensionless measure of density difference, α = a p /a b is the relative particle size and dr p /dt = u p .As in Agarwal et al. (2021), we write F = 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 r p 0 on the slow time scale T = 2 t, Agarwal et al. (2021) and where c = 1 + √ 3λ/2 and d = 1 + √ 3/(2λ) 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 G(λ), where F R = ( κλ/( κ + 1)) u • ∇u is a time-averaged force formally identical to the acoustic radiation force induced by an incident sound field with velocity field u (Bruus 2012).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 1996).As the acoustic formalism is based on the assumption of inviscid flow, G(λ) 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 G(λ → ∞) = 1, recovering the inviscid case, while the viscous limit depends on the density contrast, G(λ → 0) = −(1 + κ).
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.

Validation with DNS
We have shown that the present analytical formalism generalizes previous attempts at predicting the behaviour of particles in oscillatory flows.It is crucial to confirm the quantitative accuracy of our model.To this end, we compare our analytical predictions with independent, first-principles DNS of the full Navier-Stokes equations, previously validated in a variety of streaming flow scenarios (see Gazzola et al. 2011 In order to make quantitative comparisons, we restrict the background flow field to a spherical, oscillatory monopole.These flows are typically generated near volumetrically excited bubbles and have been shown to actuate inertial forces on particles in oscillatory microfluidics (Rogers & Neild 2011;Chen & Lee 2014;Zhang et al. 2021a), showcasing their practical utility.This specialization offers an ideal framework for validating our analytical formalism, as this radially symmetric flow by itself does not induce viscous streaming, enabling us to neatly isolate the effect of inertial forces.
Accordingly, we insert u(r, t) = (1/r 2 )e it e r into (3.4) to obtain the following time-averaged equation of motion (we drop the subscript 0): where −( κλ/r 5 p ( κ + 1)) = F R and r p is in units of the radius of the oscillating source.This simple ordinary differential equation (ODE) provides clear predictions for the particle fate that can be compared with results from DNS.Two terms in (4.1) determine the direction of particle motion: the second term involving F is always negative (Agarwal et al. 2021), representing attraction towards the source while the sign of the first term changes with κ and G. Therefore, the magnitude and sign of the net force depend on several parameters, including λ, κ, and also on r p , as the first term dominates the second at large distances.
Note that this set-up is specifically constructed such that all effects on the right-hand side of (4.1) are due to inertia.Thus, the comparison between analytical predictions and DNS solutions provides a direct and accurate test of particle-inertial effects in oscillatory microfluidics.We will focus on the key quantities of practical interest: the particle trajectories, velocities and forces.

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 ω and velocity amplitude ω (the nominal source size a b is normalized to 1).We choose = 0.01, a p = 0.05, ω = 16π throughout, and r p0 = 2 unless otherwise stated.The fluid viscosity is determined from the corresponding values of λ in each simulation.The top half of the panel in figure 3 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 ρ p /ρ f ≈ 0.95.Accordingly, the particle motion in the simulation (figure 3b) reverses direction: particles lighter than ≈ 0.95ρ f are repelled over time, while those of greater density are attracted towards the monopole (which includes the density-matched case).

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(a-d) depict examples of such averaged radial dynamics for different density ratios and different Stokes numbers λ, 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 G(λ) 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. (2018)  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 λ values, figure 4(e) condenses all results by extracting a G(λ) value from best-fitting (4.1) to the numerically simulated particle trajectories (see Appendix F for details), given the previously established accuracy of the F(λ) function (Agarwal et al. 2021).Both for heavier (ρ p /ρ f = 1.1, red) and lighter particles (ρ p /ρ 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 λ = 20, there are significant deviations of G(λ) 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 ∼24-48 core hours on a single node on the Expanse supercomputer (see Appendix E), while the theory ODE is trivial to solve.

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 2012;Friend & Yeo 2011).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 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 G(λ).That such a λ-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.
Predictions using the MR equation fail to correctly reproduce the inviscid limit (λ → ∞) 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 (2012) instead misses the opposite viscous limit (λ → 0), as it ignores viscosity completely.In previous work (Agarwal et al. 2018), 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 (1994) in both the viscously dominated (λ → 0) and the inviscid limits (λ → ∞), while quantitative discrepancies remain in the intermediate λ regime, where the G(λ) of Doinikov (1994) is larger than that of Agarwal et al. (2018).
The theory of the present work agrees with the previously established viscous and inviscid limits, and makes new predictions in the intermediate λ range, with values between those of Doinikov (1994) and Agarwal et al. (2018).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 ≈5 %-10 % across the simulation range 1 ≤ λ ≤ 20.This range is restricted by computational cost on both ends, as the boundary layer scale becomes large for λ 1 and the size of the computational domain must be increased, while for λ 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.

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 dr p /dT = 0 in (4.1) obtains the critical radial position (in units of the particle radius a p ) as In most practically relevant situations, λ O(1), and thus G > 0 (cf.figure 5a).A real r p c then exists if the particle is lighter than the surrounding medium ( κ < 0).Such an equilibrium position is necessarily unstable, as the repulsive term in (4.1) decays more slowly.Thus, for light particles and λ 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. 2016;Zhang et al. 2021a,b).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 λ and κ, 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 (ρ p /ρ f < 1) are repelled above a threshold Stokes number.Comparison with DNS data (circles in figure 5b) confirms these predictions.The sign change of G(λ) at λ 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.

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 1910) and Saffman (Saffman 1965) 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 (1993), who showed that an outer region is not present when the magnitude of oscillatory inertia in the disturbance flow ∂w (1) /∂t is much larger than that of the advective term f , i.e. the characteristic unsteady time scale ω −1 is shorter than the convective inertial time scale ν/(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(α), while for non-neutrally buoyant particles, w (0) = O( κ).With our definition of Re p , this translates into λ max(α 2 , κ 2 )Re 2 p , interpreting λ ∝ ω as a measure of unsteadiness.However, in the scenarios of oscillatory microfluidics most relevant to our analysis, Re p itself scales with ω as well as with the oscillation amplitude , so that the criterion becomes 2 λ min(α 2 / κ2 , 1). (5.1) As long as the density contrast between the particle and fluid is small, or | κ| 1, (5.1) is easily satisfied in most experimental situations, and it reverts to the criterion 2 λ 1 established for neutrally buoyant particles (Agarwal et al. 2021).An interesting point to note is that the density-dependent condition 2 λ α 2 / κ2 can be rewritten in the a p -independent form δ S /(a b κ).This is because the leading term of the background flow field expansion at the particle position contains no information about the particle length scale.

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. (2018), 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 α 0.2 and λ 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 a b and, upon solving (4.1) with these nominal parameter values, we find that our formalism predicts a time scale of ∼10 ms compared with ∼50 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 ∼500 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 a b ), the major contribution to inertial forces is due to F Γ κ , as discussed in Agarwal et al. (2021), However, since F Γ κ decays more strongly with the distance from the source than F σ Γ , 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. (2018), 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 | κ| 1, which is the case in an overwhelming majority of practical applications, the conclusions of Agarwal et al. (2018) 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.

Conclusions
We have developed a rigorous formalism to accurately describe the motion of particles in general, fast oscillatory flows.The present work systematically accounts for finite inertial forces in viscous flows that result from the interaction between the density-contrast dependent slip velocity and flow gradients.Confirmed by DNS, these forces are found to be important and often far larger than the density-contrast dependent effects present in the original formalism of MR.Our theory allows for quantitative predictions of the sign and magnitude of forces exerted on particles in many customary microfluidic settings, in particular for nearly density matched cell-sized particles -the most relevant case in medicine and health contexts.The theory encompasses special cases such as Auton's correction and acoustic radiation forces in the inviscid limit, and provides their quantitative generalization in the presence of viscous effects.other tensors are determined analogously.Using components in spherical coordinates, they read where and β = −ia 2 p /(ν/ω) = √ −3iλ is the complex oscillatory boundary layer thickness.
We emphasize that these expressions are the same for arbitrary axisymmetric oscillatory u.Accordingly, only the expansion coefficients u s , E and G contain information about the particular flow.
It is understood everywhere that physical quantities are obtained by taking real parts of these complex functions.

Appendix B. Reciprocal theorem and test flow
In order to compute the force F (1) 1 , we do not solve for the flow field w A key simplification due to oscillatory flows is that the Laplace transforms are explicitly computed.The symmetry relation employs a known test flow (denoted by primed quantities such as u ) in a chosen direction e, around an oscillating sphere such that it satisfies the following unsteady Stokes equation: where the unit vector e is chosen to coincide with the direction in which the force on the particle is desired.The solution to this problem is of the same form as (2.6), but with only In order to get explicit results for the non-trivial integration factors G 1 and G 2 , we insert f 0 (with explicitly known mobility tensors M D,Q,O ) into (2.12).Since F (1) σ Γ involves products of oscillatory terms, there are higher-order force harmonics with zero net effect on the particle dynamics which we will average out in the following to simplify the integration evaluations.
We first decompose the slip velocity into its in-phase and out-of-phase components, i.e. u s (r p , t) = u I s (r p , t) + u O s (r p , t), as noted in the main text, and time-average (2.12) over a period of oscillation to remove higher-order harmonic terms.We then perform the volume integration to obtain an explicit but rather lengthy expression that can be symbolically written as where G 1 and G 2 are explicit outcomes of the volume integration.Exploiting the orthogonality of trigonometric functions and the fact that, for fast oscillatory background flows, E is purely in-phase, we rewrite the in-phase component as u s • E and the out-of-phase component as ∂ t u s • E , where angled brackets denote time-averaging, so that Finally, we drop the time-averaging operation, producing an error in the higher-frequency force harmonics that has zero effect on net particle motion, resulting in (2.14) in the main text.
The explicit expression for the in-phase inertial force component for oscillatory flows reads The contribution due to the O( 2 ) nonlinear forcing term ∂ τ (r p 1 • ∇u osc ) is identically zero for oscillatory flows, after time-averaging.Additionally, the effect on the steady flow component is higher-order in and is, therefore, neglected.Thus, the main contributions due to the history integral appear as subdominant corrections to the Stokes drag and added mass terms, given by (D1), at O( ).
We now proceed with the formal separation of time scales of (3.3).At O(1), This equation is trivially satisfied if r p 0 = r p 0 (T); thus, the leading-order particle position r p 0 depends only on the slow-time T. At O( ), we obtain the following after explicitly evaluating the history integral:

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

Figure 3 .
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 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).

Figure 4 .
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 G and F agree with DNS (magenta) for the entire range of λ and density contrast values (all results are for r p0 = 2).Two exemplary density contrast and λ combinations are displayed, ρ p /ρ f = 1.1 (κ = 0.067) and ρ p /ρ f = 0.9 (κ = −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 (2018) (light blue) has similar, though quantitatively less severe, shortcomings.(e) Best-fits of G(λ) to (4.1) are extracted from DNS and show excellent agreement with the full theory (3.5), for both heavier (ρ p /ρ f = 1.1, red) and lighter (ρ p /ρ f = 0.9, teal) particles.

Figure 5 .
Figure 5. (a) Stokes number dependence of the overall dimensionless inertial force magnitude 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 ρ p /ρ f = 1.1 (κ = 0.067), = 0.01.The DNS values are best fits of G given the full expression for F in (4.1).The present work (red line) is in excellent agreement with all DNS data, while both the Agarwal et al. (2018) (light blue) and MR formalisms (green) significantly underestimate the forces.(b) Contour plots for steady particle velocity at r p = 2 with varying λ and ρ p /ρ 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.