Thermodynamic phases in two-dimensional active matter

Thermodynamic phases in two-dimensional active matter

Juliane U. Klamser, Sebastian C. Kapfer & Werner Krauth

Active matter has been intensely studied for its wealth of intriguing properties such as collective motion[1], motility-induced phase separation (MIPS)[2], and giant fluctuations away from criticality[3]. However, the precise connection of active materials with their equilibrium counterparts has remained unclear. For two-dimensional (2D) systems, this is also because the experimental[4, 5, 6] and theoretical[7, 8, 9, 10] understanding of the liquid, hexatic, and solid equilibrium phases and their phase transitions is very recent. Here, we use self-propelled particles with inverse-power-law repulsions (but without alignment interactions) as a minimal model for 2D active materials. A kinetic Monte Carlo (MC) algorithm allows us to map out the complete quantitative phase diagram. We demonstrate that the active system preserves all equilibrium phases, and that phase transitions are shifted to higher densities as a function of activity. The two-step melting scenario is maintained. At high activity, a critical point opens up a gas–liquid MIPS region. We expect that the independent appearance of two-step melting and of MIPS is generic for a large class of two-dimensional active systems.

  • Laboratoire de Physique Statistique, Département de physique de l’ENS, Ecole Normale Supérieure, PSL Research University, Université Paris Diderot, Sorbonne Paris Cité, Sorbonne Universités, UPMC Univ. Paris 06, CNRS, 75005 Paris, France

  • Theoretische Physik 1, FAU Erlangen-Nürnberg, Staudtstr. 7, 91058 Erlangen, Germany

  • Department of Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo, Japan

  • Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany

The kinetic discrete-step MC dynamics for active matter progresses through individual particle displacements which are accepted or rejected with the standard Metropolis criterion. The proposed displacements of a single particle are correlated in time, leading to ballistic local motion characterised by a persistence length that measures the degree of activity (see Fig. 1). Correlations decay exponentially, so that the long-time dynamics remains diffusive. Detailed balance is satisfied only for vanishing activities, that is, in the limit . Cooperative effects are introduced through a repulsive inverse-power-law pair potential (a potential is used). The Metropolis rejections are without incidence on the sequence of proposed moves, which are also uncorrelated between particles so that the active many-particle system is without alignment interactions.

We observe at all densities and activities a nonequilibrium steady state in which spatial correlation functions are well-defined. At vanishing activity , we recover the equilibrium phase diagram (see Fig. 2). In particular we observe the exponential decay of positional and orientational correlation functions in the liquid, the power-law decay of orientational correlations yet exponential decay of positional correlations in the hexatic phase, and long-range orientational correlations together with positional power-law decay in the solid. The system is particularly amenable to simulation because of its “soft” hexatic phase, characterised by small positional correlation lengths[10]. Soft hexatics have much shorter MC correlation times and reach the thermodynamic limit for smaller system sizes than the hard-disk-like hexatics[9] (that correspond to a interaction with ).

In 2D equilibrium systems, a true crystal with Bragg peaks and long-range positional order exists only in the limit, as a consequence of the Mermin–Wagner theorem[12]. We find that the solid phase of the active system also exhibits algebraic positional order, just as in equilibrium. Decreasing the density or, remarkably, increasing the activity weakens positional correlations and eventually melts the solid (see Fig. 2a). Close to the melting transition, the algebraic decay of the positional correlations can be clearly observed in our simulation data (points A and B in Fig. 2b). Together with the long-range orientational order, this explicitly identifies the solid phase (Fig. 2c).

The hexatic phase in equilibrium is characterised through a lower degree of order than the solid, namely through short-range positional correlations (no order) and algebraic orientational correlations (quasi-long-range order). We find precisely such a phase in the active system, in a narrow strip of densities below and activities above the solid phase (see Fig. 2a). Starting from the solid, we indeed observe positional correlations that change qualitatively upon a minute increase in activity while leaving the orientational correlations almost unchanged, leading to hexatic order (from point B to point C in Fig. 2b). Positional correlations in the hexatic decay exponentially beyond the correlation length but the orientational correlations remain quasi-long-ranged (points C through E in Fig. 2c and d). On moving towards the liquid at any point within the hexatic phase, the positional correlation lengths decrease, resulting in a strikingly soft hexatic close to the liquid–hexatic transition (point E in Fig. 2b-d). This soft hexatic maintains orientational quasi-order with extremely short-ranged positional correlations, even at densities for which the equilibrium system is already deep inside the solid phase. Increasing the activity thus shifts the equilibrium phase boundaries towards higher densities. The stability of the partially ordered hexatic phase is remarkable especially as it takes place for persistence lengths significantly larger than the interparticle distance . Our massive computations give no indications of a direct transition from the solid into the liquid state, even at the highest accessed densities.

