A minimal model for a hydrodynamic fingering instability in microroller suspensions

TL;DR

This paper develops a minimal continuum model to explain the hydrodynamic fingering instability in microroller suspensions, accurately capturing experimental and simulation results through a simple two-line rotlet system.

Contribution

The paper introduces a minimal two-line rotlet continuum model that reproduces the fingering instability and its lengthscale selection driven solely by hydrodynamics.

Findings

Model reproduces observed lengthscale selection

Dispersion relation matches simulations and experiments

Instability driven by combined advective and transverse flows

Abstract

We derive a minimal continuum model to investigate the hydrodynamic mechanism behind the fingering instability recently discovered in a suspension of microrollers near a floor [Driscoll et al. Nature Physics, 2016]. Our model, consisting of two continuous lines of rotlets, exhibits a linear instability driven only by hydrodynamics interactions, and reproduces the lengthscale selection observed in large scale particle simulations and in experiments. By adjusting only one parameter, the distance between the two lines, our dispersion relation exhibits quantitative agreement with the simulations and qualitative agreement with experimental measurements. Our linear stability analysis indicate that this instability is caused by the combination of the advective and transverse flows generated by the microrollers near a no-slip surface. Our simple model offers an interesting formalism to characterize other hydrodynamic instabilities that have not been yet well understood, such as size scale selection in suspensions of particles sedimenting adjacent to a wall, or the recently observed formations of traveling phonons in systems of confined driven particles.

Equations

$$v_x(x,y;h)$$

$$v_y(x,y;h)$$

$$V_x^{(1,I)}\left(x_1(y),y\right)$$

$$V_y^{(1,I)}\left(x_1(y),y\right)$$

$$V_x^{(1,S)}\left(x_1(y),y\right)$$

$$V_y^{(1,S)}\left(x_1(y),y\right)$$

$$\frac{\partial x_1(y,t)}{\partial t}=V_x^{(1,S)}+V_x^{(1,I)},$$

$$\frac{\partial x_2(y,t)}{\partial t}=V_x^{(2,S)}+V_x^{(2,I)},$$

$$\frac{\partial {\rho }_1(y,t)}{\partial t}=-\frac{\partial \left[{\rho }_1(y,t)\left(V_y^{(1,S)}+V_y^{(1,I)}\right)\right]}{\partial y},$$

$$\frac{\partial {\rho }_2(y,t)}{\partial t}=-\frac{\partial \left[{\rho }_2(y,t)\left(V_y^{(2,S)}+V_y^{(2,I)}\right)\right]}{\partial y}.$$

$$x_1(y,t)=0+\delta x_1(y,t)$$

$$x_2(y,t)=d+\delta x_2(y,t),$$

$${\rho }_1(y,t)={\rho }_0+\delta {\rho }_1(y,t)$$

$${\rho }_2(y,t)={\rho }_0+\delta {\rho }_2(y,t),$$

$$\frac{\partial \delta x_1}{\partial t}$$

$$\frac{\partial \delta x_2}{\partial t}$$

$$\frac{\partial \delta {\rho }_1}{\partial t}$$

$$\frac{\partial \delta {\rho }_2}{\partial t}$$

$${\mathcal{G}}_0$$

$${\mathcal{G}}_1(y-y^{\prime })$$

$${\mathcal{G}}_2(y-y^{\prime })$$

$${\mathcal{G}}_3(y-y^{\prime })$$

$$\mathbf{u}=(\delta x_1,\delta x_2,\delta {\rho }_1,\delta {\rho }_2)=\tilde{\mathbf{u}}{e}^{iky+\sigma t},$$

$$\tilde{\mathbf{u}}=(\delta {\tilde{x}}_1,\delta {\tilde{x}}_2,\delta {\tilde{\rho }}_1,\delta {\tilde{\rho }}_2),$$

$$A\tilde{\mathbf{u}}=\sigma \tilde{\mathbf{u}}$$

