Sheathless hydrodynamic positioning of buoyant drops and bubbles inside microchannels

,


I. INTRODUCTION
Suspensions of solid particles, drops, or bubbles in a liquid are complex systems that exhibit nonlinear flow properties [1]. One such nonlinear effect is the spatial redistribution of dispersed particles due to hydrodynamic interactions within the flow [2]. The spatial distribution of particles during the flow of suspensions is relevant to many natural phenomena and technological processes; these range from the positioning of erythrocytes during flow of blood in the circulatory system [3,4] to the pumping of slurries in industry [5,6]. Understanding and controlling the position of dispersed particles in laminar flow, for example, enables analytical techniques such as field-flow fractionation [7,8] and flow cytometry [9].
We have investigated the dynamics of buoyant drops or bubbles suspended in a liquid carrier phase that flowed inside horizontal microfluidic channels, and were separated from one another by distances larger than their * gwhitesides@gmwgroup.harvard.edu sizes. Following an initial transient, these drops and bubbles moved in straight trajectories, parallel to the centerline of the channel. The position of these trajectories relative to the centerline depended on the geometry of the channel, on the size of drops, and on the properties of the carrier and dispersed phases. We propose a technique, based on this phenomenon, for controlling the steady-state transverse position of drops and bubbles by tuning the rate of flow and the viscosity of the carrier liquid.
Existing microfluidic techniques for positioning of particles [10] can be categorized into (i) methods that use additional laminar flows-sheath flows-to focus and locate the laminar fluid stream that contains the particles and (ii) "sheathless" techniques that use either externallyapplied forces or hydrodynamic lift forces (i.e. perpendicular to the local flow velocity) to move particles between different streams of fluid. Our technique is based on the competition between buoyancy and lift forces and does not require sheath flows.
Hydrodynamic lift forces cause the migration of particles across streamlines and are interesting phenomena in fluid dynamics; several mechanisms have been shown to produce lift forces [11][12][13][14][15][16][17][18][19]. Simple interpretations of these mechanisms might lead to the conclusion that lift forces are negligible in microfluidics. For example, inertial lift forces are proportional to the Reynolds number, which quantifies the ratio of inertial to viscous effects, but the Reynolds numbers characteristic to microfluidic flows [20] are typically much smaller than 1. Nevertheless, as Di Carlo et al. [21] have shown, inertial lift forces can focus and order particles in microchannels; microfluidic positioning of particles using inertial effects is a quickly growing field of research [22]. This recent work on inertial focusing shows how studying known hydrodynamic phenomena in a microfluidic setup can lead to the discovery of new phenomenology and new applications.
Here we investigated the position of drops and bubbles (deformable objects with fluid interfaces) within the cross-section of microfluidic channels. In our experiments the inertial forces were much weaker than those investigated by Di Carlo et al., and competed with lift forces caused by the deformation of drops [11,12], by the flow of interfacial fluid due to gradients in surface tension [13,19,23], and by the hydrodynamic interaction of drops and bubbles with the walls of the channel [16,18,24,25]. Figure 1a illustrates the geometry of the system: a carrier liquid that contained buoyant drops or bubbles flowed horizontally in a rectangular microchannel, with the system imaged from the side to determine the vertical position of drops and bubbles. Figure 1b shows drops of supercooled liquid water [26] (ρ ≈ 1.0 g/cm 3 ) and a particle of ice [27] (ρ ≈ 0.92 g/cm 3 ) flowing in a microchannel in a continuous phase of perfluoromethyldecalin (ρ ≈ 2.0 g/cm 3 ).
The drops of water and ice particles shown in Figure  1b did not contact the walls; their transverse position depended on the relative magnitude of buoyancy and hydrodynamic lift forces. In the direction along the horizontal imaging axis, drops did not experience buoyant forces and were centered between the vertical walls of the channel. In this experiment, liquid drops traveled close to the center of the channel while ice particles migrated to an intermediate position between the center of the channel and its top wall. The different equilibrium position of water drops and ice particles in Figure 1b indicates the presence of a non-inertial lift mechanism at work in our system; in other microfluidic inertial focusing experiments [21] liquid drops and solid particles migrated to the same positions.
Drops and bubbles deform in viscous shear flow [28], and this deformation leads to a lift force that, in most of the cases, causes the migration of fluid particles to the regions where the gradient of flow velocity (the shear rate) is smallest [11,12], such as the center of channels in Poiseuille flow. The deformation-induced lift force could balance the buoyancy of drops but might not explain fully the positioning shown in Figure 1b. We deduce this because the drops are not visibly deformed and because existing analytical models of deformation-induced lift [11] predict a lift force that would be too weak to balance the buoyancy of drops. To understand how lift forces may cause the centering of drops, we investigated the positioning of drops experimentally and numerically, and we compared our results with analytical predictions.

A. Motivation and Significance
Our goals were to determine (i) the steady state transverse positions of buoyant drops and bubbles immersed in a carrier liquid during flow in microchannels, (ii) how these positions change with the hydrodynamic properties of the flow, and (iii) to what extent analytical models of hydrodynamic lift forces [11,14] can be used to predict the positions of drops and bubbles in our experiments. We have been initially motivated to conduct these studies by an experiment in which we studied the nucleation of ice inside drops of supercooled water [26]; in that experiment, the transverse position of the drops influenced their temperatures, which we needed to measure accurately to determine the rate of nucleation of ice in supercooled water. Accurate and reproducible positioning of drops is also critical to microfluidic applications that investigate drops optically using focused laser beams, such as fluorescence-activated droplet sorting (FADS) [29] and lasing with high-speed switching in trains of drops [30]. A more general relevance of drop positioning in microflu-idics is that the transverse position of drops in channels determines their velocity because the carrier velocity profile in pressure-driven laminar flow is not uniform. The control of transverse position is thus relevant when drops are used as miniature chemical reactors in which timeresolved processes occur [31].
One issue that might occur in technologies that use multiphase flow is that buoyant particles could sediment in channels or pipes and clog the flow. The simplest solution to this problem is to use a carrier fluid that is density-matched with the particles, but there are instances in which this solution is not applicable (e.g. if the dispersed phase is gaseous, or if it has multiple components that have different densities). Also, densitymismatched carrier fluids could have characteristics that make them a better than density-matched fluids choice in some applications. Fluorocarbon liquids, which have approximately twice the density of water, are often used in microfluidics because they are chemically inert, immiscible with both aqueous solutions and hydrocarbon oils, and compatible with polydimethylsiloxane [32] (PDMS), the most common material used for the fabrication of microfluidic chips. Thus, our present study should prove useful for understanding positioning in situations where matching the density of the dispersed and continuous phases is not preferred.

II. BACKGROUND
The relative magnitude of viscous and inertial forces acting on a fluid determines the regime (laminar or turbulent) in which it flows. The ratio of inertial to viscous forces is quantified by the Reynolds number [33] Re, where µ is the viscosity of the fluid, ρ is the density of the fluid, and V the typical variation in the velocity of the fluid over a characteristic length scale L. All equations in this paper are written for variables expressed in the International System of Units (SI).
To characterize the relative strength of inertial effects in our system, we defined two Reynolds numbers: a channel Reynolds number, Re C , which characterizes the flow of a carrier fluid with viscosity µ c and density ρ c (Eq. 2), and a particle Reynolds number, Re P , based on the local shear rate, which characterizes the motion of drops and bubbles of radius r immersed in that carrier fluid (Eq. 3).
We chose the characteristic length of Re C to be the height of the channel, H, rather than its width, w, because H in our experiment is aligned with the gravitational acceleration, and with the deviation of drops and bubbles from the center of the channel. During viscous flow inside channels, the velocity of the fluid at the walls is zero, therefore the average velocity of the fluid V avg is a good measure of the variation of the velocity of the fluid across the channel. The characteristic length scale for Re P is the radius of the particle. The characteristic velocity is V avg r/H rather than V avg because for the evaluation of the viscous forces acting on the particle, the relevant velocity is the variation of fluid velocity in the vicinity of the particle (approximately equal to V avg r/H) rather than the variation of velocity in the whole channel (approximately equal to V avg ).
Flow in microfluidic channels is characterized by small to moderate Re C (between 10 −6 and 100); the flow is laminar in this range of channel Reynolds numbers. In the limit of unbounded flow at Re C = 0, lift forces on a rigid body of revolution (e.g. spheres, disks, and rods) that does not experience external forces vanish, and such particles move with the fluid rather than transversely with respect to streamlines [34]. Hydrodynamic lift forces appear because of small deviations from this limiting case [2]: these forces appear at finite Reynolds numbers, near walls, with particles that can deform, and with carrier fluids that are not Newtonian. Lift forces could also be caused by Marangoni effects [35] when thermal or chemical gradients generate a gradient in the surface tension along the interface between drops/bubbles and the carrier fluid.

