Mode coupling of phonons in a dense one-dimensional microfluidic crystal

Hydrodynamic interactions of solute objects in driven non-equilibrium systems lead to complex phenomena such as strong correlations in sedimenting suspensions [1], shock-waves in two-dimensional disordered ensembles of droplets [2], or cluster formation and alignment of red blood cells in micro-capillaries [3]. The dynamics of non-equilibrium many-body systems can be favorably studied in microfluidic systems using ordered arrangements of flowing droplets (microfluidic crystals) [4, 5]. In pioneering studies, Beatus et al [6–8] observed collective vibrations in trains of disk like droplets flowing in a microfluidic Hele–Shaw configuration. These oscillations have been characterized as acoustic phonons with properties similar to lattice modes in dusty-plasma crystals [9]. Their dispersion relation has been obtained to linear-order in terms of a small-amplitude expansion of the ensemble flow [6, 7]. However, in contrast to solid-state phonons that can be excited in a controlled way by mechanical or optical means, microfluidic phonons have so far been observed only when excited by spontaneous asymmetries or defects in microchannels [6–8]. Possible nonlinear effects and interactions between phonon modes could not be studied systematically in these settings [8, 10].

In this work, we broaden the possibility to study hydrodynamically mediated collective phenomena by exploring dense droplet systems, where spherical droplets are only slightly flattened and nearly in contact to each other. We demonstrate that long-living transverse phonon modes can be systematically excited. The controlled excited transverse phonon modes oscillate between the confining channel walls which define the preferred oscillation mode. Moreover, the non-vanishing transverse amplitude gives rise to nonlinearly coupled longitudinal oscillations which are beyond the acoustic longitudinal modes observed in previous works [6–8]. The experimental measurements are confirmed by results from computer simulations using a mesoscale hydrodynamics approach which is able to reproduce the full nonlinear interactions. While a linearized far-field approximation describes the transverse phonon modes of the dense droplet trains surprisingly well, the observed longitudinal modes reveal a nonlinear coupling due to the influence of the confinement on the dense droplet crystal. This type of hydrodynamic coupling has impact on the self-organization in crowded non-equilibrium systems, and is thus relevant, e.g., for controlling particle flow in applications such as microfluidic cytometry [11]. In such systems, the mild vertical confinement we explore is more realistic than the extreme confinement of [7, 8].

Microfluidic devices to explore the collective droplet oscillations were fabricated in Sylgard 184 (Dow Corning) and bonded to glass microscopy slides using standard soft lithographic protocols [12, 13]. The flow rates were volume-controlled using computer controlled syringe pumps. Mono-disperse water droplets with radius $R\approx 64\;\mu {\mbox{m}}$ are generated in n-hexadecane ($\rho =773\;{\mbox{kg}}\;{{{\mbox{m}}}^{-3}}$, $\eta =3\;{\mbox{m}}\;{\mbox{Pa}}\;{\mbox{s}}$) with 2 wt% of the surfactant Span 80 using a step-emulsification geometry [13, 14], cf figure 1(a). The spacing of the droplets in the regular trains (crystal spacing) is $a\approx 140\;\mu {\mbox{m}}$, i.e. $a<3R$. The microchannel has uniform height and width of $H\times W\approx 120\times 210\;\mu {{{\mbox{m}}}^{2}}$. The size of the droplets exceeds the height H of the microfluidic channel such that the droplets are flattened, and the drag coefficient in this case scales with ${{R}^{2}}/H$ [8]. Typical flow velocities are ${{u}_{{\mbox{d}}}}\approx 250\;\mu {\mbox{m}}\;{{{\mbox{s}}}^{-1}}$ for the droplet, and ${{u}_{{\mbox{oil}}}}\approx 500\;\mu {\mbox{m}}\;{{{\mbox{s}}}^{-1}}$ for the continuous oil phase, as observed directly from microscopy images, respectively calculated from the applied volumetric flow rates. The corresponding Reynolds and Peclet number are $\operatorname{Re}=\rho {{u}_{{\mbox{oil}}}}R/\eta \approx {{10}^{-2}}$ and $Pe={{u}_{{\mbox{oil}}}}R/D\approx {{10}^{8}}$, where $D\approx {{10}^{-12}}\;{\mbox{c}}{{{\mbox{m}}}^{2}}\;{{{\mbox{s}}}^{-1}}$ is the diffusion coefficient of a droplet.

Zoom In Zoom Out Reset image size

Computer simulations were conducted using multi-particle collision dynamics (MPC) [15–18]. In this mesoscopic method, hydrodynamics are fully reproduced through the dynamics of point-like solvent particles, cf appendix B. The flattened droplets are modelled as discs in two-dimensional channel flow [19]. The friction with the top and bottom channel walls is included as a force ${{{\mathbf{F}}}_{{\mbox{friction}}}}=-\gamma {{{\mathbf{u}}}_{{\mbox{d}}}}$ where ${{{\mathbf{u}}}_{{\mbox{d}}}}$ is the velocity of the droplet relative to the microchannel. Simulation parameters were chosen to match the experimentally observed velocity ratio $K={{u}_{{\mbox{d}}}}/{{u}_{{\mbox{oil}}}}$ for an isolated droplet. The collective oscillations are modelled by an initial zigzag arrangement of droplets with a given periodicity. Droplet and flow velocities were tuned to give a Reynolds number $\operatorname{Re}\approx {{10}^{-3}}$ in the range of the experiments. Since MPC includes thermal fluctuations, the Peclet number $Pe\approx {{10}^{2}}$ for feasible droplet and flow velocities is significantly smaller than in the experiments. This choice was dictated by a trade-off between numerical stability and computational cost. Although simulating at higher Peclet numbers is expected to reveal differences in the long term behaviour, the influence of thermal diffusion is sufficiently small to expect reliable results for short simulation times, respectively short travel distances of the droplets. Moreover, a broader range of frequencies is excited in the simulations at higher Peclet number which helps to reveal the full dispersion relation of the phonon modes as will be discussed below.