$$A=\left[\begin{array}[]{cccc}\rho_{0}\mathcal{G}_{0}&\rho_{0}\tilde{\mathcal{G}}_{1}&0&\tilde{\mathcal{G}}_{x}(-d,k)\\ -\rho_{0}\tilde{\mathcal{G}}_{1}&-\rho_{0}\mathcal{G}_{0}&\tilde{\mathcal{G}}_{x}(-d,k)&0\\ -ik\rho_{0}^{2}\tilde{\mathcal{G}}_{2}&-ik\rho_{0}^{2}\tilde{\mathcal{G}}_{3}&0&-ik\rho_{0}\tilde{\mathcal{G}}_{y}(-d,k)\\ -ik\rho_{0}^{2}\tilde{\mathcal{G}}_{3}&-ik\rho_{0}^{2}\tilde{\mathcal{G}}_{2}&ik\rho_{0}\tilde{\mathcal{G}}_{y}(-d,k)&0\end{array}\right].$$

$${\left(A_{14}{A}_{31}\right)}^{1/2}$$

$$=$$

$$v_x\left(x,y;h\right)$$

$$v_y\left(x,y;h\right)$$

$$A_{11}$$

$$={\rho }_0{\mathcal{G}}_0$$

$$={\rho }_{0}Sh\frac{8d(d^2-4h^2)}{3[d^2+4h^2{]}^3}$$

$$A_{12}$$

$$={\rho }_0{\tilde{\mathcal{G}}}_1(k)$$

$$={\rho }_0\underset{-\infty }{\overset{\infty }{\int }}{\mathcal{G}}_1(y)\exp (-iky)dy$$

$$={\rho }_{0}Sh\frac{2d}{15}[\frac{-3d^2+8h^2}{(d^2+4h^2{)}^{3/2}}{k}^3{K}_3\left(\sqrt{d^2+4h^2}k\right)$$

$$+\frac{8}{(d^2+4h^2{)}^2}MG\left(\left\{(-1/2),()\right\},\left\{(0,2),(1/2)\right\},\frac{1}{4}(d^2+4h^2)k^2\right)]$$

$$A_{14}$$

$$={\tilde{\mathcal{G}}}_x(-d,k)$$

$$=Sh\frac{2d^2}{3(d^2+4h^2)}{k}^2{K}_2\left(\sqrt{d^2+4h^2}k\right)$$

Full text

instability, pattern formation, active suspensions

I Introduction

When suspended in a viscous fluid at low Reynolds number, small moving particles interact through long range hydrodynamic interactions. These many-body interactions can give rise to strong density and velocities fluctuations in the bulk and lead to instabilities Guazzelli2011 . The addition of boundaries strongly modify the hydrodynamic interactions between particles and affect the dynamics of the system. For instance, sedimenting particles between two parallel plates exhibit a transverse Rayleigh-Taylor-like instability whose wavelength strongly depends on the distance between the plates Carpen2002 .

When these particles are driven by means of an external field, or by self-propulsion mechanisms, they induce active flows that modify the interactions within the suspension and sometimes lead to strong density fluctuations and long ranged orientational correlations Saintillan2013 ; Marchetti2013 . A well-known example is the instability of the isotropic state in suspensions of elongated swimmers, which was predicted by the theory AditiSimha2002 ; Saintillan2008a ; Baskaran2009 ; Ezhilan2013 and reported by both numerics Saintillan2007 ; Saintillan2011 ; Lushi2013 and experiments Wu2000 ; Sokolov2007 ; Wensink2012 ; ryan_2013 ; Creppy2015 . Under confinement, active and driven suspensions exhibit a wide variety of behaviours such as the formation of vortices, asters or polar bands Lefauve2014 ; Brotto2013 ; Tsang2015 ; Tsang2016 .

In a recent work, we have uncovered a new hydrodynamic instability in a driven system of microrollers: suspensions of colloids rotating parallel to a floor critters . These microrollers consists of spherical polymer colloids with radius $a=0.66$ $\mu$m which have a small permanent magnetic moment ($\lvert\mathbf{m}\rvert\sim 5\cdot 10^{-16}$ Am2) due to an embedded hematite cube Sacanna2012 , see schematic in Fig. 1a. Their equilibrium gravitational height is given by the competition between gravity and thermal fluctuations: $h_{g}=a+k_{B}T/mg$, where $a$ is the particle radius, $m$ its buoyant mass, $g$ the gravitational acceleration and $k_{B}T$ the thermal agitation energy. More details about the experimental system are provided in critters . When driven by a magnetic field rotating about an axis parallel to the floor, these microrollers generate strong advective flows in their vicinity. These flows are responsible for the large-scale collective effects observed in uniform suspensions. When the particle distribution is discontinuous, the microrollers form a shock front with a well defined width (Fig. 1b), which quickly becomes unstable and generates finger-like structures with compact tips (Fig. 1c). These fingertips can detach to form compact autonomous structures, called“critters” critters .

