# Inertial Migration In Micro-Centrifuge Devices

Samuel Christensen<EMAIL_ADDRESS>Marcus Roper  UCLA Department of
Computational Medicine
UCLA Department of Mathematics

###### Abstract

Within microcentrifuge devices, a microfluidic vortex separates larger
particles from a heterogeneous suspension using inertial migration, a
phenomenon that causes particles to migrate across streamlines. The ability to
selectively capture particles based on size differences of a few microns makes
microcentrifuges useful diagnostic tools for trapping rare cells within blood
samples. However, rational design of microcentrifuges has been held back from
its full potential by a lack of quantitative modeling of particle capture
mechanics. Here we use an asymptotic method, in which particles are accurately
modeled as singularities in a linearized flow field, to rapidly calculate
particle trajectories within microcentrifuges. Our predictions for trapping
thresholds and trajectories agree well with published experimental data. Our
results clarify how capture reflects a balance between advection of particles
within a background flow and their inertial focusing and shows why the close
proximity of trapped and untrapped incoming streamlines makes it challenging
to design microcentrifuges with sharp trapping thresholds.

††preprint: APS/123-QED

##  I Introduction

Microcentrifuges are a recently developed class of microfluidic devices that
can be used to selectively trap large particles from flowing suspensions. The
devices consist of a series of chambers connected with microfluidic channels.
Within each chamber, one or more eddies may form. _Inertial migration_ causes
particles in moderate Reynolds number flows to travel across streamlines.
Larger particles migrate faster and are more likely to become trapped within
the microcentrifuge chamber [1, 2]. Size-based trapping may be used to trap
the largest particles within the suspension, and shows promise as a tool for
analyzing cell types in a patient blood sample, e.g. for isolating large
circulating cancer cells from small red blood cells [3, 4, 5].

Figure 1: A diagram of the microcentrifuge (not to scale). Microfluidic
channels lead to the microcentrifuge chamber. At sufficient Reynolds numbers,
a fluid eddy forms in the chamber and inertial migration pushes larger
particles into the chamber, where they are then trapped in the eddy, while
smaller particles are not captured and continue through the device.

The range of possible particle behavior within the complex three dimensional
flows that occur in microcentifuge chamber is not well understood, nor is
there, to our knowledge, has any straightforward mechanistic description of
how a microcentrifuge operates. In addition to the intrinsic interest of
identifying new inertial microfluidic phenomena, mechanistic explanation of
microcentrifuge function opens a door to solving the reverse problem, of
designing a chamber or channel geometry to target a particular size threshold
for trapping.

A large impediment to this understanding is how much less well developed the
theory of inertial microfluidic migration is, relative to theories modeling
the behaviors of particles in zero Reynolds number flows. Although inertial
focusing starts at arbitrarily low Reynolds numbers, in practical
implementations for useful trapping thresholds to be accessed, the devices so
far built have all operated at Reynolds numbers between 50 and 300. Inertial
migration is caused by the rigid particle disrupting the otherwise smoothly
varying fluid flow throughout the channel. Although the physics of inertial
migration in (uni-directional) pipe flow has been studied extensively, (see
e.g. [6] and [7]): the key difference between pipe flow and more complex flow
is the varying flow that the particle experiences as it advects through the
channel. In this paper, we develop an asymptotic theory that reduces the
calculation of particle migration velocity to solving a quasi-steady, linear
problem.

Direct numerical simulation of particle trajectories within inertial
microfluidic devices requires solving for the motion of particles suspended in
a fluid-filled domain whose boundaries constantly change due to the movement
of the particles. Since, in practical examples, channels operate at moderate
Reynolds number, between 5-200, flows cannot be modeled by numerical methods
that are designed for small Reynolds number particle-flows such as Stokesian
dynamics [8] or boundary integral methods [9]. The time evolving geometry of
migrating particles and nonlinear terms favors numerical methods such as
immersed boundary [10, 11] or immersed interface [12] which embed moving
boundaries within a fixed computational grid. Both methods afford a lot of
freedom in the choice of numerical method for solving the Navier-Stokes
equations; in inertial microfluidic simulations the Lattice-Boltzmann method
(LBM) is a popular method [13] and has been used for calculating particle
migration [2, 14, 15]. The Force Coupling Method[16] replaces particles with
forcing terms which enforce rigid body motion on a finite volume of fluid to
emulate particle dynamics and has been used to calculate inertial
migration[17]. Although existing numerical simulations have illuminated the
physics of inertial focusing, the high computational cost of nonlinear 3D
simulations has meant that predictive simulations are not currently used to
design or optimize inertial microfluidic devices.

In our simplification of the physics of inertial migration, a 3D linear PDE
models moderate Reynolds number fluid flow and supplies a vector field
representing inertial migration that may be added to the background flow. This
simplification allows for easier understanding of the various sorting
techniques that are used by microfluidic devices. In this study, we will first
develop our equations for particle motion and then we will compare our
theoretical results with data from a microfluidic device separating different
cell types suspended in blood along with differently sized plastic beads.

##  II Mathematical Methods

Inertial migration is caused by the disturbance the particle causes to the
background flow, the disturbance flow
$\mathbf{u}^{\prime}=\mathbf{u}-\bar{\mathbf{u}}$ is the difference between
the flow with the particle, $\mathbf{u}$, and the flow without the particle,
(called the background flow here) $\bar{\mathbf{u}}$. We get our main
equations by plugging these definitions into the Navier-Stokes equations and
non-dimensionalizing by the speed scale $U$ and length scale $L$ of the
background flow:

