Hydrodynamic irreversibility of non-Brownian suspensions in highly confined duct flow

Abstract The irreversible behaviour of a highly confined non-Brownian suspension of spherical particles at low Reynolds number in a Newtonian fluid is studied experimentally and numerically. In the experiment, the suspension is confined in a thin rectangular channel that prevents complete particle overlap in the narrow dimension and is subjected to an oscillatory pressure-driven flow. In the small cross-sectional dimension, particles rapidly separate to the walls, whereas in the large dimension, features reminiscent of shear-induced migration in bulk suspensions are recovered. Furthermore, as a consequence of the channel geometry and the development and application of a single-camera particle tracking method, three-dimensional particle trajectories are obtained that allow us to directly associate relative particle proximity with the observed migration. Companion simulations of a steadily flowing suspension highly confined between parallel plates are conducted using the force coupling method, which also show rapid migration to the walls as well as other salient features observed in the experiment. While we consider relatively low volume fractions compared to most prior work in the area, we nevertheless observe significant and rapid migration, which we attribute to the high degree of confinement.


Introduction
Suspension flows of spherical particles in a Newtonian fluid pertain to a wide range of applications, including biomedical, industrial, and geographical.Such flows represent generally chaotic dynamical systems (Drazer et al. 2002), with the particles undergoing irreversible migration from an initial configuration to a final, steady state configuration.Mechanisms for particle or droplet migration across streamlines may be exploited for passive sorting in typical duct-like microfluidic channels (Marnoto & Hashmi 2023).
The role of particle contacts and surface roughness in non-Brownian suspensions is a topic of current interest as it relates to the particles' irreversible migration at low Reynolds in a shear flow.Particle surface roughness results in an inherent irreversibility in the particles' relative motion as the particles are displaced across streamlines, and also modifies the pair distribution function between particles (Rampall et al. 1997;Blanc et al. 2011;Ingber et al. 2006).It has been shown that neither the longrange hydrodynamic interactions (Metzger & Butler 2010) nor the lubrication forces (Metzger et al. 2013) contribute significantly to the irreversibility.However, the particle surface roughness plays a key role in such dynamics (Zhang et al. 2023;Pham et al. 2015;Corte et al. 2008), and can also directly impact the bulk rheological properties of non-Brownian suspensions (More & Ardekani 2020).Pine et al. (2005) discovered that there is a critical strain amplitude for a given volume fraction below which the particle locations are reversible in oscillating flow.Above this critical strain amplitude, chaos sets in.This critical strain amplitude is quite sensitive to the volume fraction, with irreversible behavior being more difficult to achieve with low volume fractions.Pham et al. (2016) found that smoother particles lead to a higher critical strain amplitude and predicted the critical strain as a function of the particle surface roughness, although very recent work has shown the role of particle roughness on irreversibility to be more subtle (Zhang et al. 2023).The critical strain amplitude has been further explored in flows with clouds of particles in an oscillating shear flow.The clouds are shown to extend until the volume fraction is less than the critical volume fraction for irreversibilty at a given strain amplitude (Metzger & Butler 2012;Howard & Maxey 2018).The vast majority of prior experimental and numerical works have focused only on the simplest possible flow geometries, either flow between parallel plates or circular pipe flow, and typically in scenarios where the particles are much smaller than the confining geometry.Taking a first step towards extending our understanding of particle migration behaviors and irreversibility to more general geometries with non-uniform shear in multiple directions and in scenarios where confinement effects are important are the primary goals of this work.
Information about suspension flows comes from simulations and experiments.Simulations are costly, and quantities such as the particle surface roughness are represented by approximations such as an interparticle contact barrier.While simulations have been shown to accurately capture experiments (Maxey 2017), there is still a need for highly resolved experiments that can track particle locations.A number of experiments have been performed (Snook et al. 2015;Guasto et al. 2010;Pham et al. 2016) that use a laser sheet to image a slice of the suspension.This technique requires the suspending fluid to be index-matched to the particles and complete trajectories are unavailable because particles are lost from view once they leave the laser plane.Laser-Doppler velocimetry (LDV) has also been used to achieve high-resolution measurements of particle velocity and concentration (Lyon & Leal 1998a,b) but the one-dimensional nature of LDV prohibits simultaneous capture of the entire flow.Magnetic resonance imaging (MRI) can provide 3D measurements of suspensions which are optically opaque and hence is well suited to suspensions in porous media (Mirbod & Shapley 2023).However, individual particles cannot be tracked due to the relatively low resolution of MRI.Single-camera refraction-based methods have been previously developed and applied to multiphase flows in the reconstruction of the three-dimensional shape of fluid interfaces (Moisy et al. 2009;Kilbride et al. 2023); however to the best of our knowledge, similar ideas have not previously been adapted to suspension flows.Although high-resolution 2D trajectory measurements have been performed previously to study irreversible dynamics in very small numbers of particles (Pham et al. 2015), we believe our work is the first to provide such experimental data in a 3D random suspension of many particles, directly enabled by our new tracking technique.
In this work, we experimentally study the migration of spherical particles suspended in a viscous fluid subjected to an oscillatory pressure driven duct flow.Extended 3D particle trajectories are resolved with a single camera by implementing a refraction-based imaging technique capable of capturing particle motion in the out-of-plane direction.Although relatively low volume fractions are considered in this work, we nevertheless observe significant migration in both dimensions of the cross-section.We report measurements  of the migration dynamics and steady-state particle distributions for different packing fractions and strain amplitudes, and quantitatively associate particle proximity during each cycle with the measured irreversibility.Our results are also directly compared to new numerical predictions of a highly confined suspension steadily flowing between parallel plates.

