Phase behaviour of active Brownian particles: the role of dimensionality

Abstract

Recently, there has been much interest in activity-induced phase separations in concentrated suspensions of “active Brownian particles” (ABPs), self-propelled spherical particles whose direction of motion relaxes through thermal rotational diffusion. To date, almost all these studies have been restricted to 2 dimensions. In this work we study activity-induced phase separation in 3D and compare the results with previous and new 2D simulations. To this end, we performed state-of-the-art Brownian dynamics simulations of up to 40 million ABPs – such very large system sizes are unavoidable to evade finite size effects in 3D. Our results confirm the picture established for 2D systems in which an activity-induced phase separation occurs, with strong analogies to equilibrium gas–liquid spinodal decomposition, in spite of the purely non-equilibrium nature of the driving force behind the phase separation. However, we also find important differences between the 2D and 3D cases. Firstly, the shape and position of the phase boundaries is markedly different for the two cases. Secondly, for the 3D coarsening kinetics we find that the domain size grows in time according to the classical diffusive t^1/3 law, in contrast to the nonstandard subdiffusive exponent observed in 2D.

---

1 Introduction

Active materials, whose constituents possess the ability to convert chemical energy into work, have recently been shown to exhibit a range of exotic behaviours compared to what is observed in passive systems.^1–3 These include rectification of bacterial motion,^4–7 giant density fluctuations,^8–11 dynamic phase transitions,^12–17 and the ability to power microscopic motors.^5,18 From the perspective of statistical thermodynamics, all these phenomena are made possible by the fact that active matter systems are intrinsically far from equilibrium even at steady state, and thus are not constrained by the strict rules imposed on systems that obey detailed balance.

The classic example of active particles is that of “run-and-tumble” bacteria, such as Escherichia coli, whose ballistic motion (“runs”) is regularly punctuated by random reorientations of the swimming direction (“tumbles”), occurring with a frequency α. On time-scales much larger than α⁻¹, these dynamics lead to a persistent random walk with a characteristic step length v₀/α, where v₀ is the swim speed of a single bacterium. Because of this, the long-time behaviour of a run-and-tumble bacterium is similar to that of a diffusing Brownian particle, with an effective diffusion coefficient D given by

where d is the spatial dimensionality of the system; for realistic parameter values, D is hundreds of times larger than the Brownian diffusivity of a similarly sized colloidal particle.³

A second important class of active particles is that of active Brownian particles (ABPs), whose swimming direction continuously relaxes through thermal rotational diffusion rather than through discrete tumbling events. The prime experimental example of such ABPs is that of “catalytic swimmers”, half-coated spherical colloids rendered motile through a surface reaction.^19,20 As shown by Cates and Tailleur,²¹ the long-time dynamics of ABPs are equivalent to those of run-and-tumble particles upon the substitution α ↔ (d − 1)D_r, where D_r is the (thermal) rotational diffusion constant.‡

Due to the fact that active particles exhibit diffusive behaviour at long length- and timescales, mass transport in such suspensions is governed by laws similar to those operating in equilibrium Brownian systems. In particular, in a suspension of non-interacting self-propelled particles with spatially uniform swim speeds and tumble rates, transport is governed by the standard form of Fick's first law for ideal systems, relating the mass flux J to the (effective) diffusivity D and concentration gradient ∇ρ:

J = −D∇ρ. (2)

For a spatially uniform diffusivity D, eqn (2) ensures relaxation towards a state of homogeneous density (∇ρ = 0). However, as shown by Schnitzer,²² if the effective diffusivity of eqn (1) varies in space through a non-uniform swim speed (e.g. due to an inhomogeneous fuel concentration), the situation becomes markedly different. This position-dependent swimming velocity v(r) leads to the following generalization of eqn (2):

Tailleur and Cates^21,23 considered an analogue of eqn (3) for the case when particles instead slow down due to collisions and steric interactions, leading to a swim-speed that depends on the local density of bacteria, i.e. v(r) = v(ρ(r)). This assumption leads to the following form of the flux:

For swim speeds that decrease steeply enough with density, the term in square brackets in eqn (4) can become negative, enabling the emergence of non-uniform steady states. Microscopically, this corresponds to an activity-induced phase separation induced by a kinetic feedback mechanism, whereby a local (positive) density fluctuation will lead to a local slowdown of particles, which in turn causes a further accumulation of particles, eventually leading to the nucleation of a dense phase. The existence of this motility-induced phase transition has recently been confirmed for purely repulsive ABPs in both experiments¹⁷ and simulations.^24,25 Furthermore, eqn (4), together with a carefully tuned noise term and a specific form of v(ρ), was recently shown to yield predictions in quantitative agreement with the structural and kinetic properties of a phase-separating ABP fluid.²⁶

The vast majority of previous work on activity-induced phase separation (with a single, very recent, exception²⁷) has focused on phenomena in 2 dimensions. However, real suspensions of active particles are often 3-dimensional. Our goal here is therefore to investigate in detail phase separation in concentrated 3D suspensions of ABPs. From a theoretical standpoint, a phase separating ABP suspension can be viewed as a binary fluid, where the two components form the dilute and dense phase. A key difference in the topology of binary mixtures in 2D and in 3D is that in 2D fluid bicontinuity is only possible by fine tuning to a single composition (50 : 50 if the fluid is otherwise symmetric), meaning that the generic situation in 2D is that of disconnected droplets of the minority phase inside a matrix of the majority phase.^28,29 In 3D, however, both fluids can remain continuously connected through the whole sample for a wide range of compositions. In general, the interplay of such dimensionality-dependent topological differences with nonequilibrium physics is subtle and difficult to predict a priori;^30,31 the specific case of active systems is similarly difficult to gauge. For instance, self-propelled particles might move quite differently in disconnected dilute pockets, as would be expected in a concentrated 2D suspension, compared to in a percolating dilute phase as may be expected in 3D at a similar density. Another important aspect is that of fluctuations: on general grounds, the role of noise is expected to decrease as the dimensionality increases, and previous work in 2D reported important effects of fluctuations on the phenomenology of phase separation.²⁶