$\displaystyle\Delta\mathbf{u}^{\prime}-\nabla p^{\prime}$ $\displaystyle=$
$\displaystyle Re\left(\frac{\partial\mathbf{u}^{\prime}}{\partial
t}+\bar{\mathbf{u}}\cdot\nabla\mathbf{u}^{\prime}+\mathbf{u}^{\prime}\cdot\nabla\bar{\mathbf{u}}+\mathbf{u}^{\prime}\cdot\nabla\mathbf{u}^{\prime}\right)$
(1) $\displaystyle\nabla\cdot\mathbf{u}^{\prime}$ $\displaystyle=$
$\displaystyle\mathbf{0}$ $\displaystyle\mathbf{u}^{\prime}$ $\displaystyle=$
$\displaystyle\mathbf{U}_{p}-\bar{\mathbf{u}}(\mathbf{x}_{p})+\Omega_{p}\times(\mathbf{x}-\mathbf{x}_{p})\qquad\text{on
}|\mathbf{x}-\mathbf{x}_{p}|=a/L$ (2)

In addition to $\mathbf{u}^{\prime}=0$ on the walls of the channel. Neither
boundary conditions nor the incompressibility equation are altered by
following transformations of our equation, so we focus on the different
incarnations of the momentum balance equation, Eq. 1. Here
$Re=\frac{UL\rho}{\mu}$ is the channel Reynolds number, and is typically in
the range 50-300, both in real experiments and our calculations. The
disturbance velocity is generated by the boundary condition on the sphere. The
particle can translate and rotate with the fluid, but it resists the shearing
motion of the fluid near its surface, setting up the disturbance velocity
$\mathbf{u}^{\prime}$. Near the particle, the velocity profile is
approximately that of simple shear;
$\bar{\mathbf{u}}\approx\bar{\mathbf{u}}(\mathbf{x}_{p})+\bm{\gamma}\cdot(\mathbf{x}-\mathbf{x}_{p})$,
where $\gamma_{ij}\equiv\frac{\partial\bar{u}_{i}}{\partial x_{j}}$, and we
may further decompose the linearized flow field into rotational and straining
components:
$\bar{\mathbf{u}}\approx\bar{\mathbf{u}}(\mathbf{x}_{p})+\mathbf{\omega}\times(\mathbf{x}-\mathbf{x}_{p})+\mathbf{E}\cdot(\mathbf{x}-\mathbf{x}_{p})$,
where $\omega_{i}=-\frac{1}{2}\epsilon_{ijk}\gamma_{jk}$, with
$\mathbf{\epsilon}$ the unit alternating tensor, and
$E_{ij}=\frac{1}{2}(\gamma_{ij}+\gamma_{ji})$. Setting
$\mathbf{U}_{p}\approx\bar{\mathbf{u}}(\mathbf{x}_{p})$ and
$\Omega_{p}=\omega$ renders the particle force and torque free at first order.

We now know the size of all the components involved in right hand side of Eq.
1: under our non-dimensionalization, both $\mathbf{u}^{\prime}$ and
$\bar{\mathbf{u}}$ are of size $O(\frac{a}{L})$ near the particle. The size of
$\frac{\partial\mathbf{u}^{\prime}}{\partial t}$ is related to the boundary
condition Eq. 2 who’s size is the same as
$\frac{\partial\mathbf{E}(\mathbf{x}_{p}(t))}{\partial
t}=\mathbf{U}_{p}\cdot\nabla\mathbf{E}$ which is small as long as the length
scale at which $\mathbf{E}$ changes is large compared to $a$. Therefore the
right hand side is size $\frac{a^{2}}{L^{2}}Re$. The size ratio
$\alpha:=\frac{a}{L}\ll 1$ is small enough that we assume $\alpha^{2}Re<1$,
which forms the core of our asymptotic expansion. Near the particle, we assume
shear dominates and set the inertial terms to zero. Suppressing these inertial
terms, we arrive at: $\Delta\mathbf{u}^{\prime}-\nabla p^{\prime}=0$, along
with the usual incompressibility and boundary conditions. This problem of
Stokes’ flow around a sphere has solution:

$\mathbf{u}^{\prime}=-\mathbf{E}:\frac{5(\mathbf{x}-\mathbf{x}_{p})(\mathbf{x}-\mathbf{x}_{p})(\mathbf{x}-\mathbf{x}_{p})}{|\mathbf{x}-\mathbf{x}_{p}|^{5}}+O(|\mathbf{x}-\mathbf{x}_{p}|^{-4})$
(3)

[18]. However, this solution is not consistent with our complete neglect of
the inertial terms from Eq. 1, because as
$|\mathbf{x}-\mathbf{x}_{p}|\to\infty$, $\Delta\mathbf{u}^{\prime}\sim
1/r^{4}$, while $Re\,\bar{\mathbf{u}}\cdot\nabla\mathbf{u}^{\prime}\sim
Re/r^{2}$, becomes co-dominant with viscous stresses when
$|\mathbf{x}-\mathbf{x}_{p}|\sim Re^{-1/2}$. We therefore posit that the flow
contains an outer region in which inertial terms may not be neglected[6].
However, since decay of the disturbance velocity means that
$|\mathbf{u}^{\prime}|\ll|\bar{\mathbf{u}}|$ within this region, we may
linearize the inertial terms within this outer region.

