Modelling double emulsion formation in planar flow-focusing microchannels

Double emulsion formation in a hierarchical flow-focusing channel is systematically investigated using a free energy ternary lattice Boltzmann model. A three dimensional formation regime diagram is constructed based on the weber number of the inner phase $We_i$, and the capillary numbers of the middle $Ca_m$ and outer $Ca_o$ phases. The results show that the formation diagram can be classified into periodic two-step region, periodic one-step region, and non-periodic region. By varying $We_i$ and $Ca_m$ in the two-step formation region, different morphologies are obtained, including the regular double emulsions, decussate regimes with one or two alternate empty droplets, and structures with multiple inner droplets contained in the continuous middle phase thread. Bidisperse behaviors are also frequently encountered in the two-step formation region. In the periodic one-step formation region, scaling laws are proposed for the double emulsion size and for the size ratio between the inner droplet and the overall double emulsion. Furthermore, we show that the interfacial tension ratio can greatly change the morphologies of the obtained emulsion droplets, and the channel geometry plays an important role in changing the formation regimes and the double emulsion sizes. In particular, narrowing the side inlets or the distance between the two side inlets promotes the conversion from the two-step formation regime to the one-step formation regime.


I. INTRODUCTION
Double emulsions are droplets with one other droplet inside. Their core-shell structure has attracted wide attentions in various fields [60]. In pharmaceuticals, one common technique is to use double emulsion for drug encapsulation of highly hydrophilic molecules. It solves the low encapsulation efficiency problem faced in single emulsion technique due to the quick partitioning of the hydrophilic molecules into the external aqueous phase [24]. Double emulsions are also suitable containers for performing chemical and biochemical reactions under well-defined conditions. Compared to bulk reactions, the greatly reduced volume needed in double emulsion technique is beneficial for high throughput screening assays [6]. Furthermore, double emulsions can be used as templates for the synthesis of more complex colloidosomes [4,31], as well as for controlled release of the inner contents [14]. To ensure the successful applications of double emulsions, one of the key issues is to provide precise control over the double emulsion structure, size and monodispersity at a sufficient production rate [54,67].
Traditional double emulsion fabrication methods, such as the bulk emulsification and the membrane emulsification methods [60], are attractive to many industries (e.g. food and cosmetic) where scalability for large production is important [59]. However, these techniques have poor size and monodispersity control [55], which makes them inadequate for applications requiring high precision, such as in medical, pharmaceutical, and material industries. The emergence of microfluidic technology [58,63] opens up a new horizon. It provides more detailed control over the operating conditions [60] and offers great flexibility in designing multi-layer [3] or multi-core emulsions [42,47]. So far, the microfluidic devices for producing double emulsions can be roughly classified into a series of two T-junctions [49], two flow-focusing junctions [2,50,53], co-axial capillaries [42,58], and the possible combinations and variations of the aforementioned geometries [46,69].
The understanding of double emulsion formation dynamics are crucial for microfluidic control and equipment optimization. Double emulsions are commonly generated either in a two-step or one-step formation regime, depending on whether the inner part of the double emulsion is sheared simultaneously with the middle layer in the outer fluid [58]. Due to the distinct flow details in the two-step and one-step formation regimes, the influence of flow conditions, fluid properties and geometrical parameters on each regime should be analyzed respectively. For the two-step formation regime, Okushima et al. [49] have systematically showed the effect of flow rates on the double emulsion sizes and the number of internal droplet for multi-core emulsions when they are formed using a series of T-junctions. The one-step formation regime is mostly encountered in co-axial microcapillary devices. Experimental studies have been carried out on the effect of flow rates [25,31] and geometrical settings [40]. Scaling laws have also been developed for the emulsion size predictions [7,58].
Complementary to experiments, numerical studies on double emulsion formation dynamics in microfluidic channels have also garnered strong interest. For instance, great efforts were made to elucidate the effects of flow rates, interfacial tension ratios, geometry [10,41], and viscosities [18] on the double emulsion properties and the flow regime predictions [40] for co-axial flow-focusing capillary devices. Simulations are particularly advantageous for providing accurate flow details and for allowing each relevant parameter in the system to be varied systematically. In the literature, a number of ternary fluid models have been successfully developed and applied in the study of multiple emulsion formation behaviors, including using the volume of fluid (VOF) method [4,10,40,41], the front-tracking method [61], the free energy finite element method [51] and the lattice Boltzmann method [17,18].
In this work, our focus is on the planar flow-focusing cross-junctions. They are promising for integration with other devices and they can be parallelized to raise the production rate of the emulsion droplets, while still ensuring good size control [32]. Furthermore, in contrast to other microfluidic geometries, systematic parametric study is rarely reported on planar flow-focusing devices. Several works, such as Abate et al. [2] and Azarmanesh et al. [4], briefly discussed the possible conversion between the two-step and one-step formation regimes and the variation of shell thickness. However, it remains unclear in which flow rate regions monodisperse double emulsions are produced; and correspondingly, how the droplet sizes can be varied in those regions. It is likely that the droplet sizes have different dependencies on the flow rates for the two-step and one-step formation processes. There are also open questions on the role of channel geometry in the formation regime conversion, and on the effects of interfacial tension ratio in determining the morphologies of the emulsion droplets, including the possibility of complete, partial and non-engulfing shapes [8,21,45,50].
We have chosen to use the lattice Boltzmann method (LBM). LBM is highly favorable for the study of emulsion formation behaviors due to its simplicity in solving interface dynamics, including droplet break-ups and coalescences, as well as its ability to deal with complex geometries, and its high efficiency in parallel computation [26]. So far, three types of ternary lattice Boltzmann models have been developed, including the free energy model [1,35,52,64], color-fluid model [17,18,29,30,66], and the Shan-Chen type models [5,62].
Here we improve on the free energy lattice Boltzmann model developed by Semprebon et al. [52]. A major progress is that our model allows a wider range of interfacial tension ratios, such that all possible biphasic emulsion morphologies can be captured [21], including complete and non-engulfing shapes. The model by Semprebon et al. [52] only allows partial engulfing shapes. Coupling the free energy model with the advantages of the lattice Boltzmann method, we conduct a systematic study on the dynamics of double emulsion formation behaviors in planar hierarchical flowfocusing junctions. We focus on the two-dimensional (2D) case to reduce the computational time needed for parametric studies. The major physical difference in the flow dynamics between the 2D and 3D systems lies in the lack of an additional Laplace pressure induced by the out-of-plane curvature [9], which may accelerate the droplet pinch-off process [23]. However, most of the fundamental flow physics are still involved in the 2D system. Thus, we believe our 2D study still captures the qualitative trends in the formation regimes and emulsion sizes as function of the flow rates of each fluid phase.
The paper is organized as follows. In §2, we describe the improved ternary free energy model, the lattice Boltzmann method, and the boundary conditions involved. In §3, we validate the model and boundary conditions by Poiseuille flow, moving droplet and static emulsion morphology tests. In §4, our systematic parametric study allows us to construct a flow regime diagram, where we describe a wide range of formation regimes, including the periodic twostep and one-step double emulsion formation regimes, decussate regime, bidisperse regime and even the continuous structure with multiple inner droplets. Scaling laws are also proposed for the double emulsions produced in the onestep formation regime, and the effects of the interfacial tension ratios and the geometrical parameters are investigated. Finally, we summarize our main findings and forecast prospects for future work in §5.