In the present study, we perform massively parallel Brownian dynamics simulations of up to ∼4 × 10⁷ repulsive ABPs using ∼8000 CPUs. Such large systems are required to avoid finite size effects which would otherwise compromise our conclusions. We then compare our results to those obtained from large-scale 2D ABP simulations, similar to those performed previously.^24–26 Our simulation results are furthermore compared to results obtained from numerically solving a recently developed continuum model,²⁶ which we first extend from 2D to 3D. The results obtained from both particle and continuum simulations confirm that activity-induced phase separations have remarkable analogies to classical liquid–gas transitions, despite their completely different physical origins. Furthermore, we find these analogies to be even more marked in 3D than in 2D, as indicated for instance by the kinetics of domain growth, which in 3D follows the classical t^1/3 growth law³² from relatively early times. This is, we hypothesise, because fluctuations are less important in this 3D scenario. We also find other important differences between the physics of the 2D and 3D systems: for instance, the shape of the phase diagram and the region within which phase separation is observed are significantly different. As the phase separation occurs, the morphologies of the growing domains are also qualitatively different, with the 3D case featuring smoother interfaces and the near-absence of disconnected droplet phases.

2 Discrete ABP model

Our microscopic ABP model consists of spherical particles interacting through a repulsive, pairwise additive Weeks–Chandler–Andersen potential, given bywith an upper cut-off at r = 2^1/6σ, beyond which U = 0. Here σ denotes the particle diameter, ε determines the interaction strength, and r is the center-to-center separation between two particles. The model was studied by solving the fully overdamped translational and rotational Langevin equations:where F_i is the total conservative force acting on particle i, F_p is the (constant) magnitude of the self-propulsion force, whose direction is defined by p_i [where p_i = (cos θ_i, sin θ_i) in 2D], D_t and D_r = 3D_t/σ² denote the translational and rotational diffusivities, β = (k_BT)⁻¹ is the inverse thermal energy, and Λ_r, Λ_θ, and Λ_p are unit-variance random variables of appropriate dimensionality, whose Cartesian components Λ_i satisfy 〈Λ_i(r, t)Λ_j(r′, t′)〉 = δ_ijδ(r − r′)δ(t − t′). All simulations were carried out using the LAMMPS³³ molecular dynamics software package, with system sizes of up to 1000σ (N ≈ 7 × 10⁵) and 350σ (N ≈ 4 × 10⁷) in 2D and 3D, respectively. Further simulation details are presented in Section A.1.

3 Continuum model

We begin by reviewing the main aspects of the continuum theories presented in ref. 23 and 26 generalising to spatial dimensionalities d > 2.

We begin by noting that eqn (4) can be reexpressed in an equilibrium-like form as follows:

where [scr F, script letter F]₀ = ∫f₀dr is an effective free energy, anddefines the corresponding free-energy density, composed of an ideal entropy-like term and a v-dependent term that resembles an enthalpic attraction. Equations (9) and (10) clearly show how a density-dependent swim speed can induce phase separation into high- and low-density phases, and suggest that the form of this transition may show analogies with equilibrium phase transitions, in spite of the strongly non-equilibrium nature of the self-trapping mechanism leading to the instability.

As previously derived for 1D run-and-tumble particles²³ and generalized to ABPs in higher dimensions,²¹ the coarse-grained density field ρ of particles with a density-dependent swim speed obeys:

Here, D(ρ) is an effective one-body diffusivity, μ an effective chemical potential, and Λ is a random vector as discussed following eqn (8). The effective bulk chemical potential μ₀ is given by the derivative of eqn (10):

To enable a full analysis of phase-separation kinetics in the same spirit as the classical Cahn–Hilliard equation,^32,34eqn (12) needs to be complemented with an interfacial energy-like term that stabilizes domain walls between the phases. Following ref. 26 we accomplish this by assuming that a single ABP samples the local density over a length scale proportional to the density-dependent persistence length [small script l](ρ) of ABP trajectories, given by

where τ_r = D_r⁻¹ is the rotational relaxation time. As we previously showed in ref. 26, this leads to an additional term in the effective chemical potential:

μ = μ₀ − κ(ρ)∇²ρ, (14)

wherewith γ₀ a dimensionless order-unity free parameter. Unlike a traditional density-independent interfacial energy, the second term of eqn (14) cannot be written as the derivative of a free-energy functional; thus, the effective chemical potential in eqn (14) violates detailed balance at second order in a gradient expansion. However, detailed balance can be restored through the modificationwhich again renders μ integrable; this allows comparison between models that differ only in whether detailed balance applies or not.

In order to emulate the physics of excluded volume interactions between ABPs, and thus to prevent the emergence of a phase of infinite density, a repulsive term has to be added to the effective bulk free energy, i.e., f₀ → f₀ + f_rep. As previously,²⁶ we choose the following quartic form of the repulsive contribution f_rep:

f_rep = k_repΘ(ρ − ρ_t)(ρ − ρ_t)⁴. (17)

In eqn (17), Θ denotes the Heaviside step function, k_rep determines the strength of the repulsion, and ρ_t is the threshold density at which the repulsion is switched on. In practice, both these quantities are treated as free parameters, but they were previously found not to significantly affect the phase separation dynamics.²⁶

The crucial component in the continuum model is the form of the density-dependent swim speed v(ρ). To determine this, we make the assumption that single ABP trajectories consist of straight runs with speed v₀ = v(0) interrupted by collision events of duration τ_c, each leading to a complete stall of the particle motion. At low and intermediate densities this assumption leads to the following form for v(ρ):²⁶

where σ_s is a scattering cross-section and τ_MF is the mean free time between two collisions. Furthermore, the effective single-particle diffusivity D(ρ) follows as [cf.eqn (1)]where D₀ = D(0) = v₀²τ_r[d(d − 1)]⁻¹. The same functional form (i.e., a linear decrease of v with density) has recently been observed²⁴ and theoretically predicted^35,36 elsewhere. As described in Appendix B, v(ρ) and D(ρ) can be accurately measured in discrete ABP simulations performed in the one-phase region of the phase diagram and fitted to the functions v(ϕ) = v₀(1 − aϕ) and D(ϕ) = D₀(1 − bϕ)², where ϕ ∈ [0,1] is the particle packing fraction and a ≃ b are fitting parameters.§

