Inertial focusing of spherical particles in curved microfluidic ducts at moderate Dean numbers

We examine the effect of Dean number on the inertial focusing of spherical particles suspended in flow through curved microfluidic ducts. Previous modelling of particle migration in curved ducts assumed the flow rate was small enough that a leading order approximation of the background flow with respect to the Dean number produces a reasonable model. Herein, we extend our model to situations having a moderate Dean number (in the microfluidics context) while the particle Reynolds number remains small. This extension allows us to capture changes in the background flow that occur with increasing flow rate, namely a shift in local extrema towards the outside wall. The change in the axial velocity profile of the background flow has an effect on the inertial lift force, while the change in the cross-sectional components directly affects the secondary flow drag. In keeping the particle Reynolds number small we approximate the inertial lift force in a similar manner to previous studies while capturing subtle effects do to the modified background flow profile. Capturing and understanding these effects is an important step towards accurately modelling inertial migration across a wide range of practical applications. Our results reveal how the changing background flow profile modifies the inertial focusing of particles. We illustrate enhanced lateral separation of particles by size in a number of scenarios and find that focusing times can be roughly separated into two regimes. These results suggest our model might aid with parameter choices for separation of particles by size.


Introduction
Inertial focusing of particles suspended in flow through curved and spiral microfluidic ducts has been studied extensively in the experimental literature, particularly in relation to its application to size based particle/cell separation (Seo et al. 2007;Bhagat et al. 2008;Martel & Toner 2012;Warkiani et al. 2014Warkiani et al. , 2016Rafeie et al. 2019). Much is known about the nature of the inertial lift force in a variety of situations involving a particle suspended in flow between two plane parallel walls (Saffman 1965;Ho & Leal 1974;Schonberg & Hinch 1989;Hogg 1994;Asmolov 1999). However, the nature of the inertial lift force is very different for a fully enclosed duct making these studies of limited use in understanding practical applications. Sufficiently small particles suspended in flow through straight ducts with square cross-section are known to focus at one of four equilibria located a small distance from the centre of each side wall Hood et al. 2015). This is also the case for straight rectangular ducts, although stable equilibria near the shorter side walls attract relatively few particles and disappear entirely for larger particles (Martel & Toner 2013;Hood 2016). In curved ducts the migration of particles becomes complicated due to the secondary vortices that are generated as part of the Dean flow through the duct. The interaction of the drag force from these vortices with the inertial lift force leads to a wide variety of particle migration dynamics (Gossett & Di Carlo 2009;Martel & Toner 2014;Ha et al. 2022;Valani et al. 2022).
Our previous study conducted a detailed examination of the migration of neutrally buoyant spherical particles suspended in a sufficiently slow flow through curved rectangular ducts (Harding et al. 2019). An accurate model of particle migration was developed by coupling the particle motion to a Navier-Stokes model of the fluid flow. By using a carefully chosen reference frame, and applying a suitable non-dimensionalisation and perturbation expansion of both the background and disturbance flows, the individual forces primarily responsible for driving particle migration were separated and then estimated via numerical simulation. These forces were then re-assembled into a system of ordinary differential equations to model particle trajectories. It was found particles migrated towards stable equilibria whose horizontal location approximately collapsed onto a single curve when plotted against the parameter κ = 4 /(4a 3 R), with being the duct height, a being the particle radius, and R being the bend radius. It was later shown how non-neutrally buoyant particles could be modelled by adding suitable perturbations to the neutrally buoyant model (Harding & Stokes 2020).
A key part of the prior modelling was an assumption that the flow rate is small enough that both Re p = Re(a/ ) 2 and K = Re 2 are small, where = /(2R) and Re = (ρ/µ)U m ( /2) is the channel Reynolds number with ρ being the fluid density, µ the fluid viscosity and U m the maximum velocity of the background flow down the main axis. This allowed us to take the leading order contribution of each perturbation expansion as a reasonable approximation. In a typical practical setting, these assumptions only hold when Re O(10). In contrast, most microfluidics experiments in the literature correspond to Re = O(100). While we expect our existing model may still have qualitative value at these flow rates, they are of less use quantitatively.
Substantially higher flow rates, e.g. Re O(1000), are of limited practical interest for a couple of reasons. The first is that the increasing strength of the Dean flow eventually inhibits the ability of particles to focus, as seen in the decreasing sharpness factor with increasing flow rate in the experimental results of Rafeie et al. (2019). Second, there is a critical Dean number, depending on the specific cross-section, above which the secondary component of the background flow exhibits multiple vortex pairs (Winters 1987) and this seems generally undesirable for most applications.
In this paper we incorporate additional terms from the perturbation expansion of the background flow into the model. This has the effect of increasing the values of K for which our model has quantitative value. Since we continue to use a leading order approximation of the disturbance flow, this model doesn't expand the applicability in cases where the magnitude of Re p is the limiting factor (e.g. when the particle is relatively large). Nonetheless, we feel this provides valuable insights and is an important step towards producing an accurate model applicable to a wide range of physical set ups. This work also illustrates how the symmetry associated with a curved duct leads to a decoupling of axial and secondary parts in the disturbance flow at leading order which each generate inertial lift in their counterpart at first order.
The paper is organised as follows. Section 2 describes the general setup of the problem and briefly summarises the modelling of forces driving particle migration as developed in Harding et al. (2019). We also introduce some notation to support the remainder of the text and remark on the applicability of the Lorentz reciprocal theorem to torque  Figure 1: Curved duct with rectangular cross-section containing a spherical particle located at x p = x(θ p , r p , z p ). The enlarged view of the cross-section containing the particle illustrates the origin of the local r, z coordinates at the centre of the duct. The bend radius R is with respect to the centre-line of the duct. Note that we do not consider the flow near the inlet/outlet. Adapted from Harding et al. (2019).
calculations. Section 3 describes our improved approximation of the background flow which utilises multiple terms from a perturbation expansion with respect to the Dean number K. Section 4 describes how the new background flow approximation is incorporated into the inertial lift calculation and ultimately leads to a system of first order differential equations which describe particle migration. Section 5 reports a range of findings obtained from the new model: firstly, we examine how the horizontal location of stable equilibria are perturbed by increasing K; secondly, we examine how the parameters , K, κ, α(= 2a/ ) influence how long it takes particles to focus; thirdly, we examine how increasing K influences previously observed trends in the horizontal location of stable equilibria with respect to −1 and κ. Section 6 summarises our findings and remarks on avenues for future exploration.