The flattened droplets experience a friction with the top and bottom channel wall and thus move slower than the surrounding liquid phase. This leads to a quasi two-dimensional dipolar flow field around each droplet, cf figure 1(b). However, note that only the part of the dipolar field is drawn that is relevant for the effective hydrodynamic forces. The droplets produced by step emulsification self-organize into a stable zigzag configuration, where the forces from the dipolar flow fields from the neighbouring droplets cancel out [6, 13]. When this droplet arrangement is guided around a 90° bend, a defect in the ordering can be achieved, as displayed in figure 1(c). In the resulting droplet arrangement the zigzag symmetry is broken and some droplets experience a net hydrodynamic force which drives them towards the opposite channel wall, see figure 1(d).

The defect in the zigzag symmetry created in this way propagates backwards in the co-moving frame of the droplets (see supplemental movie 1, available at stacks.iop.org/njp/16/063029/mmedia). If the droplet size relative to the channel width is chosen appropriately, defects can be generated periodically. The defects then initiate global sine-waves with well defined wavelength along the droplet train that travel forward in flow direction, see figure 2(a) and supplemental movies 2 and 3. The sine waves are very stable and could be observed without any change for channel lengths up to 10 cm, i.e., for distances more than two orders of magnitude larger than a typical wavelength. The wavelength λ along the droplet train depends on both the droplet size and the droplet spacing, i.e., droplet density. Varying these parameters, various wavelengths can be specifically excited, typically containing 4–9 droplets per wavelength. However, the experimentally accessible wavelength variation is limited to a factor of about three as the size of the droplets and the droplet density have to be chosen such that the droplets are close enough to affect each other by hydrodynamic interactions but do not deform each other in lateral direction by steric contact. The hydrodynamic interactions between droplets are observed to be strongly influenced by the confinement of the droplets, i.e., the degree of flattening and the absolute size of the system. The dense droplet arrangements presented in this paper could not be produced with equally sized droplets confined into channels with just $10-30\;\mu {\mbox{m}}$ height instead of the used channel height of 120 $\mu {\mbox{m}}$, as in the shallow channels the droplets are strongly repelling each other. For an increased channel height of about $150\;\mu {\mbox{m}}$ at similar vertical drop confinement, we even managed to generate a collective droplet oscillation while keeping all droplets in steric contact for all times (see supplemental movie 4). As mentioned above, the droplet friction with the channel walls and thus the relative velocity of the droplets in the continuous phase depends on confinement ratio and linear dimension and influences the strength of the dipole forces. Numerically, we account for the different velocities by using the experimentally determined velocity ratio K.

Zoom In Zoom Out Reset image size

In the following, we analyse and discuss the specific experimental and simulation results for a wavelength of six droplets with radius of $R=64$ μm and crystal spacing $a=140$ μm. The results are equally valid for other wavelengths and the discussion applies to both experiments and simulation unless mentioned otherwise. The motion of five neighbouring droplets in the co-moving frame is plotted in figures 2(b), (c). Due to the transverse confinement in the narrow channel and the longitudinal confinement by the crystal spacing, the space trajectory of each droplet describes a figure-eight and each droplet has a constant phase shift with respect to its neighbouring droplets. The transverse displacement has a comparably large amplitude which is given by the lateral confinement, given by the channel width W and the droplet radius R, $W-2R\approx 82\;\mu {\mbox{m}}$. This confinement defines the wavelength of the experimentally observed transverse droplet oscillation. Additionally, the flow velocity of the continuous phase across the channel is not homogeneous and is slower at the channel walls as a consequence of the no slip boundary conditions. Thus, the droplets also slow down in the longitudinal direction whenever they approach the channel walls and move backward in the co-moving frame. Similarly, the droplets travel forward in the co-moving frame when crossing the channel. The amplitude of this longitudinal oscillation is given by the crystal lattice, i.e. the droplet spacing a and the droplet radius R, $a-2R\approx 12\;\mu {\mbox{m}}$, which defines the stable mode in longitudinal direction and which is smaller than the amplitude in transverse direction. Thus, due to the lateral confinement, the transverse oscillation is coupled to a longitudinal oscillation where one full cycle of the longitudinal oscillation is completed while the droplet moves from one wall to the other. However, this droplet movement from one channel wall to the other only corresponds to half a transverse cycle. Accordingly, the constant phase shift between neighbouring droplets differs by a factor of two between transverse and longitudinal modes, i.e., $55{}^\circ$ and $110{}^\circ$ for the considered wave.

To analyse the droplet oscillations quantitatively, we recorded the droplet positions $\left( {{x}_{j}}\left( t \right),{{y}_{j}}\left( t \right) \right)$ for all N droplets labelled by $j=0,...,N-1$ and for T data frames obtained every $\Delta t$. From these we extract the elongations $\delta {{x}_{j}}\left( t \right)={{x}_{j+1}}\left( t \right)-{{x}_{j}}\left( t \right)$ and $\delta {{y}_{j}}\left( t \right)={{y}_{j+1}}\left( t \right)-{{y}_{{\mbox{center}}}}$, where ${{x}_{N}}\left( t \right)={{x}_{0}}\left( t \right)+Na$ and ${{y}_{{\mbox{center}}}}$ is the center of the channel. The elongations were used to calculate the power spectra of the Fourier modes