In the outer region, we model the particle as a moving singularity.
Mathematically, this means replacing the rigid body motion of Eq. 2 with a
forcing term equal to
$F(\mathbf{x})=-\frac{20\pi}{3}E(\mathbf{x}_{p}):\nabla\delta(\mathbf{x}-\mathbf{x}_{p}(t))$.
Instead of working with a moving particle, we approximate this by looking for
a traveling wave solution of the form
$\mathbf{u}^{\prime}(\mathbf{x}-\mathbf{x}_{p}(t))$, this approximation
transforms the time derivative into
$\frac{\partial\mathbf{u}^{\prime}(\mathbf{x}-\mathbf{x}_{p}(t))}{\partial
t}=\mathbf{U}_{p}\cdot\nabla\mathbf{u}$ which we will approximate as
$\frac{\partial\mathbf{u}^{\prime}(\mathbf{x}-\mathbf{x}_{p}(t))}{\partial
t}=\bar{\mathbf{u}}(\mathbf{x}_{p})\cdot\nabla\mathbf{u}$.

Our solution in the inner region does not predict particle migration at
leading order; instead the stresslet disturbance field modeled in Eq. 3,
forces the outer region disturbance velocity. This forcing can be represented
equivalently as a boundary condition as $\mathbf{x}\to\mathbf{x}_{p}$, or,
directly, by introducing a force dipole term within the equation:

$\displaystyle\nabla\mathbf{u}^{\prime}-\nabla p^{\prime}$ $\displaystyle=$
$\displaystyle
Re\left((\bar{\mathbf{u}}-\bar{\mathbf{u}}(\mathbf{x}_{p}))\cdot\nabla\mathbf{u}^{\prime}+\mathbf{u}^{\prime}\cdot\nabla\bar{\mathbf{u}}\right)-\frac{20\pi}{3}E:\nabla\delta(\mathbf{x}-\mathbf{x}_{p})$
(4) $\displaystyle\nabla\cdot\mathbf{u}^{\prime}$ $\displaystyle=$
$\displaystyle 0$

The migration velocity is found by solving this PDE and evaluating it at the
location of the particle
$\mathcal{M}(\mathbf{x}_{p})=\mathbf{u}^{\prime}_{\mathbf{x}_{p}}(\mathbf{x}_{p})$,
where the subscript denotes that the singularity was located at
$\mathbf{x}_{p}$ according to Eq. 4. We then advect the particle according to
following differential equation

$\displaystyle\frac{d\mathbf{x}_{p}}{dt}=\bar{\mathbf{u}}(\mathbf{x}_{p})+\alpha^{3}\mathcal{M}(\mathbf{x}_{p})+\frac{\alpha^{2}}{6}\Delta\bar{\mathbf{u}}(\mathbf{x}_{p})$
(5)

This equation includes the term for Faxen’s Law, which is the finite size
particle correction for force free advection of particles in $Re=0$ flows.
This formula is exact for unbounded $Re=0$ flows, but in our use is accurate
to the order $O(Re^{0}\alpha^{4})$.

When calculating the trajectories of particles, we pre-calculated the
migration velocity at 450 points throughout the three dimensional channel by
Eq. 5. We then used linear interpolation to extend our particle velocity
calculation through the entire domain.

Eq. 5 can be used to to advect particles throughout the channel, however its
accuracy is limited to conditions where the underlying asymptotic expansion is
valid. In straight channels it has been shown that this asymptotic
approximation is accurate for $\alpha^{2}Re=O(1)$[6], however both the size
and background flow speed of the microfluidic device can change significantly
depending on if we are in a tight channel or large chamber, eg. the large
notched chamber in the microcentrifuge. if we define the Reynolds number using
the average flow rate in a cross section with length $L$, we see that
$Re=\frac{\int_{A}\bar{\mathbf{u}}dSL}{L^{2}\nu}=\frac{C}{\nu L}$ where C is
the flow rate of the device and is constant for all cross sections in the
channel. In this calculation we see that the Reynolds number decreases like
$\frac{1}{L}$ in larger channels, so the asymptotic expansion that forms Eq. 5
is more accurate in the larger testing chambers than it is in the smaller
channels leading up to it.

Another potential source of error in the asymptotic approximation is our
assumption of the traveling wave solution. While
$\frac{\partial\mathbf{u}^{\prime}(\mathbf{x}-\mathbf{x}_{p}(t))}{\partial
t}=\bar{\mathbf{u}}(\mathbf{x}_{p})\cdot\nabla\mathbf{u}$ is exact in the bulk
flow, the particle is moving with respect to the walls of the device, which
will cause additional changes to the flow as $\mathbf{x}_{p}$ moves. It can be
shown that this error is proportional to
$\mathbf{U}_{p}\cdot\nabla\mathbf{u}^{\prime}|_{\text{walls}}$ which is small
provided the particle is either far away from the wall or not headed directly
at it.

##  III Analysis of Microcentrifuge

###  III.1 Background Flow Patterns at Different Reynolds numbers