Experiments
The particles used in our experiments are optically clear acrylic (PMMA) spheres with diameter 1.598 ± 0.009 mm and root mean square surface roughness 0.08 ± 0.03 µm.They are suspended in a water-glycerol mixture with Newtonian viscosity 36.3 ± 1.8 cP (Cheng 2008) filling a 300 mm long Borosilicate glass test channel (VitroCom R0309, rectangular internal cross section 9 × 3 mm).The mixture ratio is chosen such that the particles are neutrally buoyant, each with a density of 1.179 ± 0.003 g/mL.The 3 mm channel height is just less than two particle diameters, preventing complete particle overlap.The 9 mm channel width was selected to be small enough so that a significant shear gradient was present across the entire x-direction, and wide enough to allow particles to move a few diameters away from the side walls.Experiments performed by Guasto et al. (2010) use the same channel aspect ratio, though in their case the particle size to channel depth ratio is far smaller and only a single central plane is imaged, with the intent of approximating a two-dimensional flow.The complete experimental setup is shown in figure 1 and dimensional and non-dimensional parameters summarized in table 1.The glass test channel is oriented vertically and attached to a laser cut stand by screw-adjustable mounts.3D-printed end caps with flow orifices ensure the particles remain in the channel and adapt the ends of the channel to hoses.The upper hose goes to a reservoir of suspending fluid which is open to the atmosphere and the bottom hose attaches to a syringe pump (Harvard Apparatus Pump 11 Elite).
A back light panel (3W LED tracing pad) illuminates the test section.We generate speckle patterns using the Rosta algorithm implementation (dot size 1, density 0.3 and zero smoothing) in the µDIC Python library (Olufsen et al. 2019) which are then ink-jet printed with 1200 x 2400 DPI resolution on transparency film.The pattern is fastened to the channel with a transparent acrylic clamping block.An Allied Vision Mako 503B monochrome CMOS camera with 105 mm micro Nikkor lens is used for imaging the suspension at a 1 Hz frame rate and 2592 × 610 resolution.The camera is placed at a working distance of approximately 1 m, resulting in a 40 mm by 10 mm imaged section.Hence lens distortion is minimized and the incident rays from the camera can be assumed to be parallel.The outlets of the channel are 130 mm away from the imaged section so entry length effects can be neglected.
We prepare particle suspensions with nominal bulk particle volume fractions φ ′ B = 0.07, 0.11 and 0.14.The corresponding nominal bulk area fractions φ ′ A -defined as the area of the experimental image that would be filled with particles assuming no overlap -are 0.2, 0.3 and 0.4.Similar bulk area fraction values have been previously used in simulations of a monolayer suspension (Nott & Brady 1994).A magnetic 3D-printed plow-shaped particle pusher is used to manually set the initial particle distribution in the channel by guiding it with a handheld magnet outside the channel.The geometry of the pusher is such that the particles are driven to the channel walls in the x direction and the channel center (y = h/2) in the y direction.However, we observe that the particles quickly migrate towards the walls along the y direction once the experiment starts.Since the manual particle arrangement process produces some heterogeneity in the particle distribution, the measured bulk volume fraction φ B in the field of view (averaged over all frames of the experiment) is slightly different than the nominal φ ′ B , and takes values ranging from 0.06 to 0.12.
During the experiment, the syringe pump generates oscillatory flow in a square wave, prescribing the volumetric flow rate Q to be where T is the flow oscillation period.The strain amplitude γ is the average displacement of a fluid parcel over a half cycle T /2 divided by the channel half width w/2.We perform experiments with γ = 6 up to an accumulated strain of 550 (46 oscillation cycles) and with γ = 1 to an accumulated strain of 1100 (550 oscillation cycles).The period T is chosen to maintain a particle Reynolds number of 0.016.Initial experiments performed at a Reynolds number of 0.008 produced indistinguishable suspension dynamics.Experiments were also performed with a single particle in the channel showing fully reversible motion in 3D, regardless of its initial position.
Prior to each experiment, a particle-free reference image of the speckle pattern is taken.Then, throughout the experiment, a series of photographs of the test section with particles are captured.Figure 2(a) outlines the image processing steps which include background subtraction (Zivkovic 2004), filtering with morphological operations, and determining 2D particle positions with the Euclidean distance transform and h-maxima filter (Vincent 1993).The out-of-plane y-component of each particle's position may be inferred by comparing the reference and particle images since the speckle pattern in the particle image will be distorted depending on the particles' positions.Figure 2(b) illustrates the working principle of this particle tracking technique.The incident rays from the camera refract at the fluid-particle interfaces and focus as they reach the pattern attached to the outside of the channel.As seen in the photographs corresponding to each particle height, increasing the y component of the particle position has the approximate effect of increasing the magnification of the pattern.By applying Snell's law at each optical interface in this axisymmetric geometry (Hecht 2012), an analytical expression may be constructed which predicts the observed distortion of the reference pattern through the transparent particle.For the case of partially overlapping particles, the distortion is determined numerically by recursively propagating each incident ray from the camera through the cluster of particles until it reaches the speckle pattern.The particle's position is then determined by solving the inverse problem illustrated in 2(c): given a reference image I r of the pattern and the resulting particle image I p , we find a position r for the particle that reproduces the observed distortion by maximizing the cross correlation between regions J and D p .We use the Nelder-Mead algorithm (Nelder & Mead 1965) in the NLopt library (Johnson 2021) to optimize the cross correlation.The OpenCV library (Bradski 2000) for C++ is used to implement the image transformations and cross correlations.The tracking method requires calibration in order to precisely determine the refraction indices of the particles, fluid and channel walls.A calibration device was constructed by attaching a particle on a thin wire to a high precision translation stage.The particle may then be moved to known y positions across the height of the channel while capturing images of the distorted speckle pattern.Using the algorithms described, the particle positions are measured from each image and the refraction indices are adjusted in order to minimize the sum of squared residuals between the actual particle position and the measured position.After calibration, the translation stage was again used to independently determine the accuracy of the tracking method.Using a number of different particles at various locations in the field of view, the comparison between the measured height and the actual height reveals excellent agreement (figure 2(d)).Over the full measurement range, the particle's y position can be determined to within 157 µm or 5.23% of the channel height with 95% confidence.The final processing step is to link the 3D particle positions into trajectories using Crocker and Grier's algorithm (Crocker & Grier 1996) in the open-source library TrackPy (Allan et al. 2021) so that each particle's identity is known throughout multiple frames, with the result shown in the "tracking" panel of figure 2(a).

