# Structure and Dynamics of a Phase-Separating Active Colloidal Fluid

###### Abstract

We examine a minimal model for an active colloidal fluid in the form of self-propelled Brownian spheres that interact purely through excluded volume with no aligning interaction. Using simulations and analytic modeling, we quantify the phase diagram and separation kinetics. We show that this nonequilibrium active system undergoes an analog of an equilibrium continuous phase transition, with a binodal curve beneath which the system separates into dense and dilute phases whose concentrations depend only on activity. The dense phase is a unique material that we call an active solid, which exhibits the structural signatures of a crystalline solid near the crystal-hexatic transition point, and anomalous dynamics including superdiffusive motion on intermediate timescales.

[Supplement-]supplement

Active fluids composed of self-propelled units occur in nature on many scales ranging from cytoskeletal filaments and bacterial suspensions to macroscopic entities such as insects, fish and birds Vicsek and Zafeiris (2010). These systems exhibit strange and exciting phenomena such as dynamical self regulation Gopinath et al. (2011), clustering Peruani et al. (2006), anomalous density fluctuations Ramaswamy et al. (2003), unusual rheological behavior Giomi et al. (2010); Saintillan (2010); Cates et al. (2008), and activity-dependent phase boundary changes Shen and Wolynes (2004). Motivated by these findings, recent experiments have focused on realizing active fluids in nonliving systems, using chemically propelled particles undergoing self-diffusophoresis Palacci et al. (2010); Paxton et al. (2004); Hong et al. (2007), Janus particles undergoing thermophoresis Jiang et al. (2010); Volpe et al. (2011), as well as vibrated monolayers of granular particles Narayan et al. (2007); Kudrolli et al. (2008); Deseigne et al. (2010).

In this letter we explore a minimal active fluid model: a system of self-propelled smooth spheres interacting by excluded volume alone and confined to two dimensions. Unlike self-propelled rods Baskaran and Marchetti (2008a, b); Peruani et al. (2011); Yang et al. (2010); McCandlish et al. (2012), these particles cannot interchange angular momentum and thus lack a mutual alignment mechanism. Recent simulation and experimental studies have shown that this system exhibits giant number fluctuations Fily and Marchetti (2012) and athermal phase separation Fily and Marchetti (2012); Theurkauff et al. (2012) that are characteristic of active fluids Ramaswamy et al. (2003); Mishra and Ramaswamy (2006); Cates et al. (2010). Here we employ extensive Brownian dynamics simulations to characterize the phase diagram of this system and we develop an analytic model that captures its essential features. We show that this nonequilibrium system undergoes a continuous phase transition, analogous to that of equilibrium systems with attractive interactions, and that the phase separation kinetics demonstrate equilibrium-like coarsening. These structural and dynamic signatures of phase separation and coexistence enable an unequivocal definition of phases in this nonequilibrium, active system. Finally, we find that the dense phase is a dynamic new form of material that we call an “active solid”. This material exhibits structural properties consistent with a 2D colloidal crystal near the crystal-hexatic transition point Peng et al. (2010); Han et al. (2008), but is characterized by such anomalous features as superdiffusive transport at intermediate timescales and a heterogeneous and dynamic stress distribution (see Fig. (1)).

Model and Simulation Method: Our system consists of smooth spheres immersed in a solvent and confined to a plane, similar to experimental systems of self-propelled colloids sedimented at an interface Theurkauff et al. (2012). Each particle is self-propelled with a constant force, and interactions between particles result from isotropic excluded-volume repulsion only. We include no mechanism for explicit alignment or transmission of torques between particles.

The state of the system is represented by the positions and self-propulsion directions of all particles. Their evolution is governed by the coupled overdamped Langevin equations:

(1) | ||||

(2) |

Here is an excluded-volume repulsive force given by the WCA potential if , and zero otherwise Weeks et al. (1971), with the nominal particle diameter. We use , but our results should be insensitive to the exact strength and form of the potential. is the magnitude of the self-propulsion force which in the absence of interactions will move a particle with speed , , and . and are translational and rotational diffusion constants, which in the low-Reynolds-number regime are related by . The are Gaussian white noise variables with and .

We non-dimensionalized the equations of motion using and as basic units of length and energy, and as the unit of time. Simulations employed the stochastic Runge-Kutta method Brańka and Heyes (1999) with maximum timestep . Simulations mapping the phase diagram were run with particles until time , while larger systems (up to particles) were used to explore kinetics and material properties. The simulation box was square with periodic boundaries, with its size chosen to achieve the desired density. The system is parametrized by two dimensionless values, the packing fraction and the Péclet number, which in our units is identical to the non-dimensionalized velocity (). In this work, we varied from near-zero to the hard-sphere close-packing value , and Pe from zero to .

