Spontaneous aggregation and global polar ordering in squirmer suspensions

We have developed numerical simulations of three dimensional suspensions of active particles to characterize the capabilities of the hydrodynamic stresses induced by active swimmers to promote global order and emergent structures in active suspensions. We have considered squirmer suspensions embedded in a fluid modeled under a Lattice Boltzmann scheme. We have found that active stresses play a central role to decorrelate the collective motion of squirmers and that contractile squirmers develop significant aggregates.


Introduction
Collective motion can be observed at a variety of scales, ranging from herds of large to bacteria colonies or the active motion of organelles inside cells. Despite the long standing interest of the wide implications of collective motion in biology, engineering and medicine (as for example, the ethological implications of the signals exchanged between moving animals, the evolutionary benefits of moving in groups for individuals and for species, the design of robots which can accomplish a cooperative tasks without central control, the understanding of tumor growth or wound healing to mention a few), only recently there has been a growing interest in characterizing such global behavior from a statistical mechanics perspective [1].
Although a variety ingredients and mechanisms have been reported to describe the signaling and cooperation among individuals which move collectively, it is important to understand the underlying, basic physical principles that can provide simple means of cooperation and can lead to emerging patterns and structures [2]. We want to analyze the capabilities of basic physical ingredients to generate emerging structures in active particles which self propel in an embedding fluid medium. These systems constitute an example of active fluids, systems which generate stresses by the conversion of chemical into mechanical energy. To this end, we will consider model suspensions of swimming particles (building on the squirmer model introduced by Lighthill [3]) and will analyze a hydrodynamicallycontroled route to flocking. We will use a hybrid description of an active suspension, which combines the individual dynamics of spherical swimmers with a kinetic model for * falarcon@ffn.ub.es * * ipagonabarraga@ub.edu the solvent. We can identify the emergence of global orientational order and correlate it with the formation of spontaneous structures where squirmers aggregate and form flocks of entities that swim along together. This simplified approach allows us to identify the role of active stresses and self-propulsion to lead both to global orientational order and aggregate formation. Even if in real systems other factors can also control the interaction and collective behaviors of active suspensions, the present description shows that hydrodynamics itself is enough to promote cooperation in these systems which are intrinsically out of equilibrium. This work is organized as follows. In section 2.1 we present the theoretical frame of the simulation technique that we have applied, while in section 2.2 we describe the squirmer model that we have used and introduce the relevant parameters which characterize its hydrodynamic behavior and in section 2.3 we give a detailed explanation of the simulation parameters and the systems we have studied. Section 3 is devoted to analyze the global polar order parameter and to study quantitatively the orientation that squirmer suspensions display. In Section 4 flocking is studied via generalized radial distribution functions, moreover to characterize the time evolution of the formed flocks, we calculated the time correlation function of the density fluctuations, the results are shown in this section also. We conclude in Section 5 indicating the main results and their implications.

Lattice Boltzmann Scheme
We consider a model for microswimmer suspensions composed by spherical particles embedded in a fluid. The fluid is modeled using a Lattice Boltzmann approach. Accordingly, the solvent is described in terms of a distribution function f i ( r; t) in each node of the lattice. The distribution function evolves at discrete time steps, ∆t, following the lattice Boltzmann equation (LBE): that can be regarded as the space and time discretized analog of the Boltzmann equation. It includes both the streaming to the neighbouring nodes, which corresponds to the advection of the fluid due to its own velocity, and the relaxation toward a prescribed equilibrium distribution function f eq j . This relaxation is determined by the linear collision operator Ω ij [4,5,6]. It corresponds to linearizing the collision operator of the Boltzmann equation.
If Ω ij has one single eigenvalue, the method corresponds to the kinetic model introduced by Bhatnagar-Gross-Crook (BGK) [7]. The LBE satisfies the Navier-Stokes equations at large scales. In all our simulations we use units such that the mass of the nodes, the lattice spacing and the time step ∆t are unity and the viscosity is 1/2, the lattice geometry that we have used was a cubic lattice with 19 allowed velocities, better known as D 3 Q 19 scheme [5].
The linearity and locality of LBE makes it a useful method to study the dynamic of fluids under complex geometries, as is the case when dealing with particulate suspensions. Using the distribution function as the central dynamic quantity makes it possible to express the fluid/solid boundary conditions as local rules. Hence, stick boundary conditions can be enforced through bounce-back of the distribution, f i ( r; t), on the links joining fluid nodes and lattice nodes inside the shell which defines the solid particles, also known as boundary links [8]. A microswimmer is modeled as a spherical shell larger than the lattice spacing. Following the standard procedure, the microswimmer is represented by the boundary links which define its surface. Accounting for the cumulative bounce back of all boundary links allows to extract the net force and torque acting on the suspended particle [9]. The particle dynamics can then be described individually and particles do not overlap due to a repulsive, short-range interaction among them, given by where ǫ is the energy scale, and σ the characteristic width. The steepness of the potential is set by the exponent ν 0 . In all cases we have used ǫ = 1.0, σ = 0.5 and ν 0 = 2.0.