Non-dimensionalizing eqs (11)–(12), (14)–(15) and (17) in terms of the length and time units D₀/v₀ ≡ λ and D₀/v₀² = τ_r/d(d − 1) (equivalent to fixing D₀ = v₀ = 1), respectively, finally gives

μ₀ = ln[ϕ(1 − aϕ)] (21)

μ_rep = 4k_repΘ(ϕ − ϕ_t)(ϕ − ϕ_t)³ (22)

μ = μ₀ + μ_rep − κ₀(1 − aϕ)∇²ϕ, (23)

where κ₀ = a[(d − 1)⁻¹γ₀τ_r]². Since the order-unity factor γ₀ is unknown, κ₀ will be treated as a free parameter. Furthermore, N₀ = λ^d/V_p is the number of particles in a cube of side length λ at nominal packing fraction unity, V_p being the volume of a single particle. Note that a mapping between the continuum model and discrete ABP simulations emerges automatically, including the strength of the noise, with the choice of units made here. The only free parameters in the model are thus κ₀, which controls the strength of the interfacial-like term, and k_rep and ϕ_t, that control the strength and onset of the repulsive free energy, respectively.

4 Results and discussion

4.1 Phase diagram

Fig. 1 shows 2D and 3D phase diagrams as a function of the average particle packing fraction ϕ₀ and the Péclet number, Pe, defined bywhere σ is the particle diameter. Clearly, Pe plays a role similar to that of the inverse temperature in equilibrium phase diagrams, in accordance with what has been observed previously in 2D.^25,26 The computational cost of accurately predicting phase diagrams for phase separating fluids is very large; hence, the system sizes used for determining these are, particularly in 3D, smaller than what is needed to guarantee the absence of finite-size effects. Thus, Fig. 1 should be viewed as providing only approximate phase maps. Moreover, as will be further discussed in the next subsection, our continuum theory is not capable of accurately reproducing the binodals and/or spinodals of this system; thus, no quantitative comparison between ABP simulations and continuum theory will be attempted regarding the phase diagrams (as opposed to the phase separation kinetics).

Fig. 1 Phase diagrams in the Pe–ϕ₀ plane as determined from ABP simulations of systems with size L_box = 150σ (2D, left panel) and L_box = 30σ (3D, right panel). Open symbols denote state points with a homogeneous density, and filled symbols denote points with phase-separated steady states. Solid lines show the approximate spinodals, while the dashed lines denote the approximate high-Pe asymptotes for the binodal densities, determined for Pe = 100 (2D) and Pe = 300 (3D). The locations of the spinodals were determined both by visual inspection and by examining the growth in the characteristic length-scale L(t), as further described in Section 4.3. The lower binodal densities (0.25 and 0.24 in 2 and 3 dimensions, respectively) were determined using larger systems (L_box = 1000σ and L_box = 100σ in 2D and 3D, respectively) by starting from fully phase-separated states, while the upper boundaries were determined by approximate extrapolation to the case of vanishing volume of the dilute phase.

The 2D phase diagram broadly resembles that determined previously²⁵ for an essentially identical ABP model, while the 3D phase diagram is in qualitative agreement with that determined recently²⁷ using a Yukawa-type pair potential to model the ABPs. Perhaps most strikingly, it is clear that significantly larger Péclet numbers are needed for phase separation to occur in 3D compared to in 2D. In experiments, where the value of Pe is limited and difficult to vary (for a bacterial strain it is fixed by the swim speed: for E. coli, Pe would be at most ∼100), and where concentrated suspensions are not easy to achieve, this fact may well provide practical obstacles to the observation of phase separation in 3D. Qualitatively, the difference in the critical Pe for phase separation between 2D and 3D may be explained by the fact that orientational correlations decay faster in 3D than in 2D for a given value of τ_r; thus, a higher Péclet number is needed to accomplish a given collision time τ_c, which is the parameter that determines the critical Pe for phase separation.²⁵

We also observe that there are both binodal and spinodal lines in both phase diagrams in Fig. 1, i.e., there are regions between these lines in which the fluid will not spontaneously phase separate from a homogeneous suspension, but will remain phase separated when initialized from a two-phase configuration. This metastable region is furthermore slightly larger in 3D than in 2D. Note, however, that the calculated “binodals” are actually asymptotes representing the high-Pe limit, and were obtained from simulations with Pe = 100 and Pe = 300 in 2D and 3D, respectively, by investigating the stability of the phase-separated state when starting from an initial configuration consisting of a single close-packed domain. Thus, we do not assess the behaviour of this boundary at low Pe, including whether the binodals and spinodals will merge at the critical point as is the case in equilibrium phase diagrams. Nevertheless, the existence of both binodal and spinodal lines in the phase diagrams is once again in analogy with equilibrium phase diagrams, in spite of the purely non-equilibrium nature of the phase separation studied here. The relatively large size of the metastable region (compared to the size of the binodal region) is slightly surprising considering the fact that the apparent noise level (or, equivalently, “effective temperature”) is significantly higher in the ABP systems studied here than in the corresponding thermal systems, something which should help overcoming any (effective) activation energies keeping the system in a metastable state.