$$\begin{eqnarray}&&\tilde{x}\left( k,\omega \right)=\mathop{\sum }\limits_{t=0}^{T-1}\mathop{\sum }\limits_{j=0}^{N-1}\delta {{x}_{j}}\left( t \right)\cdot {\rm exp} \left[ -i\left( jka+t\omega \Delta t \right) \right]\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\tilde{y}\left( k,\omega \right)=\mathop{\sum }\limits_{t=0}^{T-1}\mathop{\sum }\limits_{j=0}^{N-1}\delta {{y}_{j}}\left( t \right)\cdot {\rm exp} \left[ -i\left( jka+t\omega \Delta t \right) \right]\end{eqnarray} \tag{ 2 }$$

for the transverse and longitudinal oscillations both from experimental and simulation results. The main feature of the experimentally observed stable mode are distinct peak values for both modes, cf figures 3(a), (b). The maximum of spectral density at the main peaks indicates that the collective oscillation of the droplets is dominated by one specific wavelength. The corresponding phonon spectra obtained from simulations are displayed in figures 3(c), (d). They show a continuous signature revealing the full dispersion relation with a sine-like dependence of the frequency on the wave vector. This complete signature results from the small Peclet number, i.e., from the influence of thermal fluctuations that excite all the phonon modes.

Zoom In Zoom Out Reset image size

A comparison of the experimental and numerical transverse power spectra shows that the dominant frequencies are identical. Moreover, they are in excellent quantitative agreement with the dispersion relation for transverse phonons ${{\omega }_{\bot }}\left( k \right)$ predicted by a linearized far-field theory for two-dimensional channel flow [7, 8], cf details of the calculation in appendix C:

$$\begin{eqnarray}&&{{\omega }_{\bot }}\left( k \right)=2B\mathop{\sum }\limits_{j=1}^{\infty }{\rm sin} \left( jka \right){{\left[ 1+{\rm cosh} \left( j\pi \beta \right) \right]}^{2}}{\mbox{csc}}{{{\mbox{h}}}^{3}}\left( j\pi \beta \right)\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\omega }_{\parallel }}\left( k \right)=-4B\mathop{\sum }\limits_{j=1}^{\infty }{\rm sin} \left( jka \right){\rm coth} \left( j\pi \beta \right){\mbox{csc}}{{{\mbox{h}}}^{2}}\left( j\pi \beta \right),\end{eqnarray} \tag{ 4 }$$

where $B={{u}_{{\mbox{d}}}}/{{u}_{{\mbox{oil}}}}\left( {{u}_{{\mbox{oil}}}}-{{u}_{{\mbox{d}}}} \right)\left( {{\pi }^{2}}R/{{W}^{2}} \right){\rm tan} \left( \pi R/W \right)$, and $\beta =a/W$. Note that ${{\omega }_{\bot }}\left( k \right)$ in (3) has a slightly different form than reported in [7, 8], as it takes into account all linear-order terms in the expansion, cf appendix C. The agreement of the transverse power spectrum obtained both from experiment and simulation with the analytic dispersion relation confirms that the observed transverse oscillations are acoustic phonons. From the dispersion relation, we obtain the propagation speed ${{v}_{{\mbox{p}}}}\approx 312\;\mu {\mbox{m}}\;{{{\mbox{s}}}^{-1}}$ of the phonon mode along the droplet crystal in the direction of the flow. The dispersion relation is also in agreement with the short wavelength travelling backwards as observed experimentally for a single defect. The quantitative agreement also shows that the linearized far-field hydrodynamic interactions describe the transverse phonon dynamics remarkably well, even though the droplets nearly touch each other [8]. From the experimental findings it even seems that the validity of the far-field approximation for an ensemble of flattened and comparably large droplets is extended to smaller drop distances compared to an ensemble of micro-particles with similar confinement [4, 20].

The observed longitudinal modes are not actively excited and have a smaller amplitude compared to the transverse oscillation, but still show a clear signature especially in the simulation results. The dominant frequencies derived from experimental and simulation result agree quantitatively. Interestingly, however, the longitudinal modes show a dispersion relation that differs from the prediction of the linearized theory, cf figures 3(b), (d) where ${{\omega }_{\parallel }}\left( k \right)$ from (4) is plotted as a grey line. We observe that the maximum of the longitudinal dispersion relation appears shifted towards $k=\pi /a$, and for all k, the longitudinal modes propagate in the flow direction, in contrast to the analytical prediction (4). The observed shape of the longitudinal dispersion relation resembles the scaled shape of the transverse dispersion relation. This apparent resemblance is confirmed by the correlation strength

$$\begin{eqnarray}&&C\left( {{k}_{\parallel }},{{k}_{\bot }} \right)=\frac{{{\left\langle \tilde{x}\left( {{k}_{\parallel }},t \right)\tilde{y}\left( {{k}_{\bot }},t \right) \right\rangle }_{t}}}{\sqrt{{{\left\langle |\tilde{x}\left( {{k}_{\parallel }},t \right){{|}^{2}} \right\rangle }_{t}}{{\left\langle |\tilde{y}\left( {{k}_{\bot }},t \right){{|}^{2}} \right\rangle }_{t}}}}\end{eqnarray} \tag{ 5 }$$