Squirmer Model.
We follow the model proposed by Lighthill [3], subsequently improved by Blake [10], for ciliated microorganisms. In this approach, appropriate boundary conditions to the Stokes equation on the surface of the spherical particles (of radius R) are imposed to induce a slip velocity between the fluid and the particles. This slip velocity determines how the particle can displace in the embedding solvent in the absence of a net force or torque. For axisymmetric motion of a spherical swimmer, the radial, v r and tangential, v θ components of the slip velocity can be generically expressed as n-th at the squirmer spherical surface, where P n stands for the n-th order Legendre polynomial and V n is define as e 1 describes the intrinsic director, which moves rigidly with the particle and determines the direction along which a single squirmer will displace, while r 1 represents the position vector with respect to the squirmer's center, which is always pointing the particle surface and thus |r 1 | = R.
Since the squirmer is moving in an inertialess media, the velocity u and pressure p of the fluid are given by the Stokes and continuity equations The velocity field generated by a squirmer is the solution of these equations (5) under the boundary conditions specified by the slip velocity in the surface of its body, eq. (3). We will disregard the radial changes of the squirming motion, and will consider A n = 0, to focus on a simple model that captures the relevant hydrodynamic features associated to squirmer swimming. Accordingly, we will also disregard the time dependence of the coefficients B n and will focus on the mean velocity of a squirmer during a beating period [11]. Hence, from the solution of eqs. (5) using the slip velocity as a boundary condition (eq. (3)), we can write the mean fluid flow induced by a minimal squirmer as where we have taken B n = 0, n > 2, keeping only the first two terms in the general expression for the slip velocity, Eq.( 3). The two non-vanishing terms account for the leading dynamics effects associates to the squirmers. While B 1 determines the squirmer velocity, along e 1 , and controls its polarity, B 2 stands for the apolar stresses that are generated by the surface waves [12]. The dimensionless parameter β ≡ B 2 /B 1 quantifies the relative relevance of apolar stresses against squirmer polarity. The sign of β (determined by that of B 2 ) classifies contractile squirmers ( or pullers) with β > 0 and extensile squirmers (or pushers) when β < 0. The limiting case when B 1 = 0 corresponds to completely apolar squirmers (or shakers [13]) which induce fluid motion around them without self-propulsion. The opposite situation, when B 2 = 0 corresponds to completely polar, self-propelling, squirmers which do not generate active stresses around them. We will disregard thermal fluctuations; therefore B 1 and B 2 are the two parameters which completely characterize squirmer motion.

Simulation Details.
All the results that we will discuss correspond to numerical simulations consisting of N identical spherical particles in a cubic box of volume L 3 with periodic boundary conditions. In all cases we have considered N = 2000, R = 2.3 and L = 100 (expressed in terms of the lattice spacing). This corresponds to a volume fraction φ = 4πN R 3 /(3L 3 ) = 1/10, with a kinematic viscosity of ν = 1/2 (in lattice units) [14]. As we will analyze subsequently, active stresses play a significant role in the structures that squirmers develop when swimming collectively. In Fig. 1 we compare characteristic configurations of suspensions for completely polar, contractile and extensile squirmers. Apolar stresses favor fluctuations in the squirmer concentration and for contractile squirmers there is a clear tendency to form transient, but marked, aggregates. The figure also shows that one needs to distinguish between how squirmers align to swim together and how do they distribute spatially. In the following section we will analyze how active stresses interact with self-propulsion to affect both aspects of collective swimming.