In Figs. 2 and 3, representative snapshots taken along the density tieline of 2D and 3D phase-separating ABP systems well within the phase separated region are shown. These snapshots were obtained from much larger systems than those used to determine the phase diagrams, and were run for Péclet numbers of 100 and 300 in 2D and 3D, respectively. In 2D (Fig. 2), the sequence of microstructures evolves from a phase of isolated dense droplets in a dilute background at low density, via an almost bicontinuous phase at intermediate packing fractions, to a phase exhibiting dilute droplets in a dense matrix at still higher densities. In 3D (Fig. 3), the geometries of the phase separated structures are more well-defined, or, equivalently, the fluctuations at the interface between the dense and dilute phases are smaller. This effect is likely also connected with the apparently lower noise level in 3D compared to in 2D, which is clearly visible in the movies available as ESI.† Furthermore, the evolving 3D morphologies always exhibit a single domain of each phase, unlike the corresponding 2D snapshots.¶ While the 3D structures seen in Fig. 3 are clearly influenced by the presence of periodic boundaries, they still convincingly suggest the appearance of spherical, cylindrical, flat, and saddle-like interfaces between the two phases. Furthermore, and as can be seen in Fig. 2 for the 2D case, the density of the two phases remains essentially constant along the tieline in both 2D and 3D,|| an observation which is once again in accordance with the phenomenology of equilibrium phase transitions.

Fig. 2 Snapshots taken at t = 1000τ_r obtained from 2D ABP simulations with varying overall area fraction ϕ₀ as indicated. The systems are of size L_box = 1000σ (N ≈ 7 × 10⁵). The simulation with ϕ₀ = 0.25 was started from a fully phase-separated state with a single close-packed circular domain, while the remaining runs were started from equilibrated (homogeneous) suspensions of passive particles. The density field was obtained by numerical coarse-graining on a grid, as described in Section A.1.

Fig. 3 Snapshots taken at t = 120τ_r obtained from 3D ABP simulations with varying average volume fraction ϕ₀ as indicated. The systems are of size L_box = 100σ (N ≈ 10⁶). The interfaces represent isosurfaces drawn at ϕ = 0.4, where the red side faces the dense phase and the blue side faces the dilute phase. The simulations with ϕ₀ = 0.24, 0.3, and 0.7 were started from fully phase-separated states consisting of a single close-packed domain with a spherical interface for ϕ₀ = 0.24 and with flat interfaces for ϕ₀ = 0.3 and 0.7, while the remaining runs were started from equilibrated (homogeneous) suspensions of passive particles.

A subtle point which needs to be highlighted in the present context is the way in which the Péclet number is varied. In most previous studies of ABPs (with one very recent exception³⁶ employing soft particles), Pe was varied by changing the bare swim speed v₀ while keeping the rotational relaxation time τ_r constant. Viewing the problem from an energetic perspective, Pe quantifies the ratio between the “ballistic energy” F_pσ and the thermal energy k_BT. For a system of infinitely hard spheres, these two are the only energy scales present, and thus Pe uniquely determines the balance between active and thermal forces. However, for the slightly softer spheres used here and in most previous ABP studies, another energy scale arises, describing the steepness of the repulsive potential; for the WCA potential used here, this energy scale is quantified through the Lennard-Jones parameter ε. Thus, Pe is no longer sufficient to describe the system, but the ratio k_BT/ε, quantifying how “hard” the particles are, also comes into play. In previous studies on self-propelled WCA particles, k_BT/ε = 1 was used throughout, and F_p (or, equivalently, v₀) was used as the free parameter controlling Pe. However, this implies that the ratio F_pσ/ε is not constant throughout the phase diagram, potentially leading to undesired effects. For example, an increase of F_pσ/ε will lead to a decreasing effective radius of the particles, meaning that colliding particles will exhibit increasingly large overlaps as Pe increases: for two ABPs with Pe = 300 and k_BT/ε = 1 colliding at a 90 degree angle, the center-to-center distance where the active and repulsive forces balance will be ∼0.85σ, i.e., 15 percent smaller than the “thermal” diameter σ. Apart from yielding an unphysical Pe-dependent effective particle size, controlling Pe through this method will also lead to very large repulsive forces between overlapping particles. This effect may also help to explain the recently reported reentrant behaviour where the motility-induced phase transition is suppressed for large enough values of Pe.³⁵ In order to avoid these effects, we chose to keep F_pσ/ε (and thus v₀) constant at a value of 24, leading to a constant “effective diameter” of σ for two particles colliding at a 90 degree angle; Pe was then adjusted by changing k_BT (and thus τ_r). This subtle difference in how the Péclet number is varied may be important, especially in determining the region where phase separation arises in 3D: our preliminary results indicate that, for this particular ABP model, the shape of the phase diagram is dramatically changed when varying Pe by changing F_p, possibly even removing the phase coexistence region altogether. Presumably, this effect is more visible in 3D due to the higher Péclet numbers required for phase separation compared to the 2D case.

4.2 Comparison with the continuum model

Fig. 4 shows a comparison between the results of ABP and continuum simulations at equal time and system size in 2 and 3 dimensions, at average packing fraction ϕ₀ = 0.5. In accordance with our previous observations in 2D,²⁶ the agreement between domain topologies observed using the ABP and continuum models is excellent, especially considering the fact that the length and time units as well as the noise level are completely determined by the mapping detailed above (with κ₀ as the only fitting parameter). In particular, the mapping successfully captures the lower apparent noise level in 3D compared to in 2D (see further movies provided as ESI†).** We therefore conclude that Cahn–Hilliard-type continuum models like this one are successful in describing phase-ordering kinetics even in non-equilibrium systems, and at a tremendously reduced computational cost compared to direct ABP simulations: in 3D, the computational cost for solving the continuum equations for the system seen in Fig. 4 is ∼0.5 CPU hours, compared to the ∼6 × 10⁶ CPU hours needed for the corresponding ABP simulation. However, a striking feature observed in our 2D ABP simulations but not captured by our continuum model is the spontaneous formation of “gas” voids inside the dense domains that are visible in Fig. 2. This effect is clearly a far-from-equilibrium one in that it breaks time reversal symmetry: voids are usually formed in the center of a dense cluster and then diffuse towards the interface with the dilute phase (see further movies available as ESI†). This intriguing effect is neither observed in our continuum simulations, nor in 3D ABP simulations, and we currently do not have a theoretical explanation for it.