Simulations
To compare with the experiments, we also completed a limited set of simulations using the Force Coupling Method (FCM) (Yeo & Maxey 2010, 2011).Due to computational limitations, the simulations are completed with a single set of walls located at y = 0 and y = 4a representing a very narrow channel (corresponding to a channel height of h = 4a), and are periodic in the streamwise (z) and spanwise (x) directions.The simulations are scaled by the particle radius a = 1.To set up the simulations, the particles are seeded randomly in the channel to reach the proper area fraction.The particles are first seeded with radius 0.85a, and then are "inflated" to reach radius a = 1 with a molecular dynamics simulation.In doing this work, we experimented with different methods of seeding the particles to match the experimental results.With simulations, it is possible to seed the particles in a perfect monolayer, located at the channel half height y = 2a.However, this does not match the experiments, because the simulation particles are not subject to noise in their location, and therefore would remain in a monolayer for an unphysically long time.Instead, we seed the particles in a "pseudo-monolayer," or a noisy monolayer.The particles are seeded with average y location at y = 2a, and a standard deviation of 0.05a.The minimum y value for φ B = 0.14 is 1.833a, and the maximum y value is 2.221a at the beginning of the simulation.This way, the particles begin close to the center of the channel, but particle interactions allow for the particles to migrate in the cross-stream (y) direction.The channel width is w = 30a and the channel length is l = 60a, which results in a total of 115 particles for φ B = 0.07, 172 particles for φ B = 0.11, and 230 particles for φ B = 0.14.
The flow is driven by a pressure gradient ∂p/∂z = 0.5.Because the strain amplitude in the experiments are calculated using the x−direction length scales, we calculate the accumulated strain as γ a = 2 u z t/(11a), where the scaling factor of 11 is chosen to approximate the experimental strain amplitude scaling factor of w.
Figure 3. Experimental results for migration in the cross-stream (y) direction.(a) Average distance between a particle center and the nearest wall versus accumulated strain for experiments with γ = 6 and φB = 0.06 (red), 0.10 (green) and 0.12 (blue).The error bars show the standard deviation between 5 trials.(b) Initial and steady state particle position distributions over the height of the channel show that, in the narrow channel dimension, the particles prefer the walls at steady state.The initial distribution is averaged over the first 20 frames (γa < 2.74) and the steady profile over frames 601 to 656 (γa = [82, 90]).The vertical dashed lines indicate a distance of one particle radius from the wall.
The simulations use a time step of ∆t = 0.005.To approximate the particle surface roughness, we use a short range contact barrier that disrupts the symmetry in particle interactions and computationally prevents the particles from overlapping.The contact force between particles α and β, with centers r α and r β , acts along the line of centers of the particles and is given by (Yeo & Maxey 2010, 2011).The relationship between the contact force reference values and the pressure gradient driving the flow sets the minimum separation between particles.V ref = 200 is a constant chosen to so that the minimum gap between two particles is approximately 0.005a.Particles separated by a distance of less than 0.01a are considered to be in contact, so R ref = 2.01a.The contact force between a particle α with center r α and a wall is given by: where r is the vector between r α and the top or bottom wall.We take R ref,wall = 1.05a and V ref,wall = 200.