between longitudinal and transverse modes

$$\begin{eqnarray}&&{{\tilde{x}}_{\parallel }}\left( k,t \right)=\mathop{\sum }\limits_{j=1}^{N}\delta {{x}_{j}}\left( t \right){\rm exp} \left[ -2\pi i\frac{jk}{N} \right]\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\tilde{y}}_{\bot }}\left( k,t \right)=\mathop{\sum }\limits_{j=1}^{N}\delta {{y}_{j}}\left( t \right){\rm exp} \left[ -2\pi i\frac{jk}{N} \right].\end{eqnarray} \tag{ 7 }$$

The correlation strength is shown in figure 4(a). Whereas acoustic phonons are expected to have high correlation for ${{k}_{\parallel }}=-{{k}_{\bot }}$ [10], we observe high correlation for ${{k}_{\parallel }}\approx +{{k}_{\bot }}$. For small wavelengths the dominant correlation is located around ${{k}_{\parallel }}=2{{k}_{\bot }}$. The origin of the longitudinal modes can be explained in terms of the observed figure-eight motion of individual droplets (figure 2). The droplets slow down at both channel walls, i.e., twice during a transverse cycle such that ${{\omega }_{\parallel }}=2{{\omega }_{\bot }}$. Spatially, a full longitudinal wave extends over each crest or trough of the transverse wave, hence ${{k}_{\parallel }}=2{{k}_{\bot }}$. The combination leads to the dispersion relation

$$\begin{eqnarray}&&{{\omega }_{\parallel }}\left( k \right)=2{{\omega }_{\bot }}\left( k/2 \right),\end{eqnarray} \tag{ 8 }$$

which is plotted in figures 3(b), (d) and shows quantitative agreement with the observed longitudinal spectra. We conclude that the longitudinal modes do not originate from direct dipolar interactions between the droplets, but they are a consequence of the inhomogeneity of the imposed flow $\delta u\left( y \right)$ which is slower near the lateral walls. The longitudinal motion is hence coupled to the considerable amplitude of the actively excited transverse modes, because the spatial dependence of $\delta u\left( y \right)$ introduces nonlinear terms in the equation of motion.

Zoom In Zoom Out Reset image size

In the numerical power spectra in figures 3(c), (d) a weaker and apparently linear branch is present besides the main branch. The secondary branch is consistent with the faint secondary peaks in the experimental power spectra (figures 3(a), (b)). This branch has a positive slope for transverse modes and a negative slope for longitudinal modes. In the simulation results, the linear branches are relatively weak for short simulation times, as observed experimentally but become more pronounced for longer times, which indicates a relation to the long-time behaviour of the droplets. The space-time diagram of the numerically obtained longitudinal elongations in figure 4(b) reveals the formation of gaps and break-up of the crystal into smaller subunits with time. This is potentially caused by nonlinear interactions between longitudinal and transverse modes, similar to observations for unconfined crystals [8] and in dusty-plasma crystals [21].

A potential explanation for the emergence of the linear branches in the power spectra at longer times is the following: as long as the microfluidic crystal structure remains intact, the longitudinal positions of the droplets oscillate around their 'equilibrium position' $x_{j}^{\left( 0 \right)}\left( t \right)={{x}_{0}}\left( 0 \right)+ja+t{{u}_{{\mbox{d}}}}\Delta t$. In the co-moving frame, where the first droplet oscillates around $\hat{x}_{0}^{\left( 0 \right)}=0$, this equilibrium position is just $\hat{x}_{j}^{\left( 0 \right)}=ja$, which was already implicitly used in the Fourier transformations (1) and (2). In addition to the oscillatory motion, the droplets may diffuse away from this position over time. It may then be more appropriate to regard the droplet arrangement as a disordered ensemble, which suggests to replace the finite-differences in the equation of motion by their continuum approximations [22, 23]

$$\begin{eqnarray}&&\delta {{x}_{n+j}}-\delta {{x}_{n-j}}\approx 2ja\frac{\partial \delta {{x}_{n}}}{\partial x}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\delta {{y}_{n+j}}-\delta {{y}_{n-j}}\approx 2ja\frac{\partial \delta {{y}_{n}}}{\partial x},\end{eqnarray} \tag{ 10 }$$

where $\delta {{x}_{n}}$ and $\delta {{y}_{n}}$ are the elongations introduced above. Using this continuum approximation to derive the equations of motion and the dispersion relation as in appendix C, we finally arrive at

$$\begin{eqnarray}&&{{\omega }_{\parallel }}\left( k \right)=-4\;ka\;B\;\mathop{\sum }\limits_{j=1}^{\infty }j{\rm coth} \left( j\pi \beta \right){\mbox{csc}}{{{\mbox{h}}}^{2}}\left( j\pi \beta \right)\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\omega }_{\bot }}\left( k \right)=2kaB\mathop{\sum }\limits_{j=1}^{\infty }j{{\left[ 1+{\rm cosh} \left( j\pi \beta \right) \right]}^{2}}{\mbox{csc}}{{{\mbox{h}}}^{3}}\left( j\pi \beta \right).\end{eqnarray} \tag{ 12 }$$