Phase Separation: We first show that our results are consistent with prior simulations Fily and Marchetti (2012) and confirm that this system, despite the absence of aligning interactions, shows the signature behaviors of an active fluid. In particular, the active spheres undergo nonequilibrium clustering (Fig. (1)) similar to other model active systems Chaté et al. (2008); Peruani et al. (2006); McCandlish et al. (2012); Yang et al. (2010).

We establish that this clustering is indeed athermal phase separation by measuring the density in each phase at different parameter values (Fig. (2a)). We observe a binodal envelope beyond which the system separates into two phases whose densities collapse onto a single coexistence curve which is a function of activity alone. The phase diagram is thus analogous to that of an equilibrium system of mutually attracting particles undergoing phase separation, with Pe (playing the role of an attraction strength) as the control parameter. This surprising result contradicts the expectation that increased activity will destabilize aggregates and suppress phase separation (as seen in Schwarz-Linek et al. (2012)) and indicates that the effects of activity cannot be described by an “effective temperature” in this system.

Additionally, we identify a critical point at the apex of the bimodal (near , ). In the vicinity of this point, the system exhibits equilibrium-like critical phenomena which will be detailed in a future publication.

The Phase-Separated Steady State: To characterize the steady state, we measured the fraction of particles in the dense phase at time (Fig. (3)). In contrast with recent work Fily and Marchetti (2012) which placed the phase transition boundary at a constant density, we observe that this cluster fraction is a nontrivial function of the system parameters . To understand this relationship we developed a minimal model in which this function can be found analytically. Let us assume the steady state contains a macroscopic cluster which we take to be close-packed. Particles in the cluster are stationary in space but their continue to evolve diffusively. We treat the gas as homogeneous and isotropic, and assume that a particle colliding with the cluster surface is immediately absorbed.

Within this model, we can write the rate of absorption of particles of orientation from the gas phase as , where is the gas number density. Integrating yields the total incoming flux per unit length: . To estimate the rate of evaporation, note that a particle on the cluster surface will remain there so long as its self-propulsion direction remains “below the horizon”, i.e., , where is normal to the surface. When its direction moves above the horizon, it immediately escapes and joins the gas. This rate can be calculated by solving the diffusion equation in angular space with absorbing boundaries (for clusters large enough to treat the interface as flat, at ) and initial condition given by the distribution of incident particles: , with and . Further, the departure of a surface particle creates a hole through which subsurface particles (whose may point outwards) can escape. With we denote the average total number of particles lost per escape event, which we treat as a fitting parameter. The total outgoing rate is then .

Equating and yields a steady-state condition for the gas density: . can be eliminated in favor of , yielding (in terms of our dimensionless parameters):

(3) |

This function is plotted in Fig. (3) with , in good accord with our simulation results. Further, the condition allows us to deduce a criterion for the onset of clustering. Restoring dimensional quantities, this condition gives . Note that is a collision frequency; thus the system begins to cluster at parameters for which the collision time becomes shorter than the rotational diffusion time.

The mechanism we have presented here is purely kinetic and requires only an intuitive picture of local dynamics at the interface. An alternative view has been described by Tailleur and Cates Tailleur and Cates (2008); Cates and Tailleur (2012) who subsume all interactions into a density-dependent propulsion velocity which decreases with density as collisions become more frequent. From this they construct an effective free energy which shows an instability in the homogeneous phase if falls quickly enough. In a sense our kinetic model represents an extreme case of this picture in which contains a step function such that free particles are noninteracting, and particles in a cluster are completely trapped (see Fig. LABEL:Supplement-fig:d_eff in Sup ()).

Structure of the Dense Phase: Since the system is composed of monodisperse spheres, the dense phase is susceptible to crystallization Bialké et al. (2012). As shown in Fig. (1) the static structure factor of the cluster interior shows a liquid-like isotropy at low Pe, but develops strong sixfold symmetry as activity is increased. Further, the radial distribution function shows clear peaks at the sites of a hexagonal lattice (see Fig. LABEL:Supplement-fig:radial_distribution in Sup ()) which sharpen and increase in number as Pe is raised. We also measured the bond-orientational order parameter , where runs over the neighbors of particle (defined as being closer than a threshold distance), and is the angle between the - bond and an arbitrary axis (Fig. (4)). We find a structure characterized by large regions of high order with embedded defects that are predominantly 5-7 pairs (Fig. (4a) inset and LABEL:Supplement-fig:defects-movie in Sup ()). Next, we examined the correlation function (Fig. 4) which exhibits a liquid-like exponential decay for systems of low activity, while at higher activity the decay slows to a power law which is indicative of a hexatic Nelson (2002). A further transition to a crystal-like plateau is observable in larger systems (see Fig. (4) and LABEL:Supplement-fig:q6correlation-size and LABEL:Supplement-fig:defect-density in Sup ()). In all cases, this material is unique in that it is held together by active forces alone, and that the arrest of motion is due to frustration. In this sense it is similar to amorphous materials such as granular packs as reflected by the highly heterogeneous stress distribution (Fig. (1)) Blumenfeld and Edwards (2009).