In [1], Khojah et al. observed a microcentrifuge trapping differently sized
cells from blood samples. They found that at different Reynolds numbers it
preferentially captured different sizes of cells. At $Re=125$ the micro
centrifuge consistently captured larger particles and smaller particles could
pass through. At $Re=175$ they found that the microcentrifuge would
inconsistently capture particles of all sizes. At $Re=225$ the microcentrifuge
would capture smaller particles with more consistency than it would capture
larger particles.

Figure 2: The fluid flow within the microcentrifuge experiences a topological
change as Re changes from 125 to 225. The top row shows an example streamline
and the bottom row shows the Poincaré section of many streamlines and how they
progress through the channel. Blue shows where the streamlines enter and red
shows where the streamlines exit. The trajectory cartoon demonstrates how the
entrance and exit of the streamlines loop back on themselves and change their
locations in the channel, representing a topological change of how the fluid
moves through the channel.

In numerical studies using COMSOL Multiphysics [19], we found that these
Reynolds numbers correspond with major changes to the background flow.
Streamlines were calculated for Reynolds numbers 125, 175, and 225 were
calculated using P2+P1 tetrahedral elements at the ’extremely fine’ mesh size
setting.

Fig. 2 shows how the way the fluid flows through the microcentrifuge chamber
changes as we increase the Reynolds number. The top row shows 3D visualization
of example particle trajectories and the bottom row shows via Poincaré section
how the streamlines travel through the chamber once they are caught, the dots
represent where the trajectories intersect a plane at the top of each of the
orbits the streamlines performed within the chamber. Blue represents where the
streamlines entered and red represents where they exited the chamber. The x-
and y- coordinates of streamlines entering and exiting the chamber are shown
with asterisks instead of dots.

In Fig. 2, we see that at $Re=125$ the fluid flows through the microcentrifuge
chamber by having the corner streamlines flow down and into the chamber, then
the fluid streamlines rotate towards the center of the chamber and then spiral
outwards, eventually exiting the chamber. At $Re=175$ the path of the
streamlines is looping back on itself and is no longer consistently focusing
towards the center of the channel. At $Re=225$, the loop has inverted: instead
of flowing from out to in, the streamlines enter the chamber from the center
of the channel and leave out the sides. Inset diagram in the top row of Fig 2.
show the exchange of starting and ending streamline locations as $Re$ is
increased.

###  III.2 Particle Trajectories

Figure 3: The trajectories of background flow (left), $24\mu m$ diameter
particles (middle), and $28\mu m$ diameter particles. Small particles follow
the background flow and do not enter the microcentrifuge chamber, capture is
predicted to begin at approximately $22.9\mu m$, but only particles along the
bottom focusing position are captured. For larger particles, particles along
all but the top focusing position are captured. All captured particles stably
focus towards the mid plane at this Reynolds number.

We will now start looking at the microcentrifuge capture particles, but first
we need to look into the initial condition. Before the particle reaches the
microcentrifuge, it will travel through a narrow rectangular channel (in our
case $40\mu$m$\times 70\mu$m), at $Re=125$ the particles gather along 4 stable
focusing streamlines near the center of each channel wall [7, 20]. This means
the particles will be tightly grouped along 4 focusing streamlines before they
enter the channel, with the majority of particles evenly split between the two
larger walls of the channel (corresponding to the blue and purple trajectories
in Fig. 3 [21]).

In Fig 3, we place differently sized particles along the 4 focusing
streamlines and advance their position using Eq. 5. All 4 focusing streamlines
do not enter the notch, at this Reynolds number the streamlines that enter the
channel are along the bottom corners [22], however at particle size $24\mu$m
particles along the bottom streamline migrate away from the streamline enough
that they become captured in the eddy within the notch. at particle size
$28\mu$m the particles along the minor streamlines are also captured. At
particle size $35\mu$m the particles along the top streamline are captured as
well, however these particles are $87.5\%$ of the channel height and the
identification of 2 distinct focusing streamlines may no longer be relevant.
Our calculations also show that above above diameter $28\mu$m, the particles
from the minor streamlines will quickly focus towards the mid plane 3. Capture
of particles along the minor axis is driven by the strong inertial migration
towards the microvortex, inertial focusing actually fights the particle from
migrating towards the center of the channel but is overpowered by the
background flow and the particle is ultimately driven towards the mid plane,
The rapid convergence of all particles toward the symmetry plane agrees with
numerical results from [2].

Changes in the background fluid flow recorded in section A completely change
the way particle capture is achieved and explain the differences in what size
particles are captured by the micro-centrifuge. At $Re=125$, particles are
captured when inertial migration pushes the particle into the microcentrifuge
chamber, the limit cycle is created as a product of the balance between the
fluid flow spiraling the particle upwards and inertial migration pushing the
particle downwards at the top of the spiral. At $Re=175$ all particles are
predicted to enter the microcentrifuge chamber, but it is not well understood
what keeps them in the close to turbulent eddy that is spiraling within the
chamber. We can integrate the physics of capture with the description in
section A of the three dimensional paths of streamlines through the chamber.

Particles that enter the microcentrifuge chamber continue to experience
inertial migration within that chamber. At $Re=125$, the streamlines that
particles follow, absent inertial migration, spiral outwards within the mid
plane Fig. 2, while inertial migration points downward into the chamber. The
balance of these two effects caused particle trajectories to converge to a
limiting orbit, shown in Fig. 3. We report in more detail upon the limiting
orbit in section C.

