Inertial and viscoelastic forces on rigid colloids in microfluidic channels

The reliable sorting and filtration of colloidal particles based on their size, shape, and rigidity is a challenging problem which is relevant for a wide variety of applications, ranging from water purification¹ to biotechnology.^2–4 For example, colloids can be sorted by size using sedimentation, since particles of different weights settle at different rates and thereby separate. However, such an approach has several shortcomings. The buoyancy of the solute particles has to be tuned carefully.⁵ Moreover, thermal motion plays a significant role on colloidal length scales, further impeding the partition into monodisperse samples.

In the past decade, microfluidic technologies have emerged as promising candidates for manipulating colloidal particles, since they allow for high-throughput automation of physical, chemical, and biological processes by linking multiple functional elements into a pipeline to form a “lab on a chip.”⁶ In principle, it is possible to spatially separate the flowing colloids by imposing external stimuli, such as electric^7,8 or magnetic fields.⁹ In practice however, such active methods are often impractical as they rely on an intrinsic particle property or careful doping of the solute particles. Alternatively, the lateral displacement of solute particles can be controlled by exploiting non-linear flow effects.^10,11 The advantage of these passive techniques is that they are independent of solute chemistry and therefore can be applied to a wide range of materials. Furthermore, these effects can be sustained or even amplified at high flow rates, guaranteeing high throughputs.

There have been multiple experimental, theoretical, and simulation studies of particle segregation under flow. Specifically, Segré and Silberberg observed the formation of annular distributions of neutrally buoyant 0.8 mm–1.6 mm (non-colloidal) spheres in cylindrical channels, when the flow was dominated by inertial effects.¹² Under these conditions, the Navier-Stokes equation describing the fluid dynamics is not time-reversible anymore, and the parabolic velocity profile of the solvent results in a net lift on the solute particles, which pushes them away from the channel center. At the same time however, the wall-induced asymmetry in the wake vorticity field of the solute particles leads to a net repulsion away from the walls.^11–14 Recent simulations have demonstrated that this mechanism holds down to the colloidal regime,¹⁵ although thermal fluctuations can lead to strongly scattered solute distributions around the annulus.

For most practical filtering and sorting applications, it is however more desirable to achieve spatial focusing of the solute particles onto a single line or plane, since this better facilitates the extraction and purification of the target species. One promising pathway for accomplishing this goal is through viscoelastic focusing.^16–21 Early experiments^16–18 revealed the migration of non-colloidal solute particles to the region of lowest shear rate (i.e., the channel center in the case of Poiseuille flow) when the solute was dispersed in a viscoelastic medium and exposed to non-inertial flow. The lateral migration of solute particles has been described by continuum fluid mechanics as an imbalance of non-linear compressive elastic forces in the viscoelastic medium.^18,19

In inertial flow, the elastic force from the viscoelastic medium F_E is in direct competition with the outward hydrodynamic lift forces F_L near the channel centerline. Scaling laws for these forces have been formulated in the limit of macroscopic solute particles. Di Carlo et al.²² showed that the lift force on non-colloidal particles far from the walls in square channels scales as

where ρ_s is the solvent density, v_max is the maximum velocity in the channel (directly related to the Reynolds number Re), a is the particle diameter, and L_x is the channel width. The elastic force F_E has been attributed to a normal stress gradient in the viscoelastic medium across the rigid particle.¹⁹ Assuming a negligible second normal stress difference,²³ this leads to

where N₁ is the first normal stress.

For dilute solutions of high molecular weight polymer, $N 1 ∼ γ ̇ β$ where $γ ̇$ is the local shear rate and 1 < β ≤ 2.²⁴ At low shear rates, N₁ is related to the storage modulus G′ in the low frequency limit¹⁹ 

with oscillatory frequency ω. In most cases, the storage modulus of a viscoelastic medium scales as G′ ∼ ω² in the low frequency regime,²⁵ and so $N 1 ∼ γ ̇ 2$ near the centerline. For the upper convected Maxwell model, the first normal stress can be calculated analytically as $N 1 =2η τ R γ ̇ 2$⁠,²³ where τ_R is the characteristic relaxation time of the polymer. The shear rate across the channel can be characterized by $γ ̇ ∼ v max / L x$⁠, leading to