A. Free-energy model
The present model is developed based on the equal-density ternary free-energy lattice Boltzmann model proposed by Semprebon et al. [52]. The thermodynamics of the system is described by the free-energy functional which is constructed using a double-well potential form for the bulk free energy E 0 (C m ) = (κ m /2)C 2 m (1 − C m ) 2 and a square gradient term for the interface region E interf ace (∇C m ) = (α 2 κ m /2)(∇C m ) 2 . Ω is the system volume. C m (m = 1,2,3) are the concentration fractions with two minimum values at C m = 0 and C m = 1 for each component m. In the current model, all components have the same density ρ m , which we have set to be 1.0 for simplicity. Thus the total density is related to the concentration fractions by defining ρ = C 1 ρ 1 + C 2 ρ 2 + C 3 ρ 3 = C 1 + C 2 + C 3 . Three physically meaningful bulk phases termed red, green and blue here could be denoted by C = [C 1 C 2 C 3 ] = [1 0 0], [0 1 0] and [0 0 1], respectively. α is the interface width parameter. The interfacial tension between red-green phases σ rg , red-blue phases σ rb , and green-blue phases σ gb are related to κ's through To capture the interface dynamics, two order parameters φ and ψ are introduced as and the original concentration fields can be reversely obtained from ρ, φ and ψ via C 1 = (ρ+φ−ψ)/2, C 2 = (ρ−φ−ψ)/2, and C 3 = ψ. The order parameters and the hydrodynamics of the ternary fluid system are governed by two Cahn-Hilliard equations, the continuity and the Navier-Stokes equations: Here, v is the fluid velocity and η is the dynamic viscosity. M φ and M ψ are the mobility values for φ and ψ. The default value of M φ is 0.5, and M ψ = M φ /3 is considered to assign symmetric mobility for each concentration component [52]. The thermodynamic properties are related to the equations of motion via the chemical potential µ φ , µ ψ and the pressure tensor P. The chemical potential is defined as the variational derivative of the free energy µ q = δF/δq, where q = ρ, φ or ψ. The pressure tensor term in Eq. (7) is constructed by ∇ · P = ∇(ρc 2 s ) + ∇ · P 0 , where c s = 1/ √ 3 is the lattice speed of sound and ∇ · P 0 = (ρ∇µ ρ + φ∇µ φ + ψ∇µ ψ ). The first term ∇(ρc 2 s ) is the standard ideal gas term in LBM and it is simply incorporated in the equilibrium distribution function. The F = −∇ · P 0 term is treated as an external force term in the lattice Boltzmann implementation. The explicit expressions of µ ρ , µ φ , µ ψ , and P are given in Semprebon et al. [52].
Consider now a case where a red droplet is completely engulfed by a green one and they are submerged in a blue phase fluid at thermodynamic equilibrium. According to the theoretical analysis of Guzowski et al. [21], the interfacial tensions should satisfy σ rb > σ gb + σ rg . Given the relation between the σ's and κ m 's in Eq.(2), one can easily find that κ 2 should be negative while κ 1 and κ 3 are positive. In the free energy model, the negative κ 2 will invert the shape of the bulk free energy profile E 0 (C m ): the two minimum values at 0 and 1.0 become two maximum values as shown by the blue solid line in figure 1. As such, setting one of the κ m 's to be negative often leads to incorrect results or even numerical instability as the concentration value deviates significantly from [0, 1.0]. A similar situation has been encountered in other LB models, and a simple remedy has been proposed by introducing an additional free energy term, see Lee & Liu [33] and Abadi et al. [1]. Inspired by these works, to solve the problem induced by negative κ m , here we introduce an additional free energy term given by where β is an adjustable positive parameter controlling the slope of the energy profile E 0 + E a as depicted by the red dashed line in figure 1. Since we add a new free energy term in Eq. (8), additional terms should be included in the chemical potentials accordingly, which are listed in Appendix A.

C. Boundary conditions
The boundary conditions involved in the present study contain: no-slip boundary, wetting boundary and the inletoutlet boundary. No-slip boundary condition is used on the solid walls, which is realized by the half-way bounceback rule [28]. The solid walls should have a preferential affinity with the continuous phase fluid to generate stable droplets/emulsions [2]. Fu et al. [18] successfully implemented the wetting boundary condition by setting a fictive density on the walls in a LB ternary color-fluid model. Similarly for the free energy model used here, the macroscopic values of ρ, φ and ψ on the walls are designated to be the same as those of the continuous phase fluid that is assumed to completely wet the walls. For the velocity inlet, the Zou-He velocity boundary condition [70] is applied to solve the unknown density distribution functions of f i . To obtain the unknown g i and h i values at the inlet, the method used by Hao & Cheng [22] and Liu & Zhang [37] is adopted. For instance, given an inlet boundary with the inflow direction pointing to the right, g 1,5,8 and h 1,5,8 are unknown after the streaming step. We assume that one pure single fluid exists at the inlet, where the prescribed values of φ and ψ are φ in and ψ in , respectively. The sum of the unknown distribution functions can be solved according to Eq. (17), and then g 1,5,8 and h 1,5,8 are allocated by their weight factors as For the outlet boundary, the convective boundary condition (CBC) [9,38] is used for its good performance in multicomponent flow simulations. In the present model, the CBC is harnessed for two types of quantities at the right outlet. One is for the unknown distribution functions χ i = f i , g i and h i at the outlet layer (x = N ), The other is for the macroscopic quantities, such as χ = ρ, φ, ψ and P 0 at the ghost layer (x = N + 1), which are needed for computing the derivative terms at the outlet fluid layer, Here, ζ is the characteristic velocity normal to the outlet boundary. For simplicity, we have explicitly computed χ i and χ through ζ at time t. Three common choices for ζ in CBC's are the average velocity (CBC-AV), local velocity (CBC-LV) and the maximum velocity (CBC-MV) [38].
where y 1 = 0.5 and y 2 = L y − 0.5 are the locations of the bottom and top walls. All three options of the CBC mentioned above are implemented at the right outlet, and their accuracy is quantified using the relative velocity error computed by ana is the analytical velocity given by Eq. It is seen that all three outlet boundaries give satisfactory results for flow far away from the outlet. However, the accuracy at the outlet layer varies: the CBC-AV provides the highest accuracy, CBC-MV is slightly lower and CBC-LV shows the poorest performance.
In the moving droplet test, a droplet with radius R = 20 is centered at (60, 49.5) in a channel of L x × L y = 199 × 99, as illustrated in figure 2 (a). The two fluid phases have the same viscosity of 0.167 and their interfacial tension σ is 0.005. All the boundary conditions are the same as those in the single-phase Poiseuille flow simulations. The whole fluid domain is initialized with a uniform parabolic velocity profile as given by Eq. (22). Three different values of u max are tested, i.e., u max = 1.5 × 10 −3 , 3.0 × 10 −4 and 7.5 × 10 −5 . To make a quantitative comparison, the time history of the distance X d measured from the inlet to the leftmost point of the droplet is recorded and shown in figure 2 (b1-b3). The X d and time are normalized using X * d = X d /D and t * = tu max /D, where D is the droplet diameter. The X * d curve of the droplet moving in a longer channel (L x × L y = 399 × 99) computed with CBC-AV is used as the reference result for each flow condition. Note in figure 2 (b1-b3) that the sharp decrease of X * d occurs when the droplet completely moves out of the channel.
It is seen in figure 2 (b1-b3) that the X * d increases linearly with time and agrees with the reference line before the droplet interface touches the outlet boundary for each of the tested flow conditions. The option of the CBC's has little effect on the flow behaviors away from the outlet. Deviations in X * d curves from the reference lines occur at around t * = 4 when the droplet passes through the outlet. Compared to the reference lines, the case with CBC-AV slightly lags behind, and the case with CBC-MV moves a bit faster. Also, the case with CBC-LV gives the most accurate results for moderate characteristic velocities, as illustrated in figure 2 (b1)-(b2). The deviation in X * d increases as u max decreases for the cases with CBC-AV and CBC-MV. When the u max is on the same order of magnitude as the spurious velocities of the present model, i.e., u max = 7.5 × 10 −5 in (b3), numerical instability arises for the case with CBC-LV, whereas the cases with CBC-AV and CBC-MV show better robustness. Due to the low velocity often encountered in double emulsion generation, the robustness of the outlet boundary at low velocities is of great significance. On the other hand, for low velocity cases shown in (c2)-(c3), the velocity in regions close to the walls is less affected for the case with CBC-AV than that with CBC-MV. In addition, CBC-AV shows higher accuracy in the single-phase Poiseuille flow simulations. Therefore, the CBC-AV is used in the following studies.