A. Inertial Lift
Most of the inertial lift phenomena that have been investigated experimentally and theoretically belong to two classes: (a) migration of particles in the absence of external forces (e.g. buoyancy) and (b) migration of particles that are subject to an external force aligned with the direction of the flow. Most of the recent experimental work on microfluidic inertial focusing of particles is based on class (a) of lift phenomena [22].
Class (a). Inside circular pipes, neutrally buoyant spherical particles migrate to a circular annulus situated approximately halfway between the center and the walls of the pipe [36] when the channel Reynolds numbers have a moderate value (Re C ∼ 10 − 100). In channels with rectangular cross-section the inertia of the carrier fluid focuses particles to two or four positions situated near the walls of the channel [37]. These non-trivial focusing positions are caused by two opposing inertial lift mechanisms [14]: one mechanism pushes particles away from walls, and the other moves particles towards regions with larger shear rates (i.e. near the walls). The forces generated by both mechanisms are proportional to the particle Reynolds number Re P .
Class (b). Inertial lift forces push a particle towards the center of the pipe if the particle lags the flow and towards the wall if the particle leads the flow [17]. Leading and lagging the flow requires an external force that is parallel to the direction of flow to counter the viscous drag. This phenomenon was first investigated in vertical pipes using particles more dense than the fluid [38].; the particles either led or lagged behind the flow depending on the direction (up/down) of the flow. Kim et al. [39] have recently focused particles along the axis of a horizontal microchannel channel by slowing their motion relative to the channel; they used electrically charged particles, and they applied an electric field along the direction of the flow to slow down the particles.

B. Analytical and Empirical Expressions of the Inertial Lift Force
The calculation of inertial lift forces is a difficult problem, and several different formulas have been proposed. Ho and Leal [14] calculated the inertial lift force on solid spherical particles in planar Poiseuille flow in the limit of small particle Reynolds numbers (Re P 1). Matas et al. [40] used a numerical procedure described by Aslomov [41] to calculate inertial lift forces in circular pipes with Reynolds numbers up to Re C ∼ 1000; for Re C < 30 the lift force, F inertial L,analytic , is given by Eq. 4. For flow of large (i.e. comparable to the size of the channel) spherical particles in rectangular channels, Di Carlo et al. [37] calculated numerically the lift forces and presented an approximate formula (Eq. 5) for inertial lift near the center of the channel, F center L,analytic . These formulas for inertial lift forces depend on the distance, d, between center of the rigid particle and the center of the channel.
The functions g 1 (d/H) and g 2 (d/H) have a nonmonotonic dependence [37,40]  Drops and bubbles deform when they are immersed in the shear flow of another fluid [28]. The extent of deformation is characterized by the particle capillary number Ca P which quantifies the ratio of viscous stresses to the capillary forces acting on the drop or bubble. If the capillary number is small (Ca P < 0.5), the relative deformation of the drop (ratio of radius variation to radius average) is approximately equal to Ca P . For a drop or bubble with radius r and an interfacial surface tension γ, Ca P is given by The particle capillary number has the form shown in Eq. 8 rather than the simpler form µ c V avg /γ because, for reasons similar to those used in the definition of Re P (Eq. 3), the relevant fluid velocity is V avg r/H. The numerical factor of 2 ensures that the deformation of a drop is equal to Ca P at small deformations.
Drops and bubbles that are deformed by shear flows migrate towards the regions of flow that have the smallest shear rate [12]. Analytical models [11,42] that apply to small drops and bubbles predict that such migration only occurs if the ratio of viscosities of the suspended phase (drop/bubble) and of the carrier fluid is smaller than 0.5 or larger than 10; in between these values, the direction of migration is reversed. This analytical prediction was confirmed by the numerical simulations of Mortazavi and Tryggvason [43]. In pressure-driven flow, the shear rate is smallest in the center of the channel; drops and bubbles will therefore experience deformation-induced lift forces which will push them, depending on the viscosity ratio, towards or away from the center of the channel.

D. Analytical Expressions of the Deformation-Induced Lift Force
Different formulas to predict the deformation lift force in circular tubes [11,12,42] apply to small drops (r H) and are identical with each other except for a factor f (κ), where κ is the ratio of viscosities of the particle, µ particle , and the carrier fluid, µ c . Equations 9 and 10 list the deformation lift force F L,def ormation and the factor f (κ) proposed by Chan and Leal [11].

E. Wall-Effect Lift
In our experiments, drops and bubbles were largethat is, they were separated from the walls of the channel by distances comparable to their diameter. Drops and bubbles interact with the walls and experience wall-effect lift forces when their distance from the walls becomes smaller than a few times their diameter. The wall-effect forces in the case of small particles were considered both in Ho and Leal's treatment of the inertial lift force [14], and in Chan and Leal's treatment of the deformationinduced lift force [11]. When the distances of drops and bubbles from the walls become as small as those investigated in our experiment, however, the wall effects become more complex and difficult to predict [16,18,25]; for example, bubbles rising in a liquid near a wall are either attracted to or repelled from the wall, depending on the Reynolds number and on the degree of surface contamination of the bubble [24].

F. Viscoelastic Lift
In flow characterized by a non-uniform shear rate, the nonlinear relation between viscous forces and the shear rate in viscoelastic fluids can cause lift forces on particles. The viscoelastic lift acts in the direction of the regions which have the lowest shear rate. This lift force was predicted to move both solid [15] and fluid [11] particles to the center of channels in conditions of laminar flow. Leshansky et al. [44] observed viscoelastic focusing of particles along the center plane of a wide microfluidic channel.
Experimental investigations revealed more complex migration phenomena when viscoelastic and other types of lift act simultaneously on drops or bubbles. Sullivan et al. [45] observed that the transverse position of bubbles in a viscoelastic fluid became unstable when surfactant was added to the viscoelastic fluid. Yang et al. [46] explored conditions in which inertial and viscoelastic lift forces had comparable contributions in a microfluidic channel with a square cross section, and observed that particles migrated to the centerline. In cases where viscoelastic effects dominated inertial effects the particles migrated to multiple focusing positions.

G. Marangoni Lift in Systems with Thermal or Chemical Gradients
If the surface tension has a gradient along a fluid interface, the fluid at the interface will be pulled from regions with lower surface tension to regions with higher surface tension [35]. This interfacial flow entrains the bulk liquid and can lead to the migration of drops and bubbles. Gradients in the surface tension could occur due to thermal gradients or due to gradients in the concentration of surface-active solutes at the interface. Al-though several techniques have been developed to move fluids at the microscale by modulating the surface tension [47], Marangoni effects have not been used experimentally to manipulate the positioning of drops and bubbles in Poiseuille flow. Hanna et al. [13] have recently investigated theoretically the motion of a spherical drop in Poiseuille flow in the presence of surfactants, and reported that the flow-induced redistribution of surfactants causes the migration of drops to the centerline.

