Discrete fluidization of dense monodisperse emulsions in neutral wetting microchannels

The rheology of pressure-driven flows of two-dimensional dense monodisperse emulsions in neutral wetting microchannels is investigated by means of mesoscopic lattice simulations, capable of handling large collections of droplets, in the order of several hundreds. The simulations reveal that the fluidization of the emulsion proceeds through a sequence of discrete steps, characterized by yielding events whereby layers of droplets start rolling over each other, thus leading to sudden drops of the relative effective viscosity. It is shown that such discrete fluidization is robust against loss of confinement, namely it persists also in the regime of small ratios of the droplet diameter over the microchannel width. We also develop a simple phenomenological model which predicts a linear relation between the relative effective viscosity of the emulsion and the product of the confinement parameter (global size of the device over droplet radius) and the viscosity ratio between the disperse and continuous phases. The model shows excellent agreement with the numerical simulations. The present work offers new insights to enable the design of microfluidic scaffolds for tissue engineering applications and paves the way to detailed rheological studies of soft-glassy materials in complex geometries.


Introduction
Dense emulsions and suspensions of soft particles are widespread in natural as well as engineering applications, ranging from biofluids and food to pharmaceutical and cosmetic industrial processes 1,2 . A deeper understanding of the basic rheology of these soft glassy materials (SGMs) is therefore beneficial to the advancement of many fields of science and technology 3,4 .
The flow characteristics of these systems depend, often in a not trivial way, on many parameters, such as volume fraction, interfacial properties, deformability of the elementary, constituents, confinement, etc. Consequently, depending on the load conditions, SGMs display a variety of non-Newtonian features like yielding, shear thinning/thickening, rheopexy and thixotropy that make their rheology an outstanding open issue in non-equilibrium staa Center for Combustion Energy; Key Laboratory for Thermal Science and Power En-tistical physics 1 . During the past decades, significant attention has been paid to the microstructures of colloidal crystals and foams and their transformations under shear stress [5][6][7][8][9] . For example, Chen et al. studied structural changes and orientational order for dense colloidal suspensions under shear 5 . At low shear, the lattice structure will change from a crystalline state to a polycrystalline state, while at large shear, they reported the appearance of a sliding layer flow 5 . Debregeas et al. investigated the deformation and flow of 2D foam (composed of jammed bubbles) under continuous shear in a Couette geometry 7 . They found rapid decay of the average velocity from the moving inner wall to the fixed outer wall, as well as large velocity fluctuations with self-similar dynamical structures inbetween. They also developed a stochastic model relating the plastic flow to the stress fluctuations. Cicuta et al. reported rheological measurements on dense soft glass with an interfacial stress rheometer (ISR) 8 . They found an interesting response of the elastic modulus 8 : a transition from viscous liquid towards an elastic solid with the increase of concentration, and an inverse transition (from elastic to viscous) with the increase of shear frequency.
The shear rheology of a wide class of SGMs [10][11][12] can be described by the Herschel-Bulkley relation 13 between the applied stress σ and the responsive shear rateγ, In Eq. (1), σ Y is the yield stress, below which the material does not flow, instead it reacts elastically (showing solid-like behaviour) to the imposed load; for σ > σ Y the material flows, in general with a non-Newtonian constitutive law characterized by the two phenomenological parameters K (the consistency factor) and c (the flow index) that depend on the detailed microstructure of the material 2 . For c < 1 the material is said to be shearthinning, whereas c > 1 corresponds to shear-thickening. For dense suspensions of solid particles the phenomenology can be even richer: a number of rheometric studies, provided evidence of flow curve hysteresis and a discontinuous shear-thickening [14][15][16][17] (see also 18 and references therein). On the other hand, discontinuous shear-thinning has remained so far more elusive. Chen and Zukoski 19 observed, in constant-stress rheometry of concentrated suspensions with crystalline order, and well beyond yield, that a critical stress exists at which the power law index of the flow curve changes. Similar results have been obtained in numerical simulations of adhesive dispersions of soft disks 20 . It must be stressed, though, that the occurrence of such discontinuous shearthinning/thickening behaviours can be assessed in shear rheometry, where the shear stress is homogeneous across the sample, and for values of the applied load above the yield stress, but not in microchannel flows, unlike what was recently proposed in 21 . Different from the previous study within the Couette geometry 5,7,8,19,20 , the present work focuses on the deformation and flow of two-dimensional (2D) emulsions in the Poiseuille geometry. When a pressure difference is imposed along a neutral wetting microchannel, the shear stress varies linearly from wall to wall, being maximum at the sidewalls and zero in the channel center. With the increase in pressure difference, it is possible to achieve a fluid-like state near the wall, i.e., σ w > σ Y , where σ w indicates the wall stress. However, if the material develops a nonvanishing yield stress σ Y , there will always be a region, in the middle of the channel, where the emulsion is below yield and must be expected to move as a plug. For the meso-constituents near the wall, different from sliding on a non-wetting microchannel 22 , they are expected to stick on the neutral-wetting wall. In such a setup, the way a soft-glassy material fluidizes, i.e. the plug-to-fluid flow transition, still lacks a full understanding. As we mentioned, the microstructure in which the meso-constituents (droplets in the case of emulsions) are arranged is known to play a determining role; two structural indicators are crucial in this respect, namely the polydispersity, that measures the statistical distribution of sizes of the meso-constituents, and the degree of order in their spatial arrangement (crystalline vs amorphous) [23][24][25] . In addition, a high volume fraction Φ of the dispersed droplets is necessary to support such a microstructure transition. For example, recent studies on the T1 transition in complex microchannels are mainly based on emulsions/foams with Φ >∼ 0.85 [26][27][28] .
The rest of the paper is organized as follows. In Section 2, we give a brief introduction of the numerical model: the two-species Lattice Boltzmann (LB) model with suitable pseudo-potential interactions. The main features of the adopted numerical model are discussed and well validated. In Section 3, we prepare perfectly 2D monodisperse dense emulsions (constituted of soft droplets), with crystalline order, confined in microchannels and driven by a pressure gradient. The droplets are packed at the jamming vol-ume fraction, such that the emulsion behaves as an elastic solid for vanishingly low applied loads. We focus, then, on the fluidization of the material and find that the latter takes place as a discontinuous process with the increase of pressure gradient, characterized by successive local yielding events, where layers of droplets, starting from those closer to the boundaries, are more and more deformed and roll over each other. Scanning the imposed forcing, we observe that these events cause sudden drops of the relative effective viscosity between from one plateau to another; the relative effective viscosity is, therefore, "quantized", in the sense that it takes values only over a discrete set. The first and largest jump occurs for a critical pressure difference, which we find to be such that the wall stress equals the yield stress of the material. We show, moreover, that the value of the relative effective viscosity in the plug state depends, in a way that we characterize quantitatively, on the viscosity ratio between the disperse and continuous phases and confinement parameter (channel height over droplet radius). Finally, a brief summary and some research perspectives are given in Section 4. We perform numerical simulations based on a two-species lattice Boltzmann (LB) model, with suitable pseudo-potential interactions. Recently, the LB method has become a very efficient and powerful simulation method for complex flows [29][30][31][32][33][34][35][36][37][38][39][40] . The LB equation takes the form [41][42][43][44] :