###  III.3 Comparison With Data

Looking only along the mid plane, we make predictions about the critical
Reynolds number for particle capture and compare with experimental data.
Kohjoah et al.[1] performed an experiment in which a micro-centrifuge was used
to separate breast cancer cells (MDA-MB-231 cell line) based on size.
Inputting their parameters into our asymptotic method we made a prediction on
what the critical cell diameter the device sorts the cells by. Our models
predict a critical cell diameter of $22.9\mu$m and this prediction shows good
agreement with the data 4.

Figure 4: Predicted critical cell diameter for cell sorting matches well with
data from micro-centrifuge device in [1]. Top left: A continuum of differently
sized particle trajectories, smaller particles are in blue and larger
particles in red with the critical particle diameter (22.9$\mu$m) is
highlighted in black. Top right: Capture data from experiment shows difference
in cell size distribution between captured cells and the initial distribution
of cells ran through the device. Bottom left: The estimated percentage
captured by cell size shows that capture that the beginning of capture is
centered around the critical capture diameter. The estimation was necessary
because the initial cell size data was given in much coarser bins than the
captured cell size data. Bottom right: Estimated percentage captured by cell
size is shown along with the percentage captured of trajectories simulation
with Eq. 5 with added Gaussian white noise with $\sigma=12.5\%$ of expected
particle velocity.

The experiment consisted of running a blood sample containing 224 cells
through the microfluidic chamber 3 times, the experiment was performed twice
for a total of N=448 cells run through, with 166 cells captured between the
two trials. The device consists of an inlet, a $40\mu m\times 70\mu$m$\times
3$cm channel connecting to the $800\mu$m$\times 70\mu$m$\times 800\mu$m
microcentrifuge chamber, another $40\mu$m$\times 70\mu$m$\times 3$cm, and an
outlet. The distribution of inflowing and captured cells were measured
directly from images of the flowing cells, and we use the histograms reported
in [1] to estimate capture.

We approximated the initial cell size distribution by fitting a log-normal
distribution (mean: $20.2\mu$m, standard deviation: $8.2\mu$m) to available
data using the maximum likelihood method. we then approximated the percentage
of cells captured in a given size range by dividing the observed number of
cells that were captured by the expected number of cells in that size range
based on our fitted distribution. For the very largest cells (diameters
exceeding $32.75\mu$m) there was insufficient sampling sampling of inflowing
particles sizes, leading to capture probabilities that could exceed 1; we set
capture probability equal to 1 for these largest cells.

In order to estimate the magnitude of the noise in the system, we fit a
modified version of Eq. 5 with added Gaussian white noise. If we define the
RHS of Eq. 5 as
$\mathbf{V}(\alpha,\mathbf{x}_{p})=\bar{\mathbf{u}}(\mathbf{x}_{p})+\alpha^{3}\mathcal{M}(\mathbf{x}_{p})+\frac{\alpha^{2}}{6}\Delta\bar{\mathbf{u}}(\mathbf{x}_{p})$,
our SDE is equal to

$\frac{d\mathbf{x}_{p}}{dt}=\mathbf{V}(\alpha,\mathbf{x}_{p})+||\mathbf{V}(\alpha,\mathbf{x}_{p})||\mathcal{N}(0,\sigma^{2})$
(6)

We found qualitatively that $\sigma=0.125$ fit the randomness present in the
data accurately (fig. 4). Percentage capture of differently sized particles
were estimated by computing 1000 trajectories for each level of $\alpha$ and
calculating the percentage that were captured.

###  III.4 Explanation of Variation

Our model is deterministic, so all cells above the critical diameter are
predicted to be captured and below are not. However, the experimental results
reported in [1] show that capture probability increase continuously from 0 to
1 as particle radius increases over a narrow range from $10\mu$m to $30\mu$m.
The likely cause for this indeterminacy of capture is likely a combination of
cells on the top streamline not becoming captured and hydrodynamic interaction
between the cells. The micro-centrifuge is intended to be a high throughput
device and the cells are focused tightly along the streamlines of the channel
feeding the micro-centrifuge chamber. As these cells are near each other, they
affect each other’s inertial migration due to shear induced dispersion [23,
24].

Figure 5: Trajectories of different particles are displayed (purple) over the
capture region associated with each sized particle. The red area represents
starting points where the particle will not be captured by the vortex and the
blue area represents the region where the particle will be captured by the
vortex. All trajectories are fairly close to the bifurcating streamline,
however trajectories of particles near the critical diameter of $22.9\mu m$
are very close to the separating streamline and whether or not they get
captured is influenced by noise.

Inertial focusing is generally the strongest contribution in a particle’s
cross-stream velocity, so it is not immediately apparent that these unmodeled
interactions could be large enough to prevent or ensure capture. However,
particle capture trajectories are sensitive to noise. In Fig. 5 we visualize
the basin of attraction of the microcentrifuge chamber, and we show how close
the trajectories of particles approaching the chamber on the lower major
focusing streamline are to the separating manifold between the captured region
and the un-captured region. In Fig. 5 the two colored regions represent
invariant sets where a trajectory that starts inside a region will never leave
it.

