Universality in microdroplet nucleation during solvent exchange in Hele-Shaw like channels

Micro and nanodroplets have many important applications such as in drug delivery, liquid-liquid extraction, nanomaterial synthesis and cosmetics. A commonly used method to generate a large number of micro or nanodroplets in one simple step is solvent exchange (also called nanoprecipitation), in which a good solvent of the droplet phase is displaced by a poor one, generating an oversaturation pulse that leads to droplet nucleation. Despite its crucial importance, the droplet growth resulting from the oversaturation pulse in this ternary system is still poorly understood. We experimentally and theoretically study this growth in Hele-Shaw like channels by measuring the total volume of the oil droplets that nucleates out of it. In order to prevent the oversaturated oil from exiting the channel, we decorated some of the channels with a porous region in the middle. Solvent exchange is performed with various solution compositions, flow rates and channel geometries, and the measured droplets volume is found to increase with the P\'eclet number $Pe$ with an approximate effective power law $V\propto Pe^{0.50}$. A theoretical model is developed to account for this finding. With this model we can indeed explain the $V\propto Pe^{1/2}$ scaling, including the prefactor, which can collapse all data of the"porous"channels onto one universal curve, irrespective of channel geometry and composition of the mixtures. Our work provides a macroscopic approach to this bottom-up method of droplet generation and may guide further studies on oversaturation and nucleation in ternary systems.