B. Morphology diagram
Since the interfacial tension relations are crucial in determining ternary emulsion morphologies [21,50], another validation test is conducted to show the capability of the current model in simulating a wide range of interfacial tension ratios. Following the theoretical analysis of Guzowski et al. [21], two equal-sized red and green droplets are initially put next to each other and dispersed in the outer blue fluid. Three typical thermodynamic equilibrium morphologies can be obtained depending on the interfacial tension ratios of σ gb /σ rg and σ rb /σ rg , as divided by the solid lines shown in figure 3 (a): (I-A) complete engulfing with the red droplet inside the green one for σ rb /σ rg > 1 + σ gb /σ rg ; (I-B) complete engulfing with the green droplet inside the red one for σ gb /σ rg > 1 + σ rb /σ rg ; (II) non-engulfing, for σ rb /σ rg + σ gb /σ rg < 1, where the red and green droplets tend to separate from each other; (III) partial engulfing (Janus droplet), for |(σ rb /σ rg ) − (σ gb /σ rg )| < 1 and σ rb /σ rg + σ gb /σ rg > 1, where the interfacial tensions satisfy a Neumann triangle.
In our numerical test, the red and green droplets are both initialized with radius R = 60 surrounded by the blue fluid in a domain of L x × L y = 399 × 399. All the fluid viscosities are 0.167. The initial concentration fractions for three fluids are given by [66] Periodic boundary conditions are used for all boundaries. To reproduce all the possible morphologies, simulations are performed at various groups of (σ gb /σ rg , σ rb /σ rg ): (I-A) complete engulfing with red droplet inside: It is worth noting that we have investigated the optimal value of the coefficient β in the additional free energy term introduced in Eq. (8), varying β = 0, 0.0001, 0.001, 0.01, 0.1 and 1.0 for one typical double emulsion morphology at (σ gb /σ rg , σ rb /σ rg ) = (1.4, 0.35). The obtained result at β = 0 (corresponding to the model without the additional term) is shown in figure 3 (b). As highlighted by the dashed squares in figure 3 (b), two unphysical light blue regions caused by negative C 1 appear around the three-phase contact line and lead to incorrect result. The incorrect region is also observed for β = 0.0001 in figure 3 (c). For β varying from 0.001 to 0.1, the complete engulfing morphology could be successfully reproduced and invisible difference is observed for different values of β. However, further increasing β to 1.0 induces numerical instability, which indicates that the β value cannot be large enough to dominate the doublewell potential terms. Meanwhile, for the partial engulfing cases, correct morphologies could be captured even without the additional term, and they are generally unaffected by a small additional term. Based on the above findings, β = 0.001 will be used in subsequent simulations.

A. Previously observed formation regimes and grid independence test
The two-dimensional setup of the hierarchical flow-focusing device is illustrated in figure 4. The inner red fluid is injected through the leftmost inlet with a width of w 1 , and the middle green and outer blue fluids are injected by two vertical side inlets with widths of w 2 and w 3 , respectively. All the inlet widths are set equal in this section, i.e., w 1 = w 2 = w 3 . The channel connecting the two side inlets has a width of w 4 = w 1 , and the main channel width is w 5 = 1.6w 1 . The length of the first inlet is w 6 = 2w 1 , and the distance between the two side inlets is w 7 = 3w 1 . Considering the symmetry of the flow problem in the y direction, only a half of the geometry is simulated and the domain size is L x × L y = 30w 1 × 2w 1 . Zou-He velocity inlet boundary condition [70] is used for all the inlets, and the CBC-AV is applied for the outlet. In addition to the no-slip boundary condition, the wetting boundary condition is also imposed on the solid surfaces, where the first and second junctions are fully wetted by the middle and outer phase fluids, respectively.
In the following, the subscripts i, m and o are used to denote the inner, middle and outer fluids. The dimensionless parameters that characterize the double emulsion formation process could be defined as follows [2,4]: the Weber number (the ratio of inertia force to interfacial tension force) of the inner fluid W e i = ρ i u 2 i w 1 /σ im ; Capillary numbers (the ratios of viscous force to interfacial tension force) of the middle and outer fluids and interfacial tension ratios σ io /σ im and σ mo /σ im . Here, u is the average inlet velocity. Among these parameters, W e i , Ca m and Ca o are the most important parameters [2,4]. In the current study, we will change the values of W e i , Ca m and Ca o , and investigate their roles in the formation regime conversions and the obtained double emulsion sizes. The values of W e i , Ca m and Ca o are varied by adjusting the flow rate of each phase. The viscosity ratios are kept at λ im = λ om = 1.24, and the interfacial tension ratio is fixed at (σ mo /σ im , σ io /σ im ) = (1.0, 2.2).
Four basic flow regimes identified in the double emulsion preparation [2,4] are shown in figure 5 (a1-a4): (a1) the two-step formation regime, (a2) one-step formation regime, (a3) decussate regime with one empty droplet, and (a4) decussate regime with two empty droplets. Our simulations are able to successfully reproduce all four regimes. Specially, the two-step formation regime is obtained at W e i = 0.0023, Ca m = 0.011 and Ca o = 0.035 ( figure 5 (b1)). With the same Ca m and Ca o values, the one-step formation regime is observed by increasing the inner flow rate to W e i = 0.0053 (figure 5 (b2)), while the decussate regime with one empty middle phase droplet is achieved by decreasing the inner flow rate to W e i = 0.0010 ( figure 5 (b3)). Moreover, if the decussate regime happens at higher Ca o values, e.g., Ca o = 0.065 with W e i = 0.0010 and Ca m = 0.011, two empty alternate middle phase droplets are found, as shown by figure 5 (b4).
A grid independence test is conducted for the two-step formation regime mentioned in figure 5 (b1). Four different grid resolutions are tested, i.e., w 1 = 30, 40, 50 and 60. To make a quantitative comparison, the results from the highest grid resolution (w 1 = 60) is used as a reference. The relative errors (E w1 = |X w1 − X w1=60 |/X w1=60 ) of the overall emulsion size, the pinch-off location and the generation period are calculated, and their maximum value is recorded for each grid resolution. The maximum relative errors for w 1 = 30, 40 and 50 are 11.3%, 6.60%, and 2.48%, respectively. This suggests that increasing grid resolution from w 1 = 50 to 60 leads to the relative error less than 3%, and thus an inlet width of w 1 = 50 is used for the following studies.