These scaling laws have been derived in the continuum limit where the rigid solute is significantly larger than the molecular constituents of the viscoelastic medium, which allows for a mean-field approximation of the polymer solution. Interestingly, F_E exhibits a similar scaling behavior as the hydrodynamic lift force F_L [cf. Eq. (1)], which sets in at Re ≫ 1 and pushes the solute away from the channel centerline. Hence, we anticipate an intermediate flow regime, where the leading terms in F_L and F_E cancel each other out, and the dynamics are dominated by lower-order terms. In a recent publication,²⁶ D’Avino et al. employed a continuum description of the viscoelastic medium in conjunction with computational fluid dynamics (CFD) to study the effect of inertia in two-dimensional systems, and successfully identified a regime where F_L and F_E are competing.

It is unclear, however, whether such a description is valid in the colloidal limit, where the solute particles are of comparable size to the polymer. This is certainly the case in recent experiments,²⁰ where sub-micron polystyrene spheres were immersed in a poly(ethylene oxide) (PEO) solution with a radius of gyration of roughly 100 nm.²⁷ Additionally, thermal fluctuations play a significant role at these length scales and cannot be neglected anymore.¹⁵ In order to study focusing on the colloidal scale, we recently performed mesoscale molecular dynamics (MD) simulations with multiparticle collision dynamics (MPCD) where we explicitly modeled the constituent polymers of the viscoelastic medium.^28,29 We identified a broad parameter range in which the flow-induced focusing of the solute particles was possible, demonstrating the general feasibility to sort and filter colloidal particles. As in experiments, we measured the fraction of colloids focused to the channel centerline.

In the article at hand, we provide a more quantitative and rigorous analysis of the forces responsible for viscoelastic focusing in the inertial regime. From the measured forces, we extract information about the expected quality of focusing and sensitivity of focusing forces to relevant design parameters like polymer molecular weight and flow rate. Unlike previous computational efforts to study inertial viscoelastic focusing,²⁶ our simulations are three-dimensional and incorporate thermal fluctuations and microscopic detail of the viscoelastic medium. We demonstrate here that these fluctuations and details impart unique characteristics to colloidal focusing that are not present for macroscopic particles and must be considered when designing microfluidic devices.

Viscoelastic focusing has been studied experimentally for both macroscopic and colloidal particles dispersed in aqueous solvents with low polymer concentrations.^4,20,21,30 Typically, the rigid colloids are made of polystyrene, while the polymers are linear PEO or poly(vinylpyrrolidone) (PVP) chains with molecular weights in the range of 10²–10³ kg/mol. Accordingly, we modeled our systems as rigid spherical particles dispersed in a dilute solution comprised of the polymers and an underlying Newtonian solvent. As in the experiments, the polymer concentration can be varied to control the rheology of the solution. Because of the inherent disparity of length and time scales between the Newtonian solvent, the polymer, and the colloids, we employed a multiscale hybrid simulation approach based on multiparticle collision dynamics.^31–33 

We adopt a similar molecular model as in Ref. 28, which we briefly summarize here. Polymers were represented explicitly as bead-spring chains consisting of M monomers. We modeled the steric interactions between the beads through the purely repulsive Weeks-Chandler-Anderson (WCA) pair potential,³⁴ which corresponds to good solvent conditions. We defined length scales relative to the monomer diameter a_m = 1 and energy scales relative to the monomer-monomer interaction energy ϵ_m = k_BT = 1. Bonds between monomers were modeled with the finitely extensible nonlinear elastic (FENE) potential.³⁵ We adopted the standard Kremer-Grest parameters for the springs³⁶ in order to prevent unphysical bond crossing. Spherical colloids with diameter a were modeled using the WCA potential parameterized as in Ref. 28 so that they had only excluded volume interactions with the polymers. Both the colloid and polymer interactions with the channel walls were modeled through a purely repulsive 9-3 potential.³⁷ We set the wall interaction energy to ϵ_w = 0.1 to prevent the wall from excluding too much volume and narrowing the effective channel size accessible to the colloids and polymer. The equations of motion for the colloids and polymers were integrated using the standard velocity Verlet scheme with a step size of Δt_MD = 2 × 10⁻³ in MPCD reduced units.

