Thermodynamic phases in twodimensional active matter
Abstract
Active matter has been intensely studied for its wealth of intriguing properties such as collective motion[1], motilityinduced 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 twodimensional (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 selfpropelled particles with inversepowerlaw 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 twostep melting scenario is maintained. At high activity, a critical point opens up a gas–liquid MIPS region. We expect that the independent appearance of twostep melting and of MIPS is generic for a large class of twodimensional 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 ErlangenNürnberg, Staudtstr. 7, 91058 Erlangen, Germany

Department of Physics, Graduate School of Science, The University of Tokyo, 731 Hongo, Bunkyo, Tokyo, Japan

MaxPlanckInstitut für Physik komplexer Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany
The kinetic discretestep 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 longtime dynamics remains diffusive. Detailed balance is satisfied only for vanishing activities, that is, in the limit . Cooperative effects are introduced through a repulsive inversepowerlaw 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 manyparticle system is without alignment interactions.
We observe at all densities and activities a nonequilibrium steady state in which spatial correlation functions are welldefined. 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 powerlaw decay of orientational correlations yet exponential decay of positional correlations in the hexatic phase, and longrange orientational correlations together with positional powerlaw 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 harddisklike hexatics[9] (that correspond to a interaction with ).
In 2D equilibrium systems, a true crystal with Bragg peaks and longrange 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 longrange 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 shortrange positional correlations (no order) and algebraic orientational correlations (quasilongrange 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 quasilongranged (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. 2bd). This soft hexatic maintains orientational quasiorder with extremely shortranged 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 Ushaped 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 colorcoded 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 steadystate 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 singlepeaked 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 abovementioned Ushape. 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 finitesize 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 stripeshaped 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 phaseseparated 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 socalled swim pressure due to active motion[16, 17, 18].
We thus find that twostep melting and MIPS are kept separated by the homogeneous liquid phase, which is both the highdensity end of MIPS and the lowdensity 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 coarsegrained 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 twostep melting phase transitions depends on the softness of the particles, which can be tuned via the exponent of the inversepowerlaw pair potential . The twostep 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 activematter 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 equilibriumlike 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 infinitetime 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 inversepowerlaw interactions provide a tunable set of activematter models to study phase transitions and phase coexistence.
Methods
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 softsphere 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 squareshaped 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 detailedbalance 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 steadystate distribution in the infinitetime 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
(1) 
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 ,
(2) 
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 steadystate 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
bondorientational 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 directiondependent pair correlation function . Before
averaging this twodimensional 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 ensembleaveraged over 100 configurations of particles, recorded after a warmup 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.
References
 1. T. Vicsek, A. Czirók, E. BenJacob, I. Cohen, O. Shochet, Novel Type of Phase Transition in a System of SelfDriven Particles, Phys. Rev. Lett 75 6 (1995).
 2. M. E. Cates, J. Tailleur, MotilityInduced Phase Separation, Annu. Rev. Condens. Matter Phys. 6, 219244 (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, TwoStage 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, TwoDimensional Melting of Colloidal Hard Spheres, Phys. Rev. Lett. 118, 158001 (2017).
 7. B. I. Halperin and D. R. Nelson, Theory of TwoDimensional Melting, Phys. Rev. Lett. 41, 121 (1978).
 8. S. Z. Lin, B. Zheng, and S. Trimper, Computer simulations of twodimensional melting with dipoledipole interactions, Phys. Rev. E 73, 066106 (2006).
 9. E. P. Bernard and W. Krauth, TwoStep Melting in Two Dimensions: FirstOrder LiquidHexatic Transition, Phys. Rev. Lett 107, 155704 (2011).
 10. S. C. Kapfer and W. Krauth, TwoDimensional Melting: From LiquidHexatic 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 selfpropelled hard disks, Phys. Rev. E 89, 062301 (2014).
 12. N. D. Mermin and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One or TwoDimensional Isotropic Heisenberg Models, Phys. Rev. Lett. 17, 1133 (1966).
 13. L. F. Cugliandolo, P. Digregorio, G. Gonnella, and A. Suma, Phase Coexistence in TwoDimensional 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, TwoDimensional Systems, J. Chem. Phys. 42, 4268 (1965).
 15. M. Schrader, P. Virnau, and K. Binder, Simulation of vaporliquid 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, Dislocationmediated 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 PhaseSeparating Active Colloidal Fluid, Phys. Rev. Lett 110, 055701 (2013).
 23. Y. Fily, M. C. Marchetti, Athermal Phase Separation of SelfPropelled Particles with No Alignment, Phys. Rev. Lett 108, 235702 (2012).
 24. T. Speck, J. Bialké, A. M. Menzel, H. Löwen, Effective CahnHilliard 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 SelfPropelled Particles, Phys. Rev. Lett 108, 168301 (2012).
 26. W. Mickel, S. C. Kapfer, G. E. SchröderTurk, 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.

juliane_klamser@yahoo.de, sebastian.kapfer@fau.de, werner.krauth@ens.fr.