B. Effect of flow rates
In the formation of double emulsions, it is known that two-step, one-step and decussate formation regimes can be obtained by varying W e i , Ca m and Ca o values. However, the dependence of each formation regime on W e i , Ca m and Ca o values has not been systematically studied, and how the obtained emulsion sizes vary is not very clear. In figure  6, a three-dimensional phase diagram is constructed to illustrate how the formation regimes are influenced by W e i , Ca m and Ca o . The ranges for these influencing parameters are W e i = (0.0010, 0.0016, 0.0023, 0.0032, 0.0042, 0.0053, 0.0065, 0.0079, 0.0102, 0.0127, 0.0146), Ca m = (0.005, 0.011, 0.015, 0.02, 0.025, 0.03) and Ca o = (0.025, 0.035, 0.05, 0.065). It is seen that more formation regimes are obtained besides those reported in figure 5. To differentiate these regimes, each regime is represented by a unique symbol. The nomenclature for each regime is generally a combination of the breakup modes of the inner and middle phases. To distinguish the dripping and jetting breakup modes, we use the pinch-off location. It is considered as jetting if the distance between the pinch-off location of the inner (or middle) fluid and the downstream edge of the middle (or outer) fluid side inlet is longer than 3w 1 [58]. Otherwise, it is categorized as dripping. According to our nomenclature, the periodic two-step and one-step formation regimes shown in figure 5 (a1) and (a2) are therefore called as Dripping-Dripping (Regime 1) and Jetting-Dripping (Regime 8) regimes in figure 6, which distinguish them from other irregular two-step or one-step formation behaviors.
In general, the phase diagram in figure 6 can be divided into two regions by the breakup mode of the inner phase fluid on each Ca o plane, which are bounded by the red solid lines. On the left side of the red solid line, where W e i and Ca m are approximately varied from 0.001 to 0.0042 and from 0.011 to 0.03, respectively, the inner phase breaks up in the dripping mode. The shape of the red solid lines varies little with Ca o over the entire range of Ca o considered, i.e., 0.025 ≤ Ca o ≤ 0.065. This indicates that the breakup mode of the inner phase is mainly affected by W e i and Ca m , but almost independent of Ca o . All the points in this region are periodic and they could be further subdivided into seven regimes, as denoted by the Regimes 1-7.
On the right side of the red solid line, for higher W e i , the breakup mode of the inner phase shifts from dripping to jetting. This region is further divided into two parts by the breakup mode of the middle phase fluid, which are 0.065, the parallel layered flow appears (Regime 11). It is difficult to make statistical studies on the irregular regimes; thus, they are out of the scope of the present work. Detailed discussions are carried out for the two periodic regions covering Regime 1 to 8 in the following sections.

Inner phase dripping
In figure 6, two types of periodic two-step formation regimes for producing double emulsion are observed, i.e., the Dripping-Dripping regime (Regime 1) and the Dripping-Jetting regime (Regime 2). They are represented by the filled light yellow circles and the dark blue circles, respectively. We observe that Regime 1 is limited to a small range of parameters. Its formation requires a delicate balance of W e i , Ca m and Ca o . Regime 2 occupies a relatively wider region. However, as Ca o increases, the applicable range of Ca m for Regime 2 shrinks to higher Ca m , due to the appearance of decussate regime at lower Ca m .
The Let us first discuss the effects of varying W e i . It is seen in figure 7 (a-i)-(a-iii) that the thick middle phase thread cannot be emulsified by the outer fluid close to the junction for small Ca o . Instead, the middle phase front carries the pre-formed inner droplet downstream and expands in vertical and horizontal directions. Particularly, the expansion in the vertical direction leads to an increased channel blockage for the upstream outer fluid. The accumulated upstream pressure in the outer fluid then assists the shear force to pinch off the middle phase front. As W e i increases, the inner droplet size decreases but its formation frequency is increased, which agrees with the trend observed by Fu et al. [18] for double emulsions produced at high middle phase flow rates in a 2D simplified geometry of a co-axial capillary device. The increased formation frequency of the inner phase droplet shortens the time to accumulate the upstream pressure in the outer phase fluid and actuates the pinch-off of the middle phase front. Therefore, the size (area) of the middle layer of the entire double emulsion also decreases with increasing W e i . Noteworthy, if the size of the inner droplet produced is too small, such as in case (a-iv), the upstream pressure in the outer fluid could not be accumulated in the less blocked channel, and the perturbations on the middle phase fluid can not grow due to the interference from the high generation frequency of the inner droplets. Thus, the Dripping-Threading regime (Regime 3) appears, where the continuous middle phase thread contains multiple inner droplets. A similar flow pattern was reported by Nabavi et al. [42] using capillary microfluidics.
Another interesting finding in figure 7 (a) lies in case (a-ii), in which variations are observed in the obtained double emulsion sizes: a smaller double emulsion is followed by a larger one, and this pattern repeats itself. To reveal the periodicity of this behavior, the temporal evolution of the inner and middle thread tip locations are tracked from the beginning (t = 0) as denoted by X i and X m in figure 8 (a1) and (a2), respectively. The time t and locations X i,m Since the breakup mode of the inner phase fluid is rarely affected by Ca o , the occurrence of the inner phase bidisperse behavior is similar to that in a binary system. It can be attributed to the oscillations in the amount of residual liquid on the entrance side after the previous droplet is pinched off [12,20,57].
Nonetheless, compared to the binary system, the influence of inner phase bidisperse behavior brings richer dynamics in the present ternary system depending on the Ca m and Ca o values. For cases like the one shown in figure 7 (a-ii), the middle phase is easily to be pinched off and it follows the bidisperse breakup frequencies of the inner phase (Regime 5). However, when Ca m is increased for the same Ca o (here for Ca o ≤ 0.35), the thicker middle phase fluid extends its pinch-off time and it can counteract the inner bidisperse frequencies by engulfing every two inner phase droplets inside when it breaks up, as shown by one typical case at W e i = 0.0023, Ca m = 0.02 and Ca o = 0.025 in figure 8 (b1)-(b2). As such, the variation in formation frequency only happens in the inner phase fluid ( figure 8 (b1)), but not in the middle phase fluid ( figure 8 (b2)). Such regime is called the In-Bidisperse regime (Regime 6). Furthermore, even if the middle phase fluid forms a continuous thread, e.g., at W e i = 0.0032, Ca m = 0.03 and Ca o = 0.025, the inner bidisperse behavior could still happen, and it is named the Bidisperse-Threading regime (Regime 7). Figure 7 (b) shows the effect of Ca m on the two-step formation behaviors. It is seen in figure 7 (b-i)-(b-iii) that the inner droplet size decreases while the middle part of the double emulsion increases as Ca m increases. As the intermediate layer, the middle phase fluid has dual effects. With increasing Ca m , the increased viscous force of the middle phase fluid acting on the inner phase fluid leads to the size reduction of the inner droplet. At the same time, the increased middle flow rate decreases its velocity difference to that of the outer phase, which effectively lowers the outer shear stress and extends the time for the middle part of the double emulsion to grow larger. The results show that the increase in the middle part size is usually more significant than the decrease in the inner droplet size. Thus, the overall double emulsion size also increases. When Ca m is increased to 0.03 in figure 7 (b-iv), another type of Dripping-Threading regime (Regime 3) is observed for small W e i values. Compared to the case in figure 7 (a-iv), the inner droplets in figure 7 (b-iv) are formed at a lower frequency and the variation in the width of the middle phase thread is more significant. Cases (b-iii)-(b-iv) in figure 7 also remind us the varicose shapes observed in a binary system [13], where the narrow main channel suppresses the varicose deformation induced by capillary instabilities when the thread is thick. The Dripping-Threading regime detected in the current hierarchical flow-focusing device shows the capability to produce bundles of microcapsules that are promising for storing, handling and arrayed assay of small volumes [48]. To remove the Dripping-Threading regime, we can increase substantial proportion in the two-step formation regions shown in figure 6. Among them, a decussate regime with one alternate empty droplet is commonly observed, and Azarmanesh et al. [4] has numerically explained this formation process in detail. In contrast, cases with two alternate empty droplets are less common. We show two examples at (W e i , Ca m , Ca o ) = (0.0010, 0.011, 0.065) and (0.0010, 0.015, 0.065). The formation details of the two representative cases are given in figure 9.
In figure 9 (a), the two empty droplets are formed one-by-one at Ca m = 0.011, which is similar to the decussate regime with two empty droplets reported by Azarmanesh et al. [4]. A new formation process of decussate regime with two empty droplets is presented in figure 9 (b) at a higher Ca m of 0.015. A long section of the middle phase fluid is pinched off entirely, and then it breaks up into two daughter droplets during the retraction process of the stretched structure when flowing downstream. Comparing cases (a) and (b) in figure 9, the only difference lies in Ca m . A lower Ca m signifies a higher velocity difference between the middle and the outer phases, which leads to a stronger shear force exerted on the middle phase and contributed to the early pinch-off of the middle phase front around the bulb neck. For figure 9 (b), the middle phase front is not pinched off until the entrance of the inner droplet that prevents the continuous injection of the middle phase fluid to its thread tip. The formation processes of decussate regimes can be explored more deeply in future works due to its application potentials. For instance, if the downstream channel is connected to an expansion channel, the empty droplet can catch up with the double emulsion droplet ahead and merge to form a large double emulsion with thicker middle layer [4]. Another possibility is that the empty droplet and the double emulsion droplet can be viewed as two distinct inner components, engulfed by one more outer fluid, to produce more complex functional multiple emulsions [47].