Migration to the walls: cross-stream y-direction
Although we bias the particles to the channel center (y = h/2) in the narrow dimension at the start of the experiment, we observe that the particles rapidly migrate towards the walls once the suspension is subjected to a periodic strain.The rate of this migration is captured in figure 3(a) which shows the average distance between a particle center and the nearest wall as a function of accumulated strain γ a .On a timescale which is approximately an order of magnitude faster than the less confined x-direction evolution we will present in section 4.2, the particles reach a configuration in which the average particle is about 1.25 radii away from the channel wall.The corresponding FCM simulations produce a similar result in figure 4(a) where the migration to the walls is even more apparent due to the near-monolayer initial condition.The rate of migration is enhanced for higher φ B .For the simulations, the final average y-locations for the two highest values of φ B are about 1.25a, the same as for the experiments.The initial and steady state particle position distributions in the experiments and simulations are presented in figures 3(b) and 4(b), where we plot the probability that a given particle's center will lie in one of 20 evenly spaced bins distributed over the channel height.Both the experimental and simulation plots show an initial distribution that favors the channel center due to the imposed particle configuration.However, the particles strongly prefer the walls in the steady state distribution.The data points which appear closer than one particle radius to the wall in the experiments are a consequence of the y-direction tracking uncertainty, which is quantified in figure 2(d).The experimental measurements of Snook et al. (2015) in circular pipe flow show the presence of local peaks in the radial concentration distribution at the wall which they attribute to a particle wall-layering effect.This becomes more pronounced in their data for the smaller pipe radius considered, an effect they attribute to the higher confinement.In the case of extreme y-direction confinement studied here, such wall layers dominate and ultimately overwhelm the final concentration distribution.As also noted by Snook et al. (2015), continuum models such as the Suspension Balance Model are of course unable to capture such discrete particle effects.
It is possible to understand this migration to the walls by considering two-body interactions between the particles.If two particles come into contact, the force between the particles acts along the line of centers of the particles.Unless the particles have identical y−locations, the force will have a component that displaces each particle towards the walls.After many strain units and particle interactions, the particles will migrate to be close to the walls.Once a particle enters the wall layer, that is, it is on top of the wall, it is very difficult for the particle to exit because any additional interaction will force the particle closer to the wall.

Migration to the center: spanwise x-direction
Examining the evolution of the local particle volume fraction φ as a function of accumulated strain γ a gives insight into the dynamics of the migration process and its dependence on bulk packing fraction.The average particle concentrations in the center and outer quarter-width bins when γ = 6 are plotted against the accumulated strain in figure 5.The bin concentration φ bin is computed based on the 3D particle positions as the volume in the bin filled with particles divided by the total bin volume.When a particle intersects the x = w/4 or x = 3w/4 planes, its volume is partitioned accordingly between the inner and outer bins.The results in these and all subsequent experimental plots are averaged over at least 5 trials, with error bars representing one standard deviation.Due to the initial configuration of the particles, the initial outer bin concentration exceeds the bulk packing fraction for all experiments.However, for the φ B = 0.10 and φ B = 0.12 experiments, the inner and outer bin concentrations eventually cross and, at high accumulated strain, the inner bin concentration exceeds the bulk packing due to particle migration.Migration proceeds at an increased rate for φ B = 0.12.For the low packing fraction experiments with φ B = 0.06, the bin concentrations do not cross the average channel concentration up to an accumulated strain of 550, indicating a sharp reduction in the rate of migration.The oscillations in the inner bin concentration in figure 5(a) match the period of the flow and occur due to small variations in the packing fraction along the length of the channel."Anomalous" migration was not observed in the present experiments, but is predicted to occur for significantly smaller strain amplitudes than explored in the present work (Morris 2001).
The initial and steady state particle concentration profiles over the channel width are shown in figure 6(a) for experiments with γ = 6.Since the experiments begin in a wall-loaded configuration, the initial concentration profiles contain two peaks that are approximately one particle radius away from the channel walls.Additional experiments were conducted for more homogeneous initial concentrations and showed very similar results (Appendix A).For the φ B = 0.10 and φ B = 0.12 experiments, the steady state concentration profiles indicate significant migration towards the channel center since the bulk packing fraction is surpassed there.However, compared against results with a parabolic flow profile in Snook et al. (2015) and Lyon & Leal (1998a), the enhanced concentration at the channel center is less pronounced, mirroring the flatter velocity profile in the current geometry.For φ B = 0.06, the particle arrangement relaxes from its initial configuration suggesting some irreversible behavior but there is less preferential migration towards the center of the channel.The x-direction concentration profiles are constructed by computing the local y-and-z-averaged particle volume concentration φ yz in a bin which is swept over the width of the channel as shown in figure 6(b).The bin width is w/100 and the volume concentration φ yz is computed as the volume of the bin filled with particles divided by the total bin volume.For clarity, error bars are shown at discrete points rather than at every bin station.While we do observe a net migration of particles to regions of lower shear in the channel center (in the x-direction), we note that due to the relatively large particle to channel size (i.e.high confinement) it is unlikely that our results would be well described by continuum models of similar behaviors such as the Suspension Balance Model (Morris & Boulay 1999;Guazzelli & Pouliquen 2018).In contrast to the present work, prior experimental work on shear-induced migration in less confined suspensions in circular pipes have suggested an absence of migration for suspensions around φ B ≈ 10% in both oscillatory (Snook et al. 2015) and steady flows (Hampton et al. 1997).The fact that detectable migration towards the center is observed for the relatively low volume fractions in our system is thus likely a consequence of the high levels of confinement, which naturally increase the likelihood of particle contacts.Such contacts are known to be the primary source of irreversibility in these systems (Metzger et al. 2013).