Each size of particle has different invariant regions because the focusing and
curvature term that control particle off-streamline migration in Eq. 5 change
with particle size, however the separating manifold between the two invariant
regions is always close to the trajectory of a particle. Particles close to
the critical diameter, pass extremely closely to the capturing trajectory for
almost the entire length of the chamber, meaning that very small perturbations
could permit capture of a particle below the critical capture or deny capture
in particles larger than the critical diameter.

Figure 6: A: A quiver plot of the strength of inertial focusing throughout the
channel shown alongside the trajectories of particles of diameters $0\mu m$,
$10\mu m$, $20\mu m$, and $30\mu m$. The strength of migration is related to
color, with arrows not shown being $<\frac{1}{100}$ of maximum strength.
Strength of inertial migration is measured by
$\mathcal{M}(\mathbf{x})/\mathbf{u}(\mathbf{x})$ which represents the
infinitesimally small distance the inertial migration will push the particle
at a given point in space. The cutout shows how inertial focusing is strongest
in the channel leading to the chamber and near the walls. B: The particle
limit cycle does not trap stream lines, shown with the dotted line, which
slowly spiral outwards eventually exiting the microcentrifuge cavity. C: Time
lapsed limit cycles from micro-beads of various sizes is overlaid with
predictions from our asymptotic simulations. Smaller Particles have larger
limit cycles, with the larger particles ($>30\mu m$) having very small limit
cycles. While the smaller particles are not predicted to be captured via
inertial focusing, their limit cycle is well captured by our asymptotic model.

Next we demonstrate why the trajectories of the particles are so close to
capturing stream lines using a phase plane diagram showing their inertial
migration velocity vectors. this diagram can be scaled to particles of
different sizes by scaling the captured velocity vectors by $\alpha^{3}$ and
we superimpose upon the vector field plot the trajectories of three different
particle sizes. We see that the inertial migration is strongest in the channel
leading up to the micro-centrifuge and fairly weak near the vortex itself.
This is because inertial migration depends upon shear, curvature, and walls
[25]. Near the vortex the shear is low so the inertial migration is low. The
trajectories of particles near the critical particle diameter all experience
very similar inertial migration, the larger particles are captured because
they experience slightly more. We are unaware of any experimental observations
of particle trajectories during capture, or more generally, for particles
crossing the microcentrifuge and this experimental gap is likely due to the
high speed of particles within the channels feeding each cavity. By contrast,
the limiting trajectories of captured particles within the channels have been
extensively reported on. It is known that captured particles are drawn into
closed orbits (limit cycles) within the micro centrifuge; our simulations
replicate these dynamics (Fig. 6A) and agree quantitatively with measured
limit cycle paths for particles of different sizes at $Re=125$ (Fig. 6) [26,
27, 2]. In particular we find that larger particles converge to smaller
orbits. Orbit shape is well predicted, including a common center of gyration,
a long straight trajectory that is tilted and closely conforming to the
separating streamline that divides captured and non-captured particles (Fig.
5). To fit experimental data we had to compute particle sizes from the figures
of [1]. We chose diameters $15.3\mu m$, $16.2\mu m$, and $20.5\mu m$ because
they fit the sizes listed in out parts of [1]. The relative sizes we measured
in figure were 8.1, 8.6, and 12.2 pixels respectively. There are some
discrepancies, though it is not possible to determine whether they are due to
imperfect matching of the Reynolds number or particle sizes, unmodeled effects
such as the device flexing slightly when under pressure or whether they are
due to approximations made in our model.

The two smaller particles are not predicted to be captured by our simulation,
but nevertheless they are predicted to have a stable limit cycle if they enter
the cavity. This feature complicates the calculation of capture probabilities;
small particles do not migrate across the separating manifold in Fig 5, but if
perturbations, such as due to multi-particle effects, induce them to cross it
they will be stably entrained in the vortex eddy. The smallest particle that
is predicted to have a stable limit cycle in our simulation had a diameter of
$10.9\mu$m which lines up well with data where the smallest cell captured had
diameter between $10-12.5\mu$m. the complicating presence of these smaller
particles within the microcentrifuge chamber emphasizes that the vortex is
unrelated to the capture of particles, the vortex only controls the retention
of particles.

###  III.5 How to Change Critical Particle Diameter

While the calculations results above are a significant analysis of the
particles, they cannot be extended to other channels as the background flow
depends on the Reynolds number. However, if we change the channel size and
flow speed to keep our Reynolds number constant, we can reuse all of our
calculations but with a different parameter set, allowing us to select for the
critical particle diameter that our microfluidic device selects for:

$\displaystyle Re$ $\displaystyle=$
$\displaystyle\frac{UL}{\nu}=\frac{\frac{U}{k}kL}{\nu}$ $\displaystyle\alpha$
$\displaystyle=$ $\displaystyle a/kL$ $\displaystyle\alpha^{*}_{k}$
$\displaystyle=$ $\displaystyle\alpha^{*}/k$

For example, if you need to change the device so that instead of separating
particles above and below $22.9\mu m$ to selecting particles above and below
$11.45\mu m$, all you need to do is half the size of your device and double
your speed.

There is a practical limitation however, the PDMS polymer that the
microfluidic devices are etched from do not come at an arbitrary thickness.
However, there have been previous studies that if you scale the fluid speed
and channel/cavity height accordingly the change in thickness does not effect
the fluid flow [22]. We speculate, although leave the testing to future work
that this change is also unlikely to effect the inertial migration
significantly, that as the optimal Reynolds number for separation drives
particles to the symmetry plane, the furthest place from the side walls.