III. EXPERIMENTAL DESIGN
We focused our investigations on the measurement of the distance, d, between the steady state transverse position of drops and bubbles and the centerline of the channel (Figure 1a). The degree of centering of drops inside the channel depended on the hydrodynamic parameters of the system. These properties are the viscosities, densities, and flow rates of both the carrier liquid and of the dispersed phase, the surface tension, and the geometry of the system (i.e. the size of the channel and of drops). This large number of variables makes a thorough experimental investigation difficult; we therefore chose to investigate systematically only a subset of all parameters: the size of the drop, the viscosity of the carrier fluid, and the velocity of the carrier fluid. Varying the parameters in this subset had the most significant effect on the position of drops and bubbles.
Changing the surface tension and the properties (viscosity, density) of the drop or bubble can have an effect on their steady state transverse positions, but such changes were difficult to control experimentally. Surface tension can be varied by changing the fluid of the drop or bubble but this change could also affect viscosity or density; moreover, using surfactants to tune the surface tension has the disadvantage of generating Marangoni forces due to the redistribution of surfactant. Using multiple channels with different cross-sections or drops of fluids with different densities was cumbersome. We did not expect that varying the viscosity of drops a few-fold would change the lift force significantly as long as the ratio of the viscosities of the drops to the carrier fluid remained small (< 0.1); the inertial [16] and deformation-induced [11] lift forces were predicted to depend weakly on the viscosity ratio if the ratio is smaller than 0.1.
The geometry of our system is typical for microfluidics and differs from most analytical models used in the calculation of lift forces [11,14,42] in two aspects: the crosssection of the channel was rectangular (the models were derived for circular cross-sections), and the diameters of drops and bubbles were not very small compared to the height of the channel (0.35-0.6 H in our experiments). Because of these differences, analytical models of inertial [14] and deformation-induced [11,42] lift cannot be expected to predict accurately the position of drops in our system. We therefore ran numerical simulations of some of our experiments to learn whether our confined geom-etry affects the hydrodynamic lift forces qualitatively or quantitatively.

A. Setup for Experimental Measurements
The experimental setup, shown in Figure 2a, produced drops and bubbles in immiscible carrier liquids and transported these drops along a horizontal rectangular channel; this setup allowed us to vary continuously and independently the size of drops or bubbles, the viscosity of the carrier fluid, and the velocity of the carrier fluid. To produce drops and bubbles, we used soft lithography [48] to fabricate a microfluidic device composed of a flow-focusing drop generator [49][50][51][52] connected to a main straight channel with a rectangular cross-section (125µm wide, 200-µm high). The device also included a side channel connected to the long channel downstream from, and near to, the flow-focusing nozzle. The side channel allowed us to pump additional carrier fluid to increase the rate of flow in the main channel without interfering with the production of drops.
We used a main channel which was higher than it was wide because this geometry allowed us to produce spherical drops and bubbles with diameters as low as one quarter of the height of the channel using a simple design for the flow-focusing nozzle. The deviation d from the centerline could then vary over a range of ∼ 50 µm depending on experimental parameters, and we could measure d with satisfactory resolution and accuracy because our imaging system had an optical resolution of ∼ 3 µm. Along the width of the channel, the range of d was limited to at most 20 µm. In addition, we could not image drops and bubbles along this direction with good quality because the soft lithographic fabrication produced top and bottom walls with rough surfaces.
During experiments we fed fluids into the device to produce monodisperse drops or bubbles at a constant frequency of generation using syringe pumps operated in constant volumetric mode and a supply of gas at constant pressure. These drops and bubbles were initially located in the center of the channel, but then drifted to a steady state position above the centerline (Figure 2b). Drops and bubbles usually reached stable heights after they traveled several millimeters down the channel. To ensure that the vertical position of the drops and bubbles was stable, we investigated their position in a region that was near the end of the channel (∼35 mm away from the flow-focusing nozzle). We imaged the channel using a microscope that was rotated such that its imaging axis was horizontal, and we used high-speed cameras [53] to record the positions of drops and bubbles. A flow-focusing microfluidic device fabricated in PDMS using soft lithography was placed on a temperature-controlled plate and imaged along a horizontal axis to determine the steady state vertical position of drops and bubbles. The device also had a side channel that could be used to pump carrier liquid in the main channel without affecting the generation of drops and bubbles. b) The drops and bubbles were produced in the center of the channel, but drifted to higher steady state positions after traveling a few millimeters along the channel. The picture was constructed from three partially overlapping images taken at different times during the experiment. In this experiment the lift forces were weaker than in the experiment shown in Figure 1b and did not center the drops. The fluorocarbon carrier fluid had an average velocity of 9 mm/s and a viscosity of 6 mPa s.

Independent and Continuous Variation of Experimental Variables
We investigated the effect of a given experimental parameter (drop size, carrier viscosity, or carrier velocity) on the position of drops and bubbles while we held the other two constant. To change the size of the drops and bubbles without affecting the velocity of the carrier, we changed the temperature of the flow-focusing nozzle. This technique tunes the viscosity of the carrier fluid at the nozzle and has the same effect on the generation of drops and bubbles as does changing the velocity of the carrier at the nozzle [54]. We used a multi-zone temperature-controlled plate [26] to set the temperatures of the channel and of the nozzle separately. To change the velocity of the carrier, we pumped additional carrier fluid through the side channel; the additional fluid did not influence the generation of drops, but we did have to adjust the pressure of nitrogen gas to keep the diameter of bubbles constant.
We varied the temperature of the channel to tune the viscosity of the carrier fluid. The viscosity the hydrocarbon and fluorocarbon liquids that we used varied between 10 mPa s and more than 200 mPa s over the temperature range from 50 to -15 • C. Over this range of temperatures, the other hydrodynamic properties (fluid densities, surface tension, the size of the channel) changed by less than 10%. Using the variation of temperature to tune viscosity simplified measurements and expedited the acquisition of data since we did not have to change the carrier phase for each measurement. The temperature control system had a reproducibility of 0.1 • C and an accuracy of 2 • C. The relative accuracy of the values of the properties of the fluids was 0.1 for viscosity, 0.01 for density, and 0.1 for surface tension.

Choice of Fluids
In the ice nucleation experiments that inspired this work, and in the experiment shown in Figure 1b, we used drops of liquid water in a carrier phase of perfluoromethyldecalin (PFMD, 98% purity, F2 Chemicals) that contained 2% by volume 1H,1H,2H,2H-perfluorooctanol (97% purity, Sigma-Aldrich) as a surfactant. For systematic investigations, however, we chose continuous and dispersed phase fluids that did not contain surfactant to avoid the possibility of surfactant-induced Marangoni effects [13]; in addition, we selected carrier phase fluids that were compatible with PDMS [32], and whose viscosity varied approximately ten-fold as the temperature changed over a range of ∼50 • C near room temperature. Here we report data for (i) drops of liquid water (ρ ≈ 1.0 g/cm 3 ) in perfluoroperhydrophenanthrene, (PFPHP, ρ ≈ 2.0 g/cm 3 , 87.5%, F2 Chemicals Ltd.) and (ii) bubbles of nitrogen (ρ ≈ 0.001 g/cm 3 ) in Dynalene SF (DySF, ρ ≈ 0.9 g/cm 3 , Dynalene Inc.), a mixture of alkylated aromatic hydrocarbons that was designed for use in heat exchangers. We chose DySF over other hydrocarbon liquids because its temperature-dependent properties (density, viscosity, thermal conductivity) had been published by the manufacturer. We also report a few experiments performed at room temperature on bubbles of nitrogen in silicone oil (RT 500, ρ ≈ 1.0 g/cm 3 , Cannon Instrument Company).