Problem setup and theoretical background
Our curved duct setup remains identical to that in Harding et al. (2019) and is depicted in Figure 1. The (stationary) lab reference frame is x = xi + yj + zk with the duct bending around the z-axis. The duct domain is most readily described using the cylindrical coordinate system (r, θ, z) for which where R is the bend radius of the duct measured from the origin (of the lab frame) to the centre of the cross-section. The cross-section itself is described by (r, z) ∈ C (with origin (r, z) = (0, 0) in the centre of the cross-section). The duct interior is then described by D = {x(θ, r, z) | (r, z) ∈ C}. While our approach can be applied to any desired crosssection C, this study is concerned with rectangular cross-sections having width W and height H, thus We take = min{W, H} to be the characteristic length scale of the duct. Of principle interest will be ducts with W H, and thus = H, as these are most common in the experimental literature.
Steady pressure driven flow through the duct (in the absence of any particles) is (c) a spherical particle and the primary cross-sectional forces which drive its migration. Here F S is the drag from the secondary component of the background flow, and F L is the inertial lift force. The magnitude and direction of each vector are for illustration only. Gravitational and centrifugal/centripetal forces are omitted. The background flow is shown to be skewed towards the outside wall of the curved duct (here on the right), as is expected at moderate Dean numbers. Adapted from Harding & Stokes (2020).
referred to as the background flow. The fluid is assumed to be incompressible with constant density ρ and viscosity µ. The pressure and velocity fields are denotedp andū, respectively, and are modelled using the Navier-Stokes equations. We take U m to be a characteristic velocity of this flow, approximately describing the maximum axial velocity. The channel/duct Reynolds number is then Re := (ρ/µ)U m ( /2). Additionally, letting = /(2R) denote the relative curvature, we define the Dean number as K = Re 2 after Dean (1927); Dean & Hurst (1959) who studied the secondary flow which develops within curved duct flow. The interaction of drag from this secondary flow with the inertial lift force, primarily due to the axial flow, produces a variety of particle migration dynamics that are exploited in a wide range of microfluidics applications. The nature of the background flow and its approximation for the purposes of this study will be discussed further in Section 3. Figure 2 depicts the axial (a) and secondary (b) components of the background flow, and (c) depicts the competing secondary drag and inertial lift forces on a particle. Now, consider a single particle suspended in the flow. Let P := {x : |x−x p | < a} denote a spherical particle with constant density ρ p , radius a and centred at x p (t) (such that P ⊂ D). We let r p , θ p , z p denote its cylindrical coordinates, i.e. x p = x(r p , θ p , z p ), where r p , θ p , z p are each functions of t. The particle travels with a velocity u p (t) := dx p /dt and spins freely about its centre with angular velocity Ω p (t). Thus, a point x within P has instantaneous velocity u p + Ω p ×(x − x p ). The fluid domain is denoted F := D\P and its boundary ∂F consists of the duct walls ∂D and the particle surface ∂P. The fluid flow is now described by the pressure and velocity fields p, u, respectively, which are also modelled by the Navier-Stokes equations, specifically where σ(p, u) is the stress tensor σ(p, u) := −pI + µ (∇u + ∇u ) . (2.4) The gravitational body force g is only important for non-neutrally buoyant particles (ρ p = ρ). In Harding & Stokes (2020) it was demonstrated that with g = −gk the influence of non-neutral buoyancy can be separated from the migration model and subsequently treated as an additional perturbation of the force experienced by a neutrally buoyant particle. This remains true in this work so we simplify the development by considering a neutrally buoyant particle (i.e. ρ p = ρ) and dropping the gravitational acceleration from (2.3a).
The motion of the suspended neutrally buoyant particle is driven solely by the hydrodynamic force and torque exerted on it. Specifically describe the force and torque, respectively, where n is taken to be the normal with respect to the fluid domain (i.e. pointing in towards the particle centre). The development of the migration model in Harding & Stokes (2020) may be summarised as the following sequence of steps: 1. Introduce a rotating reference frame in which the particle's angular coordinate is constant. This frame rotates with angular velocity Θ = Θk where Θ := ∂θ p /∂t. Coordinates in this rotating frame are mapped to the 'stationary frame' via x (r , θ , z ) = (R + r ) cos(θ )i + (R + r ) sin(θ )j + z k (2.6a) = (R + r ) cos(θ + θ p )i + (R + r ) sin(θ + θ p )j + z k , (2.6b) where i := cos(θ p )i + sin(θ p )j and j := − sin(θ p )i + cos(θ p )j. It follows that in the rotating frame x p = x (r p , 0, z p ), u = u − Θ×x,ū =ū − Θ×x, and so on.
2. Assume the cross-sectional components of u p are sufficiently small that the flow in the rotating frame is approximately stationary. As we are concerned with crosssectional migration after any initial acceleration to terminal velocity in the axial direction, then acceleration effects (e.g. added mass and Basset force) may be neglected.
3. Introduce the disturbance flow q , v which satisfies 4. Non-dimensionalise using the velocity scale U s = (α/2)U m and length scale a, where α := 2a/ . Most other scales may be derived from these (as per usual for a viscous flow). The resulting Reynolds number Re p = (ρ/µ)U s a is often called the particle Reynolds number. The force on the particle is non-dimensionalised via the scale Re p µU s a. Hats are used to describe non-dimensionalised variables, for example v = U sv .

5.
Apply a perturbation expansion to the disturbance flow with respect to Re p , that iŝ The force on the particle is expanded in a similar fashion, although the leading term has order Re −1 p and is denoted F −1 to reflect this, that iŝ and similarly for the torque. Observe we drop both the hat and prime from variables upon applying the perturbation expansion to each ofq ,v ,F ,T . Following these steps one arrives at: the equations governing q 1 , v 1 ; and lastly, the force and torque terms of principle interest We refer the interested reader to Harding et al. (2019) for a complete derivation. Note that in (2.10d) we retain aû p contribution whereas in Harding et al. (2019) this was taken to be Θ × x p (i.e. the θ component only) with drag coefficients associated with the cross-sectional components introduced later. Additionally, here we include particle velocity and spin contributions in (2.11d). Notice these contributions in both (2.10d) and (2.11d) are expressed in terms ofû p ,Ω p as viewed in the stationary reference frame. These minor changes allow us to separate the leading and first order contributions to the particle motion, and account for inertial lift contributions arising from the secondary component of the background flow.
It is reasonably well-known that the last term of (2.12c) can be computed without needing to explicitly calculate q 1 , v 1 via an application of the Lorenz reciprocal theorem.
Thus we only need to consider the solution of (2.10). Before moving on to describe our improved background flow approximation, we introduce some additional notation and make a remark about the use of the Lorentz reciprocal theorem.

