Experimental results

Figure 1 shows a representative image of our granular flock. Despite their low concentration the rods display a high degree of orientational order. This and further results emerge from our experiments on a mixture of stepwise tapered rods and spherical beads confined between surfaces subjected to rapid vertical vibration (see Fig. 2a and Methods). The system is initialized by introducing the beads into the sample cell and distributing the rods amidst the beads, to yield a homogeneous rod–bead monolayer. We then study the nature of the statistically stationary state of the vibrated rod–bead mixture as a function of the area fractions Φ b and Φ r , defined as the fraction of the surface covered by the two-dimensional projections of the beads and rods, respectively. We monitor the centroid positions r i of all the particles and the tail-to-head unit vectors n i , where i labels the particle in question, and all vectors lie in the horizontal plane.

Figure 1: A granular flock. A monolayer of millimetre-sized tapered brass rods align spontaneously on a vibrating surface amidst a background of aluminium beads. Bead and rod area fractions are 0.51 and 0.16 respectively. Full size image

Figure 2: Structure and kinetics of the disordered and ordered states. (a) Dimensions and confinement geometry of the polar rod used in the experiment, with arrow indicating direction of self-propulsion. (b–e) Typical configuration and corresponding time-evolution of polar order parameter in the disordered phase, bead and rod area fractions Φ b =0.41, Φ r =0.05 (b,c), and in the ordered phase, Φ b =0.68, Φ r =0.05 (d,e); for the simulation in the flower geometry: (f) rod used in simulation, (g–j) as in (b–e): Φ b =0.48, Φ r =0.06 (g,h), Φ b =0.66, Φ r =0.06 (i,j); for the simulation with periodic boundary conditions, configuration and order-parameter evolution for the disordered phase, Φ b =0.20, Φ r =0.11 (k,l) and the ordered phase, Φ b =0.60, Φ r =0.11 (m,n). Full size image

With increasing Φ b (or, less surprisingly, Φ r ) the system is seen to undergo a phase transition from a disordered state to one in which the rod orientation vectors are aligned. Correspondingly, the rods move in random directions in the low-density state (Fig. 2b and Supplementary Movie 1) and coherently at high densities (Fig. 2d and Supplementary Movie 2). At long times the global movement in the ordered phase is a clockwise or anticlockwise circulation with equal probability (Supplementary Fig. 1), with the orientation vectors aligned parallel to the boundary. The direction of circulation once selected does not change for the duration of the experiment.

We have checked that neither crystalline order nor mechanical rigidity of the bead bed is an essential ingredient in promoting the flocking of the rods. For example, the flock in Fig. 1 is on a liquid-like bead background. Further, rods in an amorphous bead bed, generated by using a mixture of bead sizes, are also seen to self-organize into a flock (Supplementary Fig. 2 and Supplementary Note 1). Creating and maintaining a homogeneous state in bidisperse systems is complicated by size-segregation effects, which is why we have continued to work with beads of a single size. The rich physics that enters when the bead medium is a single- or multidomain crystal requires a separate study, now underway.

In a confined geometry the simplest flock consists of particles globally circulating around a common centre, which is the motion seen in our experiments. Accordingly, we resolve n i into local polar-coordinate components , which we use to evaluate the order parameter P≡‹p i › distinguishing isotropic and oriented states, and G(r)=‹p i ·p j › r −P2, which measures correlations of orientational fluctuations about the mean for pairs i, j of rods with separation |r i −r j |=r. The instantaneous in-plane bead-velocity unit vectors u i , transformed to the local polar frame to give , provide another measure of flocking, v≡‹v i ›. The averaging denoted by the angle brackets is carried out over all rods and over all times in the steady state for P, about 100 statistically independent images for G(r) and 25 frames for v (frame rate=40 s−1, and the plate is vibrated at 200 Hz).