Our large scale particle simulations and experiments suggest that the fingering instability is a linear instability controlled by the height of the particles above the floor (Fig. 1c,d,e). With the addition of Brownian motion, our simulations (Fig. 1d) achieve quantitative agreement with the experiments on the measured wavelength for various gravitational heights $h_{g}$ Usabiaga2017 . The fingering instability was also reproduced with a much simpler 2D system (Fig. 1e): rotlets, i.e. point torques, interacting only through hydrodynamic interactions in a plane at a fixed height above the floor critters .

These simulations showed that the instability is controlled by far-field hydrodynamic interactions in the plane parallel to the floor, but they do not provide a clear physical picture of the instability mechanism. Nonlinear phenomena such as the shock formation do not elucidate the exact mechanisms at play behind the transverse instability of the front that forms the fingers.

We recently developed a theoretical model to study density fluctuations in uniform suspensions of microrollers and to investigate the formation of the shock front Delmotte2017 . This model relies on a nonlocal description of the far-field hydrodynamic interactions between microrollers treated as rotlets. The theoretical results and comparisons with experiments and simulations showed that the front forms due to the nonlocal nature of far-field hydrodynamic interactions in the plane above the floor, and that it has a well defined finite width which is controlled by the particle height. However, the model used in this prior work to describe the front Delmotte2017 is one-dimensional and does not allow for transverse variations, and therefore cannot be used to study the two-dimensional transverse fingering instability.

In this paper, we derive a minimal two-dimensional continuum model, based on a nonlocal description, in order to: (1) confirm that the transverse instability is a linear instability, (2) study the dependence of the characteristic wavelength on the control parameters, (3) identify precisely the hydrodynamic mechanisms at play. We compare the linear stability analysis with numerical simulations and experiments and discuss the validity of the model. We take advantage of the flexibility and simplicity of the model to study the stability of the microrollers above a slip surface and finally offer promising extensions to other particulate systems within the same framework.

II Model and linear stability analysis

Our goal is to model the shock front in a simplified manner in order to carry out an analytic linear stability analysis of the system. Using a finite sheet would make the calculations untractable because of the nonlocality of the hydrodynamic interactions. Suspensions of particles have been studied with a two-fluid approach Carpen2002 ; Pan2001 ; Wysocki2009 , which requires computing the effective viscosity of the suspension. However, a two-fluid model is not applicable to the confined microroller suspension where hydrodynamic interactions are nonlocal and depend on the particle height $h$. We found that the simplest and most relevant approach to analytically model an unstable front of microrollers is to consider two lines of rotlets interacting hydrodynamically in a plane parallel to the floor at a fixed height $h$, as one line of rotlets is linearly stable to perturbations.

II.1 Governing equations

Consider two infinite lines of rotlets, labelled “1”, at the back, and “2”, at the front, respectively, with rotlet densities $\rho_{1}(y)$ and $\rho_{2}(y)$, rotating around the $y$ axis in a plane at a fixed height $h$ above a no-slip boundary. Their $x$ position is parametrized by $y$, with an initial state $x_{1}(y,t=0)=0$ and $x_{2}(y,t=0)=d$ (see Fig. 2a). Here, we neglect out of plane motion in the $z$-direction and only consider the velocities in the $xy$-plane. Although out of plane motion is seen in the experiments, our quasi-2D simulations in Fig. 1c confirm that the fingering instability still occurs when considering particles confined to a plane (with 3D hydrodynamics).

A point rotlet located at a position $(x^{\prime},y^{\prime};\,h)$ induces a horizontal fluid velocity at point $(x,y;\,h)$ given by Blake1974 ; Delmotte2017 :

[TABLE]

[TABLE]