For small $ka$ these expressions coincide with the expansions of (3) and (4), i.e., for small $ka$ a microfluidic crystal and a continuum show the same linearized dispersion. We note yet that the approximation in (9) and (10) demands small differences of the equilibrium positions, $ja$, and for the relatively strong lateral confinement between neighbouring droplets we study, hydrodynamic interactions are screened such that we can restrict the sum in the dispersion relation to the nearest neighbour interactions j = 1.

$$\begin{eqnarray}&&{{\omega }_{\parallel }}\left( k \right)=-4\;ka\;B\;{\rm coth} \left( \pi \beta \right){\mbox{csc}}{{{\mbox{h}}}^{2}}\left( \pi \beta \right)\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\omega }_{\bot }}\left( k \right)=2\;ka\;B{{\left[ 1+{\rm cosh} \left( \pi \beta \right) \right]}^{2}}{\mbox{csc}}{{{\mbox{h}}}^{3}}\left( \pi \beta \right).\end{eqnarray} \tag{ 14 }$$

These dispersion relations are linear in k and are plotted as dashed lines in figure 3. They match the secondary features reasonably well. Interestingly, this also seems to be the case at the edge of the Brillouin zone, which could be a consequence of short wavelength oscillations in the subunits of the droplet train. In the simulations, these features are more pronounced because the breakup of the crystal is promoted by thermal fluctuations which may trigger instabilities leading to fast growing modes [8]. However, in our experiments, the amplitudes of the longitudinal fluctuations are not growing on the observed length and time scales, hence the secondary peaks remain faint in the experimental power spectra. This experimental observation seems to indicate that oscillations in a dense microfluidic crystal with small droplet spacing are more stable than in a microfluidic crystal with large droplet spacing [5]. A detailed stability analysis is beyond the scope of this work, however, our simulation results suggest that the considered systems of microfluidic droplets are also suited to study the nonlinear mechanisms that lead to instabilities.

Using slightly flattened droplets in a quasi two-dimensional microfluidic channel, coupled phonon modes could be experimentally studied in a dense one-dimensional microfluidic crystal and quantitatively compared to simulations and analytical results. The experimentally excited oscillation modes can be selected by the lateral droplet confinement, i.e. the droplets radius and the droplet spacing. Consequently different oscillation modes and wavelength can be controllably excited when adjusting the droplet confinement. The specifically excited transverse oscillations can be described as acoustic phonons by a linearized far-field theory even in a dense system where droplets are in steric contact. A similar behaviour could not be observed in our experiments for disk like droplets, which demonstrates that the hydrodynamic interactions are strongly influenced by the dimensionality of the system. The considerable amplitude of the excited transverse oscillations leads to a nonlinear coupling of longitudinal to transverse modes. The coupling is induced by the lateral variation of the imposed flow across the channel. The coupled longitudinal modes are beyond the existing analytic description, but can be quantitatively described by a phenomenological dispersion relation.

The small longitudinal droplet spacing in the observed dense droplet crystals supports the coupling of the longitudinal and the transverse modes leading to very stable collective oscillations. Dense and laterally confined microfluidic crystals are therefore a promising system to study collective effects and nonlinear phenomena arising from hydrodynamic interactions in crowded environments out of equilibrium and under confinement. Our approach allows to gain insights into the self-organizing mechanisms and can be used to improve the performance of microdroplet applications in various areas [24].

The authors would like to thank Tsevi Beatus, Roy Bar-Ziv, Marisol Ripoll and Roland G Winkler for valuable discussions. RS and J-BF gratefully acknowledge financial support by the DFG grant SE 1118/4.

PIV allows to measure the flow of fluid media. The fluid is seeded with sufficiently small fluorescent tracers who follow the flow without delay. To determine the flow profile PIV does not track individual tracers but, relates on the numerical correlation of a large number of tracers whose position has been changed between two consecutive images (double images) which have been captured within a short time interval [25, 26]. Exposure time and time delay between two consecutive images have to be adjusted to freeze the motion of the tracers and to allow for a reliable correlation of the changed tracer position. In micro-PIV, as used in our experiments, the illumination and the detection is carried out through the same microscope objective. For the experiments presented in the paper we used a commercial micro-PIV setup from LaVision (Germany) with an illumination wavelength of 473 nm and a typical time delay of 10 ms. The measured flow field is presented in figure 1 of the paper as extracted by the commercial software DAVIS (LaVision, Germany) without further image processing.

As commercially available tracers are typically stabilized for an aqueous phase and we were interested in the flow field of the continuous phase, the PIV experiments presented in figure 1 were conducted using an oil in water emulsion. Consequently, the stable generation of monodisperse oil droplets in a surrounding water phase requires a sufficiently hydrophilic device material. Thus, the microfluidic device used for the PIV measurements has the same step geometry and channel dimensions as used for all experiments in the paper but was fabricated from the UV-curable resin NOA83H (Norland Products). The NOA83H device was molded from a PDMS master fabricated by standard soft lithography technologies. The UV curable resin was pre-cured by an uniform UV-illumination for 10 min (wavelength of 365 nm and intensity of 7.5 mW cm$^{-2}$) removed from the PDMS mold and bonded to glass microscopy slides by heating the device for 10 min at 120 °C. Prior to this step, holes were drilled in the glass slides and later connected with plastic tubing connector (Nanoport N-333. Upchurch).

Using this device monodisperse oil droplets (n-hexadecane) were generated in a surrounding aqueous phase containing 2wt% of the surfactant sodium dodecyl sulfate (Sigma) and green fluorescent micro-particles (Thermo Scientific - $G0100$). The typical diameter of the beads is around 1 μm. A suitable concentration of tracers was achieved by diluting 0,64 ml of the original particle solution with 1,65 ml ultra pure water. Flow conditions were similar for all experiments discussed in this article.