This limitation can be overcome in a bottom-up approach such as solvent exchange (Lou et al. 2000;Zhang et al. 2015;, 2020, where a large number of micro and nanodroplets is generated by nucleation out of an oversaturated solution. This method, also called nanoprecipitation or solvent shifting (Fessi et al. 1989;Galindo-Rodriguez et al. 2004;Aubry et al. 2009;Lepeltier et al. 2014;Hajian & Hardt 2015), though commonly used, is much less well understood.
In solvent exchange, a good solvent of the target droplet component (the solute) is replaced by a poor solvent, where the two solvents are miscible. A typical example is an oil saturated aqueous ethanol solution being replaced by oil saturated water. Upon contact of the ethanol & water solution, the two solvents start to mix with each other. Due to the addition of water, the solubility of oil is lowered, and the subsequent oversaturation leads to droplet nucleation and growth. The micro & nanodroplets can nucleate in the bulk (Vitale & Katz 2003) or on a hydrophobic surface. For solvent exchange in a microchannel, it has been found that droplets nucleated in the bulk tends to migrate to and then stay in the center in a co-flowing device (Hajian & Hardt 2015), where the droplet movement is controlled by solutal Marangoni flow and composition of the mixture. On the other hand, for droplets that nucleated on the surface, their average volume is found to increases with the Péclet number P e as ∝ h 3 P e 3/4 , where the flow rate Q is included in P e and h the channel height. Later, the effect of flow geometry ) and solution composition (Lu et al. , 2016 on the average oil droplet size were also qualitatively investigated. The mutual interaction between a multitude of surface droplets and the resulting effect on the growth dynamics was also studied (Xu et al. 2017;Dyett et al. 2018). Despite all of these studies, a thorough understanding of the oversaturation pulse -which is crucial to droplet nucleation & growth -is still lacking, because of its transient nature and a lack of means to directly measure it.
To have a quantitative understanding of the oversaturation pulse, we study the total amount of oversaturated oil inside a Hele-Shaw like microfluidic channel, by measuring the total volume of the oil droplets that nucleate out from it using confocal microscopy. In this paper, we are neither interested in the nucleation process itself nor in the droplet morphology, since they do not help to quantify the oversaturation pulse. A theoretical model for the total nucleated oil volume is developed, based on the ternary phase diagram and Taylor-Aris dispersion. The model accurately predicts the scaling behavior of the total volume V of oil with respect to the Péclet number, V ∝ P e 1/2 , including the prefactor, in which the influence of the solution composition and channel geometry is reflected. However, to compare the prefactor with the experiments, we need to prevent the oversaturated oil -especially the bulk droplets that nucleated out of it -from leaving the channel. To achieve this, a porous region consisting of circular pillars is put in the middle of the channel. For channels with such a porous region, the prefactor can collapse different groups of data onto one universal master curve for different channel geometries and mixtures. For channels without the porous region, the measured oil volume is smaller than the theoretical prediction because some of the oversaturated oil -including nucleated droplets in the bulk -leaves the channel.  . Side view of the experimental setup. Solvent exchange is performed in a square microfluidic channel with or without the porous region, which is made of circular pillar arrays. Solution A (yellow) is decane (the oil) saturated aqueous ethanol solution and is injected into the channel first. Solution B (light blue) is decane saturated water and is injected later at a fixed pressure. In the mixture between the two solutions, oil is oversaturated so that droplets (red) nucleates out of it, both in the bulk and on the hydrophobic walls. A flowmeter is used to measure the flow rate of solution B. The channel height is h, width is W (in y direction, not shown), and the channel length before the porous region is L. Solvent exchange is observed from below by confocal microscopy. The inset shows the bottom view. Pillar diameter d and pillar spacing a and b are varied to change the porosity φ = 1 − πd 2 /(2ab). See table 1 for the parameters.

Experimental procedure & methods
Solvent exchange is performed in a thin square channel with various height h, width W and length L, see figure 1 for the definitions and table 1 for the parameters. The microfluidic channel is made of a glass wafer covered on a silicon wafer which is decorated with an inlet, an outlet, and some of them a porous region in the center. The porous region is made of an array of circular pillars with pillar diameter d ranging from 6.7 µm to 8.7 µm (see figure 1, inset), and pillar spacing in the transverse and axial directions a and b are varied to change the porosity φ = 1 − πd 2 /(2ab) (see table 1 for parameters). The length between the inlet & the porous region is L. For "smooth" channels, i.e., those without the porous region, L is the entire length of the channel. The whole channel is made hydrophobic by OTS (octadecyltrichlorosilane) coating: 1 M hydrochloric acid is first pumped through the chip at 50 µL/min for 20 min by using a syringe pump (Harvard, PHD 2000). The chip is then put in a vacuum chamber at 1.8 mbar for overnight to dry. A solution of OTS dissolved in hexadecane (Sigma-Aldrich, 99 %) at 0.4 v/v% is pumped through the chip at 50 µL/min for 20 min. The chip is then sequentially cleaned by chloroform, toluene, and ethanol, and finally dried in vacuum for use.
Solution A, which is rich in oil, consists of decane (oil) saturated aqueous ethanol solution. To make a solution A with the desired concentration, its ethanol-to-water weight ratio w e,A /w w,A is first determined. The solution is prepared by first mix ∼ 400 g mixture of ethanol (Sigma-Aldrich, 99.8 %) and water (Milli-Q) with this specific weight ratio w e,A /w w,A , then add in decane (Sigma-Aldrich, 95 %) until phase separation is observed, so that the solution is saturated. The weight fractions of each species are calculated from the actual ethanol-to-water weight ratio w e,A /w w,A in the solution. The oil weight fraction w o,A of solution A is increased by increasing its ethanol-to-water weight ratio w e,A /w w,A , since the solubility of oil (decane) in ethanol is higher than in water. Finally, solution A is labelled yellow by adding a small amount of perylene (Sigma-Aldrich, 99 %) at 0.2 mg/mL. Solution B, which is poor in oil, is made of decane saturated water and is labelled light blue by Rhodamine 6G (Sigma-Aldrich, 99%) at 0.2 mg/mL. Its oil weight fraction is w o,B = 5.2 × 10 −8 at 25°C.
Solution A is first injected to fill the entire channel, then solution B is injected at a constant driving pressure to perform the solvent exchange: An oil oversaturation pulse is generated in the mixture of solutions A and B, which leads to oil droplets nucleation both in the bulk and on the hydrophobic surfaces in the channel. The contact angle of oil in water is θ = 15±3°on the same treated silicon surface and θ = 11±2°on the same treated glass (see SI for more details). The flow rate Q of solution B, measured by a flowmeter (ML120V21, Bronkhorst, Netherlands), is varied by changing the driving pressure. The Péclet number of the flow is calculated by P e h = uh/D, where u = Q/(W h) is the average velocity of solution B, and a typical diffusivity of water in ethanol D = 0.84 × 10 −9 m 2 /s is used.
After about 0.3 mL (more than 1000 times of the channel volume ≈ 0.25 µL) of solution B is injected, the injection is stopped by closing the valve, and a 3D scan of the whole channel is recorded by confocal microscopy (Nikon A1, Nikon, Japan) from below.

Experimental results
Typical mid-plane snapshots of the upstream part of the channel with the porous region (Chip No. 1, see table 1 for the geometrical parameters) are shown in figure 2(a). The channel without porous region (Chip No. 6, see table 1) is shown in figure 2(b). Black signals pillars, light blue signals water, and red signals oil (because ethanol must have been dissolved in and washed away by the excess amount of water). In general, more oil droplets (red) are found in the channel after the solvent exchange. Furthermore, oil droplets are observed before and inside the porous region (black), but not behind the porous region. However, for channels without the porous region, droplets are observed in the entire channel (see SI for the snapshots of the full channel). A closer look at the porous region is shown in figure 2(c), where the oil droplets and the circular pillars are clearly shown. A typical 3D scan of the same area is shown in figure 2(d). Note that in this work, we only focus on the total oil volume V , but not on the droplet morphology.
Solvent exchange is performed for all the 6 chips (see table 1 for the geometrical parameters) at different w e,A and flow rate Q. After the solvent exchange, 3D scans of the whole channel, similar to that shown in figure 2(d) are performed. The total volume V of these droplets is measured by counting the number of red pixels of the 3D confocal image and then multiplied by the volume of one pixel, then they are plotted against P e h W (c) . It is found that for all the chips, V increases with the Péclet number as ∝ P e α h , with α ≈ 0.50, see table 2 for details. V also increases with w o,A , i.e., the more oil is in solution A, the more oil is nucleated. In the next subsection we will develop a theoretical model to quantitatively account for these two observations.

Theoretical model
Figure 3(a) shows the schematic of the initial oversaturation pulse (red shaded region), which is a mixture of the two solutions due to advection and diffusion. The initial oversaturation pulse takes a parabolic shape in laminar flow, which is the case here. It broadens because of (Taylor-Aris) diffusion of the three components: oil, ethanol and water. The concentration of the mixture changes continuously from that of the solution A to that of solution B (Ruschak & Miller 1972). Figure 3(b) is the phase diagram of the oil-ethanol-water system, it shows the concentrations of the solutions & the mixture, with the ethanol weight fraction w e and the oil weight fraction w o being the x and y axes, respectively. Then the water weight fraction of the solution/mixture can be calculated by  Table 2. Control parameters of the solvent exchange, and the best fitted power laws: V ∼ P e α h . wo,A is the oil weight fraction of the good solvent, we,A is the corresponding ethanol weight fraction, and α is the fitted exponent of the power law. 95 % confidence bounds of the fittings are also shown. The average exponent α is 0.50 ± 0.01. this is the so called "binodal curve". The curve is fitted from the measured data points (see Appendix A for details). Then, the concentrations of solutions A and B are on the binodal curve, denoted by points A and B in figure 3(b). The concentrations of the mixture lie on a curve connecting points A and B, denoted by the red curve AB -this is the so called "diffusion path" (Ruschak & Miller 1972). All the possible concentrations of the liquid -both of the solutions and the mixture -lie on this curve. The oil concentration c o in the mixture is higher than its saturated oil concentration c o,s , thus the mixture is oversaturated with oil. The (absolute) oil oversaturation is then denoted as (4.1) We now use mass-per-volume concentrations c in the equations for simplicity but keep using weight fractions w elsewhere. Notice that the weight fractions w are the concentrations c nomalized by the density of the liquid ρ: w = c/ρ.
For any oversaturated liquid parcel M in the mixture as represented by point M on the diffusion path, when the oversaturated oil nucleates, it is considered that only oil nucleates out of the mixture, and the ethanol-to-water ratio in the mixture is kept constant (Lu et al. 2016). The concentration of the mixture moves along the so-called "dilution curve" (Lu et al. 2016) to point S on the binodal curve. The dilution curve MS is a straight line passing through point (0, 1) in figure 3(b) (this point means pure oil in the phase diagram, see SI & Appendix for more details), as shown by the black solid line. Therefore, the (absolute) oil oversaturation of this liquid parcel M is: (4.2) The exact shape of the diffusion path is determined by the diffusion speeds of the three components. In a simplified case, we assume that the diffusion of each component only depends on its own concentration gradient. To get a first order approximation of the problem, we further assume that the diffusion coefficients of oil and water are equal: D o = D w ≡ D . Then the diffusion path AB becomes a straight line (see Ruschak & Miller (1972) and Appendix B. These two assumptions are reasonable since the oil concentration is small in all the tested cases: This function of course also depends on the initial condition c e,A . Eq.(4.3) is calculated numerically (see SI for more details), and the oil oversaturation normalized by the oil density, ∆c o /ρ o , is plotted as a function of the normalized ethanol concentration w e = c e /ρ in figure 3(c), for different initial conditions w o,A = c o,A /ρ: 0.063 and 0.017, respectively. Here ρ is the density of the liquid parcel (see SI for the calculation of ρ), and ρ o = 730 kg/m 3 is the density of the oil (decane).
The oversaturation pulse shown in figure 3(a) evolves as a function of space x and time t, so that both the oil oversaturation ∆c o and c e also depend on space and time: ∆c o [c e (x, t, c e,A )]. As mentioned previously, the total oil volume V is nucleated from all the oil oversaturation generated in the liquid. The total oil volume V in the porous media can then be calculated by integrating the oil oversaturation over the entire channel volume, at the time just before the oversaturation pulse reaches the porous region: The spatial and temporal distribution of ethanol concentration c e (x, t, c e,A ) of the parabolic shaped oversaturation pulse is non-trivial. However, the channel used here is long and thin. The time scale for the axial advection is τ a ∼ L/u 4 s, which is much larger than the time scale of the vertical diffusion τ v ∼ (h/2) 2 /D = 0.12 s. In this situation, the concentration gradient in the vertical direction is smoothed out because of Taylor-Aris dispersion (Taylor 1953;Aris 1956Aris , 1959, and the oversaturation pulse becomes like a plug, as shown in figure 3(d).
The analytical solution of the ethanol concentration c e (x, t, c e,A ) for pure water entering a thin square channel which at time t = 0 contains only an aqueous ethanol solution with (uniform) concentration c e,A (no oil present) in this plugged regime is (Taylor 1953): where x 1 = x − ut is the distance to the central position of the plug, as shown in figure 2(d). κ = 1 2 k − 1 2 t − 1 2 is the wavenumber (with unit m −1 ), where k = h 2 u 2 /(210D) for a thin square channel of thickness h (Dorfman & Brenner 2001). The wavenumber κ has the temporal dependence and its inverse scales as the length of the plug: L p ∼ 1/κ.
For the solvent exchange in our case, oil is present in all the liquid, but Eq.(4.5) still holds here because of the previous assumption D o = D w ≡ D. With this assumption, any portion of the water can be replaced by oil without influencing the above equation. Substituting Eq.(4.5) into Eq(4.3), we have: (4.6) Eq.(4.6) is solved numerically, and figure 3(e) shows how the normalized oil oversaturation ∆c o /ρ o changes as a function of κx 1 at the three different initial conditions. It is worth noting that the oil oversaturation is not symmetric, but its peak position is to the right (κx 1 ≈ 0.55), or in other words, more to the downstream. This asymmetry originates from the asymmetry in ∆co ρo (w e ) shown in figure 3(c), whose asymmetry originates from the shape of the binodal curve shown in figure 3(b). Also, the oil oversaturation is higher everywhere in the channel for the case with higher initial oil concentration (larger w o,A ).
Since the oversaturation pulse becomes a plug-like shape, which means the vertical concentration gradient is negligible, Eq.(4.4) can be further simplified as the integration along the single dimension x. Substituting Eq.(4.6), we have: (4.7) Substitute t = L/u into κ as the time when the oversaturation pulse is just entering the porous region. Then the integration is performed with respect to κx 1 , which leads to: is a dimensionless prefactor which is proportional to the area under the curve ∆c o (κx 1 , c e,A ), as shown in figure 3(e). From eq.(4.8) we see that the power law dependence of V on P e h is indeed predicted to be 1/2. The prefactor s is shown in figure 3(f ) as a function of normalized initial condition w e,A = c e,A /ρ. The two red dots correspond to the initial conditions w o,A = 0.063 and 0.017, respectively. It is worth noting that s changes sharply within this range.
Eq.(4.8) with eq.(4.9) is the main theoretical result of our paper. It is universal in the sense that it includes all the factors that could influence the total amount of nucleated oil: the flow rate, the channel geometry, and the solution composition.
Note that in the model, only Taylor dispersion in the post-free region is considered, because the porous structure immediately interrupts further development of the oversaturation pulse. Further dispersion/mixing in the porous region is negligible: The length of the oversaturation pulse, which covers 99 % of the oil oversaturation (see figure 3(e)), is calculated to be L p = 3.56/κ = 0.49 √ Lh · √ P e h > 0.44 mm, much larger than the pore size ≈ 20 µm. That means that c e does not really change on the length scale of the pore size, thus the porous media does not induce further dispersion/mixing. On the other hand, though in the theoretical model the length of the post-free region L is used to calculate the total oil volume, this does not mean that the oversaturation would accumulate in the pulse and then all of a sudden nucleate into droplets at some time point, for example, at t = L/u when the center of the oversaturation pulse is at the entrance of the porous region. Instead, as the oversaturation pulse develops while moving downstream, the (newly generated) oversaturation nucleates into droplets along the way, so that some of the droplets are observed upstream of the porous region.

Comparison between universal theoretical result & experiments
To compare the universal theoretical result with the experimental results, the oil volume V measured for various chips and mixtures shown in figure 4(a)&(b) are normalized by s·W h 3/2 L 1/2 , and then plotted against P e h in figure 4(c)&(d) in log-log scale. Results for channels with the porous region (chip No. 1-5) are shown in figure 4(c), and results for channels without the porous region (chip No. 6) are shown in figure 4(d). It is found that for chips with the porous region, the calculated prefactors (s = 1.2 × 10 −2 and 3 × 10 −3 ) can indeed collapse all data onto one master curve, regardless of the channel geometry and composition of the mixture. This is in agreement with the calculated prefactor of our theoretical model. However, for channels without the porous region, the measured oil volume V is smaller than the theoretical prediction -which is not true for channels with the porous region. This confirms that the porous region can indeed prevent the oversaturated oil from leaving the channel, by which the comparison between experiment & theory is enabled.

Conclusions and outlook
We have performed solvent exchanges in thin square channels with and without a porous region in the middle, at different initial conditions w e,A , flow rates Q and channel geometries. The total volume of these oil droplets V is measured by confocal microscopy, and is found to increase with Péclet number as ∝ P e α h , with α ≈ 0.50. A theoretical model is then developed, based on the ternary phase diagram and Taylor-Aris dispersion, to predict the total amount of oil V nucleated from solvent exchange. The theory indeed predicts a power law dependence V ∝ P e 1/2 h . In addition, the influence of the channel geometry and the initial condition w o,A can also be calculated and included in the model to give a complete prediction of the oil volume V on one universal curve, thanks to the porous region which can prevent the oversaturated oil from leaving the channel. This model is found to be able to predict the total nucleated oil volume V , irrespectively of the channel geometry & initial mixture.
The findings of this work contribute to a better understanding of the solvent exchange, and could guide further design and research on this topic. First, a porous media may serve as a good tool to collect all the oversaturation, providing a macroscopic approach to study this bottom-up method.
The results of this paper encourage us to propose a three-step approach to study solvent exchange: (i) Find the concentration distribution by solving the advection-diffusion equations. The boundary conditions are defined by the flow geometries and the initial conditions are set by the initial solution concentrations. (ii) Based on the concentration distribution, the oversaturation distribution can be calculated by applying the knowledge of the phase diagram. (iii) Investigate the quantity that is of interest, such as the volume of the nucleated phase, the dynamical interaction between the nucleated phase and the oversaturation, etc.
This research can be considered as a demonstration of the above proposed approach, with the total nucleated oil volume V being the subject of interest, and the aid of Taylor-Aris dispersion in a long & thin channel to obtain the analytical solution of the concentration distribution. For more general cases where the oversaturation pulse is not a plug, analytical solutions of the concentration distribution may be difficult to get, then numerical simulations should be employed to finish the first step. It is also easier to incorporate the ternary (or multicomponent) diffusion effect in numerical simulations. Moreover, the dynamical interaction between the nucleation and the oversaturation pulse can also be incorporated in numerical simulations.
The X-and Y-coordinate in Fig.3(b) are actually w e = c e /ρ and w o = c o /ρ. The binodal w o,s (w e ) is fitted by two Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) through the data points, which are measured by titration following the procedure as described by Tan et al. (2016).
For the first part (first 10 points, lower half), the polynomial is: is valid on the 9 intervals between the 10 data points. Here w i is the ethanol weight fraction of the i th data point. For the second part (last 5 points, upper half), the polynomial is: where P i (w) has the same form as Eq.A 2. Parameters are shown in table 3.
Appendix B. Transformation from a ternary phase diagram to the phase diagram in a cartesian coordinate Figure 5(a) shows a typical ternary phase diagram, using the decane-ethanol-water system as an example. The three vertices E, W, D stand for ethanol, water and decane, respectively. Figure 5(c) shows the same phase diagram in a cartesian coordinate, with the ethanol weight fraction w e being the X-coordinate and the oil weight fraction w o being the Y-coordinate. Here we prove that any straight lines in the ternary phase diagram shown in figure 5(a) are still straight lines in figure 5(c).
To prove this, we need to find the transformation matrix between these two vector spaces. Let p and q be the basis vectors of the vector space that the ternary phase diagram is in. p is parallel to line WE in figure 5(a), and q is parallel to line WD. Let i and j be the basis vectors of the cartesian coordinate. These two basis are put together in figure 5(b), with (p, q) shown in red dashed arrows and (i, j) shown in black solid arrows. Then we have:  Figure 5. Transformation from the ternary phase diagram to a right-triangle phase diagram in a cartesian coordinate. (a) The ternary phase diagram of the decane-ethanol-water system. The three apexes are E for ethanol, W for water, and D for decane. Blue curve is the binodal curve fitted from the data points (blue diamonds) taken from Skrzecz et al. (1999). Red line AB is the diffusion path when we,A = 0.837, and black line DM is the dilution curve. (b) Basis (p, q) of the ternary phase diagram shown in red dashed arrows and basis (i, j) of the cartesian coordinate shown in black solid arrows. (c) Right-triangle phase diagram of the same system. Diffusion path AB and dilution curve DM are still straight lines. we is replaced by x and wo is replaced by y.
This is a linear transformation, with being the transformation matrix.
With the assumption of no ternary diffusion and D e = D o ≡ D, Ruschak & Miller (1972) first proved that the diffusion path in the ternary phase diagram is a straight line. By definition, dilution curve is also a straight line. Because linear transformation transforms straight lines to straight lines, diffusion path AB and dilution curve DM in figure 5(c) are still straight lines.
On the dilution curve DM in figure 5(a), ratio of water weight fraction to ethanol weight fraction is kept constant, that is: where c is a positive constant. In the cartesian coordinate, w e is the X-coordinate and w o is the Y-coordinate. For easier notation, let x = w e and y = w o , then w w = 1 − x − y, and  . Diffusion path (red curve) for the case when Dw = 0.84 × 10 −9 m 2 /s and Do = 0.43 × 10 −9 m 2 /s, with initial conditions wo,A = 0.063, ww,A = 0.119. Notice now the diffusion path (red curve AB) changes to a curve rather than a straight line as shown in figure  3(b). The blue diamonds are the measured data on the binodal, and the blue curve is the binodal. The black line is the dilution curve.