Inner phase jetting
Even though both the two-step and one-step formation regimes can be used for producing double emulsions, the one-step regime is advantageous over the two-step regime in several aspects. Firstly, the one-step formation regime shows better robustness in wetting conditions, owing to the less possibility for the inner phase fluid to touch the first channel walls when isolated by the surrounding middle phase fluid. Secondly, a thinner shell thickness can be obtained by using the one-step formation regime. Thirdly, the one-step method is capable of producing non-Newtonian emulsions that are hard to be emulsified in a two-step regime [2]. Fourthly, from the point of view of the applicable parameter ranges shown in figure 6, a wide range of one-step formation points are continuously found, different from the periodic two-step regimes (Regime 1 and 2), with interference from other flow patterns. Therefore, the periodic one-step region obtained in the present work has more statistical significance over the periodic two-step regions, which enables us to investigate the one-step formation mechanism more quantitatively and construct possible scaling laws for the double emulsion sizes.
To better understand the one-step double emulsion formation process, the typical temporal evolution of the interface dynamics at W e i = 0.0053, Ca m = 0.011 and Ca o = 0.05 is shown in figure 10. In each sub-figure, the interface shapes are depicted by the solid lines. The leading shear force component is displayed in the upper part, i.e., τ xy = η(∂u y /∂x+∂u x /∂y) for a two-dimensional system, and it is normalized using τ * xy = τ xy w 1 /σ im . The streamlines are shown in the lower part. Figure 10 (a1) corresponds to the moment just after a previous double emulsion is pinched off, where a strong shear stress region is activated to resist the retraction of the highly deformed middle-outer interface. During the evolution from figure 10 (a1) to (a3), the middle phase thread tip approximately recovers to a semicircular shape under the effect of interfacial tension [16,57]. In the meantime, the highest shear stress is lowered, and a more evenly distributed high shear stress region is formed along the inner-middle interface. The inflation of the compound inner and middle thread tip partially blocks the inflow of the outer phase fluid. Then, the outer fluid squeezes back the expanded compound thread tip and stretches it downstream. An obvious neck region is formed in figure 10 (a4) and it keeps shrinking until the pinch-off happens in the inner phase fluid as shown in figure 10 (a5). It is seen that a higher positive and a lower negative shear stress regions are induced immediately near the newly pinched inner thread tip and the generated inner droplet, respectively. The weakly connected middle thin thread is pinched off just after the configuration in figure 10 (a6).
Based on the analysis of figure 10, the double emulsion formation process in one-step regime can therefore be approximately viewed as the sum of a partial blocking period and a squeezing period, which is analogous to that of a single droplet formation in squeezing/dripping regime in binary flow-focusing systems [13,16,37]. Specifically, to construct a phenomenological scaling law for the entire double emulsion size, we take inspiration from the binary work of Liu & Zhang [37], due to the similar formation processes in action, the simple form of scaling law proposed, and the accuracy of their prediction. In their work, the scaling law for the plug length L p is given by where the plug length is normalized by the inlet width w 1 , and˜ ,γ andm are the fitting parameters that mainly depend on the channel geometry. Q dispersed /Q continuous is the flow rate ratio between the dispersed droplet phase and the continuous carrier phase. The blocking and squeezing periods for the single droplet formation are reflected by the first and second terms in the bracket of Eq. (26), respectively. Moreover, the expression of Cam o indicates the power-law dependence of the droplet size on Ca o , which is also pointed out by Christopher et al. [11] for droplets formed in squeezing and dripping regimes.
When it comes to the current ternary one-step formation regime, we also need to highlight the differences from the binary single droplet generation process. Firstly, the dispersed phase is made up of the inner and middle phases for the continuous outer phase fluid. The independent control over the inner and middle flow properties make the flow details more complex than that of a binary system. Secondly, the droplet length is a quantity suitable to measure the plug shape droplets with diameter wider than the channel width [19,37]. However, the droplet volume is a more general quantity, owing to its applicability to measure the size of each part of the emulsion with diameter wider or smaller than the main channel [7,18,56]. Therefore, for the two-dimensional studies in the present work, we monitor the evolution of the emulsion area (denoted by A).
To develop a scaling law based on Eq. (26) for the double emulsion sizes, the roles of W e i , Ca m and Ca o in the one-step formation regime are investigated individually as shown in figure 11. The comparison with periodic two-step regimes is also discussed in the following paragraphs. In figure 11, the areas of the entire double emulsion, the inner part, and the middle part are measured after the double emulsion is produced periodically. The area quantities are normalized using A * = A/(πw 2 1 ). In figure 11 (a1), the effect of W e i is investigated for W e i = 0.0042, 0.0053, 0.0065 and 0.0079, at Ca m = 0.011 and Ca o = 0.05. Compared to the two-step double emulsion formation regime shown in figure 7 (a-i)-(a-iii), as W e i increases, the inner droplet size varies little except for an initial minor increase in the one-step formation regime, while it decreases in the two-step formation regime. On the other hand, the middle part and the entire double emulsion size both decrease in the one-step formation regime, similar to their trends in the two-step formation regime. In the one-step formation regime, the breakup time of the middle phase is dominated by the formation time of the inner phase droplet (see figure 10). As W e i increases, since the size of the inner phase droplet varies little (see figure 11 (a1)-(a2)), the formation time needed to form the inner droplet is shortened due to the increased inner velocity. Accordingly, the formation time for the middle phase part is also shortened, which leads to its size decrease at the same Ca m and Ca o values. By investigating all the periodic one-step data shown in figure  6, the variations in the inner, middle and the entire double emulsion sizes caused by W e i are qualitatively the same for other Ca m and Ca o conditions.
The effect of Ca m on the size of each component and the corresponding snapshots are illustrated in figure 11 (b1) and (b2) at Ca m = 0.005, 0.011 and 0.015, with W e i = 0.0053 and Ca o = 0.05. As Ca m increases, the inner droplet size decreases, but the middle part increases. The middle part size always increases faster than the decrease in the inner droplet size. Thus, the entire double emulsion size increases monotonously with Ca m . These trends qualitatively agree with those observed in the two-step regimes shown in figure 7 (b-i)-(b-iii) regardless of the formation details. It embodies the same effect of Ca m in both formation regimes. We further verify that varying W e i and Ca o conditions in figure 6 does not change the effect of Ca m .
We have learned the effect of Ca o on two-step formation regimes in figure 7 (c): the inner droplet size is almost independent of Ca o , but the breakup frequency of the middle phase increases with increasing Ca o , which could further lead to the decussate regime. However, for the one-step formation regime, a different effect of Ca o is expected to act on the inner and middle part sizes since the two phase fluids are emulsified simultaneously. In figure 11 (c1) and (c2), Ca o is increased from 0.025, 0.035, 0.05 to 0.065 at W e i = 0.0053 and Ca m = 0.011. As Ca o increases, identical variation trends occur to the inner part, middle part and the entire double emulsion sizes: the sizes consistently increase slightly at the very beginning and then decrease monotonously. For other W e i and Ca m values investigated in figure 6, the initial increase in sizes is not common with increasing Ca o , but the decreasing trend is always obtained due to the enhanced viscous force at larger Ca o . Therefore, for the purpose of constructing the scaling law on the double emulsion sizes, the occasional increasing trend is neglected, and we will assume the size has a decreasing trend with increasing Ca o .
Based on the above analysis, two modifications are made to Eq. (26) so that it is suitable for the entire double emulsion size produced in the one-step formation regime. Firstly, the equivalent radius defined by R emulsion = A emulsion /π is used to describe the size of the simulated ellipsoid-like double emulsion, where A emulsion is the measured double emulsion area. The radius is further normalized by the inlet width w 1 . Then, the second term in the bracket Q dispersed /Q continuous is replaced by Q m /Q i to consider the positive effect of Ca m and the negative effect of W e i on the entire double emulsion size. Thus, the scaling law for R emulsion is constructed as where the parameter values 0.270, 0.0526 and -0.268 are obtained by fitting all the investigated periodic one-step data shown in figure 6 with the principle of minimum residual norm. To test the obtained scaling law, the values of the double emulsion radius (R emulsion /w 1 ) pred computed from Eq. (27) are plotted against the simulated radius values (R emulsion /w 1 ) simu in figure 12 (a). The line of parity is plotted as a reference, and the closer the scattered data points are to the line of parity, the better the agreement is between the scaling law and the simulated results. It is seen that most of the points scatter around the line of parity, and the simple formula of Eq. (27) can provide a general guidance for predicting double emulsion size. Another size of interest is the ratio between the equivalent inner droplet radius R inner = A inner /π and the overall double emulsion radius: R inner /R emulsion , where A inner is the measured inner droplet area. Chang et al. [7] experimentally proposed a scaling law for the double emulsion generated in co-axial capillaries. The inner droplet and the entire double emulsion are viewed to have the same formation time before being pinched off together in the dripping mode. According to the mass conservation law, R inner /R emulsion is predicted by R inner /R emulsion = (Q i /(Q i +Q m )) n , and the power-law exponent n is 1/3 (or 1/2) for three (or two) dimensional studies. Recently, Fu et al. [18] numerically confirmed this relation in their two-dimensional work on the co-axial capillary device. However, in reality, the inner phase fluid breaks up slightly earlier than that of the middle phase fluid, especially in the current planar hierarchical flow-focusing device (see figure 10). The difference in formation time between the two phases is moderately affected by Ca m and Ca o values. To consider these effects, two power-law relations are thus assumed between R inner /R emulsion and the capillary numbers, Ca m and Ca o . A constant scale factor is also added to the entire equation. Thus, the scaling law for R inner /R emulsion is constructed as The way to obtain the numerical coefficients in Eq. (28) is the same to that used in Eq. (27). The fitted power-law exponent of Q i /(Q i + Q m ) is 0.609, which is close to 0.5 mentioned in the work of Chang et al. [7] and Fu et al. [18]. The difference can be attributed to the inconsistency in the breakup time of the inner and middle phases. Nevertheless, the difference in the formation time is small, which is also reflected by the scale factor 0.904 that is close to 1.0, and the near zero power-law exponents for Ca −0.060 m and Ca 0.030 o . Similar to figure 12 (a), the parity plot for the computed values of (R inner /R emulsion ) pred from Eq. (28) and the measured values (R inner /R emulsion ) simu are shown in figure 12 (b). The good agreement between the scattered points and the parity line justifies the validity of the scaling law of Eq. (28) for the R inner /R emulsion values.