MIPS[2] has been frequently reported in 2D active systems but agreement on its interpretation was not reached. Recent work in an active dumbbell system proposes that MIPS continuously extends from the equilibrium liquid–hexatic transition region and that one of the separated phases preserves some degree of order[13]. This is not the case in our system: We observe MIPS as a U-shaped region of liquid–gas coexistence (see Fig. 3a) with an onset at high activity and at relatively low density. Both competing phases in the phase separated state feature exponential decay of orientational and positional correlation functions (see the color-coded configuration snapshots in Fig. 3a). MIPS is thus clearly separated from melting (see Fig. 2a). In a MC simulation, a coarsening process generally precedes macroscopic phase separation in the time evolution towards the steady state. In the active system, the transient coarsening leading up to MIPS can be expected to be overcome at earlier times than for hard disks. This makes it easier to distinguish it from the formation of a steady-state gel[11], although we do not expect the nature of the coexisting phases participating in MIPS to depend on the softness of the potential.

The analogy with the liquid–gas coexistence in equilibrium simple fluids suggests the interpretation of the onset of MIPS as a critical point. Indeed, below the onset, the system remains homogeneous at large length scales with a single-peaked local density distribution (see point H in Fig. 3b). Above the onset of MIPS, the local densities develop a bimodal distribution. The peak positions separate further as increases, quantifying the above-mentioned U-shape. Moreover, at a given , the peak local densities in the coexistence region are independent of the global density (see Fig. 3c). This is further substantiated by a finite-size scaling analysis at constant (see Fig. 4). The phenomenology thus agrees with that of an equilibrium phase coexistence where the relative proportions of the liquid and gas adapt to the global density, but where the degree of order of each of the phases and their densities remain unchanged. Inside the coexistence region, at small densities, we observe an approximately circular bubble of liquid inside the gas, followed by a stripe-shaped form that winds around the periodic simulation box, and then followed by a bubble of gas inside the liquid. In finite equilibrium systems, this complex behaviour is brought about by the interface free energy[14, 15] which vanishes at the critical point. A detailed analysis of the homogeneous phases, but also of the phase-separated systems, reveals the origin of the phase separation in the kinetic MC model. In the bulk of the coexisting liquid and gas phases, the directions of motion of the individual particles are uncorrelated beyond a very small length scale (see Fig. 4b). At the liquid–gas interface, however, a majority of the increment vectors point inwards towards the liquid phase. Even though a theoretical framework as reliable as statistical mechanics is currently lacking, the effective cohesion in nonequilibrium is often attributed to the so-called swim pressure due to active motion[16, 17, 18].

We thus find that two-step melting and MIPS are kept separated by the homogeneous liquid phase, which is both the high-density end of MIPS and the low-density end of the order–disorder transitions. Intriguingly, the subtle hexatic phase survives at considerable activities. Our massive computations allow us to reach the steady state even for strong activities and for high densities, but only theory will be able to ascertain the stability of the separating liquid phase and of the hexatic state in the limit. Another open question is how the KTHNY theory of 2D melting[7, 19, 20], built for equilibrium systems, can be extended to active systems. At the phase transitions, we find the exponents and predicted by the KTHNY theory for the positional and orientational correlations respectively (e. g. see Fig. 2b and c). This suggests the existence of a coarse-grained functional that plays the role of an equilibrium free energy. More specifically, increased activity appears to reduce a state variable in our system that corresponds to the pressure at equilibrium. This interpretation reconciles both the fact that increased activity induces melting through reduction of the effective pressure, but also that liquid–gas phase coexistence is possible only at high activities, i. e., at low effective pressure, in striking analogy to the behavior of simple fluids in equilibrium. Further work is required to test this hypothesis.

In equilibrium, the nature of the two-step melting phase transitions depends on the softness of the particles, which can be tuned via the exponent of the inverse-power-law pair potential . The two-step melting scenario for very soft particles () comprises two continuous transitions, whereas for harder particles with the liquid–hexatic transition becomes first order[10]. One may speculate that the activity plays a similar role as the hardness of the particles and that at high activities the liquid–hexatic transition becomes first order.