Migration rate
We quantify the rate of migration in our experiments by measuring the number of oscillations required for the particle concentration profile to reach steady state following a method similar to the one used in Snook et al. (2015).As seen by comparing figures 3(a) and 5, the time to steady state in our experiments is dominated by the migration in the x-direction.As such, in order to determine steady state, we compute φ yz using bins of width w/20 for each oscillation of the experiment, averaging over frames 22 to 76 of the oscillation.When a majority of the bins are within one standard deviation of the particle concentration over the remaining oscillations, the suspension is considered to have reached steady state.Steady state is reached in 29 ± 14, 11 ± 7 and 4 ± 1 oscillations for experiments with φ B = 0.06, 0.10, and 0.12, respectively.The values following the means characterize the distribution in the time to steady state between trials.This result matches the observation in Snook et al. (2015) that the mean accumulated strain to reach steady state scales approximately as φ −2 B .Despite the lower volume fractions considered in our work, the mean number of oscillations to steady state are comparable in magnitude to the numerical values reported by Snook et al. (2015).Furthermore their data indicates a decrease in time to steady state with both volume fraction and degree of confinement.It is thus likely that the increased confinement in our system promotes more particle contacts than in an otherwise equivalent lesser confined scenario, which compensates for the reduced likelihood of contact associated with smaller volume fractions.

Particle velocities
Enabled by the 3D particle tracking, figure 6(d) shows the average magnitude of the particle velocity over the channel cross section.The experimental velocity profiles are time-averaged over a full experiment and the theoretical velocity is calculated for a pure Newtonian fluid in Stokes flow.The white cells in the experimental profiles are regions where no velocity measurement could be obtained since the particle centers must be at least one particle radius away from the channel walls.For the φ B = 0.06 case, the magnitude of the measured velocity profile is close to the theoretical profile but the particle velocities are reduced for the higher packing cases.The mid-line particle velocity profiles at y = h/2 are compared in figure 6(c) against the theoretical center-line velocity.The experimental velocity profiles are constructed using particles whose centers are within 0.12h of the channel center.The profile becomes increasingly blunt with higher φ B , consistent with results in a circular pipe (Snook et al. 2015) and two-dimensional channel flow (Lyon & Leal 1998a;Guasto et al. 2010).