C. Effect of interfacial tension ratio
In figure 3, we show that a variation in the interfacial tension ratio could result in distinct equilibrium morphologies of two droplets of different fluids. To elucidate the role of interfacial tension ratios on the emulsion structure in different double emulsion formation processes, six groups of interfacial tension ratios that cover different regions of : 1 and 0.286 : 0.390 : 1. To obtain different interfacial tension ratios, σ im is fixed at 0.005 except for the case at (σ mo /σ im , σ io /σ im ) = (100, 100), where σ im = 0.00001 is used, similar to those used in figure 3. Figure 13 illustrates the snapshots of the (a) two-step and (b) one-step flow rates for each interfacial tension ratio group. Note that the first column before (a) series shows the corresponding static equilibrium morphology of each interfacial tension ratio group as shown in figure 3. Relevant experimental works are marked next to the related snapshots.
It is seen in figure 13 that the formation details and the emulsion morphologies are greatly affected by the interfacial tension ratios in both formation regimes. Firstly, compared to the double emulsions obtained at (σ mo /σ im , σ io /σ im ) = (1.0, 2.2) (figure 13 (a1) and (b1)), the inverse engulfed double emulsion is captured in figure 13 (a2) and (b2) by reversing the interfacial tension ratios to (σ mo /σ im , σ io /σ im ) = (2.2, 1.0). With the inverse interfacial tension ratios, the inner phase fluid is more favored to the outer phase fluid and tends to engulf the middle phase droplet to lower the system's interfacial energy. In the two-step formation regime shown in figure 13 (a2), as the individually generated inner droplet approaches the second cross junction, it is getting closer to the middle-outer interface. Once the inner droplet touches the middle-outer phase interface, the attraction between the inner and outer phases would prompt the pinch-off of the middle phase layer between them and actuate the formation of the middle phase droplet. Afterwards, the inner droplet itself becomes a bridge connecting the newly formed middle phase droplet and the remaining middle phase front. Soon it breaks into two parts under the viscous force of the outer fluid. The inner phase portion adhered to the middle phase droplet evolves to wrap the middle phase droplet and the inverse double emulsion morphology is finally formed. In the one-step formation regime of figure 13 (b2), the inverse double emulsion is also obtained. However, the formation details are different due to the continuous supply of the inner phase fluid in the jetting mode. A string of small middle phase droplets are formed and connected by the inner phase fluid. The compound thread tip is then emulsified by the outer fluid for every two front middle phase droplets. The detached two middle phase droplets covered by the inner phase fluid soon merge with each other and produce a pure double emulsion.
Chao et al. [8] experimentally captured the conversion from an initial double emulsion to its inverse structure using a glass-based capillary microfluidic device. Using the terminology of our work, an intermediate red-in-greenin-blue double emulsion is initially produced in their work, and the thermodynamic equilibrium green-in-red-in-blue configuration is only obtained after the external flow is stopped. However, in our work, the final configuration is formed directly without the intermediate red-in-green-in-blue configuration. This implies that the moment for interfacial tension dominating over the hydrodynamic effects in the formation behaviors is earlier in our simulations than that in the experimental work of Chao et al. [8]. This could be explained by the experimental findings of Pannacci et al. [50]. They pointed out that it is necessary for the inner droplet to touch the inner boundary of its host to evolve to thermodynamic equilibrium under the capillary forces. In other words, the sooner the three-phase contact line is formed, the faster the interfacial tension starts to dominate. For instance, if we look into the formation details in figure 13 (a2), there should be an instantaneous moment, like highlighted in the square region in figure 13 (a1), where the inner droplet is approaching the middle-outer interface due to the squeezing of the outer fluid. It allows the capillary force to act earlier. Regarding the experimental work of Chao et al. [8], a relatively thick middle layer surrounds the inner phase orifice in the co-axial glass capillaries, which could prevent the early formation of the three-phase contact line, and hence delay the interfacial tension effect.
At (σ mo /σ im , σ io /σ im ) = (0.48, 0.48), the red and green droplets tend to separate with each other at thermodynamic equilibrium. In the two-step formation regime shown in figure 13 (a3), the inner and middle phase droplets are successively formed and flow downstream without touching each other in the outer fluid, consistent with their static equilibrium morphologies. These alternately generated single droplets of two phase fluids have possible applications in being the source materials for producing multi-core emulsions [47]. In figure 13 (b3), a more complex multiple emulsion is obtained in the one-step formation regime: an inner phase droplet is seized by two middle phase droplets on both sides in the flow direction. The contact length between the components of the multiple emulsion is decreasing when flowing downstream, but the components do not completely separate from each other in the finite computational domain. It can be attributed to two possible reasons. The first one is that the sequence structure results from a transient double emulsion rather than separately produced like in the two-step formation regime. Thus, it takes longer for the sequence structure to evolve to its thermodynamic equilibrium. Secondly, once the middle phase thread tip is pinched off, the lateral outer phase fluid rapidly fills the pinch-off region. Consequently, the most upstream component in the sequence is more accelerated and the hydrodynamic effects keep the three components staying next to each other. The complete separation of the components could be expected after the inflow pumps are stopped.
For the three cases shown in figure 13 (a4-a6) and (b4-b6), since the interfacial tensions in each case satisfy the Neumann triangle relation, the partial engulfing (Janus) emulsion should be achieved at thermodynamic equilibrium. Zhang et al. [68] experimentally captured the transformation from the core-shell structure to the Janus droplets based on prefabricated double emulsions. Here, our results in figure 13 (a4) show that the Janus droplet could be produced directly in the two-step formation regime within the same device for producing double emulsions. For figures 13 (a5), (b4), and (b5), biconcave and biconvex emulsions are formed downstream. These structures are analogous to those experimentally fabricated by Nisisako et al. [44] and Nie et al. [43]. Finally, for (σ mo /σ im , σ io /σ im ) = (100, 100) shown in figure 13 (a6) and (b6), σ im is so small that the inner and middle phase fluids can be approximately viewed as the same fluid, and the high W e i induced by small σ im easily leads to the parallel layered flow behaviors for both the two-step and one-step flow conditions.