The underlying Newtonian solvent was modeled on the mesoscopic scale using the MPCD technique.^31–33 In MPCD, the solvent is modeled explicitly as an ensemble of ideal point particles with unit mass m = 1 that propagate through a series of alternating streaming and collision steps. In contrast to classical MD algorithms, the solvent-solvent and solvent-solute interactions are not governed by pairwise potentials but instead through stochastic collisions. This approach leads to a significant computational speedup, while still preserving the correct hydrodynamic interactions.^32,33

During each streaming step, all solvent particles move ballistically over a period Δt and are reflected off the channel walls using modified bounce back collision rules³⁸ to achieve no-slip boundary conditions. The solvent particles interact with the rigid colloids through momentum exchanges, where solvent particles are reflected off the colloid surfaces and impart linear and angular momenta to the colloids, following the algorithm described in Ref. 39. In the MPCD collision step, all solvent particles and monomers are first sorted into cubic cells of edge length a_m and then undergo a stochastic collision with particles in the same cell where the particle momenta are rotated around a randomly chosen unit vector with a fixed angle α.^31–33 At each collision step, the cells are shifted by a random three-dimensional vector with components drawn uniformly on [ − a_m/2, + a_m/2], to ensure Galilean invariance.^40,41 We applied cell-level velocity rescaling as a thermostat in order to maintain an isothermal system, which is necessary in non-equilibrium simulations or when viscous heating is significant.³² 

The viscosity and diffusivity of the solvent are dictated by the solvent density ρ_s, rotation angle α, and collision time Δt.^32,40,41 For our simulations, we chose ρ_s = 5, α = 130°, and Δt = 0.1, which results in a dynamic viscosity of η = 3.96 for the pure Newtonian solvent. We chose a lower solvent density and higher collision time than in Ref. 39, where ρ_s = 10 and Δt = 0.04, for computational expediency. The choice of MPCD parameters in the present work still corresponds to a liquid-like solvent. The mesoscopic character of the MPCD scheme makes it difficult to directly map the reduced viscosity back to physical units,³³ and it is instead better to characterize the fluid by dimensionless numbers. A reliable measure of the importance of hydrodynamics is the Schmidt number Sc = ν/D, which is the ratio between the kinematic viscosity ν and the diffusivity D of the fluid. We find Sc ≈ 12 for the MPCD parameters employed in this work, which is sufficiently large for achieving liquid-like behavior.^42,43

Polymers and colloids were density matched to the solvent according to the corresponding rules for our embedding schemes: $m m = ρ s a m 3$ for the monomer mass and m_c = πρ_sa³/6 for the colloid mass. We maintained a constant monomer number density ρ_m in all simulations. We define this density in terms of a volume fraction $ϕ m =π ρ m a m 3 /6$⁠, which we fixed at ϕ_m = 0.05. This leads to a polymer mass fraction of 9%, which is comparable to typical experiments with PVP as the viscoelastic agent.^4,30 As we will show in Sec. III B, this monomer concentration creates a nearly constant viscosity elastic fluid, which agrees well with experimental observations for aqueous PVP solutions.³⁰ 

All simulations were conducted in a slit channel as schematically shown in Figure 1. We set the dimensions to L_x = 50 in the gradient direction x, L_y = 40 in the vorticity direction y, and L_z = 150 in the flow direction z. The channel was periodic in y and z, and we checked that it was sufficiently large in these dimensions to avoid artificial self-interactions with the periodic images. We generated Poiseuille flow by applying a body force in the flow direction to all solvent particles, as is the case for gravity driven flow in a pipe. For the pure Newtonian solvent, this procedure gives rise to a parabolic velocity profile v_z(x), which in a slit channel takes the form