Polar Order Parameter.
In order to quantify the degree of ordering associated to collective squirmer motion, we have computed the global polar order parameter (eq. 7) [15], expressed in terms of the squirmer intrinsic orientation e, which determines the direction of swimming for isolated squirmers, In Fig. 2 we show the temporal evolution of P (t) as a function of time for completely polar, contractile and extensile suspensions. The time is normalized by t 0 which is the time that a single squirmer needs to self-propel a distance of one diameter, t 0 ≡ 2R/(2/3 B 1 ) = 3R/B 1 The three suspensions start from a completely aligned initial configuration where squirmers are homogeneously distributed spatially. This figure shows clearly that squirmers relax from the given initial configuration to the appropriate steady state and that active stresses have a profound impact on the ability of squirmers to swim together. The limiting situation of completely polar swimmers, β = 0, keeps almost perfect ordering. This is because the irrotational flow generated by the translational velocity of the particles is strong enough to maintain a symmetrical  [16] with the Normal Mode Wizard (NMWiz) plugin [17].
distortion in the fluid. Hence, a value of P (t) close to one indicates high polarity. The other two curves, corresponding to extensile (β = −1/2) and contractile squirmers (β = 1/2) , indicate that active stresses generically decorrelate squirmer motion due to the coupling of the intrinsic direction of squirmer self-propulsion with the local vorticity field induced by the active stresses generated by neighbouring squirmers. However, we do observe a clear difference because extensile squirmers have completely lost their common degree of swimming while contractile ones still conserve a partial degree of global coherence. In order to quantify in more detail the role of active stresses in the global degree of ordering in squirmer suspensions, we have computed the steady-state value of the polar order parameter, P ∞ , as a function of the relative apolar stress strength, β. Fig. 3 displays P ∞ , computed as the mean average of P (t) over the time period after the initial decay from the aligned state [15]. There are two remarkable observations of the results shown in Fig. 3. First of all, the larger |β| the smaller values of P ∞ observed, which indicate less squirmer coherence due to hydrodynamic interactions controlled by the induced active stresses, or |β|. Secondly, for a given magnitude of the apolar stress, |β|, pullers are more ordered than pushers. Hence, there is an asymmetry between pullers and pushers. This asymmetry can be explained in terms of the differences in the near-field interactions between squirmers [15,18]. Squirmer self-propulsion favors head-to-tail collisions [19] and generates an internal structure that competes with the tendency of squirmers to rotate due to local flows. In fact, head-to-head orientation is stable to rotations for pusher suspensions (as can be clearly appreciated in the last snapshot of Fig.  1, where we can see a lot of pushers interacting head-to-head). In this case, the active stresses favor head-tohead configurations, which competes with self-propulsion and decorrelates faster the comoving swimming configurations of squirmers. On the contrary, the stresslet generated by pullers destabilizes head-to-head configurations favoring the motion of squirmers along a common director. It is worth noting that puller suspensions with β > 3 will evolve to isotropic configurations, in agreement with the long-time polar order parameter displayed in Fig. 3. In order to clarify that global ordering is generic for squirmers composed of spherical particles, and hence that orientation instabilities do not require non spherical propelling particles [20], we have analyzed the collective evolution of squirmer suspensions with initial isotropic configurations. It is clear in Fig. 4.a, that both cases of puller suspensions either initially aligned or isotropic, have a similar long-time polar order; hence we can infer that puller suspensions in either an isotropic or aligned state are unstable and that the steady state is independent of the symmetry of the initial configurations.

P(t)
On fig. 4.b one can clearly appreciate that isotropic puller suspensions (red circles) are also unstable, as shown in Fig. 4.a. On the contrary, isotropic pushers suspensions are stable (black circles) for this regime of β. Similarly to the result for puller suspensions showed in Fig. 4.a, one can appreciate in Fig. 4.b that pushers are driven to the same long-time polar order parameter, and therefore that the final alignment is independent of the initial configuration.