Dynamics in the Dense Phase: Within the active solid material, self-propulsion forces continuously evolve by rotational diffusion, breaking local force balance and leading to defect formation and migration (see LABEL:Supplement-fig:defects-movie in Sup ()). A compelling way to view the motion produced by this athermal process is a simulated FRAP experiment van Blaaderen et al. (1992), in which particles within a contiguous region are tagged, making subsequent mingling of tagged and untagged particles visible (see LABEL:Supplement-fig:frap-2.avi in Sup ()). To quantify this behavior, we measured the mean square displacement (MSD) of particles in the cluster interior. As shown in Fig. (1), we observe subdiffusive motion on short timescales, followed by a superdiffusive regime, returning to diffusive motion on long timescales. The exponents of the subdiffusive and superdiffusive motion ( and , respectively) are well-conserved across a wide range of propulsion strengths. Note that an isolated self-propelled particle will exhibit diffusive, ballistic and diffusive behavior on time scales , and respectively (see Fig. LABEL:Supplement-fig:msd in Sup ()). These dynamical regimes are modified by the active solid environment; in particular, the ballistic regime is modulated by “sticking” events as the particle is localized in crystal domains, resulting in the observed Lévy-flight-like behavior Klafter et al. (1996); Pfister and Scher (1978).

Kinetics of Phase Separation: Despite the athermal origins of phase separation in this system, simulations quenched to parameters within the binodal experience familiar phase separation kinetics (Fig. (5)). Systems quenched close to the binodal exhibit a nucleation delay which can be long enough that artificial seeding is necessary for phase separation to be computationally accessible. Systems quenched more deeply undergo spinodal decomposition, leading to a coarsening regime in which the mean cluster size scales surprisingly as , with a corresponding length scale (Fig. (5) inset, also see LABEL:Supplement-fig:coarsening in Sup ()). This differs from the standard 2D coarsening exponents, but matches recent simulation results for the Vicsek model and related active systems Dey et al. (2012). This result should be viewed as preliminary due to the limited range of our data, but nevertheless this unexpected similarity between the coarsening of point-particles with polar alignment and that of spheres with no alignment suggests a deep relationship between these very different types of systems. Future work is needed to uncover the origins of these scaling exponents and their implications for universality in active fluids.

Summary: A fluid of self-propelled colloidal spheres exhibits the athermal phase separation that is intrinsic to active fluids and is a primary mechanism leading to emergent structures in diverse systems Gopinath et al. (2011); Cates et al. (2010). We have shown that the physics underlying this phase behavior can be understood in terms of microscopic parameters. From a practical perspective, our simulations show that the active solid dense phase exhibits a combination of structural and transport properties not achievable in a traditional passive material. Further development of experimental realizations of this system (e.g. Ref. Theurkauff et al. (2012)) will advance the development of materials whose phase behavior, rheology, and transport properties can be precisely controlled by activity level.

Acknowledgments: This work was supported by NSF-MRSEC-0820492 (GSR, MFH, AB), as well as NSF-DMR-1149266 and NSF-1066293 and the hospitality of the Aspen Center for Physics (AB). Computational support was provided by the Brandeis HPC.

## References