Numerical Techniques and Software
For direct numerical simulations of the flow we used TransAT (Transport Phenomena Analysis Tool, AS-COMP GmbH [55]), a commercial multi-physics finitevolume code based on the multi-fluid Navier-Stokes equations and capable of handling multiphase and microfluidic flows. This software allows to model the two-phase flow with both the Level-Set Method (LSM) and the Vol-ume of Fluid (VOF) interface tracking methods [56,57]; phase field methods [58] have been recently implemented in TransAT as well. The VOF method associates to every element of the computational domain a function equal to the volume fraction of one of the phases, while LSM associates to the two-phase system a continuous function ψ which has positive values in one phase, negative in the other, and is zero at the interface; ψ evolves temporally during the simulation and the interface is located through the condition ψ = 0. Phase field simulations use an auxiliary function that has fixed values for each phase and varies sharply but continuously at interfaces; unlike LSM and VOF, the phase field method treats surface tension as being spread over a diffuse-interface instead of a sharp interface.
We chose to use LSM over VOF because the accuracy of simulations was better. Phase field simulations parameterize surface tension differently and could produce, for our problem, more accurate results than LSM and VOF because our interfaces are characterized by large energies (i.e. near-spherical drops). However, phase field simulations would require significantly denser grids and thus more and longer computations than LSM. Computations were performed on an 8-core workstation with Intel R Xeon R processors (E5630 CPUs, 2.53 GHz clock rate) and 24 Gb of RAM.

The Geometry of the Model and its Discretization
To reduce the time needed for computation we performed most of the simulations in a two-dimensional (2D) geometry. We also performed two simulations using a three-dimensional (3D) geometry for the conditions corresponding two experiments with drops of water in PF-PHP. The steady state positions of drops were different by only a few percent between the 2D and the 3D simulations. See Appendix B, Figures 10-11, for plots of the trajectories of drops in 2D and 3D cases. Our choice of 2D simulations is also supported by the results of Mortazavi and Tryggvason [43], who simulated the flow of neutrally buoyant drops using similar 2D and 3D geometries and found that drops reached approximately the same steady state positions in both cases.
We tested several geometrical domains that were discretized using grids with different densities to find the conditions that would provide the convergence of numerical calculation while allowing the fastest computations. See Appendix B, Figures 10-11, for numericallycalculated drop trajectories for identical physical parameters but different grid densities. For drops that had a radius equal to 0.2-0.3 H (a range that covers most of our measurements), the optimal 2D domain was rectangular and had a length that was twice the height of the channel along the direction of flow; we used periodic boundary conditions along the boundaries that intersected the direction of flow. The 2D domain was discretized using a structured mesh with 54×96 points (height×length) or more, and was uniformly distributed or slightly stretched in order to have a finer grid close to the walls.

Validity Tests for Numerical Simulations
In addition to verifying that our results did not depend on the density of the grid and on the 2D/3D geometry, we ran simulations for the same parameters of some of the results for neutrally buoyant drops and moderate Reynolds number (Re P = 0.0625 to Re P = 12.5) presented in Mortazavi and Tryygvason [43]. Comparison against their results provided a good validation test because (i) it involved the dynamics of drops that experienced both inertial and deformation-induced lift forces and (ii) Mortazavi and Tryggvason used a finite difference/fronttracking method [59] that is distinct from LSM. Starting with conditions used in Mortazavi and Tryggvason's Figure 6 (curves a, b, d) and Figure 8a, we calculated trajectories that were very similar to their results. See Appendix C, Figure 13, for a graph with our trajectories.

Selection of Conditions for Numerical Simulations
The numerical modeling for a given set of conditions required significantly more time than performing experimental measurements; therefore, we could not run as many simulations as experiments. We selected for a given channel and pair of fluids several experiments that could illustrate the effect of varying only one of the experimental parameters (drop size, carrier viscosity, and carrier velocity), and we simulated numerically the flow of drops or bubbles under those conditions. The results of these simulations are listed in Table I in Appendix A. In addition, we investigated the combined effects of inertial and deformation-induced lift by varying Re P and Ca P independently of each other during the simulations. Figure 3 shows the vertical steady state position of drops of water as a function of their diameter, of the viscosity of fluorocarbon, and of its average velocity. As the fluorocarbon became more viscous, and as it flowed at a higher rate, the drops of water traveled closer to the center of the channel.

Drops of Water in Liquid Fluorocarbon
Larger drops traveled closer to the center of the channel when all other flow conditions were constant. We investigated quantitatively only drops or bubble that had a smaller diameter than the smallest dimension of the channel (its width w) to avoid deformation of drops by  confinement against the side walls of the channel. Drops that were squeezed by the walls along the width direction traveled closer to the center than spherical drops. Figure 3a illustrates our variable-temperature technique for the control of viscosity. We could increase the viscosity of PFPHP by a factor of almost 20 by changing its temperature from 50 to -20 • C; see EPAPS Document No. [number will be inserted by publisher ], Figures S1-S3, for the dependence of the relevant properties of PFPHP, water, and DySF (viscosity, surface tension, and density) on temperature. The variation of these properties with temperature (with the exception of viscosity) was small; however, the visibility of drops varied significantly because the refractive indexes of water and PHPHP are nearly matched at room temperature. These changes in the visibility of drops did not affect our measurements.
We could investigate channel temperatures below the equilibrium freezing point of water because drops of pure water can be supercooled significantly in microfluidic channels [26]. Figure 3a shows drops of pure water that that remained in liquid form at -3 • C and -17 • C.   was that the smallest bubble in Figure 4c did not follow the general trend of the data. The range of accessible bubble diameters was narrower than for liquid drops, because we could not produce, using our flow-focusing generator, bubbles as small as the smallest drops of water.

The Dependence of the Steady State Position on the Product of Viscosity and Velocity
The dependences of the steady state position on the viscosity and velocity of the carrier phase were similar for both PFPHP (Figures 3b and 3d) and DySF (Figures  4b and 4d); this fact suggests that the relevant variable in determining the steady state position was the product of velocity and viscosity, µV avg . Figure 5 shows the dependence of d on µV avg . For DySF we recorded the dependence of d on viscosity and velocity using bubbles of the same size, and the data sets for viscosity and velocity variation overlap. The corresponding data sets for PFPHP used drops of different diameters; these sets did not overlap.
The dependence of d on the product of viscosity and velocity suggests that d might be determined by the capillary number Ca P . The data in Figure 5, however, does not collapse into a single curve when it is plotted against Ca P , indicating that the positioning mechanism is more complex than that based on the deformation of the drops; one possible complication is the contribution of inertial lift forces, which are quantified by Re P .

Large Bubbles of Nitrogen in Silicone Oil
In channels with a height of 200 µm, and for the pairs of fluids that we investigated, the deformation of drops and bubbles due to shear was too small to observe. The data shown in Figures 3 and 4 are characterized by much lower capillary numbers than those tested in previous experiments [12], for which the theory of deformationinduced lift [11,12] gives accurate predictions.
To bridge the gap between the conditions of our experiments and those performed by Goldsmith and Mason [12] we pumped an emulsion of nitrogen bubbles in silicone oil (RT 500, Cannon Instrument Co.) in a large "millifluidic" channel (2 mm high, 1.17 mm wide, and 110 mm long) using the arrangement shown in Figure 2. The geometry of this large channel, built from laser-cut acrylic sheets, was similar to that of the microfluidic channel used in our other measurements. Due to the larger size of the whole "millifluidic" device (75 × 150 × 7 mm), we could not use the temperature-controlled plate; operation at room temperature reduced the range of hydrodynamic parameters that we could investigate in the large channel. Figure 6 shows the steady state vertical position of ni- trogen bubbles (∼1 mm diameter) carried by a stream of silicone oil for four different conditions. In the large channel, bubbles deformed visibly, and the extent of deformation increased for larger capillary numbers. Bubbles of approximately the same size traveled closer to the center of the channel if Ca P was larger.

