Thermodynamic phases in two-dimensional active matter

Abstract

Active matter has been much studied for its intriguing properties such as collective motion, motility-induced phase separation and giant fluctuations. However, it has remained unclear how the states of active materials connect with the equilibrium phases. For two-dimensional systems, this is also because the understanding of the liquid, hexatic, and solid equilibrium phases and their phase transitions is recent. Here we show that two-dimensional self-propelled point particles with inverse-power-law repulsions moving with a kinetic Monte Carlo algorithm without alignment interactions preserve all equilibrium phases up to very large activities. Furthermore, at high activity within the liquid phase, a critical point opens up a gas-liquid motility-induced phase separation region. In our model, two-step melting and motility-induced phase separation are thus independent phenomena. We discuss the reasons for these findings to be common to a wide class of two-dimensional active systems.

Fig. 1

Kinetic Monte Carlo 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 є(t) is sampled from a bivariate normal distribution with standard deviation σ (here $\sigma \ll \delta$), centred at the previous displacement є(t − 1). Positions sampled outside the box (black points) are folded back, implementing the reflecting boundary condition (see Methods). c Trajectory for a single particle $\mathbf{r}\left(t\right)={\sum }_{k=1}^t$ є(k), with the corresponding history of displacements from (b). The colour gradient changing with time allows to connect (b, 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–e Sampled trajectories, illustrating the transition from a passive random walk ($\sigma \gg \delta$ in e) to a persistent random walk ($\sigma \ll \delta$ in c)

Fig. 2

Full phase diagram and two-step melting. Depicted results are for N ~ 4.4 × 10⁴ particles and δ = 0.1. a Activity λ vs. density ϕ phase diagram. MIPS between a liquid and a gas, at high λ, is situated far above the solid–hexatic–liquid melting lines. The red dot indicates a possible critical point. b 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 λ. (See Methods for definition of error bars.) c–e Activity-induced two-step melting high above the equilibrium melting densities (ϕ = 2.4, A: λ = 0.987; B: λ = 0.999; C: λ = 1.018; D: λ = 1.033; E: λ = 1.045; F: λ = 1.064; G: λ = 1.079). c Positional correlation function g(x,y) along the x axis, in units of the global mean interparticle distance d = (πN/V)^−1/2. d Orientational correlation function g₆(r). e Snapshots of configurations, particles colour-coded with their local orientation parameter ψ₆. A and B are quasi-long-ranged in g(x, 0) and long-ranged in g₆(r), thus corresponding to the solid phase. C, D and E are short-ranged in g(x, 0) and quasi-long-ranged in g₆(r), thus corresponding to the hexatic phase. F and G decay exponentially in both correlation functions and thus correspond to a liquid. (For definitions see Methods)

Fig. 3

Characterisation of motility-induced phase separation. Data for N ≈ 1.1 × 10⁴, δ = 0.7. a Snapshots of configurations close to the onset of liquid–gas coexistence. (Particles represented in arbitrary size and colour-coded, as in Fig. 2e, according to their local orientation parameter.) U-shaped phase boundary is apparent. Orientational order is short-ranged in both liquid and gas phase, thus the coloured snapshots appear grey. 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 ϕ = 0.4 (A: λ = 450; B: λ = 359; C: λ = 268; D: λ = 214; E: λ = 161; F:λ = 107; G: λ = 80; H: λ = 54; For better presentation, histograms are shifted along the y 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 ϕ_loc diagram. This demonstrates independence of phase densities on global density for fixed λ and validates the phase-separation picture. For a clear presentation, we plot exemplifying error bars for ϕ = 0.4 only, which were obtained from the width of the corresponding peak in b

Fig. 4

Finite-size scaling of local density histograms. Data with λ = 214 and δ = 0.7. a, c Local density histograms (see Methods for definitions) for global densities ϕ = 0.2 and ϕ = 0.6 (identical x axis used). Local density peaks sharpen with increasing system size N, and are located at the same value of ϕ_loc, demonstrating that in the MIPS region gas and liquid densities are independent of the global density. b, d Snapshots of configurations at global densities corresponding to a, where the liquid is the minority phase, and c, where the liquid is the majority phase by volume fraction (cf. height of the peaks in a and c). Inset in b 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

Fig. 5

Convergence from two different initial conditions as a function of run-time t. Data at ϕ = 2.4 with λ = 1.0333 and δ = 0.1. a Positional correlation function g(x,y) along the x axis, in units of the global mean interparticle distance d. b Orientational correlation function g₆(r). Both panels show the evolution with run-time t (see colour code for run-time ordering in a) starting from a random initial configuration (dashed lines) and from an initial configuration with particles arranged in a perfect hexagonal lattice (solid lines). Both initial configurations approach the same hexatic steady state (black solid line)