where g is the acceleration associated with the body force and sets the flow rate. We quantify the strength of the flow in terms of the channel Reynolds number Re = v_maxL_x/ν, defined by the slit diameter and the theoretical maximum velocity $v max =g L x 2 / ( 8 ν )$ reached at the channel center. In non-equilibrium MPCD simulations, it is recommended to keep the solvent velocities well below the medium’s speed of sound c_s to avoid artificial compression effects.^15,38 Accordingly, we restricted our simulations to Mach numbers Ma = v_max/c_s < 0.5. Given the isothermal speed of sound c_s = 1 in the MPCD liquid, this leads to a maximum Reynolds number of Re ≈ 32.

In order to test the scaling of the forces on a colloidal particle, we directly computed the effective (mean) force on the particle using a method similar to Prohm et al.¹⁵ They calculated the solvent-mediated forces on a single colloid F_s in a purely Newtonian liquid using MPCD by fixing the lateral position of the colloid and measuring the time-averaged momentum exchanged from the solvent to the colloid

where Δp is the momentum exchange during each MPCD streaming step. F_s is then the average lift force on the particle during this time interval Δt. This method was shown to give a potential of mean force correctly describing the particle probability distribution (position of colloids in the channel).

We extended the method of Ref. 15 to measure the additional polymer-mediated force exerted by the pairwise repulsions between the colloid and polymer F_p

where F_i,c is the force exerted by polymer i on the colloid, and the sum is taken over all polymers. The average effective force on the particle is then given by F = F_s + F_p. The direction of F controls the direction in which the particle will migrate and so can be used to identify focusing positions in the channel.

We measured F as a function of displacement from the channel centerline. We fixed the lateral position of the colloid in the xy-plane but allowed it to flow freely along the z direction. Because the force should be antisymmetric around the channel centerline, we displaced the colloid only in the positive x direction. We ran extended simulations with the colloid fixed at each position to generate initial configurations and confirmed that a steady state had been reached by monitoring the polymer density profiles and average gyration tensors as well as the solvent momentum and drift in the average force on the colloid. We selected three statistically independent configurations at each lateral position once the steady state was reached, and restarted the MPCD simulations with a 50 000 MPCD collision warm-up period for the solvent flow. The stochastic nature of the MPCD algorithm randomizes the polymer velocities during the warm-up, and so the three simulations are statistically independent. We then calculated the effective forces on the colloid during 500 000 MPCD collision intervals. We computed the autocorrelation time of the total gradient force on the colloid and found that for our parameters this force was effectively decorrelated between MPCD collisions. Accordingly, we sampled the total forces every MPCD collision, averaged these forces over each independent simulation, and estimated uncertainties based on the averages from the three simulations.

Figure 2 shows the effective force in the gradient direction F_x vs. the lateral center of mass displacement for a colloid of diameter a = 16. We measured F_x for a purely Newtonian solvent as well as polymer solutions of increasing chain length at flow rate Re = 32. We found that a single run was sufficient to calculate the effective forces for the Newtonian solvent with reasonable uncertainty, and so error bars are based on a single run for this case. The displacement from the centerline is measured in reduced coordinates ξ = (2x − L_x)/L_x, defined so that ξ = 0 when the colloid center of mass is at the channel centerline and ξ = 1 when it is at the wall. Due to its rigid nature, the colloid cannot move all the way to ξ = 1. The contact distance for the colloid and wall is ξ = 1 − a/L_x, and so for a = 16, the colloid touches the wall at ξ = ± 0.68. Forces very close to the wall are sensitive to the choice of wall potential used to approximate the interaction between a hard sphere and a hard surface. Moreover, we found it challenging to obtain reliable statistics on the forces very near the wall, and so we restricted the colloid to |ξ| ≤ 0.6.

Because displacements from the centerline were made for ξ > 0, a positive force corresponds to an outward motion of the colloid, while a negative force pushes the colloid towards the centerline. For the Newtonian solvent (no polymer), F_x is positive near the channel centerline and then becomes strongly negative near the wall. These two regimes correspond to the lateral lift force driving the colloid outward and the inward hydrodynamic wall repulsion. The point where F_x crosses zero (ξ ≈ 0.45) gives the location of the equivalent of the Segré-Silberberg annulus¹² in slit flow. The value from our simulations is slightly smaller than the experimentally measured one of ξ ≈ 0.6, which has been determined for non-colloidal particles with a ≪ L_x in channels with circular¹² and square²² cross sections. This discrepancy is most likely due to the rather large blockage ratio in our simulations, which is a/L_x = 0.32 for a = 16, and is consistent with the results of Ref. 15 in a cylindrical channel for large blockage ratios.