Other Combinations of Dispersed and Continuous Phases
We surveyed other pairs of carrier and dispersed phases was to identify flow conditions for which drops or bubbles would center tightly in the channel (as they are in the pictures in the right side of figures 3a and 4a). We investigated the flow of drops and bubbles immersed in carrier liquids other than PFPHP, DySF, and RT 500, and we also tested channels with cross-sections half and twice as large as the channel used to acquire the data shown Figures 1-4. We surveyed (i) bubbles of nitrogen in pure water, (ii) bubbles of nitrogen in aqueous sucrose solutions (40% and 50% sucrose by weight), (iii) bubbles of nitrogen in pure ethylene glycol, (iv) bubbles of nitrogen in liquid fluorocarbons (PFPHP, perfluoromethyldecalin, perfluorodimethylcyclohexane, ρ ≈ 2 g/cm 3 , pure or mixed with 1H,1H,2H,2H-perfluorooctanol), (v) drops of water in the same fluorocarbon liquids, and (vi) drops concentrated aqueous cesium chloride solutions (9 mol/kg molality, ρ = 1.8 g/cm 3 ) in DySF (we used cesium chloride solutions instead of pure water to have a difference in density between the hydrocarbon carrier and the aqueous drop on the order of 1 g/cm 3 ). The difference in the densities of drops/bubbles and the carrier phase varied between 0.9 and 2 g/cm 3 .
The only cases in which we could not center buoyant drops and bubbles were when we used water or perfluorodimethylcyclohexane as the carrier phase. These two liquids had the smallest viscosities among the liquids we surveyed, approximately 1 mPa s and 2 mPa s, respectively, near room temperature. As we increased the rate of flow of water and perfluorodimethylcyclohexane, we reached conditions for which Re C > 10 and drops and bubbles moved away from the centerline, indicating that inertial lift forces determined their positions. The results of this survey indicate that for a given channel geometry, the viscosity of the carrier must be larger than a certain value to achieve the centering of buoyant drops and bubbles; if the viscosity is lower, inertial lift forces dominate the dynamics of drops and bubbles and prevent their centering. For our 200-µm high channel, this critical value was approximately 20 mPa s; at this viscosity, Re C ∼ 10 for an average carrier velocity of 100 mm/s.

B. Numerical Simulations
Existing analytical formulas of various lift mechanisms do not apply directly to our system. The major differences from analytical models of inertial [14] and deformation-induced [11] lift are (i) the presence of an external force (the buoyancy) acting perpendicular to the direction of flow, (ii) the relatively large size of our drops (previously reported analytical formulas apply rigorously in the limit r/H 1), (iii) the rectangular geometry of our channels, and (iv) the fact that several lift mechanisms could contribute simultaneously to the total lift force. Our numerical simulations used the same underlying assumptions about the fluid dynamics as did the analytical models (i.e. an isothermal system with homogenous fluid properties), and allowed us to investigate how rectangular geometry, confinement, and buoyancy affect the positioning of drops and bubbles.

The Dependence of Positioning of Fluid Particles on Their Deformation
The centering of drops and bubbles that we observed in experiments indicates that deformation-induced lift is a plausible mechanism for the force that balances buoyancy because drops and bubbles became more centered as the viscosity of the carrier fluid increased. Although, in most of our experiments, we were not able to observe any deformation of drops or bubbles, it is possible that drops and bubbles were nonetheless deformed, and that the deformation was too small to observe. The minimum observable change of drops and bubbles from a circular shape was limited by our imaging resolution to ∼2%. According to the theory of the deformation-induced lift, the lift force is proportional to the deformation; it is therefore possible that in the absence of other lift mechanisms (e.g. inertial), weakly deformed drops would still migrate across streamlines because of their deformation.
Such small deformations could be investigated in numerical simulations, and one of our goals was to investigate numerically the effect of small deformations on the positioning of drops. Figure 7a shows the effects of deformation (quantified by the particle capillary number Ca P ) and of buoyancy on the positioning of drops. We performed three sets of simulations at different capillary numbers (Ca P = 0.001, 0.03, 0.3), and within each set we simulated one buoyant and one neutrally buoyant drop. We kept the particle Reynolds number constant at Re P = 0.011 in all simulations. For the simulation of a buoyant drop with Ca P = 0.001 and Re P = 0.011 we used the same geometric and hydrodynamic parameters as in experiment 7 in Table I, Appendix A.
The drops reached different steady state heights after being released from a common position, which is the same response that we observed in our experiments (Figure 2b). The steady state height was reached after the drops traveled a distance along the channel equal to at most 10-20 times the height of the channel. Buoyant drops traveled closer to the top wall than neutrally buoyant drops did. Drops characterized by a larger capillary number (e.g. having a lower surface tension) were more deformable and traveled closer to the center of the channel. None of the neutrally buoyant drops migrated exactly to the center of the channel; this observation suggests that deformation-induced lift (which centers drops and bubbles) and inertial lift (which pushes them away from the center) had the same order of magnitude. Figure 7b shows the drift of drops towards their steadstate position, for several values of Re P . These drops were not buoyant, had the same diameter, and the same capillary number. Since inertial lift becomes stronger as Re P increases, we expected that the inertial lift would dominate the deformation-induced lift and that for larger values of Re P neutrally buoyant drops would settle between the centerline and the walls. The drops nevertheless migrated to the center of the channel for all but one of values of Re P we tested; counter-intuitively, the only one case in which drops did not center was the one with the smallest Re P . For very small Re P (i.e. very slow flows), the periodic boundary conditions imposed at the inlet and outlet of our computational domain led to small transverse oscillations of the drops, such as those visible in the trajectory for Re P = 0.011. Though the transverse oscillations did not change the overall shape of the trajectory, they are illustrating nevertheless the difficulty of simulating flows characterized by very small Re P . For the simulation at Re P = 0.011 the lack of resolution around the center line of the channel (the spacing of the grid was equal to 0.025 H) produced small oscillations of the shape of the drop; these oscillations might have prevented the drop from reaching the center. The slope of the curves in Figure 7b is proportional to the transverse velocity of drops, and the transverse velocity is proportional to the strength of the lift. These slopes thus provide information about the strength of the lift forces acting on drops. As Re P increased among the simulations represented in Figure 7b, the transverse drift velocity decreased, at least in the first part of the trajectory. This dependence is consistent with an enhancement in the inertial lift force relative to the deformationinduced lift force. Mortazavi and Tryggvason observed the same dependence on the Reynolds number in their simulations [43].

Comparison between Measurements and Simulations
We modeled the flow of drops and bubbles for several conditions corresponding to the data in Figures 3-5. Figure 8 summarizes the comparisons between experimental and numerical results. For a given set of fluids and channel geometry, numerical simulations overestimated the distance between the steady-state positions of drops and bubbles and the centerline. Except for the case of experiments in which drops or bubbles were approximately centered (d/H < 0.05), the simulations predicted the position of drops within a factor of two of the experimental results. See Appendix A, Table I, for the parameters and the results of these simulations, and for the corresponding numerically-calculated trajectories (Figures 14-17). Numerical simulations agreed with the experimental observation that drops and bubbles traveled closer to the center as the diameter of the drops, the velocity of the carrier, and the viscosity of the carrier increased.
The simulations predicted best the behavior of large bubbles of nitrogen in silicone oil. These experiments were characterized by the largest capillary numbers among all systems that we investigated, and thus were the least prone to numerical errors. The simulations for bubbles in DySF (Table I, experiments 9-12) and a few among the simulations for droplets of water in PFPHP (Table I, experiment 4), which were characterized by very small Ca P and Re P , were affected by significant parasitic currents (numerically generated vortices near interfaces [56]) whose amplitude became in same cases comparable to the magnitude of the flow internal to the drop.