Particle irreversibility and interactions
The mean square particle displacements (MSD) are a quantitative measure of irreversibility in the particles' motion and thus allow more detailed comparison between the experiments.The non-dimensional effective particle diffusivity D x is defined following Pine et al. (2005) according to where (∆x/(2a)) 2 is the non-dimensional MSD, 2a is the particle diameter, γ a is the total accumulated strain and ∆x is the displacement of the particle over one cycle.Given the three-dimensional nature of the flow geometry, effective diffusivities D y and D z for the other directions are defined in a similar manner.In order to ensure that the results are statistically significant, MSDs are computed up to the maximum accumulated strain for which there are at least 50 continuous particle trajectories (starting at any point in the   experiment).Since higher bulk packing fraction increases the likelihood that a particle is lost during tracking, MSDs are calculated to a lower accumulated strain for higher φ B experiments.The MSDs after an integral number of cycles in the γ = 6 experiments are graphed against the accumulated strain in figure 7(a).The y-component is omitted as it does not exhibit diffusive behavior, presumably due to the high confinement.The least squares fits for the effective particle diffusivities are shown as well.For a given packing fraction, the effective diffusivity along the flow direction D z is significantly higher than the spanwise diffusivity D x .This result is consistent with findings in a circular Couette flow where two or more components of diffusivity were measured (Pine et al. 2005;Breedveld et al. 2001Breedveld et al. , 2002)).Guasto et al. (2010) also report that the streamwise diffusivity is consistently higher than the spanwise diffusivity in a rectangular channel flow.Figure 7(b) compares the effective diffusivities as a function of φ B .For every value of φ B , the ratio between D x and D z is about 20 but their values increase approximately one order of magnitude as φ B increases from 0.06 to 0.12.The increase of effective diffusivity with bulk volume fraction is consistent with experimental results in the literature for other flow geometries (Pine et al. 2005;Guasto et al. 2010).We note that our experiments involve non-uniform strain and time-evolving non-uniform particle concentrations.As such, the particle diffusivities will depend on space (Guasto et al. 2010) and timedependencies that a single particle diffusivity value does not capture.Even in the case of a simple shear flow, it has been observed experimentally that the particle diffusivity can have distinct short and long-time values (Breedveld et al. 2001).Nevertheless, the average diffusivities computed here remain a reasonable quantitative proxy for the general level of chaos in the system (Drazer et al. 2002), and allow for a simple comparison between our experimental conditions.Direct quantitative comparison of these numerical values with simpler cases involving uniform strain and particle concentration such as Pine et al. (2005) should be made with care, however.With significantly more data, the spatial dependence and temporal evolution of the particle diffusivities could be faithfully resolved using our experimental setup.
Since the 3D position of each particle in the experiment is known, interactions between particles and the resulting irreversibility may be directly examined.Figure 8 shows the number of particles that migrate a given amount based on their proximity to their neighbors over a cycle.A particle's level of interaction is quantified along the abscissa by d min /(2a) which represents the average 3D surface-to-surface distance between a particle and its nearest neighbor over a cycle, normalized by the particle diameter.The nearest neighbor is chosen at each frame in the cycle and the data is collected over the entire experiment.While other metrics for particle contact pressure such as the radial distribution function have been used in the literature (Rintoul & Torquato 1996), we find that our metric d min /(2a) is more strongly correlated to the irreversible particle displacements in our experiments.The magnitude of the spanwise displacement over a single cycle |∆x| is plotted on the ordinate, and is a measure of the irreversibility of particle motion during the cycle.These plots demonstrate that for all experiments, the particles with more interactions (smaller d min ) have a greater tendency to migrate (larger |∆x|).The presence of close interactions does not guarantee migration however, but it is statistically more likely to occur.As the packing fraction increases, so does the number of particles which interact closely with their neighbors resulting in larger migratory excursions.For the φ B = 0.06 experiments, there are few particles with d min /(2a) < 0.2 compared to the φ B = 0.10 or φ B = 0.12 experiments and consequently little migration towards the channel center, which is corroborated in the concentration profile results.We observe no appreciable difference when the heat maps are generated using only data from before or after steady state is reached.However, the heat map results do exhibit some spatial dependence which is explored in Appendix B. Example videos of particle behavior at specific points on the heat map are provided in the supplementary materials.Heat maps are also presented for the corresponding FCM simulations in figure 9. Since the simulations feature steady rather than periodic flow, the particle migrations are measured in 12 strain unit increments so that the average particle z-displacement between snapshots is the same as in the experiments.Like in the experiments, particles with more interactions are more likely to migrate and the particle proximity tends to increase with packing fraction.

Influence of strain amplitude
Experiments are also performed at φ B = 0.10 (φ ′ B = 0.11) and γ = 1 in order to observe the influence of strain amplitude on the migration dynamics.The average particle concentrations in the center and outer quarter-width bins are plotted against the accumulated strain in figure 10(a) up to an accumulated strain of 1100.Despite the large accumulated strain, the bin concentrations do not cross and the experiments exhibit quasi-reversible behavior where the bin concentrations change very slowly, potentially due to the locally high volume fraction at the walls.Furthermore, the effects of residual buoyancy (due to distribution in particle density) also lead to irreversible behaviors that are more apparent on long timescales.The initial and steady state concentration profiles in figure 10(b) feature a similar lack of migration which is consistent with observations in the literature (Pham et al. 2016;Guasto et al. 2010) of a critical strain amplitude below which the particles can subtly reorganize into a configuration that is essentially reversible.For the migration heat map in figure 10(c), the peak intensity occurs near d min /(2a) = 0 yet almost no migration occurs.Since the suspension is quasi-reversible at this low strain amplitude, the particles remain near their initial tightly packed wallloaded configuration.Measurements of the mean squared particle displacements after integral numbers of cycles yield values that are at least an order of magnitude smaller than those reported in figure 7(a).However, given the particle tracking resolution, we cannot make a statistically significant estimate for the effective diffusivity without substantially more data.

