Swarming behavior of simple model squirmers

We have studied experimentally the collective behavior of self-propelling liquid droplets, which closely mimic the locomotion of some protozoal organisms, the so-called squirmers. For the sake of simplicity, we concentrate on quasi-two-dimensional (2D) settings, although our swimmers provide a fully 3D propulsion scheme. At an areal density of 0.46, we find strong polar correlation of the locomotion velocities of neighboring droplets, which decays over less than one droplet diameter. When the areal density is increased to 0.78, distinct peaks show up in the angular correlation function, which point to the formation of ordered rafts. This shows that pronounced textures, beyond what has been seen in simulations so far, may show up in crowds of simple model squirmers, despite the simplicity of their (purely physical) mutual interaction.


Introduction
Large-scale patterns emerging in a crowd of interacting self-driven elements are known for a wide range of biological systems, such as a host of sparrows, a school of fish, an army of ants or bacterial colonies. The complexity of the active elements thereby varies considerably, and so do their mutual interactions. The latter can be sorted in purely physical effects, such as hard core repulsion or hydrodynamic interaction, and biological signaling, such as olfactoric and visual signals, or chemotaxis, which is present in many microbial settings. The striking similarity in the swarming behavior of a wide range of active elements has significantly increased the interest in understanding the basic mechanisms of pattern and texture formation in such systems.
While theoretical works have concentrated on active elements which were greatly simplified in their shape and interactions [1,[2][3][4][5][6][7][8], most experiments to date have been performed with bacterial colonies [9][10][11][12][13]. This is due to the great importance of understanding the collective behavior of micro-organisms and also because bacterial colonies are the only relatively simple active elements that can be obtained in sufficient numbers. Even a bacterium, however, is much more complex than model elements used so far in theories. Consequently, it is not yet clear which aspects of their collective behavior are due to physical interactions and which aspects can be traced back to more complex biological signaling.
It would therefore be desirable to perform experiments with artificial active elements, the properties and interactions of which can be well controlled and adequately described by simple physical models. Furthermore, such model systems have the potential for minimizing the individual variations pertinent to biological systems. In the search for such elements, one should try to mimic actual biological organisms in order to maximize the relevance and comparability of the system. In this context, organisms such as cyanobacteria, paramecium and volvox are particularly inspiring. They belong to a class of swimmers referred to as squirmers, and are driven by tangential and/or radial deformations of the cell surface [3,14,15]. Squirming motion is particularly appealing for the study of the hydrodynamics of micro-scale swimming, since the velocities in the near and far field around such a swimming organism can be described well analytically [4,16,17], and are similar to the flow fields around moving spherical objects. Such squirming organisms may thus be modeled by self-propelling liquid droplets, which is what we pursue in the present paper.

Squirming droplets
We present a simple model squirmer consisting of an aqueous droplet moving in an oil 'background' phase. Propulsion arises due to the spontaneous bromination of mono-olein (racglycerol-1-mono-oleate) as the surfactant. The latter is abundant in the oil phase, such that the droplet interface is covered by a dense surfactant monolayer. The bromine 'fuel' is supplied from inside the droplet, such that bromination proceeds mainly at the droplet surface. It results in saturation of the C=C double bond in the alkyl chain of the surfactant, thereby rendering it a weaker surfactant. As we will see below, this results in a self-sustained bromination gradient along the drop surface, which propels the droplet due to Marangoni stresses.
We demonstrate this squirmer scheme with nanoliter droplets containing 25 mM l −1 bromine water in a continuous oil phase of squalane containing 50 mM l −1 mono-olein (MO). The critical micelle concentration is 1.5 mM l −1 . The droplets are confined by two hydrophobic glass plates to a quasi-two-dimensional (2D) space, thus simplifying droplet tracking. Although we thereby confine the droplets to a space of reduced dimension, it should be kept in mind that their propulsion mechanism is inherently 3D, thereby closely mimicking the behavior of protozoal organisms. This clearly distinguishes our system from 'striders' which are inherently bound to an interface [18][19][20][21][22]. The top panel of figure 1 shows the trajectory of a single squirmer droplet over a duration of 400 s. The velocity of the droplet is about 15 µm s −1 . The trajectory is reminiscent of a random walk, with a persistence length that is larger than the droplet size and clearly far beyond what would be expected for Brownian motion. A particularly important observation is that the trajectory crosses itself. This is in sharp contrast to other schemes, where the propulsion mechanism itself changes the surrounding medium strongly enough to prevent self-crossing of the path [18]. That nothing like this happens here is demonstrated even more convincingly in the lower panel, which shows a time lapse representation of seven droplets moving in a micro-channel. As two drops touch each other, they reverse their direction of motion and perambulate the channel again, without significant reduction in velocity. This can be attributed to the abundance of surfactant in the oil phase, which ensures that the medium is only minutely changed by the 'exhaust' of the squirmers. This property makes our system particularly well suited for the study of collective behavior.

Mechanism of locomotion
In order to gain some insights into the propulsion mechanism, let us consider a spherical droplet with radius R. The total coverage, c, of the droplet surface with the MO, either brominated or not, is assumed to be roughly constant and in equilibrium with the micellar phase in the oil. The brominated fractional coverage shall be called b. If the droplet moves, there is (in the rest frame of the droplet) an axisymmetric flow field, u(θ), along its surface. The equation of motion where b 0 is the equilibrium coverage with brominated mono-olein (brMO). It is determined by the bromine supply from inside the droplet and the rate constant, k, of the escape of brMO into the oil phase. The first term on the rhs of equation (1) describes the balance between bromination and brMO escape. The second term describes the change in bromination density of the surfactant layer due to transport along the droplet surface. D i is the diffusivity of the surfactant within the interface.
The droplet motion and the surface flow, u(θ), are accompanied by a flow pattern within the droplet as well as in the neighboring oil, which can be calculated analytically once u(θ ) is known [23,24]. The corresponding viscous tangential stress exerted on the drop surface must be balanced by the Marangoni stress, grad γ (θ) = Mgrad b(θ), where γ is the surface tension of the surfactant-laden oil/water interface, and M = dγ /db is the Marangoni coefficient of the system. Expanding the bromination density in spherical harmonics, we can then express the velocity field [23,24] at the interface as where C α n denotes Gegenbauer polynomials, and µ is a prefactor containing the liquid viscosities outside and inside the droplet. Inserting this into equation (1) and exploiting the orthogonality relations of Gegenbauer and Legendre polynomials, we obtain for all m > 0 and small b m . We see that the different modes decouple, as far as linear stability is concerned. As long as b 0 M is small enough, the expression in brackets is negative, and the resting state is stable against fluctuations. However, when b 0 M exceeds a critical value, the resting state is unstable, and the droplet spontaneously starts to move. It is straightforward to see that for k < 3D i /R 2 , this happens first for the lowest mode at m = 1. In order to determine the flow profile around a squirming drop, we performed particle image velocimetry using a standard setup (ILA GmbH, Germany). The result as displayed in the bottom panel of figure 2 shows that the flow profile indeed resembles a field as expected for the lowest-order mode, m = 1 [16].

Properties of squirming droplets
Before we turn our attention to the collective behavior of the squirmer droplets introduced above, we have to dwell some more on their properties. To discuss the steady-state droplet velocity, V , we first reconsider the total surfactant coverage, c. This adjusts itself as a balance between the molecular adsorption energy at the water/oil interface and the mutual repulsion of the adsorbed surfactant molecules. The equilibrium coverage thus represents a minimum of the interfacial energy. As a consequence, the interfacial tension will not change to first order if c is varied. Small deviations of c from its equilibrium value, which come about necessarily for any finite V , will thus be replenished by the surrounding oil phase without substantial Marangoni stresses. For a flow pattern such as the one shown in figure 2, for example, we approximately have div u ∝ cos θ . The density of the surfactant in the oil phase thus takes the form ρ ≈ ρ 0 + δρ f (r )cos θ , where δρ ∝ V . Here, f (r ) is some function of the radius. If D m is the diffusivity of the micelles in the oil, the diffusion current, D m ∂ρ/∂r , must balance the depletion rate at the drop interface, c div u = 3V c 2R cos θ [16]. As long as this holds, the only source of appreciable Marangoni stresses is the gradient in bromination density, b(θ). The drop thus keeps speeding up according to equation (4).
This comes to an end when δρ ≈ ρ 0 . The surfactant layer at the leading end of the drop surface can then not be replenished anymore, and c reaches substantially below its equilibrium value. This leads to an increase of surface tension accompanied by a 'backward' Marangoni stress, and thus eventually to a saturation of the velocity. According to the reasoning above, we expect that V ≈ 2ρ 0 D m /3c. There is no literature value for the diffusivity of MO micelles in squalane, but we can estimate it on the basis of the Stokes-Einstein relation assuming the radius of the micelles to be similar to the length of an MO molecule (2.3 nm). Using 36 mPas for the viscosity of the squalane, we obtain D m = 2.6 × 10 −12 m 2 s −1 . We thus predict V /ρ 0 ≈ 0.27 µm s −1 mM l −1 from the simple consideration above. In order to measure V (ρ 0 ), we used a reaction scheme similar to the Belousov-Zhabotinski reaction, with reactant concentrations adjusted so as to prevent chemical oscillations to occur (50 mM sulfuric acid, 28 mM sodium bromate, 400 mM malonic acid and 2.7 mM ferroin) [25,26]. This results in a rather constant bromine release rate in the aqueous phase for an extended period of time. Figure 3 shows the dependence of the droplet velocity as a function of the MO concentration in the oil phase. The initial linear increase is in quantitative agreement with the above prediction (dashed line). As the surfactant density is further increased, the complex exchange processes between the water/oil interface and the micelles will finally become the rate-limiting step. As a consequence, the velocity is expected to level off, which is again in agreement with our data.
It is therefore clear that the bromination profile across the droplet surface, b(θ), will no longer be determined by the fastest growing mode as suggested by equation (4). Instead, it will acquire a complex shape accounting for the nonlinearity of equation (1). It is still unknown how large the nonlinear contributions to b(θ) will be. As a consequence, a quantitative estimate of the driving force, b 0 M, on the basis of the steady-state velocities might easily be orders of magnitude off reality. The task of quantifying the driving force and its relationship with the other quantities characterizing our squirmer is left for future work. During an experiment, the concentration of brMO gradually builds up from the 'exhaust' of the many squirmers in a poorly controllable and spatially inhomogeneous fashion. Since the majority of the brMO becomes trapped in the MO micelles, its effect is mainly to reduce the effective concentration of MO in the oil. Figure 3 tells us that by choosing a high surfactant concentration, where there is almost no sensitivity of the squirmer velocity to the surfactant density, we can also minimize the sensitivity of the velocity to the concentration of brMO, thereby minimizing any unwanted (global) crosstalk between the squirmers. This latter property provides another key feature of our system that makes it particularly suited for experimental studies of collective motion.

Collective behavior
Let us therefore now turn to an investigation of the collective behavior of our squirming droplets. It is a long-standing debate whether physical effects, such as hydrodynamic interactions, are sufficient to explain textures observed in the swarming behavior of bacteria and other micro-organisms, without having to invoke chemotaxis or other genuinely biological effects [2,4,5,12]. Using our model squirmers instead of bacteria, we can tackle this problem from the reverse side, asking for the textures we can observe in dense populations of model squirmers, which are guaranteed to be free of biological interactions. As we will see, there is indeed considerable structure to be unveiled even for such simple systems.
We restrict ourselves to effectively 2D systems here, not only for the sake of simplicity, but also because most of the studies to date have concentrated on the 2D case. Our samples were prepared by creating shallow wells of a few millimeters diameter and a depth just slightly larger than the droplet diameter in poly-dimethyl-siloxane (PDMS) rubber by standard soft lithography techniques. Bonding the PDMS to a glass slide resulted in a flat compartment, which in addition was connected by a narrow channel to a step emulsification unit [27], where droplets were produced and subsequently transported through the channel into the sample well. Figure 4 4 shows a typical sample. The red arrows indicate the direction of motion of each droplet. The angular correlation of the droplet motion, C ϑ , as a function of the scaled distance of the droplet centers, r/d. Note that d = 2R is the droplet diameter. The red curve corresponds to an areal droplet density of 0.46 and the black curve to a density of 0.78. The black arrows are at multiples of 1.08 droplet diameters. The inset shows a semi-logarithmic plot of the decay of the correlation data. The red line corresponds to a decay length of 0.6 droplet diameter and the straight asymptote of the black line represents a decay length of 2.5 droplet diameters. An oscillation with a period of 1.08 droplet diameters (decaying exponentially over 0.9 droplet diameter) has been superimposed to fit the data.
As time proceeds, the formation of rather long-lived clusters of different sizes is observed [28]. The most striking feature, however, is the significant polar alignment of the velocities of neighboring droplets. In order to quantify this alignment, we use the angular correlation function which describes the propensity of velocities of neighboring particles to align with respect to each other. The average runs over all droplet pairs, (i j), as well as over time, t. ϑ is the angle between the velocities of droplets i and j and depends on time for each droplet pair. On the basis of earlier theoretical work, we might under certain circumstances expect a significant polar correlation of the velocities of neighboring droplets [3,4], just from the hydrodynamic interaction of the squirmers with each other. That this is indeed the case can be seen in figure 5. The red curve corresponds to a moderate density of 0.46, which we define as the fraction with respect to the density corresponding to a hexagonal close packing in the plane. We see that there is significant correlation for small distances. More specifically, the angular correlation function decays approximately exponentially away from the contact distance (which is equal to one droplet diameter, d = 2R), as can be seen from the inset. The decay constant is about 0.6 droplet diameters (red line in the inset). The angular correlation thus decays almost completely over one interparticle distance, suggesting that the correlation of the velocities is mediated by the pair interaction of the particles. In fact, it has been predicted theoretically that two adjacent droplets that are propelling themselves by means of low-order spherical harmonic flow fields may attract each other into a bound state in which they are swimming with virtually parallel velocities [3]. We provide here the first experimental corroboration of this prediction using a purely 'physical' system.
There are, however, also pronounced differences to simulation data. Ishikawa and Pedley [4] performed simulations of the collective behavior of spherical droplets with full hydrodynamic interactions, swimming in a monolayer. This system is very similar to ours, but the monolayer in the simulation was freely suspended in the 3D 'liquid', such that there was no nearby wall as in our case. This gives rise to long-range hydrodynamic interactions and provides a straightforward explanation as to why they found velocity correlations ranging up to more than five droplet diameters, in marked contrast to our results.
Correlations become more pronounced as the density of the droplets is increased to an areal density of 0.78. As the black curve in figure 5 shows, we observe two significant changes. Firstly, the range of the correlation becomes significantly larger, extending clearly beyond four droplet diameters. Secondly, we observe the appearance of distinct peaks in the correlation function. The black arrows are at multiples of 1.08 droplet diameters, which is close to what one would expect in the case of lateral 'layering' effects. Also this marked texture, which may be described as ordered rafts, has not been reported previously from simulations of similar systems. An obvious possible reason is the particularly high areal density in our experiment (it was up to 0.5 in the study by Ishikawa and Pedley [4]). As the inset shows, the decay length of the polar correlation is now about 2.5 droplet diameters (slope of the black asymptotic line) and has therefore increased by a factor of four. The amplitude of the superposed oscillations (describing the lateral layering) decays exponentially over 0.9 droplet diameters. These results suggest the presence of a phase transition occurring at a density somewhere between 0.46 and 0.78, at which a qualitative change in the correlation behavior takes place. Experiments are under way to search for this transition.

Conclusions
We have introduced a novel type of artificial squirmer that is particulary well suited for the study of collective phenomena of micro-swimmers. The collective behavior of these squirmers shows strong velocity correlations, which are qualitatively similar to recent simulation results. There are, however, distinct differences, which can be attributed in part to simplifications that had to be made in the simulations. Of particular interest is the occurrence of ordered rafts, with a lateral correlation extending over more than four droplet diameters. We have not yet been able to characterize the flow field of this swimmer precisely. Once this has been achieved, a more detailed comparison with theory and simulation will be possible, including the effects of different spherical harmonics in the flow field of the squirmer [1,4,16].