Fig. 4 Snapshots obtained from an ABP simulation (left) and by numerically solving eqn (20) (right), both at average packing fraction ϕ₀ = 0.5. Snapshots are taken at equal times of t = 500τ_r (t = 100τ_r) in 2D (3D). The continuum and ABP systems are of approximately equal size L = 60λ (2D, top panel) and L ≈ 21λ (3D, center and bottom panels). The bottom panel shows 2D projections of the 3D density field, obtained from the same configurations as the isosurfaces in the center panel.

Fig. 5 shows the probability distribution P(ϕ) associated with the local particle packing fraction ϕ, at overall packing fraction ϕ₀ = 0.5 for the continuum model and ABP simulations. In 3D as in 2D, the agreement between the general appearances of the curves is good, although the density of the dilute phase is slightly shifted to higher densities in the continuum model compared to what is observed in the ABP systems. As detailed in Fig. 6, this discrepancy can be explained by the fact that the coexisting densities predicted by a common tangent construction change when we change the parameters of the ad hoc repulsive potential of eqn (17).†† Since this potential is phenomenological, and the choice of parameters is furthermore restricted by the numerical difficulties associated with using very steep potentials, the coexistence densities (and thus the binodals and spinodals) cannot be expected to be quantitatively described by our continuum theory. Thus, the theory is not accurate when it comes to describing the location of the binodals or spinodals in the phase diagrams; for an accurate determination of phase boundaries, microscopic theories such as those developed in refs. 25 and 35 are better suited.

Fig. 5 Probability distribution P(ϕ) in 2D (top) and 3D (bottom) of the local particle packing fraction ϕ obtained from ABP simulations (black curves), from the continuum model as written (red curves) and with detailed balance (“DB”) restored as per eqn (16) (blue curves), all at average packing fraction ϕ₀ = 0.5. In 2D, P(ϕ) was sampled over quadratic coarse-graining areas of side length 0.8λ and averaged over the time window 500τ_r ≤ t ≤ 3500τ_r, while in 3D the curves were sampled from coarse-graining boxes of side length ≈0.4λ at t = 500τ_r.

Fig. 6 Bulk free energy f_bulk = f₀ + f_rep as given by eqs (10) and (17), obtained using v = v₀(1 − 1.3ϕ), and the parameter values {k_rep, ϕ_t} = {10³, 0.64} (“Soft repulsion”) and {10⁷, 0.75} (“Hard repulsion”), respectively. The dashed lines denote approximate common-tangent constructions, yielding the coexistence densities {ϕ_g, ϕ_l} = {0.23, 0.70} and {0.18, 0.75} for the “soft” and “hard” curves, respectively. Terms linear in ϕ are irrelevant for the common-tangent construction and have been subtracted for clarity.

We finally observe that the addition of the detailed balance-restoring term of eqn (16) leads to a small shift in the coexistence density of the dilute phase towards higher densities. This shift arises from the non-equilibrium nature of the system dynamics, and can be explained by theoretical arguments detailed in a recent publication.³⁹ Its magnitude, however, turns out to be small for the present parameter set.

4.3 Phase separation kinetics

Fig. 7 shows the time-evolution of the characteristic domain size L(t), obtained from the first moment of the static structure factor as described in Section A.1, for both the continuum and ABP models. For the classical models of phase-separation kinetics, L(t) usually exhibits power-law growth, i.e., L(t) ∼ t^α, where the growth exponent α depends on the transport properties of the system.

Fig. 7 Time-dependent domain length L(t) obtained from the inverse first moment of the structure factor at average packing fraction ϕ₀ = 0.5 in 2D (top) and 3D (bottom). The dashed lines indicate the fitted exponents as follows: α_2D = 0.27(9), α_3D = 0.34(1) (ABPs), α_2D = 0.28(7), α_3D = 0.33(3) (continuum model), α_2D = 0.27(9), α_3D = 0.34(8) (continuum model with detailed balance term). The latter set of curves has been vertically shifted for clarity.

For phase separation in diffusive systems with negligible hydrodynamic interactions,‡‡ one expects α = 1/3,^32,34,42 which is indeed reproduced (Fig. 7) within our numerical accuracy in 3D from both ABP and continuum simulations. In 2D, however, both our models confirm the α ≈ 0.28 observed previously in ABP simulations by Redner et al.²⁵ The fact that the 2D exponent is slightly smaller than that which is expected from traditional scaling arguments can be attributed to the relatively high noise level in the 2D system: in equilibrium systems, it is an established fact that noise will lead to a subdiffusive intermediate scaling regime, a phenomenon which is likely to be transferable to non-equilibrium phase transitions like the one studied here.^28,29,42 We therefore conjecture that the nonstandard 2D exponent found here is merely an intermediate regime that will eventually switch over to a t^1/3 scaling at later times. Nevertheless, it is interesting that the 3D kinetics show no evidence of this intermediate regime.

We also observe that the restoration of detailed balance as per eqn (16) does not seem to have any detectable effect on the phase separation kinetics for the parameters used here. However, since the detailed balance-violating term was shown in Section 4.2 to lead to a shift in coexistence densities, kinetic consequences of this violation must eventually arise, for some values of κ₀ or ϕ₀. We have in fact recently found this to be the case for a continuum model containing a similar detailed balance-breaking term, where a decrease of the growth exponent was observed as the magnitude of the detailed balance violation was increased.³⁹ Finally, it is important to note the clear crossover between superdiffusive behaviour (α > 1/3) at short times to diffusive behaviour (α ≈ 1/3) when L(t) exceeds the “persistence length” [small script l] = v₀τ_r/(d − 1) (dotted lines in Fig. 7). On a microscopic level this is explained by the fact that, as long as L < [small script l], mass transport between domains can take place ballistically, while in the region where L > [small script l], the swimming directions of ABPs are fully randomized over the typical travelling length between domains, and thus standard diffusive behaviour is recovered. In 3D, this also highlights the problem of finite-size effects: in order to measure coarsening kinetics in the diffusive regime, the box length L_box needs to greatly exceed [small script l], which is in turn fixed by the rather high Péclet number needed to obtain a deep quench in 3D (Pe = 300 leads to [small script l] = 50σ). This fact, together with the high volume fractions required for phase separation, motivates the extremely large systems (L_box = 350σ, N ≈ 4 × 10⁷) studied here.