The liquid–gas coexistence phase we observe in kinetic MC is an example of MIPS seen earlier in active-matter models using Brownian and molecular dynamics simulations[21, 22, 23, 24, 25]. As we show, MIPS can be reproduced within kinetic MC dynamics, without added equilibrium-like mixing steps[11]. Indeed, the direction of the individual persistent particle motion suffices to produce the effect: Particles may be kinetically arrested for persistence lengths larger than the mean free path, leading to density inhomogeneities, where particles in dense regions are walled in by particles coming from less dense regions. However, this does not unhinge the coarsening mechanism: The size of the inhomogeneous regions increases with time, and the infinite-time steady state in a finite system is characterized by only two coexisting regions so that, for the studied potential, a gel phase is clearly absent. MIPS appears naturally in kinetic MC, and we suggests that it is a generic feature of active matter in 2D and in higher dimensions. Our inverse-power-law interactions provide a tunable set of active-matter models to study phase transitions and phase coexistence.


0.1 Kinetic MC Dynamics for active particles.

We use a modified Metropolis algorithm that breaks detailed balance. In each MC step, a displacement by an amount is proposed for a single randomly chosen particle . If the displacement were implemented, it would cause a change in the total energy of the system . The move is accepted with the Metropolis probability . The total energy is given by , where the power–law potential does not introduce an energy scale separate from density. Thus, in equilibrium, the only relevant scale is the dimensionless density , with the effective temperature in the Metropolis probability scaling as . In nonequilibrium, activity introduces a new control parameter, expressed in terms of the persistence length . Simulations are performed in a simulation box with periodic boundary conditions and the soft-sphere potential is truncated as follows: .

We introduce activity into the dynamics by choosing the proposed displacement of particle based on the previously proposed displacement of the same particle. The correlation is introduced in two stages. First, a random vector is sampled from a bivariate normal distribution , where is the standard deviation. In the second stage, is generated from using the folding scheme

with and , with , i. e., , where is the floor function. The folding scheme limits the size of the proposed displacements and keeps the dynamics local. The folding scheme is equivalent to a random walk of the displacement variables in a box with reflecting boundary conditions, see Fig. 1b. Note that the random walk of the displacements is independent of the positions of the particles, as the new increment persists whether the displacement was accepted or not.

The square-shaped displacement box introduces a small degree of anisotropy into the dynamics for . However, we explicitly verify that the resulting steady state is unaffected with respect to the properties concerning this letter. At small densities in the gaseous state, where is much smaller than the mean free path, the kinetic MC dynamics effectively reverts to detailed-balance dynamics as interactions between particles are rare. At higher densities, all large proposed displacements have vanishing Metropolis probabilities, thus leading to effectively isotropic dynamics. In our numerical observation anisotropic effects are undetectable within other sources of noise.

0.2 Probability distribution of increments and persistence length.

In continuous time, the increment variable evolves according to a diffusion equation , with vanishing probability flux through the boundary of the displacement box. The steady-state distribution in the infinite-time limit is the uniform distribution . By a Fourier ansatz, one readily finds that the autocorrelation time of increments is dominated by the first harmonic of the displacement box, an that for large , the autocorrelation function decays as


The position of a free particle evolves as . For times shorter than the autocorrelation time, , its mean displacement is essentially given by the increment at time ,


where the drift velocity is given by the initial condition of the increment, and the subleading term contains contributions by diffusion of , including reflections, in the displacement box. Averaging over all initial conditions with the steady-state uniform distribution, we obtain the mean drift velocity

Considering Eqs. (1) and (2), we may identify the persistence length . The persistence length offers a length scale which separates ballistic from diffusive motion. The persistence length defined in this way allows to collapse data for the increment autocorrelations at widely different activities (see Fig. 1a). The prefactors of the superexponential terms are obtained from a numerical fit.

0.3 Measurements.

The orientational order is quantified by the correlation function
of the bond-orientational order parameter calculated with Voronoi weights[26]. is a measure of the correlation of the local sixfold orientational order at distance . Positional order is studied with the direction-dependent pair correlation function . Before averaging this two-dimensional histogram over configurations, each configuration is realigned[9] such that the axis points in the direction of the global orientation parameter . Correlation functions in Fig. 2 are ensemble-averaged over 100 configurations of particles, recorded after a warm-up time of MC sweeps , with each sweep containing MC steps. Configurations were recorded in time intervals of MC sweeps.