MPC is a mesoscopic simulation method that is capable of reproducing hydrodynamics of a solvent [15, 18]. The solvent is modelled explicitly by idealized point-like particles of mass m. The fluid dynamics emerges from mass, momentum and energy transport in the particle ensemble whose equation of state is that of an ideal gas.

The particle positions and momenta are updated iteratively in successive streaming and collision steps. During the streaming step the particle moves freely

$$\begin{eqnarray}&&{{{\mathbf{r}}}_{i}}={{{\mathbf{r}}}_{i}}+h{{{\mathbf{v}}}_{i}},\end{eqnarray} \tag{ B.1 }$$

where h is the time between collisions. In the collision step, the particles are sorted into cubic collision cells of size $\Delta x$. In each cell the particles then exchange momentum while the momentum of the collision cell is conserved. In this work, we use a collision rule that also conserves the angular momentum of the cell. In addition, we apply an Anderson-like thermostat to control the temperature. This collision operator is denoted as AT+a in the nomenclature of [27] and updates the particle velocities according to

$$\begin{eqnarray}&&{\mathbf{v}}_{i}^{*}={{{\mathbf{v}}}_{C}}+{\mathbf{v}}_{i}^{{\mbox{ran}}}-\mathop{\sum }\limits_{j\in C}\frac{{\mathbf{v}}_{j}^{{\mbox{ran}}}}{{{N}_{C}}}+m{{\Pi }^{-1}}\mathop{\sum }\limits_{j\in C}\left[ {{{\mathbf{r}}}_{j,C}}\times \left( {{{\mathbf{v}}}_{j}}-{\mathbf{v}}_{j}^{{\mbox{ran}}} \right) \right]\times {{{\mathbf{r}}}_{i,C}},\end{eqnarray} \tag{ B.2 }$$

where ${{{\mathbf{v}}}_{C}}$ is the centre of mass velocity of the collision cell containing ${{N}_{C}}$ particles, $\Pi$ is the moment of inertia tensor of the particles, ${{{\mathbf{r}}}_{i,C}}={{{\mathbf{r}}}_{i}}-{{{\mathbf{r}}}_{C}}$ is the relative particle position, and ${\mathbf{v}}_{i}^{{\mbox{ran}}}$ is a random velocity drawn from a Maxwell–Boltzmann distribution. The cell grid is shifted randomly before each collision step to restore Galilean invariance of the system [28]. The dynamic viscosity η of the MPC-AT+a solvent for large number density n (particles per cell) is given by

$$\begin{eqnarray}&&\eta =\frac{n{{k}_{B}}Th}{\Delta {{x}^{d}}}\left( \frac{n}{n-\left( d+2 \right)/4}-\frac{1}{2} \right)+\frac{m\left( n-7/5 \right)}{24\Delta {{x}^{d-2}}h},\end{eqnarray} \tag{ B.3 }$$

where ${{k}_{B}}T$ is the imposed temperature and d = 2 is the dimensionality of the system.

The immersed droplets are modelled as rigid discs of radius R that are moved according to Newtonʼs equation of motion. The discs are coupled to the fluid by a no-slip boundary condition. In the streaming step, a bounce-back boundary condition is applied to the solvent particles, i.e., ${\mathbf{v}}^{\prime} =-{\mathbf{v}}+2{{{\mathbf{v}}}_{{\mbox{b}}}}$ where ${{{\mathbf{v}}}_{{\mbox{b}}}}$ is the boundary velocity. In the collision step, the cell volume occupied by the boundary is filled with virtual particles that are distributed randomly within a layer of width $\sqrt{2}a$. The velocities of the virtual particles are drawn from a Maxwell–Boltzmann distribution and the velocity ${{{\mathbf{v}}}_{{\mbox{b}}}}$ of the boundary is added [16]. The momentum change of the solvent particles during streaming and collisions is accumulated and leads to the boundary force ${{{\mathbf{F}}}_{{\mbox{b}}}}$ that moves the discs [19].

Since hydrodynamics is not resolved on short distances, we add steric droplet-droplet and droplet–wall interactions. These are modelled as Weeks–Chandler–Anderson (WCA) potentials

$$\begin{eqnarray}&&{{V}_{{\mbox{WCA}}}}\left( r \right)=4\epsilon \left[ {{\left( \frac{{{r}_{0}}}{r} \right)}^{12}}-{{\left( \frac{{{r}_{0}}}{r} \right)}^{6}} \right]+\epsilon \qquad 0<r<{{2}^{\frac{1}{6}}}{{r}_{0}},\end{eqnarray} \tag{ B.4 }$$

where $r={{r}_{ij}}-2R$ for droplet–droplet interactions, and $r={{r}_{ij}}-R$ for droplet–wall interactions. ${{r}_{ij}}$ is the distance of the droplet centres, or the distance of the droplet centre from the wall surface, respectively.

The friction at the top and at the bottom of the Hele–Shaw cell is accounted for by a friction force ${{{\mathbf{F}}}_{{\mbox{friction}}}}=-\gamma {{{\mathbf{u}}}_{{\mbox{d}}}}$ where ${{{\mathbf{u}}}_{{\mbox{d}}}}$ is the velocity of the droplet relative to the microchannel. The flow is driven by an external force ${\mathbf{g}}$ corresponding to a constant pressure gradient across the channel. The initial setup of the droplets is a triangle wave of a given wavelength λ.

