Thermodynamic phases in two-dimensional active matter
Juliane U. Klamser Sebastian C. Kapfer Werner Krauth


Space-Efficient Bimachine Construction Based on the Equalizer Accumulation Principle
    Stefan Gerdjikov1,2 Stoyan Mihov2 Klaus U. Schulz3
    
====================================================================================




 
  * 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



Active matter has been intensely studied for its wealth of intriguing
properties such as collective motion<cit.>, motility-induced
phase separation (MIPS)<cit.>, and giant fluctuations away from
criticality<cit.>. 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<cit.>
and theoretical
<cit.>
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.


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 Algorithm).  
Correlations decay exponentially, so that
the long-time dynamics remains diffusive. Detailed balance is satisfied only 
for vanishing activities, that is, in the
limit λ→ 0.  Cooperative effects are introduced through a repulsive 
inverse-power-law pair
potential (a 1/r^6 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 PhaseDiagram). 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 1/r^6 system is particularly amenable to simulation because
of its “soft” hexatic phase, characterised by small positional correlation
lengths<cit.>.  Soft hexatics have much shorter MC
correlation times and reach the thermodynamic limit for smaller system
sizes than the hard-disk-like hexatics<cit.> (that correspond 
to a 1 / r^n interaction with n →∞).   

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<cit.>.  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 PhaseDiagrama).  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 PhaseDiagramb). Together with
the long-range orientational order, this explicitly identifies the solid phase
(PhaseDiagramc).

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 PhaseDiagrama).  
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 PhaseDiagramb). 
Positional correlations in the hexatic decay exponentially beyond the 
correlation length but the orientational correlations remain quasi-long-ranged 
(points C through E in PhaseDiagramc 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 PhaseDiagramb-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 d = (π N /V)^-1/2.
Our massive computations give no indications of a direct transition from the 
solid into the liquid state, even at the highest accessed densities.

MIPS<cit.> 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<cit.>. This is not the case in our
system: We observe MIPS as a U-shaped region of liquid–gas coexistence
(see MIPSa) 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
MIPSa). MIPS is thus clearly separated from melting
(see PhaseDiagrama). In a MC simulation, a
coarsening process generally precedes macroscopic phase separation in the
time evolution towards the steady state.   In the active 1/r^6 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<cit.>, 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 MIPSb).
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 MIPSc).  This is further substantiated by a finite-size scaling analysis
at constant λ (see FiniteSizeScal). 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 (N,V,T) equilibrium systems, this complex behaviour
is brought about by the interface free energy<cit.>
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 FiniteSizeScalb). 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<cit.>.

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
<cit.>, built for equilibrium
systems, can be extended to active systems. At the phase transitions, we
find the exponents ≃ 0.33 and ≃ 0.25 predicted by the KTHNY theory 
for the positional and orientational correlations respectively (e. g. see 
PhaseDiagramb 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 n of the inverse-power-law pair 
potential U(r) ∝ 1/r^n. The two-step melting scenario for 
very soft particles (n ≲ 6) comprises two continuous transitions, 
whereas for harder particles with n ≳ 6 the liquid–hexatic transition 
becomes first order<cit.>.

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
<cit.>. 
As we show, MIPS can be reproduced
within  kinetic MC dynamics, without added equilibrium-like mixing
steps<cit.>.  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
1/r^6 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.




 §.§ 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 ϵ_i is 
proposed for a single randomly chosen particle i. If the displacement were
implemented, it would cause a change in the total energy of the system E→ 
E'. The move is accepted with the Metropolis probability P = 
min[1,e^-(E'-E)]. The total energy is given by E=∑_i<j 
U(r⃗_i-r⃗_j), where the power–law potential U(r)=(γ/r)^6 does 
not introduce an energy scale separate from density.
Thus, in equilibrium, the only relevant scale is the 
dimensionless density ϕ = γ^2 N/V, with the effective 
temperature in the Metropolis probability scaling as ϕ^3. 
In nonequilibrium, activity introduces a new control 
parameter, expressed in terms of the persistence length λ.
Simulations are performed in a (7:4√(3)) simulation box with
periodic boundary conditions and the soft-sphere potential is truncated
as follows: Ũ(r) = U(min(r,1.8γ)).

We introduce activity into the dynamics by choosing the proposed displacement 
ϵ'_i of particle i based on the previously proposed 
displacement ϵ_i of the same particle.
The correlation is introduced in two stages. First, a random vector 
η is sampled from a bivariate normal distribution ∝exp[-( 
η - ϵ_i)^2/2σ^2], where σ is the 
standard deviation. In the second stage, ϵ'_i is generated from 
η using the folding scheme
        
    ϵ'_i,z→
                    q_i,z  - δ   if  q_i,z < 2δ
    
                    3δ - q_i,z   if  q_i,z≥ 2δ ,

with z ∈{x,y} and q_i,z =(η_z + δ)
mod 4δ, with
a mod b = a - b ⌊ab⌋, i. e.,
0≤ q_i,z<4δ, where ⌊ a ⌋ 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 ϵ⃗_i in
a [-δ,δ]^2 box with reflecting boundary conditions, see
Algorithmb.  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 λ > 0.  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.



 §.§ Probability distribution of increments and persistence length.

In continuous time, the increment variable ϵ⃗ evolves according
to a diffusion equation ∂_t P(ϵ⃗,t)=σ^22Δ_ϵ⃗ P(ϵ⃗, t), with vanishing probability
flux through the boundary of the displacement box.
The steady-state distribution in the infinite-time limit is the uniform distribution
P=(2δ)^-2.
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 t, the autocorrelation function decays as

    C(t) ≡⟨ϵ⃗(t_0+t)·ϵ⃗(t_0)⟩∝e^-t/τ ,  where τ = 8δ^2/π^2σ^2 .

The position of a free particle evolves as r⃗(t) ≡r⃗(t_0) + ∫_t_0^t ds ϵ⃗(s).
For times shorter than the autocorrelation time, t≲τ,
its mean displacement is essentially given by the increment at time t_0,

    ⟨ |r⃗(t)-r⃗(t_0)| ⟩ = v t + 𝒪(√(t)),

where the drift velocity v≡|ϵ⃗(t_0)| 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 ϵ⃗(t_0)
with the steady-state uniform distribution, we obtain the mean drift velocity

    v
    =
        δ/3[√(2)+sinh ^-1(1)]
         .

Considering Eqs. (<ref>) and (<ref>), we may identify the
persistence length λ≡τ v≃ 0.62 δ^3 σ^-2.
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
C(⟨ r⟩ = λ x) ∝e^-(x + c_1 x^2 + c_2 x^3)
at widely different activities (see Algorithma).
The prefactors of the superexponential terms c_1, c_2 are obtained from a numerical fit.



 §.§ Measurements.

The orientational order is quantified by the correlation
function g_6(r) ∝⟨∑_i,j^N
ψ_6(r_i)ψ_6(r_j)δ(r-r_ij)⟩ of the 
bond-orientational order parameter ψ_6(r_i) calculated with Voronoi
weights<cit.>. g_6(r) is a measure of the correlation of the local
sixfold orientational order at distance r.  Positional order is studied
with the direction-dependent pair correlation function g(x,y). Before
averaging this two-dimensional histogram over configurations, each
configuration is realigned<cit.> such that the Δ x axis points in
the direction of the global orientation parameter Ψ_6 = ∑_i^N
ψ_6(r_i). Correlation functions in
PhaseDiagram are ensemble-averaged over 100 configurations of ∼ 4.4
× 10^4 particles, recorded after a warm-up time of 5 × 10^6 MC sweeps
, with each sweep containing N MC steps. Configurations were recorded
in time intervals of 2.7 × 10^4 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 7.5γ.  The local dimensionless density of each test circle is
ϕ = γ^2 N  / V, where N are the number
of particle centres located within a circle of area V. The
detailed analysis (see Figs. <ref>, <ref>) of MIPS uses
δ = 0.7γ. 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 MIPSa were taken after 1.1 × 10^7 
MC sweeps. Data in MIPSb
consists of ensemble averages over 100 configurations recorded in 
time intervals of 3× 10^4 sweeps.
The N ∼ 4.4× 10^4 (N ∼ 1.1× 10^4, N ∼ 2.7 × 10^3) 
data in Figs. <ref> was recorded after 5.5×10^6 
(2.2×10^7, 1.4×10^7) sweeps. The average consists of 500 
configurations recorded in time intervals of 2.2 × 10^3 (8.7 × 10^3, 
1.7 × 10^4) sweeps. 








§ REFERENCES


——

Vicsek 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).

MIPS M. E. Cates, J. Tailleur, Motility-Induced Phase Separation, 
Annu. Rev. Condens. Matter Phys. 6, 219-244 (2015).

GiantNumberFluctuations S. Dey, D. Das, R. Rajesh, Spatial 
Structures and Giant Number Fluctuations in Models of Active Matter, 
Phys. Rev. Lett 108, 238001 (2012).

















ZahnMaret1999 K. Zahn, R. Lenke, and G. Maret,
Two-Stage Melting of Paramagnetic Colloidal Crystals in Two Dimensions,
Phys. Rev. Lett. 82, 2721 (1999).

KeimMaret2010 U. Gasser, C. Eisenmann, G. Maret, and P. Keim,
Melting of Crystals in Two Dimensions,
Chem. Phys. Chem. 11, 963 (2010).

ThorneyworkDullensA. 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). 













HalperinNelson1978 B. I. Halperin and D. R. Nelson,
Theory of Two-Dimensional Melting,
Phys. Rev. Lett. 41, 121 (1978).












LinTrimper2006 S. Z. Lin, B. Zheng, and S. Trimper,
Computer simulations of two-dimensional melting with dipole-dipole interactions,
Phys. Rev. E 73, 066106 (2006).












BernardKrauth2011 E. P. Bernard and W. Krauth, Two-Step Melting in Two 
Dimensions: First-Order Liquid-Hexatic Transition, 
Phys. Rev. Lett 107, 155704 (2011).

KapferKrauth2015 S. C. Kapfer and W. Krauth, Two-Dimensional 
Melting: From Liquid-Hexatic Coexistence to Continuous Transitions, 
Phys. Rev. Lett 114, 035702 (2015).

BerthierHardSpheres 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).