Stokes' flow operators
We introduce the operators Q(f ), V(f ) which map a continuous vector field f defined on ∂P to a pressure and velocity field, respectively, which satisfies the Stokes' equations Additionally, the corresponding hydrodynamic force and torque on the particle are denoted by When f is a constant unit vector then M(f ) and N (f ) may be interpreted as force and torque coefficients, respectively, in the direction of f . Note that all four of these operators implicitly depend on the location of the particle within the cross-section since the fluid domainF and particle boundary ∂P depend onr p ,ẑ p . Additionally, there is an implicit dependence on , as this relates to the bend radius of the duct, and α, as this relates to the size of the non-dimensionalised fluid domain. The implicit dependence on is expected to be very weak for 1/100. Observe that the operators Q, V, M, N are linear and therefore, for example, we may write This linearity will be further exploited in Section 4. Appendix B describes when various components of these coefficients are zero due to the symmetry of our domain (in the θ direction about θ p ).

A note on reciprocal theorems
A variant of the Lorentz reciprocal theorem can be applied to show that where I 0 is the right side of (2.11a), that is This application of the Lorentz reciprocal theorem to the calculation of the inertial lift force is reasonably well-known. Perhaps less well-known is that it can also be applied to the calculation of the torque terms. Specifically, the third term of T 0 in (2.12d) can be calculated via Notice that the M(·) and N (·) terms in (2.16) and (2.18), respectively, encapsulate the contributions from the the boundary conditions (2.11d), while the remaining volume integrals are the result of the reciprocal theorem applied to capture the contribution of I 0 . A detailed proof of the reciprocal theorems that give rise to these volume integrals is provided in Appendix A.

Improved approximation of the background flow
The background flowp,ū is a steady solution of the Navier-Stokes equations: The resulting velocity field may be decomposed into its axial componentū a and secondary componentū s . A perturbation expansion may be applied to each component with respect to K = Re 2 which converges provided K 212 (Harding 2019d ). Specificallȳ whereū θ,i ,ū r,i ,ū z,i are the (dimensionless) components ofū in the θ, r, z directions, respectively, and are each independent of θ.
Our previous model used only the leading order termsū θ,0 ,ū r,0 ,ū z,0 to model inertial migration in curved ducts based on an assumption that K is suitably small (Harding et al. 2019). In this study we extend the use of (3.2) to accurately model particle migration for values of K up to O(100). Specifically, we will use the three leading terms i = 0, 1, 2 to construct a model of particle migration which is quadratic in K. Upon nondimensionalising with respect to the shear velocity scale (U s = (α/2)U m = (a/ )U m ) our model may be expressed aŝ whereū a,i :=ū θ,i e θ andū s,i :=ū r,i e r +ū z,i e z for each i. Notice in (3.3) we have used the fact that 2α −1 Re = κRe p (recalling κ = 4 /(4a 3 R)).
We illustrate the accuracy of our improved model of the background flow in Figure 3.  We use the streamfunctionΦ, for whichū r = −(1 + r) −1 ∂Φ/∂z andū z = (1 + r) −1 ∂Φ/∂r, for the purpose of evaluating the accuracy of the secondary components. Figure 3a illustrates the accuracy of our quadratic model over the range 1 K 200 using a fixed aspect ratio W/H = 2 and curvature parameter = 0.01. Here, the notation u θ|i ,Φ |i indicates the approximation (3.2) is truncated beyond the i'th term. In each case the relative L 2 error is proportional to K i+1 . Observe that the approximationsū θ|1 ,Φ |1 andū θ|2 ,Φ |2 are a significant improvement overū θ|0 ,Φ |0 for K = O(10). Additionally, for K = O(100) the approximationū θ|2 ,Φ |2 offers a marginal improvement overū θ|1 ,Φ |1 . When K = 100 the relative errors are 0.95% forū θ|2 and 2.05% forΦ |2 . This illustrates that the first three terms is a reasonable trade-off between computational cost and achieving a reasonably accurate model for K = O(100).
For the remainder of the discussion we takeū θ =ū θ|2 and similarlyΦ =Φ |2 . Figure 3b illustrates that the model accuracy is insensitive to bend radius when < 1/10, using fixed K = 100 and aspect ratio W/H = 2. Figure 3c shows the variation in model error with respect to aspect ratio of the cross-section, using fixed K = 100 and = 0.01. While there is some dependence of the model error on the aspect ratio, the case W/H = 2 is near the peak for aspect ratios W/H 1 and is therefore a reasonable choice of reference value.
We take a moment to elaborate on the choice of characteristic velocity U m used in this study. The pressure gradient which drives flow through the curved duct increases super-linearly with the flow rate and also has a (weak) dependence on the bend radius. Given a K > 0, determining the pressure gradient which produces the required flow rate becomes a non-linear optimisation problem. The non-linearity is weak over the range of K considered in this study and therefore it will be convenient to avoid the non-linear optimisation. To this end, we specify the characteristic velocity U m to be the maximum velocity of a laminar flow through a straight duct having an identical cross-section. This is equivalent to non-dimensionalising with respect to the driving pressure gradient, since straight duct flow satisfies Stokes' equation and therefore doesn't depend on the Reynolds number, that is U m = c C G 2 /µ where G is the pressure gradient down the main flow axis and c C is a constant depending only on the cross-section C. This choice makes it straightforward to compare results across different duct bend radii and incorporate K as a new parameter in our inertial migration model. Figure 3d illustrates how the average of the non-dimensionalū θ andΦ varies as a function of K for our specific nondimensionalisation, with fixed = 0.01 and W/H = 2. The average ofū θ is approximately constant over this range of K indicating that our specific choice of characteristic velocity is easily translated to an average axial flow rate if desired. Figure 4 illustrates the difference in features of the fieldsū θ ,Φ between K = 0 to K = 100. Note that K = 0 should be interpreted as the limit K → 0 as a result of a decaying flow rate. U m While the physical magnitude of the secondary vortices decays to zero as U m → 0 our specific non-dimensionalisation ofū s , via the scaling U m Re in (3.2b), ensures that the magnitude of the secondary vortices remains finite. Observe in Figure 4 a shift of local maxima and minima towards the outside/right wall when K = 100. This is due to the increased rotational inertia of the fluid and is not captured by the leading order approximation of the background flow used in Harding et al. (2019). This shift becomes more pronounced with further increases in K. Lastly, we point out that the Dean numbers K 200 considered herein are significantly smaller than the critical Dean number after which there exist four vortex solutions (Winters 1987).