For short chains of length M = 10 and M = 20, the same general trend is observed as for the Newtonian solvent, but the magnitude of the lift force is somewhat decreased by the presence of the polymer. This decrease can be attributed to the elastic force induced by the polymer, not yet large enough to overcome the hydrodynamic lift. Hence, the M = 10 and M = 20 chains do not lead to viscoelastic focusing on the centerline. However, the outward force is considerably smaller for the M = 20 chains compared to the M = 10 case. Increasing the chain length to M = 40 leads to a significant shift in F_x(ξ), and this shift becomes more sizable for M = 80. At these longer chain lengths, the elastic force from the polymer is sufficiently large to overcome the inertial lift force, and the particle focuses on the centerline from all positions (F_x(ξ) ≤ 0). These trends qualitatively match results from two-dimensional CFD calculations of discs in viscoelastic flow, which have been conducted at comparable flow rates,²⁶ as well as our earlier simulations.^28,29

The total force on the colloid is a combination of a hydrodynamic contribution from the solvent and an elastic contribution from the polymer. The hydrodynamic force exerted by all the polymers should be approximately the same provided that the different polymer chain lengths similarly perturb the flow field in the vicinity of the colloid. Based on macroscopic models (e.g., upper convected Maxwell model), the elastic force should be proportional to the relaxation time of the chain τ_R. At rest, the relaxation time scales with the radius of gyration $τ R ∼ R g 3$⁠,⁴⁴ and R_g in turn is proportional to M^3/5.⁴⁵ Hence, the elastic force is expected to scale as M^9/5.

We can separate the elastic component of the force by subtracting the hydrodynamic contribution to the total force. The most convenient way to do this is to subtract the total force exerted on the colloid by a solution of monomers (M = 1), which is essentially a purely hydrodynamic force, from the total force in the polymer solutions. In doing so, we assume that the presence of bonds in the polymers does not significantly perturb the flow field near the colloid compared to the monomer flow field. Then, all remaining forces must be solely due to the polymer elasticity.

Figure 3 shows the elastic contribution of the force F_E as a function of colloid position for the different chain lengths. For ξ ≲ 0.5, the elastic force drives the colloid towards the center of the channel for all chain lengths, and the strength of this focusing force increases as M increases. Figure 4 shows these forces replotted against M at a fixed channel position. From these data, it is clear that the elastic force has a considerably weaker scaling than M^9/5. In fact, fitting the data through a power law suggests that F_E scales in a sublinear manner with M.

Figures 2 and 3 show that there is a crossover in the force curves near the walls (ξ ≳ 0.5), where the elastic force is starting to push the colloid away from the centerline to the channel walls. This effect becomes more pronounced with increasing chain length and could potentially lead to the surface pinning of colloids at sufficiently low Re. However, at the investigated flow rates, the hydrodynamic wall repulsion outweighs the elastic contribution close to the walls, leading to focusing irrespective of the initial colloid position. The distance ξ at which the inward and outward forces on the colloid are perfectly balanced is often referred to as the “neutral surface.” D’Avino et al. confirmed its existence for shear thinning fluids at Re ≪ 1 using CFD and microfluidic experiments but found no neutral surface for viscoelastic solvents with constant shear viscosity.³⁰ 

In order to better understand the potential role of shear thinning solvent viscosity, we measured $η ( γ ̇ )$ for all investigated polymer lengths M by applying steady shear with the Müller-Plathe reverse perturbation method.^33,46,47 In this method, momentum flux along the gradient direction of the flow is induced through momentum exchanges, and the emerging shear rate of the solution is then measured. In our implementation, momentum swaps were only applied to the MPCD solvent particles, which then exchanged momentum with the embedded polymers through the stochastic collision step in the MPCD algorithm. We verified the validity of this approach by conducting additional simulations at selected state points, where we applied the momentum swaps to both monomers and MPCD solvent particles. All simulations were conducted in the bulk to eliminate wall-induced effects.⁴⁸ The range of $γ ̇$ was chosen to cover typical values of the Poiseuille flow simulations in the slit geometry.