4.
Flocking. Fig. 1 shows that puller suspensions, (β > 0), display a cluster of the size of the box. Due to the absence of attractive forces between squirmers, these observed clusters are statistically relevant but have a dynamic character. As a function of time the observed aggregates evolve and displace; the particles they are form with change. We need then a statistical approach to analyze the formation of emergent mesoscale structures and its correlation with orientational ordering. We have computed the temporal correlation function of the density fluctuations dividing the simulation box in 1000 sub-boxes of side box l = L/10 and counted all the particles N i (t) at each i-th sub-box. This provides the particle temporal mean number, N i (t) t , from which we can determine the instantaneous density fluctuations, δN i (t) = N i (t) − N i (t) t , at each box. The average density fluctuation, δN (t), can then be derived as the mean of δN i (t) over all the sub-boxes at time t, and one can use them to study their temporal correlation. The time correlation of the squirmer density fluctuations, depicted in Fig. 5, shows that pullers have an oscillatory response, associated to the displacement of aggregates with a density markedly above average, while pushers are characterized by a more homogeneous spatial distribution. We can gain more detailed insight into the aggregation and ordering of squirmer suspensions by studying the generalized radial distribution functions [6] g n (r) ≡ P n (cos θ ij ) , where θ ij stands for the relative angle between the direction of motion of the particles i and j at a distance between r and r + dr and P n is the n-th degree Legendre polynomial. For n = 0 we recover the radial distribution function, g 0 (r). The average in eq. (8) is taken over all particle pairs and over time, once the system has reached its steady state. Fig. 6 displays g 0 (r) for three kinds of squirmers, β = {0, 1/2, −1/2}. For comparison, we also show the radial distribution function of a randomly distributed configuration, which constitutes a good approximation for the equilibrium radial distribution function for hard spheres at φ = 1/10. Fig. 6 displays also g 0 (r) for β = −1/5. This case corresponds to a pusher suspension with the same polar order value, P ∞ , than the puller suspension at β = 1/2 and will help to analyze the correlation between global polar order and the suspension structure. One can clearly appreciate that activity enhances significantly the value of the radial distribution at contact, g 0 (r = 2R), compared with the corresponding value for an equilibrium suspension. This value is larger for puller suspensions indicating the larger tendency of pullers to remain closer to each other. The radial distribution function for pullers develops a marked second maximum at r = 4.25R indicating the development of stronger short range structures for pullers. Neither pushers nor totally polar squirmers have a visible second maximum even when we compare puller and pusher suspensions with equivalent polar order parameter, P ∞ . The development of the secondary peak for pullers is consistent with their tendency to form large aggregates, or flocks, in agreement with the snapshot depicted in Fig.1. g 0 (r), β= 0.5 g 0 (r), β= -0.5 g 0 (r), β= -0.2 g 0 (r), β= 0.0 g 0 rndm (r;t=0)   7 displays the generalized radial distribution function, g 1 (r), which provides information on the degree of local correlated polar order around a given squirmer. Initially, all squirmers are parallel, and hence g 1 (r, t = 0) = 1.0 (green diamonds in the Figure). The isotropic initial condition (yellow circles), when g 1 (r, t = 0) = 0, is also shown as a reference. Completely polar squirmer suspensions, β = 0, keep g 1 (r) very close to 1 (violet triangles) showing that most of the particles swim along a common direction even if they are far away from each other; this strong correlation is easily appreciated in the first snapshot of Fig.1. We can observe a similar effect for pusher suspensions at β = −1/5 where we can see how g 1 (r) relaxes to a finite plateau for r > 3R. However, unlike completely polar squirmers, now g 1 (r > 3R) ∼ 0.6 (black diamonds) indicating a loss of coherence in the swimming suspension. The relative alignment for puller suspensions is clearly different, because g 1 (r) decays asymptotically to zero (blue squares) for separations analogous to those on which the radial distribution function decays to one. This indicates that the structure we have identified through g 0 (r) in Fig. 6 corresponds to groups of nearby particles that swim along the same direction. This behavior is consistent with the middle snapshot of Fig. 1 which shows a marked flocking formed by a significant number of particles swimming coherently in the same direction. If the apolar strength is increased, increasing the magnitude of β, for pusher suspensions, the partial coherence that we have seen in the case of β = −1/5 vanishes. The curve of g 1 (r) for β = −1/2 (red triangles) does not display any significant feature, indicating a complete decorrelation in the direction of swimmers at all length scales. The corresponding configuration in Fig. 1 shows clearly the absence of any significant correlated orientation between squirmers. g 1 (r), β= 0.5 g 1 (r), β= -0.5 g 1 (r), β= -0.2 g 1 (r), β= 0.0 g 1 rndm (r;t=0). Isotropic g 1 rndm (r;t=0). Alligned

Conclusions
We have analyzed a model system of swimming spherical particles to show the capabilities of the hydrodynamic coupling as a route to pattern formation, polar ordering and flocking in the absence of any additional interaction among the swimmers (except that swimmers cannot overlap due to excluded volume). We have shown how a numerical mesoscopic model for swimmer suspensions can develop instabilities and long-time polar order and that active stresses play a relevant role to promote flocking due to the coupling of the swimming director with the local fluid vorticity induced by the neighboring squirmers. We have identified the sign of such active stress (which distinguishes pullers from pushers) as the main element which controls squirmer flocking and swimming coherence.
We have shown that spherical squirmers, starting from aligned or isotropic state, develop a unique long-time polar order due to hydrodynamic interactions. We have found that aligned pushers suspension are unstable while isotropic suspensions are stable for β < −2/5: isotropic puller suspensions are also stable for β > 3.0.
We have seen that flocking configurations for pullers leads to large elongated structures, reminiscent of the bands observed in the Vicsek model [21]. However, in this later case hydrodynamics is absent and flocking develops at high concentrations, when the aligning interaction is strong enough to overcome decoherence induced by noise. In the systems we have explored the coherence is hydrodynamic and develops at small volume fractions. The observed elongated, spanning aggregates with internal coherent orientation, in the range 0 < β < 1, are robust and independent of the initial configuration.