Our extended particle migration model
Here we incorporate the improved approximation of the background flow (3.3) into the inertial lift model so that we may efficiently approximate particle trajectories for any desired 0 K 200. We start by decomposing the leading order disturbance flow q 0 , v 0 further than was previously done in Harding et al. (2019). This approach has the additional benefit of revealing how the axial and secondary components of the background flow influence all six components ofû p ,Ω p at both leading and first order. Additionally, the decomposition admits a decoupling which can be exploited to more efficiently compute all of the required forces and torques.
Utilising the operator V, and expanding further on (2.15), observe that v 0 can be decomposed as whereû p,0 = (u p,0,r , u p,0,θ , u p,0,z ) andΩ p,0 = (Ω p,0,r , Ω p,0,θ , Ω p,0,z ). This expansion extracts all of the key parameters from the leading disturbance solution making it a linear superposition of Stokes' flow solutions. Notice that last line of (4.1) has a factor Re p and could therefore be placed in v 1 (which would be equivalent to moving the appropriate component of the boundary condition in (2.10d) to (2.11d)). However, in many practical situations κ 1 and thus κRe p may not be small. Consequently, it will be convenient to keep these terms in v 0 . An expansion identical to (4.1) is obtained for each of q 0 , F −1 and T −1 by replacing each V with Q, M and N , respectively.
Upon setting F −1 = 0 and T −1 = 0, one obtains six equations which are linear with respect to the six components of u p,0 , Ω p,0 . Applying the symmetry properties described in Appendix B allows the 6×6 system to be decoupled into two distinct 3×3 systems. One of these describes the equilibrium attained by the particle parameters u p,0,θ , Ω p,0,r , Ω p,0,z in relation to the axial motion of the background flow: The diagonal elements of A a are negative and generally the matrix is diagonally dominant. Recall that each M(·) and N (·) implicitly depends on , α and the cross-sectional coordinates of the particle (r p ,ẑ p ). Solving (4.2) therefore yields an expression for each of u p,0,θ , Ω p,0,r , Ω p,0,z which depend only on the parametersr p ,ẑ p , , α, K. The remaining 3×3 system describes the equilibrium attained by the particle parameters u p,r , u p,z , Ω p,θ in relation to the secondary motion of the background flow: Similar to A a , the diagonal elements of A s are negative and generally the matrix is diagonally dominant. We now consider the balance of force and torque at the next order of Re p by setting F 0 = 0 and T 0 = 0. The symmetries that lead to (4.2) and (4.4) analogously lead to and similarly HereR = R/a and Θ 0 = (R +r p )u p,0,θ (such that Θ 0 k = Θ 0 ). In previous studies the corrections u p,1,θ , Ω p,1 and off-diagonal elements of A a , A s were ignored. The net contribution from the first two vectors on the right side of each of (4.6) and (4.7) is typically small and can be ignored, however they have been included here for completeness.

The completed model
Using the result of the four equations (4.2), (4.4), (4.6) and (4.7) our complete migration model is: (4.10) Appendix C describes further simplifications that can be made to (4.6) and (4.7) owing to symmetries about the plane θ = 0. These result in the alternative equations (C 4) and (C 5) which provide additional insight into the interplay between axial and secondary flow components and their influence on the first order terms u p,1 and Ω p,1 . Given the spherical symmetry of the particle we don't need to calculate the specific angle of rotation which has occurred about its centre. However, the rate of rotation Ω p is needed as it influences the disturbance of the surrounding fluid and the Ω p,0 components are needed in order to calculate the u p,1 components. Observe that the Ω p,1 components are weakly coupled to u p,1 through the typically small, but non-zero, off-diagonal components of the matrices A a , A s .
This first order ODE model (4.8) assumes that the particle motion is such that the net hydrodynamic force and torque on the particle is zero at all times. In other words, the particle instantly attains terminal velocity with respect to its current position as it moves about within the cross-section. Obviously this fails to capture effects due to the acceleration of the particle and the immediately surrounding fluid. However, it is not an unreasonable assumption in regimes where the particle slowly migrates across streamlines towards a stable equilibrium. Moreover, we have found this to be a reasonably effective model in previous work so think it is worth exploring in the context of moderate Dean numbers. The development of an efficient second order model which captures acceleration effects is the subject of ongoing work and is beyond the scope of this study.