where $S=6T_{y}/(8\pi\eta)$, $T_{y}>0$ is the applied constant torque around the $y$-axis (i.e. the rotlet magnitude). Changing the sign of $T_{y}$ would change the direction of motion of the system and reverse the role of the lines. $\mathcal{G}_{x}$ and $\mathcal{G}_{y}$ are the hydrodynamic interaction kernels in the $xy$-plane. Note that these kernels are not singular because of the finite constant height in the denominator $h>0$. The velocity field around a rotating particle above a no-slip wall is shown in Fig. 2b. In an unbounded fluid, this velocity field would be zero in the whole plane. The transverse flow $v_{y}$ observed in Fig. 2b arises from the image system of the rotlet which includes a stresslet Blake1974 . Since $h$ is fixed, we will henceforth omit it in the argument of the kernels: $\mathcal{G}(x,y)\equiv\mathcal{G}(x,y;h)$.

The velocity induced by the front line “2” at a point $(x_{1}(y),y)$ on the back line “1” is given by the functionals

[TABLE]

where $\rho_{2}(y^{\prime})=\rho(x_{2}(y^{\prime}),y^{\prime})$. In $V^{(1,I)}_{x,y}$ the first superscript “1” indicates the line considered and the second superscript “I” stands for “induced” by the other line. The self-induced velocity, “S”, is given by

[TABLE]

Similar terms are derived for line 2. The governing equations for the two lines, in the limit of small deviations from a straight line, are the equations of motion in the direction perpendicular to the lines

[TABLE]

and mass conservation coming from the motion parallel to the lines,

[TABLE]

These four nonlinear, nonlocal, coupled equations (7)-(10) can be linearized about the homogeneous state.

II.2 Linear stability analysis

We perturb the $x$ positions of each line about their initial position

[TABLE]

where $\delta x_{1},\delta x_{2}\ll\min(d,h)$, and their rotlet densities about a constant value $\rho_{0}$,

[TABLE]

where $\delta\rho_{1},\delta\rho_{2}\ll\rho_{0}$. After Taylor expanding the functionals in Eqs. (7)-(10) we obtain the linearized governing equations

[TABLE]

where we have used the symmetries of $\mathcal{G}_{x}$ and $\mathcal{G}_{y}$. The star “$*$” denotes the one dimensional convolution product. The interaction kernels are given by

[TABLE]

Note that $\mathcal{G}_{0}$ is a constant which depends only on the geometric parameters of the problem $d$, $h$, and on $S$. Its physical meaning will be explained below.

Looking for periodic solutions of Eqs. (15)-(18) of the form

[TABLE]

where the wavenumber $k=2\pi/\lambda$ and

[TABLE]

we obtain the following eigenvalue problem

[TABLE]

where

[TABLE]

The complete expression for the entries of $A$ can be found in Appendix A. Note that a model with one line of rotlets would be linearly stable since all self-induced terms are at least quadratic in the magnitude of the perturbations.

II.2.1 Structure of the instability

Due to the particular structure of $A$, Eq. (19) can be solved analytically. The four solutions $\sigma_{1}(k),..,\sigma_{4}(k)$ are written in Appendix A. The first two eigenvalues are real and of opposite sign $\sigma_{1}=-\sigma_{2}$. The two other are imaginary and conjugate $\sigma_{3}=\bar{\sigma_{4}}$. All eigenvalues depend linearly on $S\propto T_{y}$. Figure 3a shows the real eigenvalues $\sigma_{1,2}$ for $d=10h$. First, one can see that the eigenvalues exhibit a well defined peak at $\lambda=\lambda_{m}$, meaning that the two-line model selects a fastest growing mode, which is characteristic of a fingering instability. Second, all modes, except the zero mode, are unstable and a clear plateau is visible for short wavelengths. The value of this plateau is exactly $\rho_{0}\mathcal{G}_{0}$, which is the only constant entry in the matrix $A$. To better understand the presence of the plateau, we study the structure of the amplified eigenmode $\mathbf{u}_{1}$, associated to the real positive eigenvalue $\sigma_{1}$, in Fig. 3c. For this eigenmode, the density and position disturbances of the front line, $\delta\tilde{\rho}_{2}$ and $\delta\tilde{x}_{2}$, vanish for $\lambda\approx\lambda_{c}$, where $\lambda_{c}$ is a critical value, while the back line disturbances, $\delta\tilde{\rho}_{1}$ and $\delta\tilde{x}_{1}$, are nonzero for all $\lambda$’s.