The basic simulation units are chosen such that the size of an MPC collision is $\Delta x=R/5$ and the MPC time step is $h=0.005\;\tau$, where the time scale is $\tau ={{\left( {{k}_{B}}T/m \right)}^{-1/2}}\Delta x$, and m is the mass of one MPC particle. The parameters for the MPC simulation are the fluid density $\rho =40\;m/\Delta {{x}^{2}}$, the driving force $g=0.1\;m\Delta x/\tau$, and the friction $\gamma =2\times {{10}^{4}}\;m/\tau$. These parameters correspond to a fluid viscosity $\eta =321.77\;m/\tau$. The mass of the droplets is given by $M=\rho \pi {{R}^{2}}\approx 3140\;m$. The parameters for the WCA potential are $\epsilon ={{k}_{B}}T$ and ${{r}_{0}}=\Delta x$. We varied the channel width W, the droplet spacing a and the initial wavelength λ. The simulation results presented in the main article are for 90 droplets in a channel of length $L=1000\Delta x$ and width $W=15\Delta x$. The initial zigzag pattern was set up with a length of six droplets.

In order to compare the simulation results to the linearized theory and to the experiments, we determined the value $K={{u}_{{\mbox{d}}}}/{{u}_{{\mbox{oil}}}}$ from independent simulation runs with a single droplet under the same confinement. For a channel width of $W=3R$ we obtained a value of $K\approx 0.62$ which is on the order of the experimental parameters. The oil velocity is in the range ${{u}_{{\mbox{oil}}}}\approx {{10}^{-2}}\Delta x/\tau$, such that the Reynolds and Peclet number of the simulations are $\operatorname{Re}\approx {{10}^{-3}}$ and $Pe\approx {{10}^{2}}$, respectively.

The complex flow potential for a point-like dipole at position ${{z}_{j}}={{x}_{j}}+i{{y}_{j}}$ in the channel evaluated at position $z=x+iy$ is [7, 8]

$$\begin{eqnarray}\begin{array}{rcl} {{w}_{j}}\left( z \right) & = & \frac{\pi {{R}^{2}}\delta {{u}_{0}}}{2W}\left\{ {\rm coth} \left[ \frac{\pi \left( z-{{z}_{j}}-i{{y}_{j}} \right)}{2W} \right] \right. \\ {} & {} & +\left. {\rm coth} \left[ \frac{\pi \left( z-{{z}_{j}}-i\left( W-{{y}_{j}} \right) \right)}{2W} \right] \right\} \\ {} & = & \frac{\pi {{R}^{2}}\delta {{u}_{0}}}{2W}\left\{ {\rm coth} \left[ \frac{\pi \left( x-{{x}_{j}}+i\left( y-2{{y}_{j}} \right) \right.}{2W} \right] \right. \\ {} & {} & +\left. {\rm coth} \left[ \frac{\pi \left( x-{{x}_{j}}+i\left( y-W \right) \right.}{2W} \right] \right\}. \\ \end{array}\end{eqnarray} \tag{ C.1 }$$

Differentiating ${{w}_{j}}\left( z \right)$ with respect to z we obtain

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\mbox{d}}{{w}_{j}}}{{\mbox{d}}z} & = & -\frac{{{\pi }^{2}}{{R}^{2}}\delta {{u}_{0}}}{4{{W}^{2}}}\left\{ {\mbox{csc}}{{{\mbox{h}}}^{2}}\left[ \frac{\pi \left( z-{{z}_{j}}-i{{y}_{j}} \right)}{2W} \right] \right. \\ {} & {} & +\left. {\mbox{csc}}{{{\mbox{h}}}^{2}}\left[ \frac{\pi \left( z-{{z}_{j}}-i\left( W-{{y}_{j}} \right) \right)}{2W} \right] \right\} \\ {} & = & -\frac{{{\pi }^{2}}{{R}^{2}}\delta {{u}_{0}}}{4{{W}^{2}}}\left\{ {\mbox{csc}}{{{\mbox{h}}}^{2}}\left[ \frac{\pi \left( x-{{x}_{j}}+i\left( y-2{{y}_{j}} \right) \right)}{2W} \right] \right. \\ {} & {} & +\left. {\mbox{csc}}{{{\mbox{h}}}^{2}}\left[ \frac{\pi \left( x-{{x}_{j}}+i\left( y-W \right) \right)}{2W} \right] \right\}. \\ \end{array}\end{eqnarray} \tag{ C.2 }$$

Note that ${\mbox{d}}{{w}_{j}}/{\mbox{d}}z$ is related to the force exerted on other droplets, and it does not contain terms involving the velocity ${{\dot{x}}_{j}}$ or ${{\dot{y}}_{j}}$ of the droplets.