Model description and validation
where f k,i is the probability of finding a particle of species k (k = 1, 2) at the space-time point (x,t), moving along the ith (i = 0, 1, ..., 8, see Fig. 1) direction e i of a regular D2Q9 lattice 43 . The right-hand side, computed at (x,t), represents the time relation towards local equilibrium f eq k,i on a time scale τ k and F k,i is a forcing term incorporating the effects of the total force F k acting upon each species 22 . The kinematic viscosity for each component is related to the relaxation time by ν k = (τ k − 0.5∆t)c 2 s , where ∆t = 1 and c s = 1/3 are the time step and lattice sound speed, respectively. During the evolution of Eq. (2), the density and momentum for each species are calculated as, and the total velocity of the fluid mixture is updated by u = For the two-species system, a repulsive force between species, promoting a positive surface tension, is give as usual 45 , where G kk =Gk k is the strength coefficient for the inter-component interaction, and the weights are w(0) = 4/9, w(1) = 1/9, and w(2) = 1/36. A crucial characteristic of our adopted two species LB model is the competing mechanism between a short range attraction and a mid-range repulsion within each species 46 . The short-range attraction acts between the nearest-neighbor D2Q9 lattice nodes, while the mid-range repulsion acts between the next-to-nearest-neighbor lattice nodes extending up to the D2Q25 lattice e j ( j = 0, 1, ..., 24, see Fig. 1). The competing interaction is explicitly written as 46 where G k,1 and G k,2 are introduced to tune the strengths of the two interactions, respectively, and the weights for the D2Q25 lattice are p(0) = 247/420, p(1) = 4/63, p(2) = 4/135, p(4) = 1/180, p(5) = 2/945 and p(8) = 1/15120. The pseudo-potential function originally suggested in 45,47 , ψ k (ρ k ) = ρ 0 (1 − e −ρ k /ρ 0 ) (with an identical reference density ρ 0 = 1.0 for each component) is adopted. Supplemented with the body force F b k , the total force imposed on each species is The above definition of F c k is to mimic the spatially complex (non-monotonic) interactions among molecules within each species. For example, the interplay between the competing interactions can give rise to a rich density-field configuration, which has proved fairly successful in reproducing many features of soft flowing systems, such as structural frustration, aging, elastoplastic rheology, in confined and unbounded flows 9,25,48,49 . In addition, it has been shown that by appropriately tuning the strengths of the two interactions, a positive disjoining pressure Π , which emerges as a repulsive force per unit area between opposing interfaces, can be obtained at sufficiently low film thickness 48 . Moreover, such a positive disjoining pressure can be tuned independently of the viscosity ratio between the two species, as recently proposed in 22 .
To verify the adopted LB model, we first consider two flat interfaces among three fluid domains (occupied by species 1, 2 and 1, respectively), separated by a distance h (the width of the mid-dle fluid domain). According to 48,50 , the disjoining pressure Π is defined as where γ is the surface tension and γ f is the overall line tension, defined as the integral of the mismatch between the normal (P xx ) and tangential (P yy ) components for the pressure tensor 51 , i.e., Π at a given h can be calculated using standard numerical integration methods 22 . The obtained disjoining pressure based on the above theory is shown in Fig. 2. It can be seen that for the original pseudo-potential model 45 (G k,1 = 0, G k,2 = 0, i.e., without the competing interaction), the disjoining pressure is always negative and decreases with decreasing h. By tuning the competing interaction appropriately (G k,1 = −10 and G k,2 = 8), the disjoining pressure first increases with the decrease of h and then goes down, achieving positive values in between. Moreover, the disjoining pressure profile is independent of the viscosity ratio χ between the two species within a moderate range.
To test the above theory, we then preform a numerical measurement based on the method suggested in 48 : two hemispherical droplets (species 1) are pushed against each other to form a thin film (species 2) inbetween. According to the mechanical equilibrium, the disjoining pressure in the thin film must be equal to the capillary pressure at the curved interface, defined as the pressure difference from both sides of the curved interface. By changing the radii of the droplets, we can obtain the capillary pressure at different values of the film width h. A good agreement between the theoretical disjoining pressure and the measured capillary pressure can be clearly seen in Fig. 2. It is noted that the capillary pressure curves at varying viscosity ratio χ are essentially identical and only the case of χ = 1 is shown in Fig. 2.  6) for different cases; Pink filled circle: capillary pressure curve based on the two hemispherical droplets system at G k,1 = −10, G k,2 = 8 and χ = 1; Inset: a stable two hemispherical droplets system. The positive disjoining pressure between the close-contact interfaces stabilizes the thin film inbetween. Therefore, the novelty of our improved numerical model, recently proposed in 22 , resides in a number of original features which prove key to probing the high packing regime typical of dense emulsions, in particular: (i) the emergence of a positive disjoining pressure between close-contact interfaces, which stabilizes the thin films between colliding droplets and prevents their coalescence, (ii) the possibility to tune the dynamic viscosity ratio between the two species, independently on surface tension and disjoining pressure, within a moderate range, and (iii) the highly efficient implementation of dense emulsions with up to hundreds of dispersed droplets, due to the explicit and simple algorithm for only two species probabilities f k,i (k = 1, 2).