This difference in behavior between the two lines can be understood by looking at the flow induced by one line on the other. Fig. 3d shows that a straight line of rotlets at $x=0$ (the back line), with uniform density $\rho_{0}$, damps the disturbances of the other line at $x=d$, regardless of their wavelength $\lambda$. Conversely, a straight line at $x=d$ (the front line) increases the disturbances of the line at $x=0$ for all $\lambda$. This phenomenon can be explained with a simple analogy: consider a force that dies off quickly with distance from a line. Take a parallel line a distance away and introduce a distortion so that parts are closer and others further away from the driving force. If the forces are attractive then the closer parts are drawn faster than the further parts and the distortion is amplified. If the forces are repulsive the closer parts are pushed away faster than the far ones and the distortion is reduced. For two lines of rollers the leading line advects the trailing one toward it while the tailing line pushes the leader away (Fig. 3d). Thus, the front line is stable while the rear goes unstable. The term which relates a straight line with density $\rho_{0}$ and a line with perturbed position is precisely $\rho_{0}\mathcal{G}_{0}$ (see Eqs. (15)-(16)). Thus the presence of the plateau is due to the growth of each mode on the back line. The term $\rho_{0}\mathcal{G}_{0}$ is compensated by $\rho_{0}\tilde{\mathcal{G}}_{1}(k)$ when $k\rightarrow 0$. Indeed $\lim_{k\rightarrow 0}\rho_{0}\tilde{\mathcal{G}}_{1}(k)=-\rho_{0}Sh\frac{8d(d^{2}-4h^{2})}{3(d^{2}+4h^{2})^{3}}=-\rho_{0}\mathcal{G}_{0}$.

II.2.2 Comparison with nonlinear simulations

We simulate the two-line system with discrete rotlet simulations. Each line is discretized with 2000 rotlets interacting in the $xy$ plane. The two initially straight lines are perturbed by adding a small random increment to the particle positions. The hydroynamic interactions are calculated using Eqs. (1)-(2), and pseudo-periodic boundary conditions are used Usabiaga2017 . Positions are updated with an Adams-Bashforth-Moulton predictor-corrector scheme (see Methods section in critters ). The full development of the fingering instability in the simulations is shown in Fig. 4a. To characterize the instability, we measure the position disturbance of each line $\delta x_{1,2}(y,t)$ that we bin along the $y$ direction with 1024 points. The resulting two vectors are Fourier transformed to obtain the power spectra for each line shown in Fig. 4b at two different times. At $t=8.65$ s, both lines have a uniformly distributed spectrum. At $t=21.15$ s, both spectra exhibit a clear peak at long wavelengths. The front line (light blue) has damped all the modes with short wavelength while the back line (dark blue) has amplified all modes. The time evolution of three Fourier modes $\lambda=20,91,146$ $\mu$m is shown in Fig. 4c. One can see that the shortest wavelength, $\lambda=20$ $\mu$m, decreases for the front line and grows for the back line. We use an exponential fit $\delta\tilde{x}(t)\sim\exp(\sigma t)$ (see inset) to extract the growth rate $\sigma$ of each mode for each line at the onset of the instability.

Fig. 3a-b compares the extracted growth rate, $\sigma(k)$, of each line, averaged over 10 simulations, to the theoretical predictions, with no adjustable parameters. The theory is in excellent agreement with the particle simulations, demonstrating that the linear stability analysis correctly captures the early-time dynamics of the fingering instability. Consistently with the theoretical predictions in Fig. 3c, all modes of the back line are unstable ($\sigma>0$), while short wavelength disturbances on the front line are damped at a rate $-\rho_{0}\mathcal{G}_{0}$ when $\lambda$ is below the critical value $\lambda_{c}$.

II.2.3 Dependence on the control parameters

The fastest growing mode of the two line model $\lambda_{m}$ is controlled by two geometric parameters: the height $h$ and the distance between the two lines $d$. Since $\sigma_{1}=Sf(h,d,k=2\pi/\lambda)$ has a complex analytic form (see Eq. (LABEL:eq:EVS)), it is not possible to express $\lambda_{m}$ as a function of $h$ and $d$ explicitly. Instead we use numerical evaluations to determine $\lambda_{m}$ as a function of $h$ and $d$. Figure 5 shows the graph of $\lambda_{m}/h$ vs. $d/h$ for various heights $h=1-10$ $\mu$m. The collapse of the curves show that $\lambda_{m}$ depends linearly on the height $h$: $\lambda_{m}=hf(d/h)$, where $f(d/h)$ is a function of $d/h$. As shown in this Figure, $f(d/h)$, becomes approximately linear when $d/h>5$: $f(d/h)\approx 2.02(d/h)+2.85$.