D. Effect of geometry
Geometrical parameters in microfluidics are usually the key factors in single or double emulsion preparations [37,41,65]. In this section, we focus on the effect of the geometrical parameters in changing the double emulsion formation regimes and the obtained double emulsion sizes. For the geometry shown in figure 4, six normalized geometrical parameters can be defined as w 2 /w 1 , w 3 /w 1 , w 4 /w 1 , w 5 /w 1 , w 6 /w 1 and w 7 /w 1 . Among them, the inlet length w 6 /w 1 can be neglected, since the fully developed velocity distribution is always provided at the inlet, and the inner phase flow profile varies little before it reaches the middle phase inlet junction. Then, for simplicity, we make two assumptions to reduce the governing geometrical parameters, i.e., the side inlets for the middle and outer phase fluids have equal widths (w 3 = w 2 ), and the width of the channel connected the side inlets is set equal to that of the inner phase inlet (w 4 = w 1 ). Therefore, the main geometrical factors are reduced to the side inlet width (w 2 /w 1 ), the main channel width (w 5 /w 1 ) and the distance between the side inlets of the middle and outer phase fluids (w 7 /w 1 ). Those geometrical factors are all investigated at two flow rates that lead to two-step and one-step formation regimes, respectively, for the original geometry. Different from the flow conditions used in the interfacial tension effect section, two closer W e i values of 0.0032 and 0.0042 are used in this section at Ca m = 0.011 and Ca o = 0.035, to show the geometrical effect more obviously in changing the formation regimes.
The effects of w 2 /w 1 , w 5 /w 1 and w 7 /w 1 on double emulsion formation behaviors are illustrated in figure 14. The (a) and (b) series correspond to the two-step and one-step flow rate conditions, and each parameter of concern increases from top to bottom in each sub-column. The results for the original geometry used in previous sections are marked with an inverted triangle. In figure 14 (a1) and (b1), w 2 /w 1 is increased from 0.8, 1.0, 1.2 to 1.4 at w 5 /w 1 = 1.6 and w 7 /w 1 = 3.0. It is seen that the breakup mode of the inner phase apparently changes from the jetting mode to the dripping mode with increasing w 2 /w 1 at both flow rates. A larger w 2 /w 1 is required to induce the inner breakup mode transition at higher W e i values. Theoretically, increasing the side inlet width increases both the viscous force and the inertia force of the side-injected fluids, and the two enhanced forces act together to overcome the unaltered interfacial tension force, which leads to the breakup mode transition of the inner phase fluid. Figure 14 (a1-ii)-(a1-iv) and figure 14 (b1-i)-(b1-iii) illustrate the effect of increasing w 2 /w 1 on emulsion sizes in the two-step and the one-step formation regimes. The size of the middle part increases in both formation regimes. However, the inner droplet size varies little in the two-step formation regime but decreases in the one-step formation regime.
The effect of the main channel width is studied for w 5 /w 1 = 1.  . It is seen that decreasing w 5 /w 1 does not affect the formation regime of the inner phase fluid, but it could increase the breakup frequency of the middle phase fluid and induce the decussate regime, as observed in figure 14 (a2-i, b2-i). For the flow-focusing geometry, all three inflow fluids converge to the main channel. Thus, narrowing the width of the main channel (w 5 ) increases the fluid velocity in the axial central region of channel, which creates a larger velocity gradient in the direction perpendicular to the main flow. During the expansion of the middle phase thread tip, it is subject to a higher shear stress, and as such the middle phase is more more likely to break up. With increasing w 5 /w 1 , the inner part and the entire double emulsion size vary little in the two-step formation regime (see figure 14 (a2-ii)-(a2-iv)), but they both increase in the one-step formation regime (see figure 14 (b2-ii)-(b2-v)). It indicates that the main channel width has a more obvious effect on the size of double emulsions generated in the one-step regime. Additionally, emulsions with two inner droplets are regularly obtained in the two-step formation regime at a wider collection channel, i.e., w 5 /w 1 = 2.0 (see figure 14 (a2-v)), similar to those experimentally captured in a double cross-junction device [15] and capillary devices [34,42].
At last, the distance between the two side inlets is investigated at w 7 /w 1 = 1.0, 2.0, 3.0, 4.0 and 5.0, w 2 /w 1 = 1.0 and w 5 /w 1 = 1.6, as shown in figure 14 (a3) and (b3). The two-step formation regime shifts to the one-step formation regime at w 7 /w 1 = 1.0 ( figure 14 (a3-i)), where the inner phase front reaches the second junction before it breaks up in the dripping mode. However, more generally, the breakup modes and double emulsion sizes vary little with w 7 /w 1 in both formation regimes, similar to the findings in binary systems using flow-focusing type geometries [57,65]. Even though the lengthening of the connection channel increases the flow resistance through it, the flow behaviors inside vary little due to the slightly affected viscous force. As such, the velocities of the inner and middle phases are almost unaffected when they flow into the outer phase junction. Therefore, the overall flow behaviors are almost unchanged. Noteworthy, satellite droplets appear at w 7 /w 1 ≥ 4.0 in the two-step flow regime, due to the highly stretched middle phase thread tip during the emulsification process. It suggests that narrowing the distance between the side inlets could be a possible solution to avoid satellite droplets in producing double emulsions.