To determine the degree of shear thinning, we considered the polymer solutions as power-law fluids and fit the shear stress through $K γ ̇ n$⁠. In this model, the viscosity of the fluid is quantified by the flow consistency index K, while the dimensionless index n characterizes the type of the fluid (n < 1 pseudoplastic; n = 1 Newtonian; n > 1 dilatant). For the pure solvent, we computed a value of K = 3.98 and n = 1, which is in excellent agreement with theoretical predictions^32,41 and viscosity measurements from our Poiseuille flow simulations. In the case of the polymer solutions, we found K = 4.80 and n = 0.98 for the shortest investigated polymers (M = 10) and K = 4.93 and n = 0.97 for the longest ones (M = 80). As expected, the addition of polymers leads to an increase of the overall viscosity, where longer polymers have a slightly larger effect. Furthermore, we can conclude from these data that the polymer solutions exhibit no significant shear thinning in the investigated $γ ̇$ regime.

Hence, our simulations suggest that a neutral surface can in principle exist for viscoelastic solvents with constant shear viscosity. Such a neutral surface exists in an intermediate flow regime where chains are somewhat shear-aligned (based on the trends in Fig. 3), but the flow rate is not high enough for significant wall repulsions. This finding is in contrast to previous predictions from CFD calculations³⁰ and might be due to microscopic details of the polymer-colloid interactions which cannot be captured in mean-field treatments of the polymer solution.

To investigate the effect of flow rate on focusing, we took the M = 40 polymers that focused for a = 16 and Re = 32, and systematically decreased the flow rate to essentially non-inertial conditions (Re = 1). Figure 5 shows the effective force F_x as a function of channel position for various flow rates; it is immediately apparent that the driving force F_x steadily decreases with decreasing flow rate. For Re = 1, there is essentially no focusing force near the channel centerline, but instead there is a noticeable attractive force close to the channel walls. We surmise that this effective wall attraction is the result of polymer depletion effects, akin to the situation at equilibrium: the conformational degrees of freedom of the polymers are significantly restricted in the immediate vicinity of the walls and therefore prefer to move to the channel center. This purely entropic effect leads to the displacement of the colloids from the channel center, which then manifests itself in the observed attraction.⁴⁹ 

Varying the flow rate affects the forces acting on the colloid in two ways. First, the outward hydrodynamic lift force F_L on the colloid decreases with decreasing flow rate. For a particle Reynolds number Re_p = Re a/L_x ≪ 1, there is essentially no hydrodynamic lift force and the colloid is confined on a streamline,⁵⁰ whereas in the inertial regime, F_L increases with the square of the flow rate near the channel centerline (see Eq. (1)). Second, while the outward lift force is weakening, the magnitude of the inward elastic force is also decreasing at the same time. Macroscopic scaling arguments (see Eq. (4)) suggest that the elastic force should also be proportional to the square of the flow rate due to the dependence on the first normal stress gradient. From a microscopic point of view, sufficient shear rates are required to deform and shear align the polymers in order to create this normal stress difference. The threshold shear rate $γ ̇ c$ is characterized by the dimensionless Weissenberg number $Wi= γ ̇ c τ R$⁠, where shear alignment is expected when Wi ≳ 1.

For M = 40, we estimate Wi < 1 for Re = 1 and Wi ≥ 1 for all other investigated flow rates, based on the average shear rate in the channel and a relaxation time estimated within the Zimm model using the equilibrium radius of gyration.^44,51 Hence, the M = 40 polymers should be isotropic at Re = 1 and stretched otherwise. We confirmed this trend quantitatively by evaluating the time averaged gyration tensor for the polymers under flow. We believe that because our M = 40 polymer chains do not stretch at Re = 1, no viscoelastic focusing is observed. For higher Re where the chains are stretched, F_x still does not show a strong dependence on Re, and we conclude that the leading scaling terms for F_L and F_E roughly cancel each other out for this state point. The weak Re dependence observed in Fig. 5 can then be attributed to prefactors, lower order terms, and the increased role of wall hydrodynamic repulsions.