Conclusions
In our experiments, suspensions of neutrally buoyant spherical particles at bulk volume fractions φ B from 0.06 to 0.12 are subjected to oscillatory pipe flow in a confined rectangular channel and the resulting particle migration is characterized through direct imaging.We study a channel of aspect ratio 3:1 (x:y) with particles that are just larger than half the channel height in y to prevent complete overlap.While numerous qualitative comparisons to the existing state-of-the-art measurements have been made throughout, quantitative comparisons are far more challenging, as it is not straight forward to unambiguously map between highly disparate flow geometries.
In experiments performed at strain amplitude γ = 6, the particles are observed to migrate preferentially in the spanwise x-direction towards the channel center, with higher bulk volume fraction hastening the dynamics.These findings are consistent with results in the literature for other, unconfined flow geometries exhibiting shear-induced migration (Snook et al. 2015;Lyon & Leal 1998a;Guasto et al. 2010), although the enhanced concentration at the center of the channel is less pronounced in our case.Other measurements of the collective particle behavior such as local particle concentration and velocity are similarly consistent with previous results (Lyon & Leal 1998a;Guasto et al. 2010;Snook et al. 2015).The relative particle concentration at the channel center increases as the bulk volume fraction is increased.Conversely, the velocity profile is increasingly blunted compared to the theoretical φ B = 0 case as φ B is increased.Measurements of particle trajectories show that individual particles exhibit diffusive-like behaviors, with streamwise diffusivities about an order of magnitude larger than their spanwise counterparts, again consistent with prior work (Pine et al. 2005;Guasto et al. 2010).Additional experiments at a lower strain amplitude γ = 1 produce quasi-reversible particle behavior consistent with the existence of a critical strain amplitude below which shear migration does not occur (Pham et al. 2016;Guasto et al. 2010).In the small crossstream (y) cross-sectional dimension, particles rapidly migrate to the walls in all cases and reach a steady-state distribution more rapidly than in the x direction.This rapid segregation is also observed in companion Force Coupling Method simulations where the spheres are not confined in the x-direction.
Furthermore, other comparisons to the literature indicate that confinement plays an impactful role in our system.In particular, previous results indicate a lack of migration for low volume fractions around 10% (Hampton et al. 1997;Snook et al. 2015), whereas we find clear evidence for migratory behavior at these low densities.Additionally, the time to steady state in our system is comparable to that of more dense (but less confined) suspensions (Snook et al. 2015).Both of these observations suggest that confinement tends to enhance irreversibility.Geometrical confinement restricts possible particle motions, and thus naturally provides a mechanism for increasing the number of particle contacts when sheared.Lastly, for the case of extreme confinement (i.e. in the small cross-stream direction of our channel), the wall layers observed in prior work (Snook et al. 2015) entirely dominate the resultant distribution.Future experiments that more continuously span from high to low levels of confinement using a single flow geometry may provide additional valuable insight into the role of confinement on non-uniformly sheared non-Brownian suspensions.
In order to conduct the experimental measurements, we also implement a new, singlecamera 3D particle tracking method.The transparent particles are imaged above a speckle pattern and the magnitude of the resulting optical distortions are used to infer the out-of-plane particle positions with high accuracy.This method trades the high optical component cost of typical particle suspension setups for computational complexity, utilizing a refraction model to predict the distortion of the speckle patterns and iterating the particle position to reproduce the observed distortion.
The 3D tracking method permits the characterization of not only collective, but also individual particle behavior, enabling direct observations of interparticle spacing and trajectories over long times.Consequently, it is possible to build a quantitative connection between particle interactions and displacement over a cycle to explain the source of irreversibility.In both experiments and simulation it is observed that, as bulk volume fraction increases, a greater number of particles experience close interactions with their neighbors and these same particles are responsible for the strongest individual migration events.While similar measurements have been performed for very few numbers of particles in 2D (Pham et al. 2015), this work represents the first realization of similar temporally resolved, high-resolution trajectories for suspensions of large numbers of particles in 3D.
In addition to allowing measurements of otherwise inaccessible quantities in the monodisperse case, our 3D tracking method has the advantage of convenient extension to the case of bi-disperse or poly-disperse suspensions as particle size can be easily determined during the optical reconstruction.Due to the presence of the channel walls which prevent particle overlap, the present experimental findings are also relevant to natural problems involving particle suspension in high confinement such as thin films (Hooshanginejad et al. 2019) or porous media (Mirbod & Shapley 2023).Furthermore, as a consequence of the standard manufacturing methods, duct-like channels are extremely common in microfluidics.It has previously been shown that inertial effects in such geometries lead to focusing of suspensions in certain regions of the channel cross-section and can be exploited for passive particle manipulation (Di Carlo 2009).Irreversible particle-particle interactions in duct-like geometries can also result in particle migration across streamlines, even in the absence of inertia, with confinement hastening the process.The initial concentration profile is averaged over the first 20 frames (γa < 2.74) and the steady profile over frames 3500 to 3700 (γa = [480, 508]).The shaded region indicates the measured bulk packing fraction φB.(c) Initial (dotted) and steady (solid) distributions of particle centers in the y-direction.The likelihood that a given particle's center will lie in each region of the channel is plotted.The initial distribution is averaged over the first 20 frames (γa < 2.74) and the steady profile over frames 601 to 656 (γa = [82, 90]).The vertical dashed lines indicate a distance of one particle radius from the wall.loaded initial configuration, as can be seen by comparing with figures 3(b) and 6(a).An approximately equal number of particles start in the inner and outer bins in this case, with particles preferentially migrating to the inner bin as time progresses.Similarly, the bin concentrations recover a state close to the equilibrium values in figure 5  the wall.This results in a consistent and stable monolayer forming on the channel walls, with very few particles leaving over long times (Howard 2018).