##  IV Perspectives

By replacing particle’s with singularities, we were able to accurately model
inertial migration without solving the nonlinear Navier-Stokes equations or
explicitly meshing the 3D and time-varying fluid domain between particles and
channel walls. This method allowed us carefully analyze all possible particle
trajectories within a microcentrifuge device. Our analysis shows two clear
takeaways for microfluidic devices designed to separate particles by size:

  1. 1.

The preformance of the microcentrifuge is deeply dependent on Reynolds number.
Lower Reynolds numbers advect particles towards the symmetry plane of the
device and allow for more stable separation of particles by size. Higher
Reynolds numbers push particles away from the symmetry plane of the device and
introduce more complicated 3D effects, which are both difficult to anatomize
and have been found empirically to cause devices to be have wider, less
consistent capture thresholds.

  2. 2.

Error in sorting is due to lack of separation between particle trajectories
and the capture separatrix. Particles within $3\mu$m of the critical capture
diameter are straddling the capture separatrix for 75% the length of the
microcentrifuge, making their capture susceptible to particle-particle
interaction and other forms of random noise. Inertial migration velocity scale
as $O(\alpha^{3})$, the high order size scaling leads to sufficient separation
for sorting particle accurately by size, but is not enough to have low error
rates by itself.

Experimental studies that have optimized for capture accuracy have
independently agreed with this design principles. [1, 2, 28] have all found
optimal Reynolds numbers between 100-150 and [28] were able to achieve 90%
capture accuracy at 90% purity by optimizing over Reynolds number. We identify
the effect of this optimization as finding maximum separation between focusing
position in the channel and separating manifold.

Our singularity method models only the first order inertial migration, further
expansions in terms of particle size and Reynolds number are possible, and
would add either higher order forcing terms within our Oseen equation, or
would require that we impose additional matching conditions to model the
particle. Such an extension could be useful for predicting differential
focusing of different particles within the same channel. Unpacking the role of
particle interactions and extending the calculation to the largest particles
on which separation is performed. The simulations here took approximately 10
minutes per solve of inertial migration, however computes could be sped up
significantly by using discontinuity decompositions from [21] to . The method
can be extended to particles of other shapes, singularity modeling of the
particle needs only the stresslet strength associated with the Stokes solution
around the particle. This stresslet strength is already known for ellipsoids
in shear flow[29] and for other particles can be found by solving for the
motion of the particle in Stokes flow, for which there are many approximate or
numerical methods [30, 18].

Changes to channel geometry can also lead to significant increases in
separation accuracy. Paié et al. [31] found a 20% increase in accuracy of
capture by creating a separate reservoir for captured particles and optimizing
over 6 parameter combinations. Numerically driven parameter optimization is
still extremely computationally expensive, here we have done our best to
analyze indepth the basic microcentrifuge geometry and describe the physical
principles that that fuel its strengths and weaknesses. Our results provide a
platform on which accelerated studies may be used to further optimize the
accuracy, controllability, or throughput of future microcentrifuges.