III Comparison with large scale simulations and experiments

In this section we compare the two-line continuum model with numerical simulations and experimental measurements. We want to evaluate how relevant our model is to the more realistic microroller system.

The quasi-2D large scale simulations are performed with the method described in critters . From the particle positions at a given time $t$ we can compute the empirical number density $n(x,y,t)$ in the $xy$-plane. We compute the Fourier transform of $n(x,y,t)$ at $k_{x}=0$, i.e., the Fourier transform $\tilde{n}(k_{y}=2\pi/\lambda,t)$ of the number density along the direction of the front. We only use particles in the shock front, specifically, we only include the 70$\%$ of the particles with the largest $x$-coordinates. This ensures that the Fourier modes are not affected by the particles left behind the shock front. We have confirmed that essentially the same results are obtained when including between 50$\%$ and 90$\%$ of the particles.

Figure 6a compares the growth rate of the two line model with our quasi-2D particle simulations, where the microrollers are restricted to the $xy$-plane critters , at several heights $h=1.97-4.92$ $\mu$m. Two parameters must be adjusted in the two-line model to match these simulations: the distance between the two lines $d$ to match $\lambda_{m}$, and the strength $S$ (or $\rho_{0}$) to match the magnitude of the maximum growth rate $\sigma_{m}$. Setting $d/h=9.5$ matches $\lambda_{m}$ quantitatively for all the simulated heights, which shows that the instability wavelength depends linearly on $h$. It also suggest that, as shown by our previous work on the shock front Delmotte2017 , the width of the front is proportional to the height above the wall. The experimental stability diagrams in Fig. 6b are qualitatively similar: they have a well-defined fastest growing mode $\lambda_{m}$ which increases with the particle height.

Thus the two-line model contains the essential physical ingredients to capture the fingering instability observed in the simulations and in the experiments. Its simplicity allows us to study the physical meaning of each term in the governing equations in order to understand the mechanisms at play.

IV Mode selection by transverse flows

In this section we examine the role of each term from the matrix $A$ (Eq. (20)) that appear in the full expression of the real positive eigenvalue $\sigma_{1}(k)$ in Eq. (LABEL:eq:EVS). We find that the location and the shape of the peak around $\lambda=\lambda_{m}$ in $\sigma_{1}$ are mainly controlled by the following term in Eq. (LABEL:eq:EVS)

[TABLE]

where $K_{n}(x)$ is the modified Bessel function of the second kind. Without this term there is no peak and therefore no lengthscale selection. The product $\tilde{\mathcal{G}}_{2}(k)\tilde{\mathcal{G}}_{x}(-d,k)$ indicates that both terms $\mathcal{G}_{x}$ and $\mathcal{G}_{2}$ must be nonzero to generate the fingering instability. $\mathcal{G}_{x}$ corresponds to the displacement of one line $\delta x_{i}$ induced by the density perturbations of the other line $\delta\rho_{j}$. $\mathcal{G}_{2}$ corresponds to the transverse velocity in the $y$ direction that drives the density perturbations $\delta\rho_{i}$ on each line due to the displacement on that line $\delta x_{i}$. A sketch of these mechanisms is shown in Figure 7. The physical meaning of this product is that both advection by the other line and self-induced transverse motion must be present to generate the fingering instability. We therefore conclude that the transverse instability is due to both the advective and transverse flows induced by a particle rotating parallel to a no-slip boundary (see Fig. 2b). Without both these advective and transverse flows, we would not observe this type of linear fingering instability.

To demonstrate this we now consider a flat stress-free surface, i.e. a free-slip boundary such as an air-water interface with small curvature. The image system of a rotlet near a flat slip surface is a counter rotating rotlet with the same magnitude DiLeonardo2011 , which ensures zero tangential stress at the interface. The corresponding velocity kernels replacing (1)-(2) are

[TABLE]

The absence of transverse flow in (24) leads to the absence of peak in the growth rate, i.e. no fastest growing mode, and thus no clear lengthscale selection for the two-line model, see Appendix B for a detailed analysis. Therefore, microrollers near a slip surface should not be subject to a linear fingering instability with a well-defined finger width.