The disordered and ordered phases are most clearly distinguished by growth kinetics. Figure 2c shows a typical disordered state with the magnitude P(t)≡|P(t)| of the instantaneous order parameter, averaged over rods but not over time t, fluctuating just above zero, and Fig. 2e presents an ordered state where P(t) grows and saturates to a value close to 1, which implies that polar rods are now moving parallel to the boundary. We construct a phase diagram based on the steady-state value of P as shown in Fig. 3a. For Φ r >0.025 we find that P increases abruptly across the line Φ b +kΦ r ≃0.65, k≃2. We identify this line, provisionally, as the location of the nonequilibrium phase transition between the disordered and ordered states, bearing in mind the effects of finite size and the sample boundary. Correspondingly, Fig. 3c shows the experimentally measured polar order parameter as a function of Φ b for Φ r =0.08. At the smallest values of Φ r the phase boundary bends sharply upwards, possibly implying a true lower limit to the Φ r required to produce a flock regardless of Φ b . At the highest densities arrested states are observed, where rods either jam against the particles of the bead medium or cluster along the boundary (Supplementary Fig. 3a). Figure 3e reports the variation of the steady-state orientational correlations G(r) (scaled by its value at the smallest r) of the rods as the system approaches and then crosses the phase boundary. For Φ b =0.45, spatial correlations consistent with a power-law are seen, while for larger or smaller Φ b , safely within the ordered or the disordered phase, correlations decay more rapidly. These observations can consistently be interpreted as pretransitional fluctuations in the vicinity of a nonequilibrium phase transition, although the system sizes explored are too small to permit a claim about the order of the transition. As effects of finite size limit the accuracy of our estimates of the value of the order parameter and the location of the phase transition, we present in Supplementary Fig. 4a the correlation functions of the orientation without the mean subtracted. The flocking of the polar rods is accompanied by coherent motion of the beads (Fig. 4a) as well, as measured by the bead-velocity order parameter v, which shows the same growth kinetics as P (Fig. 4b).

Figure 3: Flocking phase transition and growth of correlations. (a) Experimental phase diagram in the Φ b −Φ r plane, showing the phase boundary between the isotropic phase and the flock. Clustering at the sample boundary sets in at high concentration. (b) Simulation phase diagram in the flower geometry reproduces the flocking transition. We have not explored concentrations at which boundary clustering is expected. Increase in polar order parameter P as function of Φ b for experiments (c) and simulations (d). The error bars correspond to the standard deviation of P. Orientational correlation function from experiments (e) and simulations in a 55.1-rod-length periodic box (f) at various Φ b at a fixed Φ r . A considerable range of power-law correlations is seen as Φ b increases. Full size image

Figure 4: Correlation between the dynamics of the polar rods and the bead medium. (a) Arrows in the zoom show the coherent velocity field of the bead medium in the ordered phase with Φ b =0.66 and Φ r =0.06. Blank regions are occupied by polar rods. (b) The ordering kinetics of the bead velocity field tracks that of the polar orientation. Full size image

Simulation results

To understand the physics underlying this novel realization of a flocking transition, and to explore which of the many possible control parameters are essential, we create a computer model of tapered rods and spherical beads. For convenience in simulation the rods are constructed by joining overlapping spheres in a straight line (Fig. 2f), but their shape is chosen to resemble as closely as possible the stepwise tapered rods in the experiment. We have also tried other shapes such as a uniform taper and find qualitatively similar results which we do not present here. We carry out two types of studies, one in an enclosure constructed to imitate the experimental geometry, the other with periodic boundary conditions in the horizontal plane. In both cases the simulations, like the experiments, are three-dimensional, with a base and a lid permitting tightly confined motion in the third dimension. The energy input, precisely as in the experiment, is through vertical agitation of the container. The particles move according to the laws of Newtonian mechanics. Collisions of the particles with each other and with the bounding surfaces are governed by inelasticity and static Coulomb friction. We add a small rotational noise to the dynamics of individual rods. This improves the fidelity of our imitation of the experimental system, in which minute imperfections in particle shape and base or lid topography endow each rod with rotational diffusion even when no beads or other rods are present.

Our simulation results can be seen in Figs 2g–j,k–n. In the ‘flower’ geometry we find precisely the behaviour seen in the experiment. Random motion at low area fraction (Supplementary Movie 3) gives way to an ordered, circulating flock at high area fraction (Supplementary Movie 4), and a phase diagram broadly corresponding to that in the experiment is obtained. The numerical experiment, moreover, can be performed in a periodic simulation box, eliminating the effects of sample boundaries. Figure 2k,l (Supplementary Movie 5) and Fig. 2m,n (Supplementary Movie 6) show that the order parameter, now measured in a global Cartesian frame, remains at values consistent with zero for all times at low Φ b , while at high Φ b it grows to saturation. Clear evidence of a flocking transition around Φ b ≃0.35 is seen in Fig. 3d showing polar order parameter as a function of bead area fraction. Consistent with these findings and with the experiments, Fig. 3f shows evidence of power-law correlations for Φ b =0.4, but appreciable curvature on a log–log plot for Φ b well in the ordered or the disordered regime. Increasing Φ b serves to promote orientational correlations for Φ r as low as 0.03 (Supplementary Fig. 5). As in the experiments, we present the ‘unsubtracted’ correlation functions in Supplementary Fig. 4b. For good measure, we also investigate and rule out, in Supplementary Fig. 6 (and Supplementary Note 2), the possibility that the only role played by the beads is to reduce the rotational diffusion of rods. These simulation results establish beyond reasonable doubt that the phenomenon we observe is a phase transition, with the beads playing a direct role in mediating an aligning interaction among the rods.