- Vicsek and Zafeiris (2010) T. Vicsek and A. Zafeiris, (2010), arXiv:1010.5017 .
- Gopinath et al. (2011) A. Gopinath, M. F. Hagan, M. C. Marchetti, and A. Baskaran, (2011), arXiv:1112.6011 .
- Peruani et al. (2006) F. Peruani, A. Deutsch, and M. Bär, Phys. Rev. E 74, 030904 (2006).
- Ramaswamy et al. (2003) S. Ramaswamy, R. A. Simha, and J. Toner, Europhys. Lett. 62, 196 (2003).
- Giomi et al. (2010) L. Giomi, T. B. Liverpool, and M. C. Marchetti, Phys. Rev. E 81, 051908 (2010).
- Saintillan (2010) D. Saintillan, Phys. Rev. E 81, 056307 (2010).
- Cates et al. (2008) M. E. Cates, S. M. Fielding, D. Marenduzzo, E. Orlandini, and J. M. Yeomans, Phys. Rev. Lett. 101, 068102 (2008).
- Shen and Wolynes (2004) T. Shen and P. G. Wolynes, Proc. Natl. Acad. Sci. USA 101, 8547 (2004).
- Palacci et al. (2010) J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet, Phys. Rev. Lett. 105, 088304 (2010).
- Paxton et al. (2004) W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. St. Angelo, Y. Cao, T. E. Mallouk, P. E. Lammert, and V. H. Crespi, J. Am. Chem. Soc. 126, 13424 (2004).
- Hong et al. (2007) Y. Hong, N. M. K. Blackman, N. D. Kopp, A. Sen, and D. Velegol, Phys. Rev. Lett. 99, 178103 (2007).
- Jiang et al. (2010) H.-R. Jiang, N. Yoshinaga, and M. Sano, Phys. Rev. Lett. 105, 268302 (2010).
- Volpe et al. (2011) G. Volpe, I. Buttinoni, D. Vogt, H.-J. Kummerer, and C. Bechinger, Soft Matter 7, (2011).
- Narayan et al. (2007) V. Narayan, S. Ramaswamy, and N. Menon, Science 317, 105 (2007).
- Kudrolli et al. (2008) A. Kudrolli, G. Lumay, D. Volfson, and L. S. Tsimring, Phys. Rev. Lett. 100, 058001 (2008).
- Deseigne et al. (2010) J. Deseigne, O. Dauchot, and H. Chaté, Phys. Rev. Lett. 105, 098001 (2010).
- (17) See supplemental material at [url will be inserted by publisher] for movies and additional figures.
- Baskaran and Marchetti (2008a) A. Baskaran and M. C. Marchetti, Phys. Rev. E 77, 011920 (2008a).
- Baskaran and Marchetti (2008b) A. Baskaran and M. C. Marchetti, Phys. Rev. Lett. 101, 268101 (2008b).
- Peruani et al. (2011) F. Peruani, T. Klauss, A. Deutsch, and A. Voss-Boehme, Phys. Rev. Lett. 106, 128101 (2011).
- Yang et al. (2010) Y. Yang, V. Marceau, and G. Gompper, Phys. Rev. E 82, 031904 (2010).
- McCandlish et al. (2012) S. R. McCandlish, A. Baskaran, and M. F. Hagan, Soft Matter 8, (2012).
- Fily and Marchetti (2012) Y. Fily and M. C. Marchetti, Phys. Rev. Lett. 108, 235702 (2012).
- Theurkauff et al. (2012) I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet, (2012), arXiv:1202.6264 .
- Mishra and Ramaswamy (2006) S. Mishra and S. Ramaswamy, Phys. Rev. Lett. 97, 090602 (2006).
- Cates et al. (2010) M. E. Cates, D. Marenduzzo, I. Pagonabarraga, and J. Tailleur, Proc. Natl. Acad. Sci. USA 107, 11715 (2010).
- Peng et al. (2010) Y. Peng, Z. Wang, A. M. Alsayed, A. G. Yodh, and Y. Han, Phys. Rev. Lett. 104, 205703 (2010).
- Han et al. (2008) Y. Han, N. Y. Ha, A. M. Alsayed, and A. G. Yodh, Phys. Rev. E 77, 041406 (2008).
- Weeks et al. (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).
- Brańka and Heyes (1999) A. C. Brańka and D. M. Heyes, Phys. Rev. E 60, 2381 (1999).
- Chaté et al. (2008) H. Chaté, F. Ginelli, G. Grégoire, F. Peruani, and F. Raynaud, Euro. Phys. J. B 64, 451 (2008).
- Schwarz-Linek et al. (2012) J. Schwarz-Linek, C. Valeriani, a. Cacciuto, M. E. Cates, D. Marenduzzo, a. N. Morozov, and W. C. K. Poon, Proc. Natl. Acad. Sci. USA 109, 4052 (2012).
- Tailleur and Cates (2008) J. Tailleur and M. Cates, Phys. Rev. Lett. 100, 3 (2008).
- Cates and Tailleur (2012) M. E. Cates and J. Tailleur, (2012), arXiv:1206.1805 .
- Bialké et al. (2012) J. Bialké, T. Speck, and H. Löwen, Phys. Rev. Lett. 108, 168301 (2012).
- Nelson (2002) D. R. Nelson, Defects and Geometry in Condensed Matter Physics (Cambridge University Press, 2002).
- Blumenfeld and Edwards (2009) R. Blumenfeld and S. F. Edwards, The Journal of Physical Chemistry B 113, 3981 (2009).
- van Blaaderen et al. (1992) A. van Blaaderen, J. Peetermans, G. Maret, and J. K. G. Dhont, J. Chem. Phys. 96, 4591 (1992).
- Klafter et al. (1996) J. Klafter, M. F. Shlesinger, and G. Zumofen, Physics Today 49, 33 (1996).
- Pfister and Scher (1978) G. Pfister and H. Scher, Advances in Physics 27, 747 (1978).
- Dey et al. (2012) S. Dey, D. Das, and R. Rajesh, Phys. Rev. Lett. 108, 238001 (2012).