Figure 1 .
Figure 1.(a) A diagram of the channel assembly, along with a magnified view of the rectangular glass channel and speckle pattern positioning.Channel cross-section is 3 mm (y-direction) × 9 mm (x-direction).(b) A photograph of the complete experimental setup.

Figure 2 .
Figure 2. (a) Example images of the channel test section at each step of the processing procedure.The tails in the "tracking" panel illustrate how each particle's position has evolved from the initial state in the previous panels, with the coloring indicating the y position over time.(b) Cross-sectional diagram of the channel labelled with relevant dimensions.Incident rays from the camera at p are propagated based on the refraction model to R(p; r) for particles at different heights.Experimental images of particles over a regular 200 µm grid pattern at the different heights shown below.(c) The reference image Ir is transformed based on the refraction model and the guessed particle position r.The result in J is compared to the corresponding region Dp in the particle image Ip using the cross correlation.(d) Histogram of particle height measurement errors from 120 data points taken at 12 evenly spaced known positions in y.

Figure 4 .
Figure4.Simulation results for migration in the cross-stream (y) direction.(a) Average distance between a particle center and the nearest wall versus accumulated strain for simulations with steady flow and φB = 0.07 (red), 0.11 (green) and 0.14 (blue).The particles begin in a "pseudo-monolayer" but quickly migrate to the walls in the narrow channel dimension.(b) Initial and steady state particle position distributions over the height of the channel.The initial distribution is averaged over the first 20 frames (γa < 3) and the steady profile over frames 550 to 600(γa = [82, 90]).The vertical dashed lines indicate a distance of one particle radius from the wall.

Figure 5 .
Figure 5. Experimental results for migration in the spanwise (x) direction.The evolution of particle concentration in the inner and outer bins is shown for (a) φB = 0.06, (b) φB = 0.10 and (c) φB = 0.12; γ = 6.The dashed line indicates the measured bulk volume fraction φB in the field of view averaged over the experiment.(d) The outer bins are one quarter of the channel width and the concentrations are calculated from the reconstructed 3D particle positions as the volume in the bin filled with particles divided by the total bin volume.

Figure 6 .
Figure 6.(a) Initial and steady state concentration profiles for γ = 6.Initial profiles were averaged over the first 20 frames (γa < 2.74) and steady profiles over frames 3500 to 3700 (γa = [480, 508]).The shaded regions indicate the measured bulk packing fraction φB in each experiment.(b) A bin of width w/100 is swept over the reconstructed 3D particle positions in x to measure local volume concentration φ yz.(c) Center-line (y = h/2) particle velocity profiles compared to theoretical Newtonian case.The curves are normalized by the theoretical velocity at the channel center Ucenter = 1.15 mm/s.(d) Cross-sectional particle velocity profiles.

Figure 7 .
Figure 7. (a) Mean squared particle displacements in x and z for γ = 6 experiments.The dashed lines show linear fits to the data used to extract the effective diffusivities via equation (4.1).(b) Comparison of effective diffusivities Dx (squares) and Dz (triangles) at different bulk volume fractions.

Figure 8 .
Figure 8. Experimental heat maps of x-particle migration over one cycle versus average distance to nearest neighbor at γ = 6 for (a) φB = 0.06, (b) φB = 0.10 and (c) φB = 0.12.The color map is normalized by the maximum bin count which is 56, 57 and 32, respectively.

Figure 9 .
Figure 9. Simulation heat maps of x-particle migration after 12 strain units versus average distance to nearest neighbor for (a) φB = 0.07, (b) φB = 0.11 and (c) φB = 0.14.The color map is normalized by the maximum bin count which is 22, 42 and 54, respectively.

Figure 10 .
Figure 10.Experiments are performed at φB = 0.10 and γ = 1 and the (a) bin concentration evolution, (b) concentration profiles and (c) migration heat map are plotted.The initial concentration profile is averaged over the first 20 frames (γa < 2.74) and the steady profile over frames 7000 to 7400 (γa = [960, 1016]).The shaded region in (b) indicates the measured bulk packing fraction φB.

Figure 11 .
Figure 11.Experiments are performed at φB = 0.11 and γ = 6 with a homogeneous initial concentration yet the steady state of the suspension is similar to the case with wall-loaded initial conditions.(a) The bin concentration evolution is plotted, with the error bars showing the standard deviation of 5 trials.(b) x-concentration initial (dotted) and steady (solid) profiles.The initial concentration profile is averaged over the first 20 frames (γa < 2.74) and the steady profile over frames 3500 to 3700(γa = [480, 508]).The shaded region indicates the measured bulk packing fraction φB.(c) Initial (dotted) and steady (solid) distributions of particle centers in the y-direction.The likelihood that a given particle's center will lie in each region of the channel is plotted.The initial distribution is averaged over the first 20 frames (γa < 2.74) and the steady profile over frames 601 to 656(γa = [82, 90]).The vertical dashed lines indicate a distance of one particle radius from the wall. (b).

Table 1 .
Relevant parameters and their range of values in our experimental study.