MIPS is quantified by histograms of local densities. We compute local densities by covering the system with randomly placed test circles of radius . The local dimensionless density of each test circle is , where are the number of particle centres located within a circle of area . The detailed analysis (see Figs. 3, 4) of MIPS uses . The larger shifts the liquid–gas coexistence phase boundaries, in particular the critical point, to smaller densities and persistence lengths, which drastically shortens the time to reach the steady state. Configuration snapshots in Fig. 3a were taken after MC sweeps. Data in Fig. 3b consists of ensemble averages over 100 configurations recorded in time intervals of sweeps. The (, ) data in Figs. 4 was recorded after (, ) sweeps. The average consists of 500 configurations recorded in time intervals of (, ) sweeps.


  • 1. T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, O. Shochet, Novel Type of Phase Transition in a System of Self-Driven Particles, Phys. Rev. Lett 75 6 (1995).
  • 2. M. E. Cates, J. Tailleur, Motility-Induced Phase Separation, Annu. Rev. Condens. Matter Phys. 6, 219-244 (2015).
  • 3. S. Dey, D. Das, R. Rajesh, Spatial Structures and Giant Number Fluctuations in Models of Active Matter, Phys. Rev. Lett 108, 238001 (2012).
  • 4. K. Zahn, R. Lenke, and G. Maret, Two-Stage Melting of Paramagnetic Colloidal Crystals in Two Dimensions, Phys. Rev. Lett. 82, 2721 (1999).
  • 5. U. Gasser, C. Eisenmann, G. Maret, and P. Keim, Melting of Crystals in Two Dimensions, Chem. Phys. Chem. 11, 963 (2010).
  • 6. A. L. Thorneywork, J. L. Abbott, D. G. A. L. Aarts, and R. P. A. Dullens, Two-Dimensional Melting of Colloidal Hard Spheres, Phys. Rev. Lett. 118, 158001 (2017).
  • 7. B. I. Halperin and D. R. Nelson, Theory of Two-Dimensional Melting, Phys. Rev. Lett. 41, 121 (1978).
  • 8. S. Z. Lin, B. Zheng, and S. Trimper, Computer simulations of two-dimensional melting with dipole-dipole interactions, Phys. Rev. E 73, 066106 (2006).
  • 9. E. P. Bernard and W. Krauth, Two-Step Melting in Two Dimensions: First-Order Liquid-Hexatic Transition, Phys. Rev. Lett 107, 155704 (2011).
  • 10. S. C. Kapfer and W. Krauth, Two-Dimensional Melting: From Liquid-Hexatic Coexistence to Continuous Transitions, Phys. Rev. Lett 114, 035702 (2015).
  • 11. D. Levis and L. Berthier, Clustering and heterogeneous dynamics in a kinetic Monte Carlo model of self-propelled hard disks, Phys. Rev. E 89, 062301 (2014).
  • 12. N. D. Mermin and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models, Phys. Rev. Lett. 17, 1133 (1966).
  • 13. L. F. Cugliandolo, P. Digregorio, G. Gonnella, and A. Suma, Phase Coexistence in Two-Dimensional Passive and Active Dumbbell Systems, Phys. Rev. Lett. 119, 268002 (2017).
  • 14. J. E. Mayer, and Wm. W. Wood, Interfacial Tension Effects in Finite, Periodic, Two-Dimensional Systems, J. Chem. Phys. 42, 4268 (1965).
  • 15. M. Schrader, P. Virnau, and K. Binder, Simulation of vapor-liquid coexistence in finite volumes: A method to compute the surface free energy of droplets, Phys. Rev. E 79, 061104 (2009).
  • 16. S. C. Takatori, W. Yan, and J. F. Brady, Swim Pressure: Stress Generation in Active Matter, Phys. Rev. Lett 113, 028103 (2014).
  • 17. A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. E. Cates, and J. Tailleur, Pressure and Phase Equilibria in Interacting Active Brownian Spheres, Phys. Rev. Lett 114, 198301 (2015).
  • 18. A. P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar and J. Tailleur, Pressure is not a state function for generic active fluids, Nat. Phys. 11, 673 (2015).
  • 19. D. R. Nelson and B. I. Halperin, Dislocation-mediated melting in two dimensions, Phys. Rev. B 19, 2457 (1979).
  • 20. A. P. Young, Melting and the vector Coulomb gas in two dimensions, Phys. Rev. B 19, 1855 (1979).
  • 21. M. E. Cates, Diffusive transport without detailed balance in motile bacteria: does microbiology need statistical physics?, Rep. Prog. Phys. 75 042601 (2012).
  • 22. G. S. Redner, M. F. Hagan, A. Baskaran, Structure and Dynamics of a Phase-Separating Active Colloidal Fluid, Phys. Rev. Lett 110, 055701 (2013).
  • 23. Y. Fily, M. C. Marchetti, Athermal Phase Separation of Self-Propelled Particles with No Alignment, Phys. Rev. Lett 108, 235702 (2012).
  • 24. T. Speck, J. Bialké, A. M. Menzel, H. Löwen, Effective Cahn-Hilliard Equation for the Phase Separation of Active Brownian Particles, Phys. Rev. Lett 112, 218304 (2014).
  • 25. J. Bialké, T. Speck, H. Löwen, Crystallization in a Dense Suspension of Self-Propelled Particles, Phys. Rev. Lett 108, 168301 (2012).
  • 26. W. Mickel, S. C. Kapfer, G. E. Schröder-Turk, and K. Mecke, Shortcomings of the bond orientational order parameters for the analysis of disordered particulate matter, J. Chem. Phys. 138, 044501 (2013).
  • We thank Hughes Chaté and Hartmut Löwen for helpful discussions. W.K. acknowledges support from the Alexander von Humboldt Foundation.