MerminWagner 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).

Dumbbell 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).

MayerWood J. E. Mayer, and Wm. W. Wood,  Interfacial Tension 
Effects 
in Finite, Periodic, Two-Dimensional Systems, J. Chem. Phys. 
42, 4268 (1965).

SchraderBinder 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).

















Brady2014 S. C. Takatori, W. Yan, and J. F. Brady,
Swim Pressure: Stress Generation in Active Matter,
Phys. Rev. Lett 113, 028103 (2014).

SolonCatesTailleur2015 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).

SolonCatesTailleur2015_Nat 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).

NelsonHalperin1979 D. R. Nelson and B. I. Halperin,
Dislocation-mediated melting in two dimensions,
Phys. Rev. B 19, 2457 (1979).

Young1979 A. P. Young,
Melting and the vector Coulomb gas in two dimensions,
Phys. Rev. B 19, 1855 (1979).


CatesReviw2012 M. E. Cates,
Diffusive transport without detailed balance in motile bacteria: does 
microbiology need statistical physics?, Rep. Prog. Phys. 75 
042601 (2012).

KineticModel G. S. Redner, M. F. Hagan, A. Baskaran, Structure 
and Dynamics of a Phase-Separating Active Colloidal Fluid, Phys. Rev. 
Lett 110, 055701 (2013).

FlyMarchetti2012 Y. Fily, M. C. Marchetti, Athermal Phase 
Separation of Self-Propelled Particles with No Alignment, Phys. Rev. 
Lett 108, 235702 (2012).

CahnHilliard 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).

LoewenCrystalization J. Bialké, T. Speck, H. Löwen, 
Crystallization in a Dense Suspension of Self-Propelled Particles, Phys. 
Rev. Lett 108, 168301 (2012).

VoronoiWeights 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.





  Email addresses: julianeklamser@yahoo.de, sebastian.kapfer@fau.de, 
werner.krauth@ens.fr.