Results
For rectangular ducts having aspect ratios 2 and 4 it was noted in Harding et al. (2019) that there exists a pair of stable equilibria (with vertical symmetry about z = 0) for all of the values α ∈ {1/20, 2/20, 3/20, 4/20} over the range of −1 ∈ [40, 1280] that was considered. Herein, as we increase the Dean flow parameter K, this continues to be the case. Thus, it will be useful to restrict our study to the horizontal location of the stable equilibria pair and consider how this changes with respect to different parameters, especially K. In this study we consider two larger particle sizes α = 5/20, 6/20 in addition to those of the previous study. We first illustrate how the flow rate modifies the horizontal focusing location for two fixed duct bend radii corresponding to −1 = 80, 160. We then examine how long (in terms of both time and distance) it takes particles to focus and any influence K has on this. Lastly, we examine whether the parameter κ continues to characterise the horizontal focusing location of particles for increasing Dean number.
The computations required to compute the coefficient fields described in Section 4 were conducted using the open source finite element software FEniCS (Alnaes et al. 2015;Logg et al. 2012). The background flow was estimated over the appropriate twodimensional cross-section by solving the equations described in Harding (2019d ). Thirdorder Lagrange elements were used for this calculation over a triangular mesh having roughly 500,000 elements so that the relative error is expected to be much less than 10 −6 . A finite difference code was used for validation of the background flow approximation and is available at Harding (2019b). For the three-dimensional computations of the disturbance flow, Taylor-Hood elements were used with tetrahedral meshes which were refined around the particle and typically had from 500,000 to 1,000,000 cells over the duct section. A convergence analysis of the code was conducted in Harding (2019a) from which we are confident that the relative error in the computed coefficients is typically less than 10 −3 . The required coefficient fields were estimates via 1000 to 2000 samples within the cross-section (depending on aspect ratio and particle size) and bivariate cubic interpolation was used for the purpose of sampling these fields within standard ODE solvers. A Python class for accessing the computed coefficients in the case of small Dean number is available at Harding (2019c) and will be updated with moderate Dean number data in the near future. changes as K increases for curved ducts with rectangular cross-sections having aspect ratio 2 and 4, respectively.
Considering first a duct with aspect ratio 2 and −1 = 80. Figure 5a illustrates a small degree of separation of the different particle sizes when K = 0 (interpreted as the behaviour in the limit K → 0 as a result of a decaying flow rate). As K increases, the particles with α = 0.05, 0.15 become increasingly separated from the others. The stable equilibria for α = 0.15 initially has r * p ≈ −1, and initially shifts slightly towards the centre as K increases before swinging back towards the left/inside wall. The stable equilibria for α = 0.05 begins with r ≈ −0.16 and steadily shifts towards the right (outside wall), up to r ≈ 0.67. In contrast, the remaining particle sizes begin with centre between the other two and shift towards slightly right of centre without spreading significantly. Interestingly, the α = 0.10 particle experiences a relatively rapid shift toward the outside wall around K = 150, enough to give it some small separation from the α = 0.20 particle. Figure 5b illustrates the case with the larger bend radius −1 = 160. The α = 0.10, 0.15 particles begin separated from the others when K = 0. As K increases, so does the separation between the two groups. Moreover, the α = 0.10 particle becomes slightly separated from that with α = 0.15 for K > 150 as it shifts towards the inside wall, Similar to the −1 = 80 case, the α = 0.05 particle becomes slightly separated from the other group as it more readily migrates towards the outside wall. From a practical perspective it is useful that the clear separation between the initial two groups of particles is maintained for increasing K in this particular case.
We now examine the results for a duct with aspect ratio 4 shown in Figure 6. With −1 = 80, the α = 0.10, 0.15, 0.20 particles are closer to the inside wall than the others at K = 0. This separation increases for increasing K, largely due to r * p for those α = 0.05, 0.25, 0.30 particles shifting towards the outside wall. Additionally, each of the α = 0.05, 0.25, 0.30 particles attains a small degree of separation from each other for K approaching 200. The α = 0.15 particle also achieves some separation from the α = 0.10, 0.20 particles over the entire range of K, with maximum separation when K ≈ 125. For α = 0.05, as K approaches 200 there is a bifurcation that leads to a second pair of stable equilibria in the half of the duct nearer to the outside wall.
With −1 = 160, there are three distinct groups when K = 0, namely α = 0.10, 0.15, α = 0.05, 0.20 and α = 0.25, 0.30 ordered from the inside wall to the outside wall. The separation between the three groups is maintained for increasing K, with the first group becomes increasingly separated from the latter two.
All together, these results suggest several significant trends. First, clear separation between the stable equilibria of groups of different particle sizes at a low flow rate is not adversely affected as the Dean number increases. Second, particle sizes which focus near the left wall when K = 0 don't move significantly for increasing K and, as a consequence, become increasingly separated from the other particle sizes which shift towards the outside wall. Third, the α = 0.05 particle tends to focus near the larger particles but can achieve some separation for large enough K. The α = 0.15 particle appears to be well separated from the remainder of the particles most consistently across these results.

Trajectory illustrations
We provide illustrations of particle trajectories towards stable equilibria for some selected values of α, −1 , K. These help to demonstrate the effect of increasing Dean number on the migration dynamics of particles and assist with the interpretation of the results from Section 5.1. Figure 7 illustrates trajectories of two distinct particle sizes at three different Dean numbers in a curved duct with aspect ratio 2 and dimensionless bend radius −1 = 80. The stable equilibria pair for the α = 0.10 particle shifts to the right as K increases (as also evident in Figure 5a) and the trajectories deform accordingly. Conversely, the equilibria pair of the slightly larger α = 0.15 particle shifts very little over this range of K. For increasing K, particles starting near the inside (left) wall migrate further to the right before making their way back to a stable equilibrium along the slow manifold.  The marker size reflects the size of the particle.
Additionally, when K = 200 we observe a slight deformation of the slow manifold, evident as a pinch near its centre. Observe that the horizontal location of the stable equilibria for the two particle sizes has only a small degree of separation when K = 1 (compare 7a and 7d) but increases with increasing K (compare 7c and 7f). Figure 8 similarly illustrates trajectories of two distinct particle sizes at three different Dean numbers in a curved duct with aspect ratio 2, but with the dimensionless bend radius −1 = 160. The stable equilibria pair for the α = 0.10 particle doesn't shift much as K increases (as also evident in Figure 5b), but particles starting near the inside (left) wall migrate further to the right before making their way back to a stable equilibrium along the slow manifold. Additionally, when K = 200 we again observe a slight pinching deformation of the slow manifold. This is much like the behaviour when α = 0.15 and −1 = 80 from Figure 7. For the larger α = 0.20 particle, the stable equilibria pair shifts to the right with increasing K but with only subtle changes to the trajectories as a whole. Again we observe that the horizontal separation of stable equilibria between the two particle sizes increases with increasing K, however in this case it is the larger particle that shifts towards the outer wall, unlike the case in Figure 7. Figure 9 illustrates trajectories of two distinct particle sizes at three different Dean numbers in a curved duct with aspect ratio 4 and dimensionless bend radius −1 = 80. The stable equilibria for the smaller α = 0.05 particle shifts to the right with increasing K and the spiral trajectories become increasingly elongated. For K = 200 there is an additional pair of stable equilibria located nearer to the outer wall (as also seen in Figure 6a) and, although not easily seen in the figure, this pair captures more of the particles than the stable pair located closer to the centre. Given that our secondary flow velocity approximation has an accuracy of O(10)% when K = 200 it is unclear if a two pairs of stable equilibria occurs in practice (and is difficult to verify since fluorescent streak images from experiments typically show that small particles have not converged (Rafeie et al. 2019)). For the larger α = 0.20 particle there is almost no perceptible change in the migration dynamics over this range of K values. Figure 10 illustrates trajectories of two distinct particle sizes at three different Dean numbers in a curved duct with aspect ratio 4, but with dimensionless bend radius −1 = 160. We observe the dynamics of the smaller α = 0.10 particle changes very little over this range of K. The stable equilibria pair for the larger α = 0.20 particle shifts towards the centre along the slow manifold for increasing K but with almost no perceptible change in the dynamics of trajectories otherwise. Note that the horizontal separation between the larger and smaller particle again increases with K.