Figure 1: Kinetic MC algorithm. (a): Autocorrelation of proposed displacements as a function of the covered distance for a single 2D active particle. The data collapse for widely different activities allows the definition of a persistence length (see Methods) (See inset for raw data without rescaling). (b): Time evolution of proposed displacements in a box of size . The displacement is sampled from a bivariate normal distribution with standard deviation (here ), centred at the previous displacement . Positions sampled outside the box (black points) are folded back, implementing the reflecting boundary condition (see Methods). (c): Trajectory for a single particle , with the corresponding history of displacements from (b). The color gradient changing with time allows to connect (b) and (c), e. g. the last displacements in (b) are in the third quadrant thus the particle in (c) moves to the lower left etc. (c), (d) and (e): Sampled trajectories, illustrating the transition from a passive random walk ( in (e)) to a persistent random walk ( in (c)).
Figure 2: Complete phase diagram and two-step melting. Depicted results are for particles and , with the particle diameter (see Methods). (a): Activity vs. density phase diagram. MIPS between a liquid and a gas, at high , is situated far above the solid–hexatic-melting lines. The red dot indicates a possible critical point. Inset in (a): Two-step melting for small with shift of transition densities to higher values with increasing , preserving the equilibrium phases. Two-step melting from the solid is induced by density reduction (as in equilibrium) but also by an increase in . (b), (c), and (d): Activity-induced two-step melting high above the equilibrium melting densities (, A: ; B: ; C: ; D: ; E: ; F: ; G: ). (b): Positional correlation function along the axis, in units of the interparticle distance . (c): Orientational correlation function along the axis. (d): Snapshots of configurations, particles colour-coded with their local orientation parameter . A and B are quasi-long-ranged in and long-ranged in , thus corresponding to the solid phase. C, D, and E are short-ranged in and quasi-long-ranged in , thus corresponding to the hexatic phase. F and G decay exponentially in both correlation functions and thus correspond to a liquid.
Figure 3: Characterisation of MIPS. Data for , . (a): Snapshots of configurations close to the onset of liquid–gas coexistence. (Particles represented in arbitrary size and colour-coded, as in Fig. 2d, according to their local orientation parameter.) U-shaped phase boundary is apparent. Orientational order is short-ranged in both liquid and gas phase. At constant activity , the liquid volume fraction grows with increasing density until the liquid entirely fills the system. The location of the critical point depends on and it appears at a smaller density and than in Fig. 2. (b): Histograms of local densities (see Methods for definition) for a variation of activities at constant (A: ; B: ; C: ; D: ; E: ; F:; G: ; H: ; For better presentation, histograms are shifted along the -axis with increasing ). Transition from a single-peaked to a double-peaked distribution and increasing separation of peaks with increasing . (c): Densities of liquid and gas (identified through peak position as in (b)) in an activity vs. local density diagram. This demonstrates independence of phase densities on global density for fixed and validates the phase-separation picture.
Figure 4: Finite-size scaling of local density histograms Data with and . (a) and (b): Local density histograms (see Methods for definitions) for global densities and (identical -axis used). Local density peaks sharpen with increasing system size , and are located at the same value of , demonstrating that in the MIPS region gas and liquid densities are independent of the global density. (c) and (d): Snapshots of configurations at global densities corresponding to (a), where the liquid is the minority phase, and (b), where the liquid is the majority phase by volume fraction (cf. height of the peaks in (a) and (b)). Inset in (c): Direction of motion of the individual particles indicated by arrows, illustrating the origin of MIPS. The directions of motion are uncorrelated inside the homogeneous liquid and gas. Only particles at the interface move towards the interior of the liquid patch, enclosing particles of the high density region.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description