V Conclusions and discussion

We derived the simplest possible model that can capture the fingering instability observed both experimentally and numerically in suspensions of microrollers above a floor. Our model directly accounts for the nonlocal hydrodynamic interactions between the particles. Our analytic linear stability analysis confirmed that the fingering instability is linear, that it happens in the plane parallel to the floor, and is purely hydrodynamic in origin 111We also carried 3D simulations of rotating particles with no external force above a no-slip surface and still observed the fingering instability, thus confirming the hydrodynamic origin of the instability.. Our comparisons with quasi-2D particle simulations showed quantitative agreement for a range of heights by adjusting only the distance between the two lines as a function of the height to $d=9.5h$. These results showed that the instability wavelength and the front width depend linearly on $h$. Thanks to the simplicity of the model, we identified each term separately and showed that the transverse flows due to the nearby no-slip surface are responsible for the lengthscale selection. In particular, for a free-slip bottom wall, there are no transverse flows and there is no linear instability. The two-line model and the linear stability analysis have improved our understanding of the fingering instability and shed more light on the genesis of the autonomous motile colloidal structures called “critters”.

Our two-line model considers only far-field hydrodynamic interactions and does not account for finite particle size effects. One can test the effect of near-field hydrodynamics on the qualitative behaviour of the system with numerical simulations. In the case of the microroller fingering instability, we found that far-field hydrodynamics alone achieved qualitative agreement with the experiments.

More importantly, our two-line model is not limited to rotating particles and finds interesting applications in other particulate systems. It could be used to study the sedimentation of particles parallel to a wall or to a slip surface, but also to investigate the formation of phonons in microfluidic systems of confined driven particles (e.g. droplets Beatus2006 or microswimmers Tsang2016b ). Indeed, both kind of particles generate transverse flows that destabilize the systems and could be analyzed within the same framework. One would only have to change the velocity kernels (and regularize them when necessary): stokeslets and their image system for sedimenting particles Blake1974 , source dipoles for confined driven particles Beatus2012 ; Brotto2013 . Surprisingly, the sedimentation of particles adjacent to a single wall has received little attention in the literature. Our preliminary large scale particle simulations suggest that the parallel wall selects a characteristic wavelength in the unstable sedimenting suspensions. We will study this phenomenon experimentally, numerically and theoretically in more detail in a near future.

VI Acknowledgments

This work was supported primarily by the Materials Research Science and Engineering Center (MRSEC) program of the National Science Foundation under Award Number DMR-1420073 and the Gordon and Betty Moore Foundation through Grant GBMF3849. A. Donev and B. Delmotte were supported in part by the National Science Foundation under award DMS-1418706. P. Chaikin was partially supported by the Center for Bio-Inspired Energy Science, A DOE BES EFRC under Award DE-SC0000989. We gratefully acknowledge the support of the NVIDIA Corporation with the donation of GPU hardware for performing some of the simulations reported here.

Appendix A Entries and eigenvalues of the matrix A

The matrix $A$ contains only six distinct coefficients listed below:

[TABLE]

[TABLE]

where $K_{n}(x)$ is the modified Bessel function of the second kind, and $MG$ is the Meijer G-function.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The structure of the matrix $A$ permits an analytic calculation of the eigenvalues:

[TABLE]

Appendix B Microrollers above a free-slip surface

Figure 8a shows the flow field around a particle above a slip surface. Note that the streamlines in the plane of rotation are parallel, and that the induced translational velocity is in the opposite direction to the no-slip case (see Fig. 2b). After carrying the same linear stability analysis we obtain the following matrix:

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The nonzero eigenvalues of the new matrix $A_{\mbox{slip}}$ are

[TABLE]

and

[TABLE]

Figure 8b compares the stability diagram between a slip and a no-slip surface for $h=1.97-4.92$ $\mu$m. The peak around $\lambda_{m}$ disappears in the absence of transverse flows, meaning that there is no clear lengthscale selection above a slip surface. The constant term $\rho_{0}\mathcal{H}_{0}$ in the matrix $A_{\mbox{slip}}$ plays the same role as $\rho_{0}\mathcal{G}_{0}$, as shown by the plateaus at short wavelengths.