Analytical Predictions of the Steady-State Position of Drops
We calculated the steady-state position of drops by requiring that the sum of buoyant, deformation-induced lift (Eq. 9), and inertial lift (Eq. 4) forces equals zero. We digitized and then fitted with a third-order polynomial the inertial lift force for Re = 30 in Figure 14 from Matas et al. [40] to determine the value of the function g 1 (d/H) (Eq. 6).
The analytical approach gave good predictions in the case of the 2-mm high channel. Formulas fared significantly worse than simulations in predicting our measurements in the 200-µm channel, and failed to predict the placement of drops and bubbles away from walls due to the flowing carrier fluid. Several of the steady-state positions calculated analytically, including all those for which d/H > 0.5, lie outside the channel; physically they correspond to drops and bubbles that flow in contact with the walls of the channel. Although steady-state positions with d/H > 0.5 are unphysical, we plotted them to underline that confinement and wall effects played an important role in determining the equilibrium position of droplets. Figure 9 shows the velocity fields obtained from numer-   Table I, Appendix A). To investigate experimentally the case of a neutrally buoyant drop, we rotated the channel and the imaging system such that the imaging axis was vertical; in this experimental geometry, there is no component of the buoyant force in the plane of the image. Figures 9c and 9d show the equilibrium configuration for a bubble of nitrogen in Dynalene SF (experiment 11 in Table I) and for a bubble of nitrogen in RT 500 (experiment 14 in Table I). The simulations shown in Figure  9 modeled qualitatively well the shape of the drops, but the agreement between the numeric and the experimental transverse steady state positions was not perfect.

V. DISCUSSION
The hydrodynamic lift on transversally buoyant (i.e. gravitational acceleration perpendicular to the direction of flow) drops and bubbles has been relatively unexplored compared to the neutrally buoyant and longitudinally   [38,60], rather than perpendicular to this direction (as in our experiments). Hogg investigated theoretically the inertial forces acting on settling particles in horizontal flow [61], and found that inertial forces are in general weak compared to buoyant forces. With the exception of studies on bubbles rising near walls in a stagnant fluid [16,18,24] and in vertical shear flow [25], experimental, analytical, and computational studies of deformationinduced lift [11,12,43] have investigated only the case of neutrally buoyant drops.

A. Positioning of Drops as an Outcome of the Competition between Mechanisms of Hydrodynamic Lift
Theoretical investigations of inertial and deformationinduced lift have treated these phenomena separately, addressing the regimes in which either Re P or Ca P was separately dominant, and the other negligible. Numerical simulations by Mortazavi and Tryggvason [43] of the motion of neutrally buoyant deformable drops for Re C > 1 have shown that the steady state position of drops depends on the relative strength of deformation and inertial effects; very deformable drops migrate to the centerline, while drops that had a near-spherical shape stabilize at intermediate positions between the center and the periphery of the channel. They have also investigated the movement of large drops with diameters larger than half of the height of the channel, and found that large drops became centered even when the Reynolds numbers were large. Our numerical simulations agree with Mortazavi and Tryggvasons results, and indicate that the position of drops and bubbles is determined by the combined contributions of buoyancy and lift forces.

Expected Steady State Positions of Drops and Bubbles
Assuming that the only lift mechanisms are inertial and deformation-induced, and that these mechanisms are not coupled, the total lift force can be calculated as the sum of inertial (Eq. 4) and deformation-induced (Eq. 9) lift forces. This total lift force depends nonmonotonically on the distance from the center of the channel and vanishes at the steady state position, d, of neutrally buoyant drops and bubbles. The possible values of d range from zero (if the deformation-induced lift is stronger than the inertial lift near the center of the channel) to the inertial focusing position (if the deformationinduced lift is negligible compared to the inertial lift). For large drops and bubbles this description of the total lift force is likely to be complicated by the deformation of drops and bubbles near walls [16,18,25]. Deformationinduced wall effects should generate a centering lift force which, when added to the total lift force, will bring the steady-state position closer to the center of the channel.
At the steady-state position of buoyant drops and bubbles, the sum of all lift forces and of the buoyant force vanishes. Since the buoyant force is directed away from the center of the channel, the steady-state position of buoyant drops and bubbles should be further away from the center than in the non-buoyant case. To avoid con-tact with the top wall, drops and bubbles must experience a total lift force that is not smaller than the buoyant force. Our experiments and simulations show that buoyant drops can be supported by lift forces, and that buoyancy makes them travel at a further distance from the center than non-buoyant drops.

Explanation of the Dependency of Positioning on Experimental Parameters
Experimental dependencies on drop size, carrier viscosity, and carrier velocity can all be explained qualitatively assuming that the forces acting on drops and bubbles are the inertial and deformation-induced lift, and buoyancy. (i) Size of drops or bubbles. The inertial (Eq. 4) and deformation-induced (Eq. 9) lift forces are proportional to r 4 , while buoyant forces are proportional to r 3 ; therefore as the size of drops and bubbles increase, the hydrodynamic lift becomes stronger relative to buoyancy, and drops and bubbles are centered more tightly as long as the deformation-induced lift dominates inertial lift (which was the case of our measurements). (ii) Viscosity of the carrier fluid. The deformation-induced lift force is proportional to with µ 2 c , while the inertial lift force and the buoyant force do not depend on the viscosity of the carrier[62]. If the viscosity increases, fluid particles move closer to the center because the deformationinduced force increases while the inertial lift and the buoyancy remain constant. (iii) Velocity of the carrier fluid. Both lift forces are proportional to V 2 avg , therefore the total lift force is proportional to with V 2 avg as well; because the total lift force becomes stronger relative to the buoyant force, the steady state position of drops and bubbles moves closer to the center.

B. The Possibility of Other Lift Mechanisms in Our Experiments
The combined effects of buoyancy, deformationinduced lift, and inertial lift explain qualitatively our experimental and numerical results. The large quantitative differences between analytical predictions and measurements ( Figure 8) can be explained by the large size of drops and bubbles [16,43], and possibly by a coupling between inertial and deformation-induced lift when they have contributions of similar magnitudes [16].
Our numerical simulations predicted a smaller degree of centering than that observed in experiments. We believe that our numerical simulations did account for the finite size of drops and bubbles, for their proximity to the walls of the channel, and for the coupling of inertial effects and deformation. This consistent discrepancy between simulations and experiments might be explained by the presence of other lift mechanisms; we therefore investigated whether other lift phenomena could be present.

Variations in the Size of the Channel, Viscoelastic Effects, and Hydrodynamic Interactions with other Drops or Bubbles
We have chosen our fluids and experimental conditions to minimize possible positioning effects due to geometric variations in the cross-section of the channel, and due to viscoelastic effects. These local geometric variations could be an artifact of the soft lithographic fabrication, or be caused by the deformation of the channel under the pressure of the carrier liquid. For the channel that we used to record the data presented in Figures 3-5, the total height variation along the channel was ∼0.01 of the height of the channel. Such a small variation could not have caused the migration of drops and bubbles; we made measurements in a channel with the same geometry but with smaller variations in height and we observed the same positioning of drops and bubbles. The larger channel used for experiments with bubbles in silicone oil had uniform height, being fabricated from flat sheets of acrylic glass.
Hydrostatic pressure could deform a channel made in a soft material (such as PDMS) and thus create trenches; such trenches can guide drops and bubbles [63]. The microfluidic devices that we used, though molded from soft PDMS, had glass side walls which were practically rigid. One of these hard walls was the glass slide to which we bonded the PDMS slab. The slab contained another glass slide, embedded in PDMS during the molding of the device.
We verified using a rheometer [64] that the carrier fluids which we used were not viscoelastic at shear rates from 0 to 3000 s −1 (a range of shear rates that covers all our experiments).
The typical separation between consecutive drops and bubbles in most experiments was on the order of ten times the diameter of drops or bubbles. In experiments with bubbles of nitrogen in DySF, we verified that the transverse position of bubbles did not depend on their density; the shortest distance between bubbles that we investigated was approximately equal to one bubble diameter. We have also simulated numerically the flow of nitrogen bubbles in RT 500 for different densities of bubbles, and we did not observe any changes in the steadystate transverse position of bubbles. Based on this insensitivity of the transverse position to the distance between bubbles, we concluded that drop-drop or bubble-bubble interactions did not play a significant role in the transverse positioning of drops and bubbles.