We also computed F_x for varying colloid diameters at fixed flow rate Re = 32 and polymer length M = 40, shown in Fig. 6. F_x depends only weakly on a near the centerline, with a stronger inward force for larger particles, suggesting that near the centerline the leading scaling terms of F_L and F_E (both predicted to be proportional to a³) have similar prefactors for this choice of parameters. The biggest effect on colloid size appears close to the walls (ξ ≳ 0.4), where the apparent repulsion is much stronger for the largest colloid than for the smaller ones. This behavior is mainly due to the fact that the larger colloid has a much smaller wall separation compared to the smaller ones at a fixed position: for a = 8, wall contact does not occur until ξ = 0.84, whereas the largest colloid is already touching the wall at ξ = 0.68. Di Carlo et al. have shown that the strength of the hydrodynamic wall repulsion scales as ∼a⁶,²² which is qualitatively consistent with our simulations (i.e., at a fixed position, larger colloids have considerably stronger repulsions).

In practice, it is not only important to know whether a colloidal particle will focus but also to characterize the quality of focusing. Important characteristics are, for example, the speed with which focusing can be achieved, the width of the focused stream, and its stability. Due to the inherent Brownian motion of colloidal particles, there will be a distribution of particle positions at steady state, which intuitively should be more narrow for stronger focusing forces. Such particle distributions have been measured experimentally to qualitatively characterize focusing of colloidal particles.²⁰ Prohm et al. characterized the particle distribution of colloidal particles for tube Poiseuille flow in a Newtonian solvent.¹⁵ Here, we apply their methods to characterize particle distributions for colloidal particles in dilute polymer solutions.

The motion of a colloid in the channel gradient direction can be described by a one-dimensional random walker in an external potential. Proceeding as in Ref. 15, the Langevin equation describing the position of such a walker is

where γ = 6πηa is the Stokes drag coefficient and ζ(t) is white noise with zero mean and variance satisfying the fluctuation-dissipation relation 〈ζ(t) ζ(t′)〉 = 2γk_BTδ(t − t′). The steady state distribution function f(ξ) for such a random walker is a standard Boltzmann distribution

where U(ξ) is the potential of mean force obtained by integrating the calculated forces

with U(ξ = 0) = 0. In a recent publication,²¹ De Santo et al. proposed to model the forces exerted by the viscoelastic medium as an effective harmonic potential by considerations of the first normal stress on a Brownian particle. However, their model is restricted to the inertialess limit (Re ∼ 10⁻⁸) and is therefore not applicable to our simulations where inertia has a significant contribution to the total effective force.

In order to compute the particle distribution directly, we repeated our simulations starting from the same initial configurations as the previous force calculations in Fig. 2 but this time allowed the colloids to move freely in all dimensions. Due to the symmetry of the channel geometry, we consider only the probability of finding the colloid at an absolute displacement |ξ| to improve sampling.

Figure 7 shows the probability distributions calculated from our simulations along with the expected distributions based on Eq. (9). The inset of Fig. 7 shows the potentials corresponding to the forces calculated in Fig. 2. As expected, the simulated particle distributions have a maximum around the minimum of the potential. Brownian particles can fluctuate within a few k_BT of the minimum, and so for the Newtonian solvent, we expect and obtain a peak between 0.4 ≲ ξ ≲ 0.5. The distribution for the Newtonian solvent is not perfectly symmetric but slightly skewed towards the center of the channel because of the strong hydrodynamic repulsion near the wall.

In the case of the polymer solutions, we can discern a sizable discrepancy between the simulated particle distribution and the distribution predicted by Eq. (9), which becomes more pronounced as the polymer chain length is increased. For the shortest polymers (M = 10), the two distributions are still in excellent agreement, whereas for the longer chain lengths, the simulated distributions are much broader than the calculated potentials suggest.