## References

  * Khojah  _et al._ [2017] R. Khojah, R. Stoutamore, and D. Di Carlo, Size-tunable microvortex capture of rare cells, Lab on a Chip 17, 2542 (2017).
  * Haddadi and Di Carlo [2017] H. Haddadi and D. Di Carlo, Inertial flow of a dilute suspension over cavities in a microchannel, Journal of Fluid Mechanics 811, 436 (2017).
  * Che  _et al._ [2016] J. Che, V. Yu, M. Dhar, C. Renier, M. Matsumoto, K. Heirich, E. B. Garon, J. Goldman, J. Rao, G. W. Sledge, _et al._ , Classification of large circulating tumor cells isolated with ultra-high throughput microfluidic vortex technology, Oncotarget 7, 12748 (2016).
  * Dhar  _et al._ [2018] M. Dhar, J. Wong, J. Che, M. Matsumoto, T. Grogan, D. Elashoff, E. B. Garon, J. W. Goldman, E. S. Christen, D. Di Carlo, _et al._ , Evaluation of pd-l1 expression on vortex-isolated circulating tumor cells in metastatic lung cancer, Scientific reports 8, 1 (2018).
  * Chen  _et al._ [2012] J. Chen, J. Li, and Y. Sun, Microfluidic approaches for cancer cell detection, characterization, and separation, Lab on a Chip 12, 1753 (2012).
  * Schonberg and Hinch [1989] J. A. Schonberg and E. Hinch, Inertial migration of a sphere in poiseuille flow, Journal of Fluid Mechanics 203, 517 (1989).
  * Hood  _et al._ [2015] K. Hood, S. Lee, and M. Roper, Inertial migration of a rigid sphere in three-dimensional poiseuille flow, Journal of Fluid Mechanics 765, 452–479 (2015).
  * Sierou and Brady [2001] A. Sierou and J. F. Brady, Accelerated stokesian dynamics simulations, Journal of fluid mechanics 448, 115 (2001).
  * Liron and Barta [1992] N. Liron and E. Barta, Motion of a rigid particle in stokes flow: a new second-kind boundary-integral equation formulation, Journal of Fluid Mechanics 238, 579 (1992).
  * Kazerooni  _et al._ [2017] H. T. Kazerooni, W. Fornari, J. Hussong, and L. Brandt, Inertial migration in dilute and semidilute suspensions of rigid particles in laminar square duct flow, Physical Review Fluids 2, 084301 (2017).
  * Nakagawa  _et al._ [2015] N. Nakagawa, T. Yabu, R. Otomo, A. Kase, M. Makino, T. Itano, and M. Sugihara-Seki, Inertial migration of a spherical particle in laminar square channel flows from low to high reynolds numbers, Journal of Fluid Mechanics 779, 776 (2015).
  * Xu  _et al._ [2020] S. Xu, B. Gao, A. Lofquist, M. Fernando, M.-C. Hsu, H. Sundar, and B. Ganapathysubramanian, An octree-based immersogeometric approach for modeling inertial migration of particles in channels, Computers & Fluids , 104764 (2020).
  * Bazaz  _et al._ [2020] S. R. Bazaz, A. Mashhadian, A. Ehsani, S. C. Saha, T. Krüger, and M. E. Warkiani, Computational inertial microfluidics: a review, Lab on a Chip 20, 1023 (2020).
  * Sun  _et al._ [2016] D.-K. Sun, Y. Wang, A.-P. Dong, and B.-D. Sun, A three-dimensional quantitative study on the hydrodynamic focusing of particles with the immersed boundary–lattice boltzmann method, International Journal of Heat and Mass Transfer 94, 306 (2016).
  * Chun and Ladd [2006] B. Chun and A. Ladd, Inertial migration of neutrally buoyant particles in a square duct: An investigation of multiple equilibrium positions, Physics of Fluids 18, 031704 (2006).
  * Loisel  _et al._ [2013] V. Loisel, M. Abbas, O. Masbernat, and E. Climent, The effect of neutrally buoyant finite-size particles on channel flows in the laminar-turbulent transition regime, Physics of Fluids 25, 123304 (2013).
  * Abbas  _et al._ [2014] M. Abbas, P. Magaud, Y. Gao, and S. Geoffroy, Migration of finite sized particles in a laminar square channel flow from low to high reynolds numbers, Physics of Fluids 26, 123301 (2014).
  * Kim and Karrila [2013] S. Kim and S. J. Karrila,  _Microhydrodynamics: principles and selected applications_ (Courier Corporation, 2013).
  * Multiphysics [1998] C. Multiphysics, Introduction to comsol multiphysics®, COMSOL Multiphysics, Burlington, MA, accessed Feb 9, 2018 (1998).
  * Di Carlo  _et al._ [2007] D. Di Carlo, D. Irimia, R. G. Tompkins, and M. Toner, Continuous inertial focusing, ordering, and separation of particles in microchannels, Proceedings of the National Academy of Sciences 104, 18892 (2007).
  * Christensen  _et al._ [2022] S. Christensen, R. Chu, C. Anderson, and M. Roper, Fast asymptotic-numerical method for coarse mesh particle simulation in channels of arbitrary cross section, Journal of Computational Physics 471, 111622 (2022).
  * Khojah  _et al._ [2021] R. Khojah, D. Lo, F. Tang, and D. Di Carlo, The evolution of flow and mass transport in 3d confined cavities (2021).
  * Da Cunha and Hinch [1996] F. Da Cunha and E. Hinch, Shear-induced dispersion in a dilute suspension of rough spheres, Journal of fluid mechanics 309, 211 (1996).
  * Hood and Roper [2018] K. Hood and M. Roper, Pairwise interactions in inertially driven one-dimensional microfluidic crystals, Physical Review Fluids 3, 094201 (2018).
  * Di Carlo [2009] D. Di Carlo, Inertial microfluidics, Lab on a Chip 9, 3038 (2009).
  * Mach  _et al._ [2011] A. J. Mach, J. H. Kim, A. Arshi, S. C. Hur, and D. Di Carlo, Automated cellular sample preparation using a centrifuge-on-a-chip, Lab on a Chip 11, 2827 (2011).
  * Zhou  _et al._ [2013] J. Zhou, S. Kasper, and I. Papautsky, Enhanced size-dependent trapping of particles using microvortices, Microfluidics and nanofluidics 15, 611 (2013).
  * Wang  _et al._ [2013] X. Wang, J. Zhou, and I. Papautsky, Vortex-aided inertial microfluidic device for continuous particle separation with high size-selectivity, efficiency, and purity, Biomicrofluidics 7, 044119 (2013).
  * Chwang and Wu [1975] A. T. Chwang and T. Wu, Hydromechanics of low-reynolds-number flow. part 2. singularity method for stokes flows, Journal of Fluid mechanics 67, 787 (1975).
  * Happel and Brenner [2012] J. Happel and H. Brenner,  _Low Reynolds number hydrodynamics: with special applications to particulate media_ , Vol. 1 (Springer Science & Business Media, 2012).
  * Paiè  _et al._ [2017] P. Paiè, J. Che, and D. Di Carlo, Effect of reservoir geometry on vortex trapping of cancer cells, Microfluidics and Nanofluidics 21, 1 (2017).
  * Haddadi  _et al._ [2018] H. Haddadi, H. Naghsh-Nilchi, and D. Di Carlo, Separation of cancer cells using vortical microfluidic flows, Biomicrofluidics 12, 014112 (2018).