During coarsening, we find that the static structure factor S(k) typically exhibits a peak whose magnitude grows with time while its position gradually shifts towards lower values of k, corresponding to a growth in the typical length scale L(t). In classical phase-separation theory, this behaviour is usually expressed through the so-called dynamical scaling hypothesis^32,42 as

S(k, t) = [L(t)]^df(kL), (25)

where f is a time-independent scaling function and d is the spatial dimensionality of the system. Fig. 8 shows rescaled plots of S(k), sampled at different times during coarsening for both the continuum and the ABP models. The good data collapse shows that dynamical scaling is fulfilled in all cases, as has previously been observed in molecular dynamics^40,41 as well as lattice Boltzmann⁴³ simulations of equilibrium phase separations. Furthermore, the dotted lines in Fig. 8 show the prediction of Porod's law that S(k) ∼ k^−(d+1) for large k, as long as the interface is reasonably flat.§§⁴² In 2D, the correspondence is clearly satisfactory within a reasonable range of k values (10 ≤ kL ≤ 30). In 3D, however, the predicted k⁻⁴ behaviour is poorly reproduced by the ABP simulation, although it can be observed for 10 ≤ kL ≤ 20 in the continuum simulation. Our explanation for this is the comparingly small spatial dimensions (L_box = 350σ) of the 3D ABP simulation box compared to the 2D one (L_box = 1000σ), which presumably does not allow the k⁻⁴ behaviour to develop fully before boundary effects start to take over.

Fig. 8 Static structure factor S(k) obtained at different times as indicated during ABP (left) and continuum (right) simulations in 2D (top) and 3D (bottom). The curves have been rescaled in accordance with the dynamical scaling hypothesis, as described in the text. Dotted lines indicate the large-k behaviour predicted by Porod's law. All results were obtained from the same systems as the results in Fig. 7. The oscillations at high values of k in the continuum plots correspond to very small lengthscales and are discretisation artifacts.

5 Conclusions and outlook

In this study, we have performed very large scale Brownian dynamics simulations of dense suspensions containing up to 40 million active Brownian particles (ABPs), focusing on their phase behaviour and phase separation kinetics. We have also compared our results to those of a continuum model recently developed by us,²⁶ which uses an effective free-energy mapping based on a particle swim-speed that decreases with density due to collisions.

Our results show some remarkable similarities to phenomena seen in thermodynamic (attraction-induced) gas–liquid phase coexistence:

(1) Both the 2D and 3D phase diagrams exhibit binodal and spinodal lines, with the role of (inverse) temperature being played by the Péclet number.

(2) The coexistence densities remain essentially constant along the tieline: only the amount of each phase is affected by changing the average density ϕ₀, in accordance with what is seen in equilibrium phase diagrams.

(3) The sequences of microstructures observed during coarsening closely resemble those observed in classical spinodal decompositions; this is most apparent in our largest 3D simulation (left panel in Fig. 4).

(4) The phase separation kinetics exhibit traditional power-law growth of the characteristic length-scale L(t), with the 3D systems showing the classical t^1/3 behaviour.

(5) Phase separation kinetics can be accurately described by a classical continuum theory, analogous to the Cahn–Hilliard equation but using an effective free energy density that depends on v(ϕ). Although the effective interfacial free energy weakly violates detailed balance, the effect of this violation remains small, and only leads to a moderate shift of the coexistence densities for the parameters used here.

(6) The time-evolution of structure shows dynamical scaling behaviour.

The points listed above are valid in both 2 and 3 dimensions, but there are also significant differences between the two cases:

(1) The “critical Péclet number” needed for phase separation is significantly higher in 3D than in 2D (see Fig. 1). This may have practical implications for future 3D experiments, as it is not easy to experimentally control the Péclet number in active suspensions.

(2) The “metastable” region between the binodal and the spinodal is larger in 3D than in 2D; in particular, such a region does not appear to exist at all at the high-ϕ₀ end of the 2D phase diagram.

(3) There is an apparently lower level of noise in 3D than in 2D. This lower noise level in 3D leads to patterns with more well-defined geometries than in 2D (see Figs. 2 and 3, and movies available as ESI†).

(4) The phase transition in 3D always appears to take place via a single domain of each phase, whereas in 2D the microstructures usually consist of a continuous matrix of the majority phase containing separated domains of the minority phase which eventually coarsen through an Ostwald-like process.

(5) The kinetics of domain growth differs between 2D and 3D, with the 2D case exhibiting a non-standard scaling exponent α ≈ 0.28, while in 3D the standard α ≈ 1/3 exponent is reproduced beyond the ballistic regime.

While our ABP model clearly relates only approximately to experimental systems, such as suspensions of bacteria or catalytic swimmers, the equilibrium-like behaviours listed above indicate surprising analogies between non-equilibrium and equilibrium phase transitions. However, several important issues within the field remain to be resolved: perhaps most crucially, the impact of hydrodynamic interactions on the complex collective behaviour of active suspensions. While the kinetic consequences of far-field hydrodynamic flows on phase-separation kinetics in classical fluids have been extensively studied,^32,42 no such studies exist for active systems (although the role of near-field hydrodynamics in the clustering of self-propelled disks has recently been investigated^44,45). Furthermore, the subtle interplay between thermodynamic interparticle attractions and motility-induced phase transitions has only recently begun to be studied,¹⁶ and is likely to generate many more interesting non-equilibrium phenomena.

A Simulation details

A.1 ABP model