V. CONCLUSIONS
In this work, a two-dimensional ternary free energy lattice Boltzmann model is developed and used to systematically study the double emulsion formation behaviors in a planar hierarchical flow-focusing channel under variations of the flow rate, interfacial tension ratio and geometrical settings.
The periodic two-step, one-step and decussate double emulsion formation regimes previously reported in the literature are qualitatively reproduced. A three-dimensional phase diagram is then constructed to show the distribution of each formation regime governed by W e i , Ca m and Ca o values. Depending on the breakup mode of the inner and middle phases, three distinct domains are classified as the periodic two-step, periodic one-step and non-periodic regions. The range for the periodic two-step region is almost unaffected by Ca o , and it can be subdivided into seven formation regimes according to the pinch-off locations and the uniqueness of formation frequencies. Among them, periodic double emulsions are produced in Regime 1 and 2. In these two regimes, the entire double emulsion size decreases with W e i , increases with Ca m , and varies little with Ca o . Dripping-Threading regime (Regime 3) occurs when the middle phase fluid forms a continuous protective layer and carries multiple inner droplets. Decussate regimes (Regime 4) with one or two alternate empty droplets are both obtained. Noteworthy, the two empty droplets in the decussate regime could be produced either in a one-by-one sequence, or by breaking an initially formed large empty droplet into two daughter droplets. The bidisperse behaviors in double emulsion size and formation frequency are captured in a certain range of W e i values in the two-step formation regime. The bidispersity could exist simultaneously for both the inner and middle phase fluids (Regime5), or only occur to the inner phase fluid (Regime 6 and 7). In the periodic one-step region for double emulsions (Regime 8), the entire double emulsion size is found to decrease with W e i and Ca o , but increases with Ca m . Compared to the two-step formation regime, Ca o has a more obvious effect on the size of double emulsions formed in the one-step regime. Based on the one-step data (Regime 8), two empirical scaling laws are constructed for the size of the entire double emulsion and the proportion of the inner droplet. The good predictions of both scaling laws justify that the one-step formation process of double emulsions can be analogously viewed as a sum of a blocking period and a squeezing period.
Another contribution of this work is that the presented free energy model is capable of dealing with a wide range of interfacial tension ratios, and provides accurate results for predicting complete engulfing double emulsions, partial engulfing Janus droplets and non-engulfing separate droplets. In particular, it was necessary to include an additional free energy term to capture the complete engulfing double emulsions. In the current microfluidic device, a variation in the interfacial tension ratios leads to distinct emulsion morphologies, including the inverse engulfing double emulsions [8], non-engulfing single droplets [47], Janus droplets [68], biconcave and biconvex emulsions [43,44], and even parallel flows.
Regarding channel geometrical parameters, the breakup mode of the inner phase fluid is changed from dripping to jetting by decreasing the side inlet width w 2 /w 1 , or by narrowing the distance between the two phase side inlets w 7 /w 1 . This leads to the conversion from the two-step formation regime to the one-step formation regime. The main channel width w 5 /w 1 should not be too small in order to avoid the decussate regime. Moreover, narrowing w 7 /w 1 is a possible solution to get rid of the satellite droplets for double emulsions generated in the two-step regime. The entire double emulsion size increases with w 2 /w 1 , but is rarely affected by w 5 /w 1 or w 7 /w 1 in the two-step formation regime. For the one-step formation regime, the double emulsion size increases with w 2 /w 1 and w 5 /w 1 , but is independent of w 7 /w 1 .
Finally, we would like to point out that the above work is carried out in a two-dimensional scheme. Based on the fundamental knowledge achieved in the present work, a three-dimensional study can be expected in the near future. Furthermore, equal density fluids are used at present. Our newly developed high-density ternary free energy model [64] could be applied to investigate double emulsion formation behaviors with other fluid types, such as gasin-oil-in-water. It is also worth extending the current ternary free energy model to deal with multiple emulsions with more components (N > 3), or introducing variable interfacial tensions governed by the surfactants [36] to study more complex fluid systems.