We expand around ${{x}_{n}}=na$, ${{x}_{n+j}}=\left( n+j \right)a$, ${{y}_{n}}=0$, and ${{y}_{n+j}}=0$ and keep all linear terms

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\mbox{d}}{{w}_{j}}}{{\mbox{d}}z} & \approx & -\frac{{{\pi }^{2}}{{R}^{2}}\delta {{u}_{0}}}{4{{W}^{2}}}\left\{ 4{\mbox{csc}}{{{\mbox{h}}}^{2}}\left( \frac{aj\pi }{W} \right) \right. \\ {} & {} & -\frac{4\pi i}{W}{\mbox{csc}}{{{\mbox{h}}}^{3}}\left( \frac{aj\pi }{W} \right){{\left[ 1+{\rm cosh} \left( \frac{aj\pi }{W} \right) \right]}^{2}}\delta {{y}_{n+j}} \\ {} & {} & +\frac{8\pi i}{W}{\rm coth} \left( \frac{aj\pi }{W} \right){\mbox{csc}}{{{\mbox{h}}}^{2}}\left( \frac{aj\pi }{W} \right)\delta {{y}_{n}} \\ {} & {} & +\left. \frac{8\pi }{W}{\rm coth} \left( \frac{aj\pi }{W} \right){\mbox{csc}}{{{\mbox{h}}}^{2}}\left( \frac{aj\pi }{W} \right)\left( \delta {{x}_{n}}-\delta {{x}_{n+j}} \right) \right\} \\ \end{array}\end{eqnarray} \tag{ C.3 }$$

Now we add the expansion of ${\mbox{d}}{{w}_{-j}}/{\mbox{d}}z$ around ${{x}_{n-j}}=\left( n-j \right)a$ and arrive at

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\mbox{d}}{{w}_{j}}}{{\mbox{d}}z}+\frac{{\mbox{d}}{{w}_{-j}}}{{\mbox{d}}z} & \approx & -\frac{{{\pi }^{2}}{{R}^{2}}\delta {{u}_{0}}}{4{{W}^{2}}}\left\{ 8{\mbox{csc}}{{{\mbox{h}}}^{2}}\left( \frac{aj\pi }{W} \right) \right. \\ {} & {} & -\frac{4\pi i}{W}{\mbox{csc}}{{{\mbox{h}}}^{3}}\left( \frac{aj\pi }{W} \right){{\left[ 1+{\rm cosh} \left( \frac{aj\pi }{W} \right) \right]}^{2}}\left( \delta {{y}_{n+j}}-\delta {{y}_{n-j}} \right) \\ {} & {} & -\left. \frac{8\pi }{W}{\rm coth} \left( \frac{aj\pi }{W} \right){\mbox{csc}}{{{\mbox{h}}}^{2}}\left( \frac{aj\pi }{W} \right)\left( \delta {{x}_{n+j}}-\delta {{x}_{n-j}} \right) \right\}. \\ \end{array}\end{eqnarray} \tag{ C.4 }$$

Since the expansion of the compressibility factor C around ${{y}_{j}}=0$ has a vanishing linear term, we can readily multiply with $C=\frac{W}{\pi R}{\rm tan} \left( \frac{\pi R}{W} \right)$ and K, and sum over j to obtain the equations of motion for the nth droplet

$$\begin{eqnarray}\begin{array}{rcl} \mathop{\delta {{x}_{n}}}\limits^{} & = & 2\;K\delta {{u}_{0}}\frac{{{\pi }^{2}}R}{{{W}^{2}}}{\rm tan} \left( \frac{\pi R}{W} \right) \\ {} & {} & \cdot \mathop{\sum }\limits_{j=1}^{\infty }{\rm coth} \left( \frac{aj\pi }{W} \right){\mbox{csc}}{{{\mbox{h}}}^{2}}\left( \frac{aj\pi }{W} \right)\left( \delta {{x}_{n+j}}-\delta {{x}_{n-j}} \right) \\ \end{array}\end{eqnarray} \tag{ C.5 }$$

$$\begin{eqnarray}\begin{array}{rcl} \mathop{\delta {{y}_{n}}}\limits^{} & = & -K\delta {{u}_{0}}\frac{{{\pi }^{2}}R}{{{W}^{2}}}{\rm tan} \left( \frac{\pi R}{W} \right) \\ {} & {} & \cdot \mathop{\sum }\limits_{j=1}^{\infty }{\mbox{csc}}{{{\mbox{h}}}^{3}}\left( \frac{aj\pi }{W} \right){{\left[ 1+{\rm cosh} \left( \frac{aj\pi }{W} \right) \right]}^{2}}\left( \delta {{y}_{n+j}}-\delta {{y}_{n-j}} \right). \\ \end{array}\end{eqnarray} \tag{ C.6 }$$

Plugging in a plane wave solution we arrive at the dispersion relations

$$\begin{eqnarray}&&{{\omega }_{\parallel }}\left( k \right)=-4B\mathop{\sum }\limits_{j=1}^{\infty }{\rm sin} \left( jka \right){\rm coth} \left( j\pi \beta \right){\mbox{csc}}{{{\mbox{h}}}^{2}}\left( j\pi \beta \right)\end{eqnarray} \tag{ C.7 }$$

$$\begin{eqnarray}&&{{\omega }_{\bot }}\left( k \right)=2B\mathop{\sum }\limits_{j=1}^{\infty }{\rm sin} \left( jka \right){{\left[ 1+{\rm cosh} \left( j\pi \beta \right) \right]}^{2}}{\mbox{csc}}{{{\mbox{h}}}^{3}}\left( j\pi \beta \right),\end{eqnarray} \tag{ C.8 }$$

where we introduced the prefactors $B={{c}_{0}}\left( {{\pi }^{2}}R/{{W}^{2}} \right){\rm tan} \left( \pi R/W \right)$, ${{c}_{0}}=K\left( {{u}_{{\mbox{oil}}}}-{{u}_{{\mbox{d}}}} \right)$, $\beta =a/W$, and $K={{u}_{{\mbox{d}}}}/{{u}_{{\mbox{oil}}}}$ as in [7, 8].