Discrete fluidization
We prepare the monodisperse soft suspensions by packing equal sized droplets, of radius R = 18 lattice units, in a 2D microchannel, with the length and width varied in the range L ∈ [200, 800] and H ∈ [201, 804], respectively; no-slip boundary conditions for the velocity and neutral wetting conditions for the two fluids (contact angle of 90 • ) are imposed on both walls. Large values of L × H, with up to ≈ 500 droplets, are necessary to verify that the discontinuous rheological response is not an artifact due to small sizes, and let us unravel the impact of confinement. We would like to stress that such sizes would be computationally unfeasible with other LB methods where coalescence is prevented by introducing a distinct species for each droplet 21 , even equipped with advanced acceleration algorithm based on the exclusion principle 52 .
We introduce the two non-dimensional parameters that play a key role in the dynamics: the confinement parameter ξ = H/R, where R is the droplet radius (large ξ means low confinement and vice versa), and the viscosity ratio χ = µ D /µ C , where µ D and µ C are the dynamic viscosity of the droplets and continuous phase, respectively. The viscosity ratio is tuned by changing µ D with fixed µ C = 0.1. The pressure difference ∆p between inlet and outlet is imposed via a homogeneous body force ∆p/L; the lattice in-teraction parameters are chosen as, G kk =Gk k = 2.33, G k,1 = −10, G k,2 = 8, to realize a surface tension (γ ≈ 0.02), such that the corresponding Laplace pressure for droplets is p L ≡ γ/R ≈ 10 −3 (LB units) and to achieve a positive disjoining pressure, sufficient to prevent droplet coalescence (see Fig. 2). We performed a set of runs spanning the 2D parameter space (ξ , χ), at effective volume fraction of the dispersed phase Φ ≈ 0.9 53,54 .
In Fig. 3 [(a) and (b)], we show two typical snapshots of the simulated system (the density field of the dispersed phase fluid ρ D is depicted, the red/blue colours indicating high/low values; see also the Supplemental Movies S1 and S2): an emulsion with viscosity ratio χ = 4 and under low confinement, ξ = 44.7, is subject to low [panel (a)] and high forcing [panel (b)]. The corresponding streamwise velocity profilesū(y) (averaged along the x-direction and in time, over the steady state) are plotted in Fig.  3c, with blue/red lines standing for low/high pressure gradients, respectively. In both cases the profiles flatten (as compared to the parabolic profile expected for a Poiseuille flow in a simple liquid), signalling that a portion of bulk remains below yield. However, while for the lower forcing this region extends over the whole volume, such that the material moves coherently as a solid block (a "plug"), for the higher forcing the droplets in the two layers next to the walls are more deformed and the droplets in the nextto-nearest layers roll over those in direct contact with the wall periodically: in other words, the emulsion starts to be "fluidized" (see Movie S2). This is confirmed in the inset of the same figure, where we plot the velocities of the first two layers of droplets (counting from the bottom wall), called u 1 and u 2 : in the plug state these two values coincide (no deformation within the material, as it would be for a solid), whereas, at increasing the pressure difference, they bifurcate, with u 2 > u 1 and periodic velocity oscillations.
One may wonder whether other activation events of this kind will involve more and more droplet layers, for larger and larger imposed pressure gradients. This is, indeed, what happens, as we can grasp from the main panel of Fig. 4, by inspection of the velocity profiles (same system as in Fig. 3  activation of the third and fourth droplet layers (and the symmetric ones on the top wall side) that extends the region of shear localization towards the bulk.
To characterize the different flow regimes, at changing the imposed pressure drop, we define the relative effective viscosity, where Q 0 is the flow rate for the pure continuous phase for a given ∆p (defined later) and Q = H 0ū (y)dy is the emulsion flow rate. The dependence of µ r on ∆p for various viscosity ratios and confinement parameters is shown in Fig. 5: for fixed confinement and varying viscosity ratio (5(a), ξ = 11.2, and 5(b), ξ = 44.7), and for fixed viscosity ratio and varying confinement (5(c), χ = 0.5, and 5(d), χ = 4), respectively. In every data set, we observe an initial steep decrease of µ r for very low pressure differences (∆p ≤ 10 −3 , not shown in the figure), followed by a clear-cut step-wise behaviour: µ r is roughly constant, µ r ≈ µ P r (the superscript P standing for "plug state"), over a certain range of ∆p values, and then suddenly jumps to a lower plateau µ F r ("fluidized state"), as the pressure difference is increased beyond a "critical" threshold ∆p c . The applied pressure drops (on the xaxis) are made dimensionless, hereafter, as follows: where the latter expression represents the ratio between the wall stress and the Laplace pressure of the droplets. As ∆p is increased further, the relative effective viscosity develops more additional jumps (see Figs. 5 (b)and (d), as well as the inset of Fig. 4), and related after-jump plateaux. Essentially, the material responds elastically to the increasing applied load (whence the µ r -plateaux), up to a point where the stress cannot be sustained any more and the system yields (corresponding to the µ r -jumps), with successive boundary trains of droplets that start to roll over each other (see Movies S3 and S4). This phenomenology is what we dub discrete fluidization, with an accompanying discrete relative effective viscosity. Let us underline that jumps at higher ∆p are smaller and smaller, suggesting that the fluidization process is recovering continuity and eventually the material starts flowing as a viscous, non-Newtonian, fluid.  Recently, Foglino et al. also found the first and largest jump for emulsions with unit viscosity ratio in a small channel (χ = 1 and ξ = 12) 21 . In the present work, such a behaviour is further confirmed and significantly extended, and its robustness is well verified, in the sense that we find even more jumps in a larger parameter space (χ, ξ ). In addition, the first jump was interpreted as a discontinuous shear thinning behaviour 21 . We would like to stress that the shear-thinning process happens in a material which is well beyond yield, i.e., σ > σ Y , as we have discussed in the Introduction. In our simulations, the first discontinuity in the µ r vs ∆p curve, instead, coincides for all ξ 's, upon the proper rescaling Eq. (7), at roughly the same critical pressure drop ∆p ≈ 0.1. Remarkably, the latter condition is equivalent to where σ w is the wall stress. The expression at the rightmost hand side is comparable with the theoretical value, obtained by Princen 56 , of the yield stress of a perfectly monodisperse twodimensional, highly concentrated emulsion, at a volume fraction Φ ≈ 0.9. These considerations suggest that the fluidization transition occurs when the imposed pressure difference is such that σ w ∼ σ Y , i.e. when the maximum imposed stress across the material equals its yield value, lending to the idea that the origin of the discontinuity can be ascribed to local yielding events, rather than the discontinuous shear-thinning events proposed in 21 .
Taking the channel height H as the characteristic length, V = Q/H as the characteristic velocity, the Reynolds number and Capillary number can be defined, respectively, as Re = ρ D V H/µ D and Ca = µ D V /γ. For all the cases considered in the present work (i.e., Fig. 5), the Re and Ca span approximately from 0.015 to 500 and 0.0006 to 0.7, respectively. It is noted that the maximum Capillary numbers in the plug state, indicated as Ca * , for all the considered cases, almost collapse around 0.01 ∼ 0.03, as shown in Fig. 6, while we do not find a similar collapse for critical Reynolds number. This criterion is consistent with the previous study for the dislocation dynamics in 2D emulsion crystal in a tapered channel that horizontal slip planes parallel to the x axis are observed at Ca >∼ 10 −2 26 . In the experiments, the extrusion speed of the crystal needs to be sufficiently small (Ca < 10 −2 ), so that the flowing foam is in the solid-like plug state (below yield) and the spatial distributions of T1 evens remain localized in distinct rearrangement zones 26,28 . In other words, the collapse of the critical Capillary number, Ca * ∼ 10 −2 , further supports our previous argument that the plug-to-fluid transitions indicate local yielding events.

Dependence on confinement and viscosity ratio
We then move to address the confinement and viscosity ratio dependence of the rheological response more quantitatively, we focus next on the value of µ r at the first plateau, denoted as µ P r . In Fig. 7(a) we plot µ P r , averaged over the low ∆p plateau, as a function of the viscosity ratio χ, finding a linear scaling, µ P r ∝ χ, highlighted by the dashed line. This growth of the µ P r with χ is a consequence of the lower deformability of droplets, as their intrinsic viscosity exceeds that of the continuous phase. The linear scaling, which is less obvious and to which we will return shortly, is preserved also for larger system sizes ( Fig. 7(b)), albeit with a larger prefactor. Analogous plots, showing the effective plug viscosity as a function of the confinement parameter ξ , for two different viscosity ratios, are reported in Figs. 7(c-d). In either case, µ P r increases linearly with ξ , µ P r ∝ ξ , with a χ-dependent prefactor.
We now proceed to rationalize the results on the plug viscosity under the light of simple phenomenological arguments grounded on mechanical considerations. Let us first notice that, since in the plug state the whole material moves as a solid, the mass flux is Q P ≈ UH, where U is the travelling speed shared by all droplet layers and whose value is dictated by the sliding velocity of the droplets on the walls, that we need to evaluate. To this purpose, we recall that, when a pressure gradient ∆p/L is applied to a continuum system limited by solid walls (at a distance H far apart), the force balance provides a wall stress σ w = ∆pH/(2L). A droplet in contact with the wall will be, then, subjected to a force F ∝ ∆pHR/(2L), delivering a power input P in ∝ ∆pHRU/(2L); such input has to be balanced by the total dissipated power P diss inside the droplet, which can be estimated as the product of the viscous dissipation rate per unit volume ε µ D ≈ µ D (U/R) 2 and the droplet "volume" , ∝ R 2 i.e. 57,58 : Although this estimate and the one for the force acting on a droplet are based on a 2D system, it is easy to notice that the scaling relation for U remains the same in 3D. Solving the balance equation P in = P diss for U yields: delivering for the mass flux, Since the mass flux in pure continuous phase Poiseuille flow is Q 0 = (2H/3)U 0 = ∆pH 3 /(12Lµ C ), we readily obtain for the relative effective plug viscosity µ P r : Eq. (11) informs us that µ P r grows linearly with both the viscosity ratio and the confinement parameter, in agreement with the numerical results shown so far.
To check the validity of our prediction in Eq. (11), we plot in Fig. 7(e) the values of µ P r as a function of the product χξ , for all combinations of (χ, ξ ) considered in this work, which span the ranges χ ∈ [0.5; 4] and ξ ∈ [11.2; 44.7]. We observe that all points tend to collapse, over a comparatively wide range of parameter values, onto a unique master curve which is well fitted by a single linear relation µ P r = Aξ χ, in agreement with Eq. (11), with a prefactor A ≈ 0.92.

Effect of polydispersity
Although the monodisperse emulsions have been used in many microfluidic systems 26,28,59 , the realization of perfectly monodisperse emulsions are essentially impossible in experimental systems. Nevertheless, using the advanced experimental techniques 60,61 , monodisperse droplets can be generated with very small volume polydispersity, e.g., less than 3% 27,59 . To link the present work with future experiments, it is necessary to discuss the effect of small volume polydispersity on the discrete fluidization process. To this end, we initialize the monodisperse emulsions, with some perturbations for each droplet's density (ρ = ρ 0 + δ ρ). As a result, a emulsion system with finite volume polydispersity can be produced in the steady state, and a snapshot for the case with 8% polydispersity is shown in the top panel of  and χ = 4. Bottom panel: relative effective viscosity, µ r , as a function of the rescaled pressure difference, ∆p, for three different volume polydispersities. A metastable state for the 8% polydispersity case is marked by the red arrow. Inset: time evolution of the first two-layer droplets' velocities in the metastable state.
We then perform the same pressure driven simulations as before and obtain the relative effect viscosity curves for the produced polydisperse emulsions. We find that, for small and experimentally achievable volume polydispersity (4%), the presented phenomenology is well preserved, in the sense that we still observe a steep change in the relative effective viscosity around a characteristic pressure difference, but the "plug"-"fluidized" transition is gradually smoothed out with the increase of polydispersity (see Fig. 8). For example, we can see a metastable state for the 8% polydispersity case (marked by the red arrow in Fig. 8), in the sense that the first and second layer droplets coincide and bifurcate periodically (in the inset of the bottom panel) and as a result the relative effective viscosity at the specified pressure gradient switches between the "plug" and "fluidized" state intermittently.

Conclusions
In summary, monodisperse crystalline dense emulsions in a pressure-driven microchannel flow have been investigated numerically, using a new mesoscopic lattice kinetic model. Our study unveils the discrete nature of the relative effective viscosity of this kind of soft materials under confinement and rationalizes their discrete fluidizations. We showed, for the first time, that such phenomenology is preserved in large (low confined) systems and even at changing the viscosity ratio between the two phases. We also provide physical interpretation of the phenomenon that the discrete fluidization is ascribed to local yielding evens. Furthermore, we proposed a theoretical argument for the scaling relation of the relative effective viscosity of the "plug" state with confinement and viscosity ratio, and verified its validity over a wide parameter range.
We would like to stress that our results might be of more general applicability. We expect that our model of SGMs, within suitable constraints on the imposed forcing, can capture the physics of a broad family of monodisperse suspension of deformable particles (droplets and bubbles, as in foams, vesicles, etc) densely packed and confined in a microchannel.
Due to the simplicity of the configuration, it would be completely realistic to recreate the discrete fluidization in the lab, via suspensions of soft and non-coalescing droplets. In the future, it would also be informative to study the fluidization of dense emulsions or foams in more realistic geometries 27,28 .

Conflicts of interest
There are no conflicts to declare.