Equations (6)–(8) were solved using a slightly modified version of the LAMMPS³³ molecular dynamics software package. For determining the phase diagram, a cubic system with side length 150σ (2D) or 35σ (3D) was used, both amounting to N ≈ 20 000 particles. Periodic boundary conditions were applied in all directions. For studying the phase separation kinetics, significantly larger systems with side lengths 1000σ (2D) and 350σ (3D) were used; these system sizes amount to N ≈ 7 × 10⁵ and N ≈ 4 × 10⁷ particles in 2 and 3 dimensions, respectively. In terms of Lennard-Jones time units τ_LJ = σ²/(εβD_t), a time step of 5 × 10⁻⁵τ_LJ was used throughout, and each simulation was run for ∼10⁸ timesteps, with the largest 3D system amounting to ∼700 hours of computing time on a 8192-core IBM Blue Gene/Q node. As discussed in Section 4.1, the Péclet number was varied by varying D_t while keeping F_p constant at a value of F_p = 24ε/σ. Unless otherwise stated, simulations were started from an initial configuration of equilibrated passive colloids (F_p = 0) with k_BT = ε, and were quenched by switching on F_p at t = 0. In terms of Lennard-Jones units, the reduced length and time scales (see Section 3 for definitions) in the ABP systems with Pe = 100 (2D) and Pe = 300 (3D) are given by λ_2D = λ_3D = 16.67σ and τ^(2D)_r = 1.389τ_LJ, τ^(3D)_r = 4.167τ_LJ. The length scale mapping furthermore leads to N₀ = 4λ²/(σπ²) ≈ 354 (2D) and N₀ = 6λ³/(πσ³) ≈ 8842 (3D), which was used to set the noise strength when numerically solving eqn (20).

The characteristic domain length L(t) was computed as the inverse of the first moment of the static structure factor S(k, t), i.e.,

where L is the length of the simulation box and the upper cut-off k_cut was taken to be the first minimum in S(k).

In 2D, snapshots from ABP simulations were obtained by coarse-graining the local density on a grid using a weighting function w(r) ∼ exp[−r_cut²/(r_cut² − r²)], where r is the distance of a particle from the a particular lattice point, and r_cut is a cut-off distance which was taken to be slightly larger than the size of a lattice site. In 3D, no such weighting function was employed. The coarse-grained density field was further analyzed by constructing a density isosurface using the ParaView visualization software.

A.2 Continuum model

Equations (20)–(23) were solved numerically employing a standard Euler finite-difference scheme using the parameter values reported in Table 1. For comparison with ABP simulations, lattice sizes of 150² (2D) and 52³ (3D) were employed, whereas lattice sizes of 512² (2D) and 256³ (3D) were used for the calculation of growth exponents. The initial condition was taken to be one of uniform density, with a local random offset of ≈5%.

Table 1 Parameters used to solve eqn (20). ΔL indicates the lattice spacing and Δt the size of the time step; all other parameters are defined in Section 3

	2D	3D
a	1.0	1.3
k rep	2500	1000
ϕ t	0.88	0.64
κ 0	0.3	0.9
ΔL/λ	0.4	0.4
d(d − 1)Δt/τr	1 × 10^−2	3 × 10^−2

---

B Density dependence of swim speed

As was discussed in Supplemental Material of ref. 26, to obtain an estimate for the effective free energy, v(ϕ) should be sampled in the one-phase region of the phase diagrams in Fig. 1, i.e., for sufficiently low Péclet numbers that the system does not phase separate. This might, in principle, lead to problems since the fitting parameters a and b will themselves depend on Pe. However, as we showed in ref. 26, this dependence is small enough to be negligible. We therefore measured v(ϕ) and D(ϕ) right outside the spinodal region (for Pe = 40 and Pe = 100 in 2D and 3D, respectively).

As our operational microscopic definition of v(ϕ), we sampled the average of the instantaneous velocities of all particles projected onto their self propulsion direction p_i:

v = βD_t〈(F_pp_i + F_i)·p_i〉 (27)

In the limit ϕ → 0, where interparticle forces vanish, this definition further leads to v → βD_tF_p = v₀. The density-dependent effective diffusivity D(ϕ) was independently measured by linear fitting to the long-time limit of the mean-square displacement 〈R²〉 = 2(d − 1)D. In Fig. 9, plots of v(ϕ) and D(ϕ) together with the best fits to the functions v = v₀(1 − aϕ) and D = D₀(1 − bϕ)² are shown. It can clearly be seen that the predicted linear decrease of v(ϕ) is accurately reproduced over a wide range of packing fractions. Furthermore, the predicted relation between D and v is satisfactorily fulfilled, with the fitting coefficients a and b approximately equal in both 2D and 3D. Based on the fitted parameter values, we used a = b = 1.0 (2D) and a = b = 1.3 (3D) when solving the continuum model (see further Section A.2).

Fig. 9 Density-dependent swim speed v(ϕ) (black symbols) and diffusivity D(ϕ) (red symbols), and the ratio D/v² (blue curve and symbols) obtained from ABP simulations at Pe = 40 (2D, top panel) and Pe = 100 (3D, bottom panel). The black and red curves show the best fits to the functions v₀(1 − aϕ) and D₀(1 − bϕ)², respectively, with the optimized values a = 1.05, b = 1.04 (2D) and a = 1.33, b = 1.28 (3D). Dashed lines show the predicted zero-density values. Plotted quantities are in Lennard-Jones units σ and τ_LJ, as defined in Section A.1.

Acknowledgements

Helpful discussions with Alan Bray, Aidan Brown, Tom Lion, Ignacio Pagonabarraga, Julien Tailleur, Adriano Tiribocchi, and Raphael Wittkowski are gratefully acknowledged, and JS would like to thank Fred Farrell for assistance with the numerical coarse-graining and Kevin Stratford for providing code for the calculation of structure factors. We thank EPSRC EP/J007404 for funding. JS is supported by the Swedish Research Council (350-2012-274), RJA by a Royal Society University Research Fellowship, and MEC by a Royal Society Research Professorship.

References