We believe that this discrepancy originates from the method of calculating the force on the colloid and the microscopic relationship between the colloid and nearby polymers: in order to compute the force, we fixed the lateral position of the colloid so that it is unable to move in the gradient direction regardless of the nearby polymer solvent environment. This allows us to sample the average potential of mean force experienced by the colloid. However, when the colloid is allowed to freely move, it may experience and also induce fluctuations in the nearby polymer environment, e.g., gaps where there are no polymers. We hypothesize that these additional degrees of freedom could lead to colloid trajectories that are inconsistent with the average potential of mean force computed for colloids at fixed position.

This line of reasoning suggests that the discrepancy in the predicted distribution should be the largest near the channel centerline where small movements of the polymers can induce dramatic changes in their conformations. When a polymer is perturbed from the channel centerline (where there is no shear) to a region of finite shear rate, the previously coiled polymer may suddenly stretch if the shear rate exceeds the threshold rate $γ ̇ c$⁠. This effect is more pronounced for longer polymers because they deform at smaller $γ ̇ c$⁠. Away from the center of the channel, the polymer conformation is less sensitive to position because the chains are already stretched. We confirmed this trend from the gyration tensor of the polymers, for which we found a significantly stronger position dependence near the channel centerline, especially for longer chains (M = 40 and M = 80).

Figure 8 demonstrates this effect for single colloid trajectories in polymer solutions with chain lengths of M = 10 and M = 80, respectively. In both systems, the colloid positions fluctuate about the expected distribution peak. However, the fluctuations for the M = 10 case are much smaller than for the M = 80 case, although the opposite is predicted by Eq. (9). Moreover, for M = 80, several trajectories exhibit rare but noticeable veerings away from the centerline position, before returning to the middle in a similar fashion (e.g., between t = 60 000 and t = 80 000). We believe that these events correspond to fluctuations in the calculated average potential of mean force, which are not captured appropriately in F_x.

The fluctuation of the local polymer environment in response to the motion of the colloid suggests that the colloid and polymer solution cannot be decoupled from each other. At these length scales, a mean-field treatment of the polymer as a density or stress field is no longer adequate, and a microscopic description of the polymers is necessary. Although the calculated forces do not quantitatively reproduce the exact particle distributions, we emphasize that our measured force curves still predict the qualitatively correct behavior: polymer chains with M = 80 provide better focusing than M = 40, whereas chains with M ≤ 20 do not focus at all. Furthermore, we expect that many body effects, such as hydrodynamic wakes, may play an important role when focusing systems with multiple colloids. These additional terms may have a sizable effect in applications with multiple particles.

We studied the viscoelastic focusing of rigid colloids in slit-like channels via a hybrid molecular dynamics approach that faithfully takes into account hydrodynamic interactions. A bead-spring model was used to model the constituent polymers of the viscoelastic medium in a way consistent with the colloidal nature of the solute particles. We systematically measured the forces exerted on the colloid through the surrounding solvent and decomposed them into inertial and elastic contributions to better understand the driving force behind the cross-stream migration of colloidal particles at high Reynolds numbers Re ≥ 1.

For short polymer chains, the elastic forces are not sufficiently strong to outweigh the hydrodynamic lift forces at the channel center, leading to unfocused states. However, as the polymer length increases, elastic contributions become dominant, resulting in the focusing of the colloids. Further, we found that the total forces become stronger with increasing colloid size and flow rate. However, these parameters exhibited a significantly weaker scaling compared to the predictions within a continuum approximation of the viscoelastic medium. In part, this discrepancy stems from microscopic details of the polymer solution, which become increasingly important at these length scales and which manifest as wall attractions and pronounced fluctuations in the colloid trajectories. These microscopic effects cannot be captured in mean-field descriptions of viscoelastic focusing and potentially need to be taken into account when designing and conducting microfluidic experiments on the colloidal scale.

We acknowledge financial support for this work from the Princeton Center for Complex Materials (PCCM), a U.S. National Science Foundation Materials Research Science and Engineering Center (Grant Nos. DMR-0819860 and DMR-1420541). In addition, M.P.H. received government support under Contract No. FA9550-11-C-0028 and awarded by the Department of Defense, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a.