Focusing time and distance
Here we investigate the time and axial distance required for particles to focus and whether the Dean number has a significant effect on this. The migration time t m (r, z) of an individual particle starting with centre at (r p (0), z p (0)) = (r, z) ∈ C (such that the particle does not intersect the duct walls), is taken to be the earliest time after which the particle remains within 0.01 of a stable equilibrium, specifically t m (r, z) := min t 0 (r p (0), z p (0)) = (r, z) and (r p (τ ), z p (τ )) − (r * p , z * p ) /100 for all τ t . (5.1) Here (r * p , z * p ) is taken to be the specific equilibrium which (r p (τ ), z p (τ )) converges towards as τ → ∞. At the focusing time t m = t m (r, z), we also record the total angle around the curved duct travelled by the particle, that is θ p (t m ), and from this we calculate a projected axial distance d m = Rθ p (t m ).
Observe that (5.1) is a slight modification of the definition in Ha et al. (2022). We make the distinction of ensuring the particle remains within 0.01 of a stable equilibrium for all τ t because a small particle, for which the secondary flow drag is dominant, may migrate within 0.01 of (r * p , z * p ) and then away again as it follows an elliptical shaped spiral towards it. This wasn't an issue in the prior square duct study because the background flow vortices are not as elongated as they are within the rectangular ducts considered here.
For a given triple α, , K we determine the migration time as the time at which 90% of particles are determined to be focused. This is estimated by calculating t m for a large number of individual particles with starting locations defining a grid of equidistant points in C. From these times we then select the 90th percentile. Similarly, the 90th percentile is selected from the migration distances d m . Here we present results using the  Figure 11: Focusing timet m (a,c) and distanced m (b,d) versus particle size α (a,b) and Dean number K (c,d). The cross-section has aspect ratio 2. In (a,b) the colours and markers differentiate between −1 and K values respectively. In (c,d) the colours and markers differentiate between −1 and α values respectively. Below the black dotted lines all data points satisfy κ > 10, while above all but one satisfies κ < 10. dimensionless migration timet m and distanced m which are related to their physical quantities via t m =t m ( /U m )( /a)Re −1 p and d m =d m ( /a)Re −1 p , respectively. The results of these focusing time and distance measurements for a cross-section with aspect ratio 2 are depicted in Figure 11. The focusing times and distances are qualitatively similar, so we will focus our discussion on the focusing times.
Two distinct regimes are identified in the figures either side of the black dotted line: data for which κ > 10 is below; while data with κ < 10 is above; excepting one data point with α = 0.20, −1 = 40 and K = 150. The boundary κ = 10 roughly distinguishes between when the inertial lift is dominant (κ 10) and the secondary flow drag is dominant (κ 10). We observe that the κ > 10 regime generally appears to focus much quicker. This is most likely explained by the absence of a slow manifold in the cross-sectional plane when the secondary flow drag is dominant. The triple α, −1 , K = 0.2, 40, 150, an outlier for which κ > 10 but is located above the dotted line, appears to migrate slower than anticipated as it is near a bifurcation (with respect to K) which gives rise to two stable pairs within the cross-section. We note that samples for K = 200 were also calculated but ultimately excluded from these figures due to many similar outliers.
The results also demonstrate that the focusing time and distance typically appears to increase slightly with increasing α. This can be explained by the (dimensionless) drag coefficient increasing with α (as predicted by Faxen's law), and also the inertial lift force decreasing slightly in magnitude as α increases. The combination of these two effects contributes to a decreasing inertial migration velocity with respect to α, thereby leading to larger focusing times.
Additionally, we observe that increasing K typically results in an increasing spread  Figure 12: Focusing timet m (a,c) and distanced m (b,d) versus particle size α (a,b) and Dean number K (c,d). The cross-section has aspect ratio 4. Above the black dotted line all data points satisfy κ 10, while below they satisfy κ > 10. In (a,b) the colours and markers differentiate between −1 and K values respectively. In (c,d) the colours and markers differentiate between −1 and α values respectively. of focusing times. Moreover, in the regime κ < 10, the focusing times generally appear to increase with increasing K. A possible explanation is that for sufficiently large K we sometimes observe a bifurcation resulting in multiple stable focusing pairs. As K approaches this bifurcation points there is a location on the slow manifold where the migration velocity tends towards zero thereby greatly increasing the focusing time of particles that must pass through this region. Figure 12 illustrates the focusing times and distances within a cross-section having aspect ratio 4. The observations are generally the same as the aspect ratio 2 case. We again use a black dotted line to differentiate data points with κ < 10 above and κ > 10 below (with no outliers in this instance). Although the gap between the two subsets of the data is smaller in this case there is a small degree of separation. Observe that the focusing times for data with κ < 10 are much larger than those seen in the duct with aspect ratio 2. This occurs because the particles have further to migrate due to the increased duct width and, moreover, much of this extra distance is covered along the slow manifold. Analogous observations are made with respect to the focusing distance.

Comparison with an experimental study
Rafeie et al. (2019) studied particle focusing in a spiral duct with rectangular crosssection having aspect ratio 4 (designated design R3). For backward flow (going from the innermost to the outermost turn) −1 ≈ 180 near the outlet, while for forward flow (from the outermost to the innermost turn) −1 ≈ 74 near the outlet. These two values are reasonably close to the two values −1 = 160, 80 presented in our study. They experimented with three particle sizes α ≈ 0.04, 0.067, 0.1 and used flow rates from 0.5 to 9.0 mL/min in steps of 0.5 mL/min (corresponding to Dean numbers from K ≈ 2 up to 840 at the innermost turn). Their results suggest the smallest particle remains unfocused for both flow directions, and appears to be distributed more uniformly for larger flow rates. On the other hand, the larger two particle sizes start focused near the inside wall at low flow rates and shift significantly towards the outside wall as the flow rate increases.
Bearing in mind that K = 200 only covers up to a flow rate of approximately 5.0 mL/min for this specific device (with K measured on the innermost turn), we observe that the shifting of the two larger particle sizes towards the outside wall is qualitatively similar behaviour to our predictions for a α = 0.05 particle. Curiously, our α = 0.10 results suggest this particle stays near the inside wall, something that is not observed in the experiments for this particle size. In the case of the smallest particle in the experiment (α ≈ 0.04), our model predicts that the duct is much too short, by roughly a factor of 10, to achieve full focusing (see section 5.3), and this is supported by the wide spread of streaks. Moreover, the amount of spread in the experimental results for α ≈ 0.067 suggests the duct is too short for complete focusing of this particle size as well. This too is supported by our model which predicts the duct needs to be roughly 3 times longer.

Approximate collapse of horizontal focusing location
In Harding et al. (2019) it was demonstrated that at low flow rates the horizontal location of stable equilibria pairs approximately collapses when plotted against κ. Here we investigate whether this remains the case when the Dean number increases.
In Figure 13 we plot the horizontal focusing location of stable equilibria, denotedr * p , for the three Dean numbers K = 50, 100, 150 in the case of a duct cross-section with aspect ratio 2. To facilitate comparison between different particle sizes, herer * p has been nondimensionalised with respect to half of the duct height. The left plots show the change in r * p versus −1 which is of practical interest when different size particles are suspended in flow through a duct having a specific bend radius. The right plots show the change in r * p versus κ and illustrates how this parameter characterises the focusing behaviour. The top row illustrates the case K = 50 which is similar to the low flow rate results of Harding et al. (2019), with just a slight shift towards the outside wall (r = 2) of stable equilibria which focus near the centre (r = 0). We also include two additional particle sizes α = 0.25, 0.30 here. The middle row shows the case of K = 100 and, qualitatively, the results are very similar. The main difference is that the increased flow rate shifts stable equlibria located near the centre (most noticeably those with |r| 1/2) further towards the outside wall (r = 2). Additionally, for a narrow range of −1 , we observe that α = 0.05 particles no longer focus to stable equilibria but, instead, onto small stable orbits around unstable equilibria. The bottom row shows the case of K = 150. The more pronounced shift in centrally focused equilibria is seen and there are subtle qualitative changes, particularly in relation to how sharplyr * p changes for larger κ values. Additionally, the stable orbits that occur for the particle size α = 0.05 now exist over a larger range of −1 and cover a greater portion of the cross-section.
The approximate collapse of focusing behaviour against κ is preserved remarkably well as the Dean number K is increased to a moderate size. This result is significant in that is indicates that sized based separation of particles via inertial focusing should be reasonably robust to changes in flow rate. Moreover we observe that particles focused closer to the centre of the duct shift towards the outer wall as the flow rate increases whereas those that are focused closer to the inside wall don't move quite so much. This results in a slightly deeper trough in the plots of r * p against κ and further support the previous observation that separation is typically improved with increased flow rate. The duct crosssection has aspect ratio 2 andr * p is non-dimensionalised with respect to /2. The light shaded area illustrates the region occupied by a stable orbit which occurs only when α = 0.05 for K 100.
As discussed in Harding et al. (2019), the tail of data points where α = 0.05 and κ < 10 is due to stable equilibria appearing near the centre of the inside and outside walls. These occur for large −1 as the duct becomes straighter and secondary flow drag becomes negligible in comparison to the inertial lift force such that the equilibria of a straight rectangular duct are attained. Similar stable equilibria near the centre of the inside wall exist for α = 0.10 and κ < 10, although we expect these to ultimately disappear again in the limit κ → 0 (as a straight duct has no such equilibria). The basin of attraction of these particular equilibria is very small and so it is reasonable to ignore them from a practical viewpoint (although were included here for completeness).
We briefly remark on the existence of stable orbits when α = 0.05 for sufficiently large K. For these parameters there are no stable equilibria but instead a vertically symmetric pair of unstable (spiral) equilibria and around each of these there is a reasonably tight stable orbit. These stable orbits occur over a relatively small parameter space for rectangular ducts in comparison to the case of square ducts whose dynamics are studied in Ha et al. (2022). In this instance the required bend radius is O(1000) times the duct height which is unlikely to occur in a practical setting. Moreover, we expect that experimental noise would impact the stability of such orbits thereby making them difficult to observe experimentally. Figure 14 shows r * p against −1 and κ for a cross-section with aspect ratio 4 and with Dean numbers K = 50, 100, 150. The main observation is that again centrally located equilibria shift closer to the outside wall with increasing K which ultimately means separation can occur across a larger portion of the duct width. The approximate collapse of the focusing curves against κ for κ < 10 is incredibly well preserved as K increases. The wider aspect ratio duct doesn't feature any stable equilibria appearing near the centre of the inside or outside wall, and stable orbits are very small where they do occur, thereby making the wider aspect ratio duct more attractive for size based separation.

Conclusions
We have extended our model of inertial migration of neutrally buoyant spherical particles in curved microfluidic ducts to moderate Dean number by incorporating further terms of a suitable perturbation approximation of the background flow. This approach is particularly powerful as the computed coefficient arrays can be cheaply re-assembled to model particle trajectories at any desired flow rate for which the model is applicable. We have demonstrated a smooth shift in stable equilibria with respect to the Dean number K and observed that increasing K can enhance the separation of particles by size. Moreover, we demonstrated that the approximate collapse of the horizontal location of equilibria when plotted against κ is robust to increasing K.
We also investigated how long it takes particles to focus with respect to a time scale that accounts for the expected scaling with flow rate. On this time scale, we found that results could be loosely separated into two distinct categories, faster focusing for κ > 10, and slower focusing otherwise. This can be explained by the existence of a slow manifold most particles must travel along when inertial lift is the dominant force.
Non-neutrally buoyant particles could easily be incorporated into this model by adding the appropriate perturbations as described in Harding & Stokes (2020). We chose to focus on neutrally buoyant particles in this work to simplify the presentation and because our previous study concluded that small variations in buoyancy have negligible effect on particle dynamics.
A natural question going forwards is to understand how particle migration is perturbed by higher flow rates. This poses several challenges to the approach used herein. First, the perturbation expansion of the background flow (3.3) fails to converge for K values larger than roughly 212 (Harding 2019d ). One solution could be to expand around a non-zero K value, or instead interpolate between solutions computed for several large K values. Second, the assumption that Re p remains suitable small becomes questionable and therefore additional terms in the perturbation expansion of the disturbance flow with respect to Re p may need to be considered (although it is unclear if the radius of convergence for this expansion is large enough for this to be useful). Third, it may no longer be appropriate to neglect the effects of particle acceleration, particularly for smaller particles which migrate rapidly around the secondary vortices. All of these challenges will be considered going forwards as we continue to expand and improve our model of inertial migration in curved ducts.
Appendix A. Application of the reciprocal theorem to calculating hydrodynamic force and torque For completeness we begin by recalling the Lorentz reciprocal theorem for Stokes flow in its most general form and then revisit its application to the calculation of the force on a portion of the boundary of the fluid domain. Following this we illustrate how it can also be applied to the calculation of the torque about a desired reference point over a portion of the boundary of the fluid domain.
Theroem. Let p, u and q, v be the solution of (A 1c) where a, b, f , g are suitably smooth vector fields over their respective domains, then where n denotes the outward pointing normal vector with respect to Ω.
over Ω. In general, given a differentiable tensor field S, one has ∇·(v·S) = v·(∇·S)+∇v : S where : is the tensor double dot product. Therefore (A 3) can be written as Subtracting the second equation from the first, then noting ∇v : pI = p(∇ · v) = 0 and ∇u : qI = q(∇ · u) = 0, and additionally ∇u : ∇v = ∇v : ∇u and ∇u : ∇v = ∇u : ∇v = ∇v : ∇u , then It remains to integrate this equation over Ω and apply the divergence theorem to the left hand side. The result is then straightforward to re-arrange to obtain (A 2). Now we consider the reciprocal theorem applied to the calculation of the hydrodynamic force from a flow p, u on a subset of ∂Ω, e.g. denoting the surface of a particle. In particular, we consider calculating the component of the force on the surface Γ ⊂ ∂Ω projected onto a constant unit vector e * . Corollary. Let p, u and q, v satisfy (A 1) with a = 0 over ∂Ω, g = 0 over Ω, b = e * is a constant vector on Γ ⊂ ∂Ω, and b = 0 over the remainder of ∂Ω\Γ . Then Proof. Immediately, due to a = 0 and g = 0, one obtains from (A 2) the result Since v = b = e * is constant over Γ (and is zero elsewhere) it follows that One simply needs to multiply both sides by −1 to obtain the desired hydrodynamic force in the direction e * .
Observe that in the corollary v could be written as V(e * ) using the operator V introduced in Section 2.1, taking Γ = ∂P . Now, we wish to apply the reciprocal theorem again to the calculation of the torque Γ (x − x c )× (−n) · −pI + µ ∇u + ∇u dS , (A 10) about a (fixed) reference point x c , e.g. the centre of mass of a particle. Corollary. Let p, u and q, v satisfy (A 1) with a = 0 over ∂Ω, g = 0 over Ω, b = e * ×(x − x c ) on Γ ⊂ ∂Ω where e * is constant, and b = 0 over the remainder ∂Ω\Γ . Then Proof. This time (A 2) yields the result Γ e * ×(x − x c ) · n · −pI + µ ∇u + ∇u dS = Ω v · f dV . (A 12) Applying the triple product shift rule (x×y) · z = (y×z) · x to the integrand on the left hand side produces The dot product with e * can be safely pulled outside the integral (as e * is constant). Multiplying both sides by −1 then yields the desired component of torque in the e * direction, i.e. (A 11).
Observe that using the operators introduced in Section 2.1 we can write v = V(e * ×(x− x c )) within this corollary, taking Γ = ∂P . Appendix B. Symmetries associated with Stokes flow around a spherical particle in a curved duct.
Here we describe a variety of quantities which are zero due to symmetries associated with our problem. Specifically, in the rotating frame the fluid domain F is symmetric about θ = 0 and many fields of interest relating to solutions of (2.13) are either odd or even with respect to θ . Specifically, the following are odd with respect to θ while the following are even with respect to θ These properties can be used to show: Moreover, given a sufficiently smooth function f (r , z ), i.e. which is independent of θ , then each of are odd with respect to θ , and each of Of particular interest is when f (r , z ) described appropriate components of the background flow velocity fieldū, as described in Section 3.
Using the symmetries described in Appendix B it can be shown that I 0,a · V(j ) is odd with respect to θ such that F V(j ) · I 0 dV = F V(j ) · I 0,s dV .

(C 3)
Analogous simplifications can be made for a number of other volume integrals over a dot product between I 0 and a V(. . . ) term. Upon applying these symmetry properties to (4.7) and (4.7) we arrive at j · P ū a · ∇ū s +ū s · ∇ū a dV i · P (x −x p )×(ū a · ∇ū s +ū s · ∇ū a ) dV k · P (x −x p )×(ū a · ∇ū s +ū s · ∇ū a ) dV    +    F V(j ) · I 0,s dV F V(i ×(x −x p )) · I 0,s dV F V(k×(x −x p )) · I 0,s dV i · P ū a · ∇ū a +ū s · ∇ū s dV k · P ū a · ∇ū a +ū s · ∇ū s dV j · P (x −x p )×(ū a · ∇ū a +ū s · ∇ū s ) dV This decomposition illustrates a clear separation in how different parts of the leading order disturbance flow solution contribute to the hydrodynamic force generated by the first order disturbance flow solution. In (C 4) and (C 5)ū a should be interpreted as u a,0 + Kū a,1 + K 2ū a,2 and similarly forū s . Both v 0,a , v 0,s from (C 1) can be substituted into I 0,a , I 0,s in (C 2), and subsequently the volume integrals overF in (C 4) and (C 5) can be expressed as a polynomial function of u p,0 , Ω p,0 , Θ 0 , κ, Re p , K with variable coefficients depending on r p , z p , , α. By sampling the coefficient fields over a suitable range of r p , z p , , α and interpolating appropriately, we can efficiently estimate the values of u p,1 , Ω p,1 given any value of u p,0 , Ω p,0 , Θ 0 , κ, Re p , K.