Finally, on theoretical grounds, in particular because the suspending medium of beads is compressible16, we expect large number fluctuations in our system as in standard flocking models of the Toner-Tu type2. In such systems, large-scale correlated particle currents arise as a result of the coupling to easily excited and slowly decaying variations in the orientation field, both in the ordered phase2 and in the isotropic phase near onset22. In Methods we show that the separate conservation of rods and beads leaves unaltered the arguments2 for anomalously large number fluctuations. We have carried out a limited set of numerical studies of fluctuations in the number density of rods in the ordered and disordered phases, and find an excess in both cases (Fig. 5). Our simulation measurements of the standard deviation ΔN in the number of polar rods in regions containing ‹N› particles on average show that grows with ‹N›. However, more extensive simulations or experiments are required for a useful fit to a power law or other dependence on ‹N›, and to probe the relation of this phenomenon to known sources2,5,22,23,24 of enhanced density fluctuations in active systems.

Figure 5: Enhanced density fluctuations. (a) The scaled standard deviation in the number of polar rods, for regions containing ‹N› rods on average, is seen to grown with ‹N› for area fractions in the ordered (labelled O) and disordered (D) phases. (b) and (c) Representative images from the two phases. Full size image

Theory

We now integrate the insights gained from experiment and simulation to construct a theory. The sequence of events in Fig. 6 also seen in simulation studies, shows that coherently moving polar rods are able to entrain a stray rod initially moving in the opposite direction, through the bead flow they generate, as seen in Supplementary Movies 2, 4 and 6. These observations suggest that the coupling of flow and orientation lies at the heart of the phenomenon of flocking at a distance, as we now show through a minimal, universal hydrodynamic theory whose parameters and their dependence on Φ b are measured in our simulation. We also discuss the connection to the particular case of Bricard et al.16

Figure 6: Flowing beads reorient polar rods. A series of experimental images showing a stray polar rod brought into alignment with the rest of the flock as a consequence of bead flow. Full size image

As we are concerned with large-scale ordering, we work with coarse-grained fields. Our system is a suspension of polar rods in a compressible two-dimensional fluid of beads on a substrate. The variables of interest are therefore the local order parameter P(r, t), given by the average of the orientation vectors p i of the polar rods, and the number densities ρ and σ of beads and rods, respectively, in a small neighbourhood around point r at time t. Below, for simplicity, we consider the limit of small rod concentration, and so treat only the bead density ρ(r, t). A treatment with both densities is in Methods. In addition, we include v(r, t), the average bead-velocity vector in the neighbourhood, as it has a crucial role in mediating the interactions that lead to order. In the dilute-rod limit v can be viewed as the total velocity field of both species.

The equations of motion follow from a few general principles and some key elements of the dynamics of the system. Conservation of particles is embodied in the continuity equation

Newton’s second law locally reads

where we present only the leading-order terms in an expansion in powers of gradients and fields. In (2), Γ>0 damps motion with respect to the substrate and lid, the term in B can be viewed as describing pressure forces, driving flow downhill in density for B>0, and α measures the degree to which orientation of the polar rods gives rise to propulsive flow. We note in passing here that such a term was argued to arise25 in the two-dimensional projection of the three-dimensional hydrodynamics of a suspension of motile organisms confined to a thin layer of fluid. We discuss below the somewhat different origin of the term in the case of our monolayer system. The viscosity η describes transmission of momentum in the plane of the system. The local order parameter P, again to leading order in gradients and fields, obeys

where a measures the rate at which an initial polarization relaxes, presumably through rotational diffusion, and the magnitude and sign of A characterize the tendency of rods to orient parallel or antiparallel to concentration gradients2,26. The parameter λ determines the rate at which a uniform velocity (with respect to the confining walls) aligns the polarization. The parameter K governs spatial variations in the direction and magnitude of the order parameter. The ellipsis in (2) and (3) denotes contributions2 at higher orders in p, v and ∇, discussed in Methods.

The terms with coefficients α and λ are the key players in (2) and (3). No symmetry rules out their existence; a phenomenological theory must therefore include them. However, the physics underlying these terms, especially as regards their signs, merits some discussion. Recall that p i for the ith rod points from the thick to the thin end. Depending on the detailed contact mechanics with the vibrating base, the rod could propel itself on average parallel or antiparallel to p i . Friction would cause the bead medium in the vicinity to be dragged in the direction of motion of the rod. A collection of such rods thus generates a force proportional to the mean polarization P. The parameter α measures the strength of this forcing. The magnitude of α should grow with the mean rod concentration and propulsion speed, and should further depend on the surface characteristics of a rod, which influence its ability to carry the ambient medium with it, by dragging or pushing. A rod with spontaneous motion parallel (antiparallel) to p i has α>0 (α<0). By suitably engineering the geometry, mass distribution, surface properties and contact mechanics of a rod, and exploring a range of vibration parameters, it should be possible to make polar rods with a range of magnitudes and either sign of α; the rods we work with here have α>0.