Marangoni Effects
Drops and bubbles are hydrodynamically different from solid particles because (i) they are deformable, (ii) they are made from a fluid material, and (iii) they have fluid interfaces. Our simulations, which took into account deformability and the flow inside drops, showed that these characteristics of drops alone might not account for the full strength of the lift force that we observed experimentally. It is possible that the dynamics of the fluid-fluid interface played a significant role in the generation of lift forces, through a Marangoni effect.
a. Thermal Marangoni Effects. Gradients of temperature in the carrier fluid cause gradients in the surface tension along the interface of drops and bubbles and could lead to their migration. In many experiments, the temperatures of the nozzle and of the channel were different; the temperature of the fluids varied along the channel, and most of the variation in temperatures occurred near the nozzle. This longitudinal temperature gradient generated a transverse gradient of temperatures in the cross-section of the channel because of the finite thermal diffusivity of the carrier fluid [26]. To minimize the transverse thermal gradients, we investigated the position of drops as far as possible (∼35 mm) from the nozzle. Using a model for thermal transfer in the carrier fluid that we described previously [26], we determined that the maximum transverse temperature gradients at the location where we investigated the transverse positions of drops and bubbles were ∼0.5 • C/mm in PF-PHP and ∼0.05 • C/mm in DySF. These gradients were much smaller than the temperature gradients that produce thermal Marangoni forces which are comparable to the buoyant forces acting on drops and bubbles [19] (on the order of 10-100 • C/mm for the fluids and sizes of drops that we investigated). Since thermal Marangoni forces are proportional to the gradient in temperature, we inferred that thermal Marangoni forces had a negligible contribution to the steady-state position of drops and bubbles.
b. Chemical Marangoni Effects. We use here the term chemical Marangoni effects to refer to the movement of a fluid interface in an isothermal system due to variations in the surface tension; these variations can occur due to inhomogeneities in the chemical properties of the interface, such as the composition of a carrier liquid that is not chemically pure. Lift forces caused by chemical Marangoni effects can act on bubbles immersed in a shear flow when the carrier fluid contains surfactants [13,65]; in this case, the shear causes an inhomogeneous distribution of surfactant along the interface.
We did not use any surfactants for data shown in Figures 3-6 to reduce the magnitude of possible chemical Marangoni lifts. We cannot exclude the possibility that chemical Marangoni lift forces acted on our drops and bubbles, however, because the carrier liquids that we used were not pure chemical compounds. PFPHP was only 87.5% pure (the remainder being perfluorinated molecules with structures similar to PFPHP), DySF is a trademarked mixture of alkylated aromatic hydrocarbon molecules, and RT 500 silicone oil is composed of a mixture of molecules (mostly polydimethylsiloxane) that have a range of molecular weights. It is therefore possible that under shear flow the interfaces of drops and bubbles with PFPHP, DySF, and RT 500 are inhomogeneous.
Shear-induced inhomogeneities might be nevertheless small because the pure components of a carrier fluid were likely to have similar properties (chemical structure, boiling points, and surface tensions). We did not expect PF-PHP and DySF to contain significant amounts of surfaceactive molecules because the values of surface tension (at room temperature, ∼50 mN/m for the interface of water with PFPHP, and ∼30 mN/m for the interface of between nitrogen and DySF) were similar to those of pure hydrocarbon [66] and pure fluorocarbon [67] liquids.

C. Guidelines for the Control of Transverse Position of Buoyant Drops and Bubbles in Microfluidic Channels
The balancing of buoyant and lift forces provides a method for controlling the steady state height at which drops and bubbles flow in a horizontal microfluidic channel. The highest degree of control possible by this method involves flowing drops and bubbles anywhere between either the top or the bottom walls (depending on the sign of the buoyant force) and the center of a channel. We found that the carrier fluid must have a viscosity of at least 20 mPa·s to center drops or bubbles in a channel with cross-sectional dimensions on the order of 100 µm.
Based on experiments with a channel with a cross-section of approximately 1 mm, we inferred that the threshold of viscosity is a function of the size of the channel: larger channels will have a larger threshold.
If the viscosity of the carrier is large enough to allow centering of drops and bubbles, it is possible to control the height of their trajectories as well. Starting from centering conditions, a reduction in the rate of flow or in the viscosity of the continuous phase will change predictably and monotonously the position of drops and bubbles; bubbles and drops lighter than the carrier will travel higher than the centerline, and heavy drops will travel below the centerline.