1. M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 11431189.
2. P. Romanczuk, M. Bär, W. Ebeling, B. Lindner and L. Schimansky-Geier, Eur. Phys. J.: Spec. Top., 2012, 202, 1–162.
3. M. E. Cates, Rep. Prog. Phys., 2012, 75, 042601.
4. P. Galajda, J. Keymer, P. Chaikin and R. Austin, J. Bacteriol., 2007, 189, 8704–8707.
5. L. Angelani, R. Di Leonardo and G. Ruocco, Phys. Rev. Lett., 2009, 102, 048104.
6. G. Lambert, D. Liao and R. H. Austin, Phys. Rev. Lett., 2010, 104, 168102.
7. A. Pototsky, A. M. Hahn and H. Stark, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 042124.
8. V. Narayan, S. Ramaswamy and N. Menon, Science, 2007, 317, 105–108.
9. J. Deseigne, O. Dauchot and H. Chaté, Phys. Rev. Lett., 2010, 105, 098001.
10. H. P. Zhang, A. Be'er, E.-L. Florin and H. L. Swinney, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 13626–13630.
11. H. H. Wensink and H. Löwen, J. Phys.: Condens. Matter, 2012, 24, 464130.
12. M. E. Cates, D. Marenduzzo, I. Pagonabarraga and J. Tailleur, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 11715–11720.
13. I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2012, 108, 268303.
14. S. R. McCandlish, A. Baskaran and M. F. Hagan, Soft Matter, 2012, 8, 2527.
15. F. D. C. Farrell, M. C. Marchetti, D. Marenduzzo and J. Tailleur, Phys. Rev. Lett., 2012, 108, 248101.
16. G. S. Redner, A. Baskaran and M. F. Hagan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 012305.
17. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
18. R. Di Leonardo, L. Angelani, D. Dell'Arciprete, G. Ruocco, V. Iebba, S. Schippa, M. P. Conte, F. Mecarini, F. D. Angelis and E. D. Fabrizio, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 9541–9545.
19. J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh and R. Golestanian, Phys. Rev. Lett., 2007, 99, 048102.
20. S. J. Ebbens and J. R. Howse, Soft Matter, 2010, 6, 726–738.
21. M. E. Cates and J. Tailleur, EPL, 2013, 101, 20010.
22. M. J. Schnitzer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1993, 48, 2553–2568.
23. J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2008, 100, 218103.
24. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 235702.
25. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 055701.
26. J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2013, 111, 145702.
27. A. Wysocki, R. G. Winkler and G. Gompper, arXiv:1308.6423, 2013.
28. A. J. Wagner and J. M. Yeomans, Phys. Rev. Lett., 1998, 80, 1429–1432.
29. A. J. Wagner and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1999, 59, 4366–4373.
30. P. Stansell, K. Stratford, J. C. Desplat, R. Adhikari and M. E. Cates, Phys. Rev. Lett., 2006, 96, 085701.
31. K. Stratford, J. C. Desplat, P. Stansell and M. E. Cates, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 030501.
32. P. Chaikin and T. C. Lubensky, Principles of condensed matter physics, Cambridge University Press, 1995.
33. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
34. Kinetics of Phase Transitions, ed. S. Puri and V. Wadhawan, CRC Press, 2009.
35. J. Bialké, H. Löwen and T. Speck, EPL, 2013, 103, 30008.
36. Y. Fily, S. Henkes and M. C. Marchetti, Soft Matter, 2013 10.1039/c3sm52469h.
37. J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936–940.
38. M. Abkenar, K. Marx, T. Auth and G. Gompper, Phys. Rev. E, 2013, 88, 062314.
39. R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo and M. E. Cates, arXiv:1311.1256, 2013.
40. R. Yamamoto and K. Nakanishi, Mol. Simul., 1996, 16, 119–126.
41. S. Majumder and S. K. Das, EPL, 2011, 95, 46002.
42. A. J. Bray, Adv. Phys., 2002, 51, 481–587.
43. V. M. Kendon, M. E. Cates, I. Pagonabarraga, J.-C. Desplat and P. Bladon, J. Fluid Mech., 2001, 440, 147–203.
44. A. Zöttl and H. Stark, arXiv:1309.4352, 2013.
45. S. M. Fielding, arXiv:1210.5464, 2012.

---

Footnotes

† Electronic supplementary information (ESI) available: Movies showing the time evolution of the density field obtained from Brownian dynamics simulations of ABPs and solutions of the continuum model in 2D and 3D. See DOI: 10.1039/c3sm52813h

‡ Note that ABPs, unlike run-and-tumble particles, cannot be realized in 1D.

§ From the kinetic argument above, a and b are predicted to be identical, which is also what we found in our previous 2D study where a = 1.05 and b = 1.04.²⁶ In the rest of the article the two parameters will therefore be assumed identical and both denoted by a.

¶ However, our results indicate that the 2D systems will also eventually coarsen until a single domain of each phase remains. Thus, we do not see any sign of the finite cluster phases observed in experiments on active colloids.^13,37

|| Note, however, that the density ϕ_g of the dilute phase in similar systems has been predicted and observed^25,38 to change with Pe as ϕ_g ∼ Pe⁻¹.

** A more detailed analysis of the mapping between ABPs and the continuum model shows that the lower apparent noise level in 3D compared to in 2D is simply a consequence of the fact that the number of particles N₀ present in a coarse-graining volume λ^d is an increasing function of d. In other words, the “granularity” of the system becomes smaller in higher dimensions, leading to smaller fluctuations. Thus, this effect is not specific to non-equilibrium systems like the one studied here.

†† Note, however, that the common-tangent construction is only expected to yield the observed coexistence densities if the detailed balance-restoring term of eqn (16) is also added to the gradient term of the chemical potential (blue curves in Fig. 5).³⁹

‡‡ Note that our simulations are distinctly different from classical molecular dynamics studies of gas–liquid coexistence in Lennard-Jones fluids (e.g.refs. 40 and 41), where momentum is conserved. While momentum conservation will lead to several hydrodynamic regimes with exponents α > 1/3, our overdamped Langevin dynamics will suppress these super-diffusive transport mechanisms.

§§ Strictly speaking, the structure factor should scale as S(k) ∼ k^{−(2d−[scr D, script letter D])} for large k, where [scr D, script letter D] is the so-called fractal dimensionality of the interface. For a flat interface, however, [scr D, script letter D] = d − 1, where d is the spatial dimensionality of the system.

---