The mechanics underlying λ is reminiscent of the interaction of wind with a weathercock. Imagine a single tapered rod lying athwart a uniform flow of particles over a surface. Although the momentum transfer is distributed uniformly across the rod, its mass is concentrated, and hence its pivot point shifted, towards one end. This leads to a net torque due to the flow, and a consequent rotation as shown in Fig. 7a. Geometrical arguments alone cannot determine the sign of λ and hence the sense of rotation, and self-propelling activity could influence its value. It should be possible to engineer λ<0 by using intrinsically higher-density material at the thin end, or a hollow interior at the thick end.

Figure 7: Weathercock effect and location of the phase transition. (a) The direction in which the polar rod is rotated by a flow depends on its natural pivot point. (b) α/Γ is seen to grow and a/λ to decrease with increasing Φ b . Thus the flow velocity generated by a local polarization as well as the aligning power of bead flow grow with increasing bead concentration. The mean-field estimate of the flocking transition corresponds to the intersection of the two curves, marked by stars, where ā in Equation (4) changes sign from positive to negative as Φ r +Φ b is increased. Full size image

The availability of a numerical simulation model provides independent measurements of all phenomenological parameters of importance to our theory. To measure λ, we place a single polar rod initially pointing in the x direction, and impose a flow of beads with velocity v in the y direction (see Supplementary Movie 7). By measuring the initial growth rate we estimate the parameter λ in (3), which we average over 50 such trials. Independently, by measuring the rotational relaxation of a single rod, we infer a. In a separate simulation, we impose a nonzero mean P by applying an orienting field on a collection of rods. Measuring the saturation value of the resulting macroscopic velocity yields α/Γ in (2). The damping Γ is readily estimated from the relaxation of an imposed initial velocity. Figure 7b shows the dependence of these key parameters on area fraction ρ 0 . We comment in Methods on the parameter A, which has only a minor role in our analysis.

We ignore density variations in (3) and (2), replacing ρ by its mean value ρ 0 , which we treat as a parameter on which the coefficients depend. Working with the spatial Fourier components v q and P q of the velocity and order parameter at wavevector q, (2) then tells us that, on timescales long compared to ρ 0 /Γ, v q →α P q /(Γ+ηq2). Substituting this value for v in (3) we obtain, to order q2,

where ā=a−λα/Γ and . The mechanism for spontaneous ordering now becomes clear. As λ is an increasing function of ρ 0 , ā can turn negative even if a is positive. The state P=0 is then linearly unstable, and a state with nonzero P on macroscopic scales will set in. Let us now test our theory using our simulation measurements of the parameters. First, for the proposed mechanism to work, ā must be less than a, that is, we predict λα/Γ>0. Our simulations show that a rod drags beads in the direction of its own motion and that the flow of beads rotates neighbouring rods to point in, and thus move in, the direction of the flow. Thus α and λ have the same sign, confirming the prediction. Second, Fig. 7b shows that as Φ r +Φ b is increased, α/Γ (for three values of Φ r ) increases and a/λ decreases, and the two quantities cross (that is, ā crosses zero) at a value of Φ r +Φ b , which decreases with increasing Φ r . Therefore the theory predicts a phase transition to an ordered phase with increasing area fraction, and negative slope to the phase boundary in the Φ r −Φ b plane, both of which are borne out by experiment. The intersection points in Fig. 7b, corresponding to the mean-field estimate of the threshold area fraction, are, not surprisingly, somewhat lower than the measured values in the simulation, Fig. 3b.

Within our mean-field treatment ((4) or the parent equations (3) and (2)) it is straightforward to see that the transition has a continuous onset and a diverging correlation length . Our limited observations in experiment and simulation show an increasing correlation length upon approach to the transition. Experience with standard flocking models would suggest27 instead a discontinuous transition, albeit with appreciable pretransitional fluctuations. Larger-scale simulation studies including finite-size scaling analyses, now in progress, will determine the character and spatial structure of the onset of order for our rod–bead flock. Note that orientability, motility and coupling to flow, rather than specifics of shape, were the essential ingredients of our theory. On symmetry grounds one cannot therefore rule out the possibility that the self-propelled disks of ref. 10, in a background of non-motile particles, could flock through the mechanism reported here—provided flow can reorient the disks in the manner required.