VI. CONCLUSION
We have investigated the movement of buoyant drops and bubbles flowing in a liquid carrier phase in horizontal microfluidic channels. Despite their buoyancy (the densities of drops and bubbles were ∼ 1 g/cm 3 smaller than that of the carrier phase), drops and bubbles flowed without touching the top wall of the channels because hydrodynamic lift forces balanced their buoyant force. The vertical position of drops and bubbles depended mainly on the size and buoyancy of drops and bubbles, and on the viscosity and velocity of the carrier liquid. In carrier fluids with viscosities larger than 20 mPa s, we could always make drops and bubbles travel within several micrometers of the centerline of a 200-µm channel by increasing the rate of flow of the carrier fluid up to average carrier velocities of 100 mm/s. The positioning of buoyant deformable particles in horizontal Poiseuille flow is a relatively unexplored problem. We hypothesized that the behavior of such a system is determined by the sum of the buoyant force and all lift forces acting in the system; this hypothesis thus ignores the possibility of a strong nonlinear coupling between lift mechanisms and buoyancy. Our experiments and numerical simulations were designed to separate the effects of inertial and deformation-induced lift, and our results show that the positioning of buoyant drops and bubbles can be understood qualitatively as an outcome of the competition between buoyancy and lift mechanisms.
We found that analytical theories of deformationinduced lift [11] and inertial lift [14,40] do not provide a satisfactory quantitative prediction of lift forces. These formulas apply to drops and bubbles much smaller than the cross-section of the channel. We infer from this disagreement that in microfluidic channels the lift forces acting on large bubbles and drops are significantly enhanced by proximity to the walls of the channel. New analytical theories of hydrodynamic lift in microchannels must therefore take into account confinement effects [37] such as the proximity of walls.
The numerical simulations predicted that drops and bubbles would be somewhat less centered than we observed experimentally. The discrepancy between numerical simulations and experimental results suggests that an additional type of hydrodynamic lift force (not inertial and not deformation-based) might contribute to the centering of drops and bubbles.
Out of all lift mechanisms that we found in the literature, and of the ones that we could imagine, the only one that could act during our experiments is a chemical Marangoni effect. It is possible that shear-induced inhomogeneities in the carrier fluid near the surface of drops and bubbles generate interfacial flows that produce lift forces. We could not evaluate the magnitude of this mechanism; such an effect is likely to be weaker than the Marangoni lift in carrier fluids that contain surfactants.
In our experiments, the viscosity of drops and bubbles was typically lower than one tenth of that of the carrier fluid, but it is possible that our conclusions are applicable to more viscous drops. Recently, Hur et al. [68] reported that neutrally buoyant, near-spherical drops of oil ∼1000 times more viscous than the carrier traveled closer to the center of the channel than solid particles. Their viscous drops had spherical shapes, therefore the deformation of their drops could not generate significant lift forces; these spherical viscous drops experienced nevertheless a stronger centering lift force than solid particles of same size. Further experimental and theoretical investigations are needed to understand the lift mechanisms acting on spherical drop and bubbles in microfluidic channels, and to predict accurately their magnitude.
Our measurements show that buoyant drops and bubbles can be reproducibly positioned and centered in microfluidic channels without using sheath fluids. The phenomenon that enables such positioning is the presence of centering hydrodynamic lift force mechanisms that act on fluid particles (drops and bubbles). Sheathless positioning of drops and bubbles using buoyancy and hydrodynamic lift forces is a simple and reliable method that could become a useful tool for microfluidic analytical applications. Table I lists the conditions of the experiments that we selected for comparison with simulations, and includes experimental results along with numerical and analytical predictions. We selected from the experiments on drops of water in PFPHP pairs that illustrated the effect of varying one hydrodynamic parameter while the others were constant. Experiments 1 and 2 vary the size of the drop; 3 and 4, the viscosity of PFPHP; 5 and 6, the velocity of PFPHP; and 7 and 8 are distinguished by the presence or absence, respectively, of a buoyant force.
Simulations of nitrogen bubbles in DySF (experiments 9-12) required the most intensive computations and we could not simulate pairs of experiments the same way as we did in PFPHP. The four experiments we selected for simulations were characterized by some of the largest Re P and Ca P among experiments with bubbles in DySF; we estimated that for experiments with smaller Re P and Ca P the calculation of bubble trajectories would require more than a month of computation time on a single-core workstation.
Simulations of nitrogen bubbles in RT 500 were the fastest because the Ca P values were largest among all pairs of fluids that we investigated. We were limited in our choice of experiments, however, because the parameters of these experiments could not be controlled independently, and because the number of experiments we could perform was small. These experiments were timeconsuming because they required large rates of flow of RT 500; we could only make a few measurements before needing to reload with RT 500 the syringes that fed the "millifluidic" device. In principle, the most accurate numerical simulations of our experiments should be three-dimensional and be performed on a computational grid that is as dense as possible so as to capture variations in velocity and shape. We therefore attempted to identify a setup for numerical simulations that would provide the results we presented here within days or weeks.
The most radical simplification in our simulations was the use of a 2D computational domain to simulate a 3D problem. We justified this choice empirically by comparing trajectories calculated on 2D and 3D domains (Figures 10-11). For 2D and 3D simulations of experiments 7 and 8, the difference in the steady state positions was at most a few percent of the height of the channel. Mortazavi and Tryggvason [43] reported a similar degree of agreement between 2D and 3D simulations of two-phase flow, for larger Re P and Ca P than those characteristic of our experiments. Figures 10 and 11 also illustrate the effect of using grids with different densities and aspect ratios (here we define the aspect ratio as the ratio of the length of the simulated section of the channel to its height). The differences in the steady state positions between simulations calculated using different grids had the same magnitude as the differences between 2D and 3D simulations. Another example of simulations performed on different grid sizes is shown in Figures 12 and 13 (these simulations used the hydrodynamic conditions of experiments [13][14][15][16]. Due to the long duration of numerical simulations, we could not perform a systematic search for the optimal grid size and geometry. Instead, we tested a few different types of grids and we found that using these grids in our simulations provided consistent steady state transverse positions. The aspect ratio of our computational domain determines the distance between consecutive drops or bubbles. This distance was equal to the length of the computational domain and in most cases it was smaller than the experimental values. Figures 10, 12, and 13 show that when the distance between drops or bubbles varied from one and two times the height of the channel, the steadystate transverse of drops or bubbles did not change sig- FIG. 11. The influence of grid density, aspect ratio, and dimensionality of the computational domain on the results of numerical simulations. The data represents the trajectory of drops for conditions identical to those of experiment 8 from Table I, and the numerical parameters are the same as those used in Figure 10. In the experiment, the drops stabilized at y/H = 0.50.  Table I. The numerical parameters are the same as those used in Figure  10. In experiments, bubbles stabilized at y/H = 0.59 for experiment 13 and y/H = 0.63 for experiment 14. nificantly. Figure 14 shows the flow field calculated in a 3D domain for experiment 7; the 2D simulation is shown in Figure 9a of the main text. In Figure 14, we show the flow field in orthogonal planes that intersect in the center of the drop. The color-scale images (vertical and transverse to the flow) show the velocity of the carrier fluid in the reference frame of the channel, normalized to the average velocity of the carrier. The streamline plots (vertical and horizontal) represent the flow of the carrier in the reference frame of the moving drop.   Table I. The numerical parameters are the same as those used in Figure  10.  Table I. This simulation used a grid with 32×48×24 points and an aspect ratio of 2.  [43] performed their simulations using a finite difference/front-tracking method. Figure 15a corresponds to their Figure 6a, and Figure 15b to their Figure  8a, which we reprinted here [69] as Figure 15c. Figure  15a-b uses the same conventions (e.g. the z-axis is vertical) for notation and graphics layout as those used by Mortazavi and Tryggvason. The channel Reynolds number Re b listed in Figure 15b-c is almost identical to our Re C (Eq. 2); the only difference is that Re b uses the centerline velocity in its formula, rather than the average velocity.

Comparison with a Finite Difference/Front-Tracking Method
Our trajectories, though not identical, are similar to those reported by Mortazavi and Tryggvason. Trajecto-ries (b) and (d) in Figure 15a exhibit more noise in the position of drops than those calculated by them. We attribute this noise to a lower density of our grid than that of Mortazavi and Tryggvason; trajectory (a), which is smoother, was calculated using a grid as dense as theirs.

Appendix C: Simulations of Experiments with
Drops and Bubbles Figure 16 shows drop trajectories calculated for experiments 1-8 on drops of water in PFPHP. We chose to report as the steady state position of drops their position at the end of their respective trajectories, even though some of these trajectories were not perfectly stabilized. We argue that this was a satisfactory choice because we estimated that the transverse distance that the drops would drift until complete stabilization was of the same magnitude as the position differences between 2D/3D dimensionality and different grid densities (i.e. a few percent of the height of the channel). With the exception of experiment 2, we launched all drops from the same position (z/H = 0.6). We attempted to launch the drop from experiment 2 at z/H = 0.6 as well, but the calculations were too slow to be practical; we therefore launched this drop from a position that we expected would be close to the steady state position.
Instead of simulating trajectories such as those shown in Figure 16 for bubbles of nitrogen in DySF, our simulations launched bubbles from three different heights to search for a launch position from which the trajectory would be perfectly horizontal ( Figure 17); such a launch position would equivalent to the steady state position of bubbles. Since we could not always guess such position in three attempts, we evaluated the steady state position by interpolating the initial slope of the trajectory as a function of the launch position. Out of the three trajectories for each simulation, the lowest launch positions used a grid with 54×96 points and an aspect ratio of 2, and the upper launch positions used a with 72×108 points and an aspect ratio of 1.5.
Simulations of nitrogen bubbles in RT 500 ( Figures  12 and 13) were characterized a transient upward drift of bubbles just after their launch. This behavior was caused by the delayed deformation of bubbles. In these simulations the bubbles had a spherical shape when launched, the deformation-induced lift force was initially negligible, and the buoyant bubbles drifted up; after a short distance of travel, they started to deform due to shear; as the deformation-induced lift force became stronger, and the bubbles reversed the direction of their drift. Comparison of our LSM numerical approach with the finite difference/front-tracking method used in Ref. [43].
Here α is the ratio of densities, and γ the ratio of viscosities of drop fluid to carrier fluid. Graphs 15a-b show our results; graph 15b corresponds to graph 15c (reproduced from Ref. [43]) and graph 15a corresponds to Figure 6a in Ref. [43].