A kinetic-Monte Carlo perspective on active matter

TL;DR

This paper uses a kinetic Monte Carlo approach to map out the phase diagram of active matter, revealing motility-induced phase separation and two-step melting with an intermediate hexatic phase in two-dimensional self-propelled particles.

Contribution

It introduces a kinetic Monte Carlo method to study active matter, providing a detailed phase diagram including MIPS and two-step melting phenomena.

Findings

Identification of MIPS in various potential forms.

Observation of two-step melting with an intermediate hexatic phase.

Phase diagram mapped for different potential stiffnesses, including the hard-disk limit.

Abstract

We study non-equilibrium phases for interacting two-dimensional self-propelled particles with isotropic pair-wise interactions using a persistent kinetic Monte Carlo (MC) approach. We establish the quantitative phase diagram, including the motility-induced phase separation (MIPS) that is a commonly observed collective phenomena in active matter. In addition, we demonstrate for several different potential forms the presence of two-step melting, with an intermediate hexatic phase, in regions far from equilibrium. Increased activity can melt a two-dimensional solid and the melting lines remain disjoint from MIPS. We establish this phase diagram for a range of the inter-particle potential stiffnesses, and identify the MIPS phase even in the hard-disk limit. We establish that the full description of the phase behavior requires three independent control parameters.

Figures

Tables

Table 1. Table 1: Decay of correlation functions in the liquid, hexatic, and solid phases in two-dimensional particle systems.

Order	Liquid phase	Hexatic phase	Solid phase
Positional	short-range111$\propto \exp \left(-r/\xi \right)$	short-range1	quasi-long-range222$\propto r^{-\alpha }$
Orientational	short-range1	quasi-long-range2	long-range333$\propto$ constant

Equations

$$\lambda \simeq 0.62\frac{{\delta }^3}{{\sigma }^2},$$

$$\tau =\frac{8}{{\pi }^2}\frac{{\delta }^2}{{\sigma }^2}.$$

$$P(E^{\prime }\to E)=\min \left[1,e^{-\beta \Delta E)}\right],$$

$$U(r)=u_0{\left(\frac{\gamma }{r}\right)}^n,$$

$${\psi }_6({\mathbf{r}}_i)=\frac{1}{\text{number of neighbors j of i }}\sum _j\exp (6ı{\theta }_{ij}),$$

$$g_6(r)\propto ⟨\sum _{i,j}^N{\psi }_6^{\star }({\mathbf{r}}_i){\psi }_6({\mathbf{r}}_j)\delta (r-r_{ij})⟩$$

$${ϵ}_{k+1}$$

$$x_{k+1}$$

$$f(x,ϵ)=\min \left\{1,\exp \left(-\frac{U(x+ϵ)-U(x)}{k_{\text{B}}T}\right)\right\}$$

$$t=k\delta ,v(t)=\frac{{ϵ}_k}{\delta },\xi (t)=\frac{r_k}{{\delta }^2},x(t)=x_k,$$

$$\dot{v}(t)$$

$$\dot{x}(t)$$

$$\frac{U(x+v\delta )-U(x)}{k_{\text{B}}T}={\Gamma }_{1}v{U}^{\prime }(x)\left[1+\frac{\delta v}{2}\frac{U^{\prime \prime }(x)}{U^{\prime }(x)}+\cdots \right],$$

$$h[x(t),v(t)]=\min \left\{1,\exp \left(-{\Gamma }_{1}v{U}^{\prime }(x)\right)\right\}.$$

$$\dot{v}(t)$$

$$\dot{x}(t)$$

$${\Gamma }_1=\frac{\delta }{k_{\text{B}}T},\text{ and }\lambda =\frac{{\delta }^3}{{\sigma }^2}.$$

$$M(x,ϵ|y,{ϵ}^{\prime })=g(ϵ,{ϵ}^{\prime })W_{ϵ}(x,y),$$

$$g(ϵ,{ϵ}^{\prime })=\frac{1}{2\delta }+$$

$$\cos$$

$$W_{ϵ}(x,y)=f(y,ϵ)\delta (x-y-ϵ)+[1-f(y,ϵ)]\delta (x-y),$$

$$P_{k+1}(x,ϵ)=\int dy{\int }_{-\delta }^{\delta }d{ϵ}^{\prime }M(x,ϵ|y,{ϵ}^{\prime })P_k(y,{ϵ}^{\prime })$$

$$P_{k+1}(x,ϵ)=\int d{ϵ}^{\prime }g(ϵ,{ϵ}^{\prime })P_k(x,{ϵ}^{\prime })+\int d{ϵ}^{\prime }g(ϵ,{ϵ}^{\prime })$$

$$\left\{f(x-ϵ,ϵ)P_k(x-ϵ,{ϵ}^{\prime })-f(x,ϵ)P_k(x,{ϵ}^{\prime })\right\}.$$

$${\tilde{P}}_{t+\delta }(x,v)=$$

$$+\int d{v}^{\prime }\tilde{g}(v,v^{\prime })\{f(x-\delta v,\delta v){\tilde{P}}_t(x-\delta v,\delta v^{\prime })$$

$$-f(x,\delta v){\tilde{P}}_t(x,\delta v^{\prime })\},$$

$$\frac{\partial }{\partial t}{\tilde{P}}_t(x,v)=$$

$$-\int d{v}^{\prime }\tilde{g}(v,v^{\prime })v\frac{\partial }{\partial x}\left\{h(x,v){\tilde{P}}_t(x,v^{\prime })\right\}+\cdots ,$$

$$a_m(v)=\frac{1}{\delta }\int d{v}^{\prime }(v^{\prime }-v{)}^m\tilde{g}(v,v^{\prime }).$$

$$\tilde{g}(v,v^{\prime })=\frac{1}{\sqrt{2\pi \delta /\lambda }}\exp \left[-\frac{{\left(v-v^{\prime }\right)}^2}{2\delta /\lambda }\right]$$

$$\frac{\partial }{\partial t}{\tilde{P}}_t(x,v)=\frac{1}{2\lambda }{\partial }_v^2{\tilde{P}}_t(x,v)-\frac{\partial }{\partial x}\left\{vh(x,v){\tilde{P}}_t(x,v)\right\}$$

$$\frac{\partial }{\partial v}{\tilde{P}}_t(x,v)=0\text{for v=±1}.$$

$$t=k{\delta }^2,v=\frac{ϵ}{\delta },x=x.$$

$$\frac{\partial {\tilde{P}}_t(x)}{\partial t}=\frac{1}{k_{\text{B}}T}\frac{\partial }{\partial x}[U^{\prime }(x){\tilde{P}}_t(x)]+\frac{{\partial }^2}{\partial x^2}[{\tilde{P}}_t(x)].$$

Full text

I Introduction

Active matter is an important field of research that considers particle systems whose microscopic components are characterized by systematic persistent dynamic rules and by various types of mutual interactions. On a microscopic scale, the persistent dynamics breaks the detailed-balance condition (underlying all of equilibrium physics) and defines active matter as out-of-equilibrium systems.

Many different models of active matter have been proposed. They feature a wide range of self-propelled dynamics and of mutual interactions. A great many theoretical studies employ either Langevin or molecular dynamicsCates (2012); Redner, Hagan, and Baskaran (2013); Fily and Marchetti (2012); Speck et al. (2014); Bialké, Speck, and Löwen (2012); Cugliandolo et al. (2017) in order to model persistent motion. Recently, a kinetic Monte Carlo (MC) approach was proposedLevis and Berthier (2014); Klamser, Kapfer, and Krauth (2018) as a minimal model for active matter in two dimensions.

Although the primary interest in theoretical models of active matter comes from the non-equilibrium nature, their properties can often be connected to their equilibrium counterparts that are realized in the zero-persistence limit. This limit is of particular interest in two dimensions, where equilibrium particle systems with short-range interactions cannot crystallizeMermin (1968). Nevertheless, it was established that, from the high-temperature (low-density) liquid regime towards the low-temperature (high-density) solid regime, two-dimensional equilibrium particle systems normally undergo two phase transitionsHalperin and Nelson (1978); Nelson and Halperin (1979); Young (1979). These transitions describe the passage into and out of a hexatic phase that is sandwiched between the liquid and the solid phase (see Table 1). In this work, we are concerned with two-dimensional models with repulsive inverse-power-law interactions, for which the two-step melting scenario in equilibrium is firmly establishedBernard and Krauth (2011); Kapfer and Krauth (2015).

In this work, we extend our earlier findings for a special caseKlamser, Kapfer, and Krauth (2018) and show that kinetic MC generically reproduces motility-induced phase separation (MIPS) in two-dimensional active-particle systems with a wide range of inverse-power-law interactions including the hard-sphere limit. We also confirm the stability of the hexatic phase up to considerable values of the activity, and conjecture that it is indeed stable at any value of the activity. We finally confirm that MIPS is generically a liquid–gas transition under the kinetic MC dynamics. Moreover, it is decoupled from the melting transitions. This separation can be understood in the limit of infinitesimal MC steps from the scaling behavior of the MIPS phase transition and the melting transitions.

The work is organized in the following order. In Section II, the essential elements of the kinetic MC algorithm are described, the interplay of persistence with interactions is illustrated on a simple case of two particles in one dimension (1D), and possible anisotropy effects in two dimensions (2D) are analyzed. In Section III, we discuss the effect of the stiffness of the interparticle potential on the full quantitative phase diagram of two-dimensional particle systems on the activity–density plane. The continuous-time limit of the kinetic MC dynamics is discussed in Section IV, with a focus on the number of relevant parameters.

II Model: kinetic MC

The kinetic MC algorithm with which we model active dynamics consists of the standard Metropolis filter combined with a memory term for the proposed moves. The memory term is characterized by a time scale $\tau$, which allows for a smooth tuning from a passive motion (described by equilibrium statistical mechanics) to a self-propelled/persistent particle motion, where a single particle moves ballistically (mean square displacement $\propto t^{2}$) for times $t\ll\tau$ and moves diffusively (mean square displacement $\propto t$) for times $t\gg\tau$ (see Fig. 1(a)).

In contrast to active Brownian particles Romanczuk et al. (2012), the velocity amplitude fluctuates in the MC dynamics. The dynamics is comparable with the active Ornstein–Uhlenbeck process (see e.g. Ref. Fodor et al., 2016), where the velocity performs a random walk in a harmonic potential. Similarly, in the discrete-time kinetic MC approach the increment performs a random walk in a box with reflecting boundary conditions. For a single particle in two dimensions, the $x$- and $y$-components of the increment $(\epsilon_{x}(t),\epsilon_{y}(t))$ at time $t$ are sampled from two independent Gaussian distributions of standard deviation $\sigma$ and the mean corresponds to the previously sampled increment ($\epsilon_{x}(t-1)$ or $\epsilon_{y}(t-1)$, respectively). This walk is confined in $x$ and $y$ by reflecting boundaries at $x=\pm\delta$ and $y=\pm\delta$. In two dimensions the average distance a single particle covers before changing the direction is given byKlamser, Kapfer, and Krauth (2018) the persistence length

[TABLE]

and the persistence time is given by

[TABLE]

This characteristic length and time scales are confirmed in numerical simulations. For example, the crossover from ballistic to diffusive motion in Fig. 1(b) appears in the rescaled time-dependence of the mean absolute displacement around the point $(t/\tau=1,\langle r(t)\rangle/\lambda=1)$.

In the many-particle case, the increment for each particle performs its proper random $\bm{\epsilon}$-walk, independent of the other particles. Interactions between particles are introduced by the Metropolis filter. In one kinetic MC step, a particle $i$ is chosen at random. The change of its position ($\bm{r}_{i}(t+1)=\bm{r}_{i}(t)+\bm{\epsilon}_{i}(t)$) is accepted with probability

[TABLE]

where $\Delta E=E-E^{\prime}$ is the total energy change caused by the particle displacement $\bm{r}_{i}(t)\rightarrow\bm{r}_{i}(t+1)$. The parameter $1/\beta=k_{\text{B}}T$ is now an energy scale rather than a temperature. The random $\bm{\epsilon}$-walk persists whether the resulting displacement is accepted or not. For a detailed description of the kinetic MC approach, see Ref. Klamser, Kapfer, and Krauth, 2018.

II.1 1D Model for Persistence

Kinetic MC, via the memory present in its displacements, generates persistence in a manner that differs from equilibrium systems. This has far-ranging consequences for many-particle systems, but the effects of a memory term in the equations of motion can already be studied for $N=1$ or $N=2$ particles. Here we study the case of two particles on a ring (a line of length $L$ with periodic boundary conditions), interacting with an inverse-power-law pair potential that we will use throughout this work

[TABLE]

where $\gamma$ reduces to the particle diameter in the hard-disk limit $n\to\infty$. The 1D inter-particle distance $r$ in eq. (4) is easily generalized to higher dimensions.

Fig. 2 shows the probability distribution $P(r)$ of the inter-particle distance. The maximum jump length $\delta$ is kept constant and the persistence length $\lambda$ is varied by changing $\sigma$. As the interaction is repulsive, the two particles repel each other for small $\lambda$, and the Boltzmann weight is maximal at $r=L/2$ for $\lambda=0$ (see Fig. 2(c)). For increased $\lambda$, a peak appears at small $r$, and its position decreases with increasing $\lambda$. This means that particles are more probable to be near each other, which is due to the case, where the particles try to move against each other (see Fig. 2(a)). The higher the persistence length, the stronger the force the particles push against the interparticle potential barrier and therefore, the peak position shifts. This shows that the self-propulsion force increases with $\lambda$ or equivalently with the persistence time $\tau$, and the shift is a result of the larger number of attempts in the Metropolis filter to increase the total energy.

At finite $\lambda$, the original (Boltzmann) peak at $r=L/2$ shifts to smaller $r$ and takes on a position that is independent of $\lambda$. This peak appears due to arrangements where one particle “hunts” the other circling around the ring with $|\epsilon_{\text{hunter}}|>|\epsilon_{\text{hunted}}|$ (see Fig. 2(b)). For this to have a non-negligible probability, the role of the slower and faster particle has to be stable for a sufficient time span, which explains that the peak appears only after overcoming a certain threshold in $\lambda$. The independence of the peak position on $\lambda$ is understood from the following argument. Moves of the slower particle are always accepted by the Metropolis filter, as they decrease the total energy. In contrast, moves of the faster particle are often rejected as they increase the total energy. This leads to a competition where the slower particle increases $r$ in every attempt, and the faster particle tries to decreases $r$, but does not succeed in every attempt, thereby leading to an average "hunting distance" $r$, which is independent of $\lambda$ (see Fig. 2(c)).

The bimodal probability distribution of $P(r)$ (see Fig. 2) is a consequence of the non-constant velocity amplitude and thus in contrast to e.g. active Brownian particles. The result clearly shows that the self-propulsion force in the kinetic MC dynamics depends on the persistence time $\tau$. This is not the case in the Langevin approaches of active Brownian particles or the active Ornstein–Uhlenbeck process. However, as discussed in Section IV, the persistent length $\lambda$ is the relevant measure for activity is this dynamics.

II.2 Anisotropic effects

In our kinetic MC algorithm, displacements $\bm{\varepsilon}=(\epsilon_{x},\epsilon_{y})\in[-\delta,\delta]^{2}$ are confined to a square box rather than being sampled from an isotropic distribution (as for example a circle of radius $\delta$). (The 2D Gaussian distribution of the Ornstein–Uhlenbeck process is also isotropic.) Although the square box is chosen for simplicity only, it is useful to check that it does not induce anisotropies. This trivially follows in the passive limit for the steady-state probability distribution because of the detailed-balance condition.

At small times ($t\ll\tau$), the reflecting boundary conditions for the sampling scheme of the displacements introduces some degree of anisotropy in the two-dimensional single-particle dynamics (see Fig. 3(a)). The probability distribution $P_{\Delta t}(|\Delta x|,|\Delta y|)$ of the absolute particle displacement $\Delta x$ (or $\Delta y$) in the $x$ (or $y$) component in a time $\Delta t$ is anisotropic for $\Delta t<\tau$ (see Fig. 3(a)) whereas, without the square box, this displacement (which is constructed from identically distributed, independent Gaussians in both dimensions) would be isotropic. However, the isotropy is reinstalled for $t\gg\tau$ (see Fig. 3(b)). As $\tau=0$ in the passive limit, the particle dynamics is isotropic for all timesKrauth (2006). The anisotropy for $t<\tau$ could also be avoided by choosing a circular sampling box of radius $\delta$ with reflective boundary conditions for the displacements (see Fig. 3(b)). However, such a choice would be more costly to implement.

We conjecture that the square box in fact renders anisotropic the long-time dynamics neither for $N=1$ nor for the many-body case. In the dilute case (where $\lambda$ is much smaller than the mean free path) the kinetic MC dynamics effectively reverts to the detailed-balance dynamics as interactions between particles happen at the diffusive time scale. At higher densities, anisotropy in the many-body properties might arise if the probability distribution of the accepted displacements is itself anisotropic. However, this is not the case (see Fig. 4). At higher densities, all large proposed displacements have a vanishing probability to be accepted by the Metropolis filter, thus leading to an effectively isotropic dynamics. Therefore, the Metropolis filter effectively realizes a circular reflecting $\bm{\varepsilon}$-sampling box without additional computational cost. Pair-correlation functions are also found to be perfectly isotropic at high density, both in the motility-induced liquid phase (at a density deep inside the equilibrium solid phase) and in the MIPS region (see Fig. 5).

III Two-dimensional simulation results

Our simulations are performed in an ensemble of $N$ particles confined to a rectangularKlamser, Kapfer, and Krauth (2018) box of volume $V$ with periodic boundary conditions. The density $\phi=\gamma^{2}N/V$ (with $\gamma$ from eq. (4)) is varied by changing $V$. In the simulations $\delta$ is kept constant and the activity is varied by changing $\sigma$.

III.1 Full phase diagram and the effect of

stiffness

Fig. 6(a) shows the full phase diagram for the potential in eq. (4) with $n=16$ on the $\phi$–$\lambda$ plane. At all $\lambda$, the equilibrium two-step melting transitionKapfer and Krauth (2015) is recovered. For increasing $\lambda$, the melting lines shift to higher densities. Remarkably, the hexatic phase separating the liquid and solid phases is stable even far from equilibrium. This non-equilibrium two-step melting can be induced either by reducing the density (just as in equilibrium) or by increasing the persistence length. In addition to these melting transitions, at low $\phi$ but high $\lambda$, a motility-induced liquid–gas coexistence region opens up. It is separated form the melting transitions by a disordered fluid phase. This generalizes the phase diagram under the same dynamics, but for $n=6$, found previouslyKlamser, Kapfer, and Krauth (2018). A change of $n$ (see eq. (4)) only shifts the positions of the phase boundaries (see Fig. 6(b)). As for the equilibrium melting transitions Kapfer and Krauth (2015), both the liquid–hexatic and the hexatic–solid phase boundaries shift at constant persistence length $\lambda$ to smaller densities with increasing $n$.

Increasing $\lambda$ shifts the melting transitions to higher densities. However, the shift is smaller for larger $n$, resulting in steeper transition lines (see Fig. 6(b)). At the same time, the onset of the motility-induced liquid–gas coexistence shifts to smaller $\phi$ and $\lambda$ and the coexisting region shrinks. This ensures that the liquid–gas coexistence and the melting transitions remain disjoint. In Section III.3, we argue that there is no singular change in the phase diagram even in the hard-disk limit $n\to\infty$.

III.2 Non-equilibrium two-step melting

In equilibrium, the Mermin–Wagner theorem forbids long-range translational (i.e. crystalline) order in a 2D particle system with short-range interactions Mermin and Wagner (1966); Mermin (1968). However, at large densities, particles can arrange in locally hexagonal configurations. This can lead to two different high-density phasesNelson and Halperin (1979); Halperin and Nelson (1978); Young (1979), which are characterized by different degrees of orientational and positional order (see Table 1). In this section, we define these measures of order and use them to quantify the two-step melting far from equilibrium.

The local bond-orientational order parameter $\psi_{6}(\bm{r}_{i})$ measures the six-fold orientation near a particle $i$. It is defined as

[TABLE]

where $\imath$ is the imaginary unit and $\theta_{ij}$ it the angle enclosed by the $x$-axis and the connection line between particle $i$ and its neighbor $j$. Here, we use the Voronoi construction to identify neighbors and $\psi_{6}(\mathbf{r}_{i})$ is calculated with Voronoi weightsMickel et al. (2013). Then, the correlation function

[TABLE]

is a measure of the correlation of the local six-fold orientational order at distance $r$ and its decay is used to quantify the degree of orientational order in the system (see Table 1).

The direction-dependent pair-correlation function $g(x,y)$ provides a measure for the positional order. This two-dimensional histogram is averaged over different configurations $C$ after re-aligning Bernard and Krauth (2011) $g_{C}(x,y)$ such that the $\Delta$x-axis points in the direction of the global orientation parameter $\Psi_{6}(C)=\sum_{i}^{N}\psi_{6}(\bm{r}_{i})$ of $C$. Then, the decay of, e.g. $g(x,0)$ determines the degree of positional order.

The correlation functions ($g(x,0)$ and $g_{6}(r)$), allow one to identify the two-dimensional phases by the properties summarized in Table 1. In equilibrium, the Kosterlitz–Thouless–Halperin–Nelson–Young theory provides an additional selection criterionNelson and Halperin (1979); Halperin and Nelson (1978); Young (1979), where the exponent (defined in Table 1) $\alpha\leq 1/4$ for the orientational order and $\alpha\leq 1/3$ for the positional order give theoretical bounds for the hexatic phase. However, these bounds are not shown to apply outside equilibrium. We thus identify the phases by the characteristic decay of $g(x,0)$ and $g_{6}(r)$.

Orientational and positional correlation functions change as the system melts from solid to liquid passing through the hexatic phase (see Fig. 7, at a density far above the equilibrium melting point). Our simulations clearly identify solid state points with power-law decay in $g(x,0)$ and constant $g_{6}(r)$, hexatic state points with exponential decay of $g(x,0)$ and quasi-long-range order in $g_{6}(r)$, and liquid state points, where both correlations decay exponentially. Snapshots illustrate these different phases (see Fig. 7 and Table 1). The power-law exponent in $g_{6}(r)$ grows when approaching the liquid–hexatic transition and therefore, weakening the order in the hexatic phase. This behavior of $g_{6}(r)$ is in agreement with equilibrium studiesKapfer and Krauth (2015) and with our earlier resultsKlamser, Kapfer, and Krauth (2018) obtained with the same dynamics for $n=6$.

We check that our simulations indeed reach the steady state (as in earlier workKlamser, Kapfer, and Krauth (2018)) by verifying the convergence of the spatial correlation functions to the same steady state starting from a crystalline and a liquid initial particle arrangement. These very time-consuming computations assure that the defining decays of the correlation functions obtained in the hexatic phase reflect the physical system and not a bias introduced by the initial condition.

III.3 Motility-induced phase separation

At sufficiently low densities and high activities, a liquid–gas coexistence region opens up with a roughly U-shaped phase boundary (see Fig. 6(a)). Following an analysis of local densities, which was applied earlierKlamser, Kapfer, and Krauth (2018) for $n=6$, we confirm that also for $n=16$ the densities ($\phi_{\text{Liquid}}$ and $\phi_{\text{Gas}}$) of the two coexisting phases depend on $\lambda$ but not on the global density $\phi$. Therefore, the low-density boundary of the MIPS region in Fig. 6(a) is given byKlamser, Kapfer, and Krauth (2018) $\phi=\phi_{\text{Gas}}(\lambda)$ and the high-density boundary by $\phi=\phi_{\text{Liquid}}(\lambda)$, respectively.

A much discussed questionFily and Marchetti (2012); Redner, Hagan, and Baskaran (2013); Buttinoni et al. (2013); Cugliandolo et al. (2017); Digregorio et al. (2018) concerns the nature of the two phases at coexistence. We can clearly identify the high-density phase as liquid, and MIPS as a liquid–gas coexistence. The particles in the snapshots in Fig. 8 are $\psi_{6}$-color-coded (see Fig. 7 for definition), illustrating short-ranged orientational order in the liquid phase. We do not observe that the orientational correlation in the liquid phase increases with increasing $\lambda$. Even at very high activities (e.g. at $\lambda=4\times 10^{3}$ in Fig. 8(a)), the local orientational order changes upon length scales of the order of the interparticle distance (also compare with Fig. 7(c) point E).

We do not observe any quantitative difference with the orientational order in MIPS previously observedKlamser, Kapfer, and Krauth (2018) for $n=6$ for the same kinetic MC dynamics. Furthermore, we also recover MIPS in the hard-disk system in the form of a liquid–gas coexistence (see Fig. 8(b)). We conjecture from these findings that the separation of the MIPS region from the melting transitions is a generic feature of self-propelled particles, at least within kinetic MC dynamics.

IV Relevant parameters

So far (as our previous work Klamser, Kapfer, and Krauth (2018)), we have considered the phase diagram as a function of the persistence length and the density, keeping the maximum step size $\delta$ constant (see Fig. 6). However, $\delta$ has a profound influence on the phase boundaries (see Fig. 9). Keeping $\phi$ constant, the melting lines shift to smaller activities as $\delta$ is decreased, while in contrast the MIPS phase boundary shifts to larger activities. In this section, we study this $\delta$-dependence of melting and of MIPS.

The kinetic MC dynamics depends on three parameters ($\delta$, $\lambda$, $\phi$). We now show that MIPS (seen at high $\lambda$) and the melting close to equilibrium are described by a different reduced set of relevant parameters in the $\delta\to 0$ limit. The single relevant parameterKapfer and Krauth (2015), which describes the melting transitions for inverse-power-law potentials is not commensurate with the reduced parameters for MIPS, as it does not capture the critical melting density in the passive limit. Therefore, there are separate descriptions of MIPS and the melting transitions.

Here we introduce the Master equation as a stochastic descriptions of our kinetic MC dynamics, in addition to a Langevin description and the associated Fokker–Planck equation. We also compare the dimensional reduction of the relevant parameters with other stochastic models of active matter.

IV.1 A simple argument

We first address the question of relevant parameters with a simple argument for a single particle. (A more detailed analysis is presented in the following Section IV.2 and Section IV.3.) We consider the small-$\delta$ limit, which is justified for the choice of parameters used in the simulations. For a single particle in a 1D confining potential $U(x)$, the kinetic MC rule is approximated by the following discrete-time ($k=0,1,2,\dots$) dynamics:

[TABLE]

where $r_{k}$ is a Gaussian random number with $\langle r_{k}\rangle=0$, $\langle r_{k}r_{k^{\prime}}\rangle=\sigma^{2}\delta_{k,k^{\prime}}$ and

[TABLE]

is the acceptance rate of the Metropolis filter in eq. (3). The reflecting boundary at $\epsilon=\pm\delta$ is denoted by $R$ (without specifying it rigorously) and $\delta_{k,k^{\prime}}$ is the Kronecker delta.

Defining a set of rescaled coordinates,

[TABLE]

and taking the small-$\delta$ limit, we get

[TABLE]

where $\langle\xi(t)\rangle=0$, $\langle\xi(t)\xi(t^{\prime})\rangle=\lambda^{-1}\delta(t-t^{\prime})+\mathcal{O}(\delta/\lambda)$, with444The persistence length used for the numerical results eq. (1) differs by a numerical constant in both one and two dimensions. $\lambda=\delta^{3}/\sigma^{2}$ and we use a continuous limit of the Kronecker delta $\delta_{k,k^{\prime}}\simeq\delta\,\delta(t-t^{\prime})$, with $\delta(x)$ being the Dirac delta function. Moreover, in the small-$\delta$ limit, using Taylor expansion, we get

[TABLE]

where $\Gamma_{1}=\delta/k_{\text{B}}T$. As long as $\delta U^{\prime\prime}(x)/U^{\prime}(x)\ll 1$, from eq. (7), we get $f[x(t),v(t)\delta]\simeq h[x(t),v(t)]$, with

[TABLE]

This gives a continuous dynamics for $\delta\to 0$

[TABLE]

Clearly, this rescaled dynamics depends on only two relevant parameters, namely

[TABLE]

However, for this description to be valid, the subleading terms have to be negligible, leading to the following range of validity for the scaling in eq. (13):

• a)

$\delta\ll\lambda$ (from $\langle\xi(t)\xi(t^{\prime})\rangle$),

• b)

$\delta\ll\lambda^{-1/2}$ (from eq. (12a)), and

• c)

$\delta\ll 1$ and $\delta\ll U^{\prime}(x)/U^{\prime\prime}(x)$ (from eq. (12b)).

Therefore, the scaling expressed in the two-parameter reduction in eq. (13) breaks down in a) the passive limit, b) at very high persistence lengths (for constant $\phi$), and c) at high density (for constant $\lambda$).

IV.2 Stochastic description of the MC dynamics

We begin the rigorous analysis by considering a single particle in a 1D confining potential $U(x)$. The kinetic MC dynamics is Markovian in the $(x,\epsilon)$ space. The discrete kinetic MC time is denoted by $k=0,1,2,\dots$. The conditional probability for a transition $(y,\epsilon^{\prime})\rightarrow(x,\epsilon)$ in one time step is given by the Markov matrix

[TABLE]

where, $\epsilon$ is sampled with probability $g(\epsilon,\epsilon^{\prime})$ and $W_{\epsilon}(x,y)$ is due to the Metropolis filter. It can be shownKlamser (2018) that

[TABLE]

whereas the Metropolis filter in eq. (3) yields

[TABLE]

with $f(x,\epsilon)$ as defined in eq. (7). Using this in the corresponding discrete-time Master equation

[TABLE]

gives

[TABLE]

This describes the exact time evolution of the probability $P_{k}(x,\epsilon)$ in the kinetic MC dynamics555The approximate dynamics in eq. (6) has a different Master equation. Nevertheless, both describe the same dynamics in the small-$\delta$ limit.. In the passive limit, this satisfies the standard detailed-balance condition with respect to the Boltzmann distribution (see Appendix A).

Rescaled coordinates

To determine the relevant number of control parameters, we use the scaled coordinates defined in eq. (8) in the Master equation (eq. (19)):

[TABLE]

with $P_{k}(x,\epsilon)=\tilde{P}_{t}(x,v)/\delta$ and $g(\epsilon,\epsilon^{\prime})=\tilde{g}(v,v^{\prime})/\delta$.

We have shown in Section IV.1 that the effective number of control parameters can be reduced by taking the small-$\delta$ limit. In this limit $f(x,\delta v)\simeq h(x,v)$ (see eq. (11)). Using a Taylor expansion, this leads to

[TABLE]

where we used $\int dv^{\prime}\,\tilde{g}(v,v^{\prime})=1$ and

[TABLE]

$a_{m}(v)$ can be computed using eq. (16). Alternatively, we can use the free case

[TABLE]

in combination with a zero-current condition on $\tilde{P}_{t}(x,v)$ for the reflecting boundary. This simplifies the calculation of $a_{m}(v)$, giving $a_{1}(v)=0$, and $a_{2}(v)=1/\lambda$, which leads to the Fokker–Planck equation at small $\delta$,

[TABLE]

This Fokker–Planck equation is equivalent to the coupled Langevin equation in eq. (12) of the approximate dynamics. In consistence with the analysis in Section IV.1, this Fokker–Planck equation depends on two parameters ($\lambda$ and $\Gamma_{1}$) given in eq. (13). However, as noted earlier, this is only true in a certain range of parameters where the subleading terms are negligible. In particular, the Fokker–Planck equation in eq. (24b) fails to describe the passive limit ($\lambda=0$), where the relevant parameter is different Kapfer and Krauth (2015) from $\Gamma_{1}$.

To describe the passive limit, a diffusive scaling with $\delta$ is required

[TABLE]

Starting from eq. (19) and following a similar procedure as for the derivation of eq. (24a), leads to the well-known Fokker–Planck equation for a passive particle in a potential

[TABLE]

A detailed derivationKlamser (2018) finds back that in the passive limit the inverse-power-law potential is described by a single control parameter, as is well knownKapfer and Krauth (2015).

In conclusion, the number of relevant parameters can be reduced when $\delta$ is small: the persistent limit has two relevant parameters ($\Gamma_{1}$ and $\lambda$) and the passive limit ($\lambda=0$) has a single parameter $\Gamma_{0}$ (see eq. (40)). However, $\Gamma_{1}$ does not converge to $\Gamma_{0}$ for $\lambda\to 0$. Therefore, the order of the limits $\delta\to 0$ and $\lambda\to 0$ cannot be exchanged.

IV.3 Multi-particle case

The discussion in Section IV.2 can be generalized to the multi-particle case. The single-particle Fokker–Planck equation in eq. (24a) generalizes to the $N$-particle Fokker–Planck equation

[TABLE]

where particles interact via the inter-particle potential $U(\mathbf{x})$, with $\mathbf{x}=\left\{x_{1},\cdots,x_{N}\right\}$ and

[TABLE]

with the force on particle $i$

[TABLE]

The reflecting boundary condition in the velocity space is

[TABLE]

The Fokker-Planck (26) in the passive limit has a very similar multi-particle generalization.

Power-law interaction potential

The numerical studies presented in this work are for the inverse-power-law potential in eq. (4). The force on particle $i$ is

[TABLE]

For this specific choice, two dimensionless parameters characterize the persistent many-particle behavior. These can be obtained from the mean inter-particle distance

[TABLE]

with $L$ being the system size or, equivalently, from the dimensionless density

[TABLE]

Then, using the scaled coordinates

[TABLE]

in eq. (27), we obtain the rescaled Fokker–Planck equation

[TABLE]

with

[TABLE]

and

[TABLE]

These two parameters ($\tilde{\lambda}$ and $\Gamma$) govern the scaled many-particle probability distribution.

It is useful to express these parameters in terms $\phi$, which leads to

[TABLE]

The generalization to higher dimensions is straightforward. For example, in 2D, where $\phi=N\gamma^{2}/V$ and $d=\sqrt{V/(N\pi)}=\gamma/\sqrt{\pi\phi}$, the two relevant parameters are

[TABLE]

In this form, the dimensionality primarily enters into the $\phi$-dependence of these two parameters (see eq. (36)). The validity of these reduced parameter sets is confirmed by numerical simulations (see Fig. 10).

As discussed earlier in the single-particle case, this parameter reduction does not extend to the passive regime. In this case a diffusive scaling is required

[TABLE]

in terms of which the Fokker–Planck equation becomesKlamser (2018)

[TABLE]

with

[TABLE]

Therefore, the probability only dependsKapfer and Krauth (2015) on $\Gamma_{0}$, which differs from $\Gamma$.

In summary, the small-$\delta$ limit of the systems described in this work is described by two relevant parameters for the active regime and a single relevant parameter for the passive regime. There is no smooth transition from one regime to the other when $\delta$ is small. This is because the limits $\delta\to 0$ and $\lambda\to 0$ are not exchangeable. Therefore, a phase diagram covering both passive and active is not possible on a two parameter space. A complete phase diagram is three dimensional on a rescaled parameter space $\delta/d$, $\lambda/d$ and $\Gamma_{0}$ in eq. (40). Expressed in terms of $\phi$, they yield

[TABLE]

Continuous-time description

In the analysis presented so far, the Fokker–Planck description was obtained in the small-$\delta$ limit. However, it is possibleKlamser (2018) to define a continuous-time description of the MC dynamics for arbitrary values of $\delta$. This can be obtained using the Kramers–Moyal expansionKampen (2007) of the Master equation eq. (19), which leads to the Fokker–Planck equation

[TABLE]

with the reflecting boundary at $\epsilon=\pm\delta$ as zero-current condition

[TABLE]

Taking the small-$\delta$ limit of eq. (42), we recover the scaling of eq. (37) that was obtained from the discrete-time Master equation eq. (19).

IV.4 Active random acceleration process

An intriguing aspect of the dimensional reduction of the phase space is the singular nature of the passive case on a two-parameter plane for $\delta\to 0$. This singular behavior of the rescaled parameters is not due to the discrete nature of the MC dynamics as it also appears in the continuous-time description discussed above. For the kinetic MC, the singularity appears in the small-$\delta$ limit. However, it can also appear in other stochastic models of active matter, even without taking a vanishing limit of a parameter similar to $\delta$. To illustrate this point, we define a continuous model which closely resembles the kinetic MC dynamics. This model, which we refer as the active random acceleration processMajumdar, Rosso, and Zoia (2010) is defined for continuous time by a coupled Langevin equation

[TABLE]

where $x_{t}$ is the position of the self-propelled particle at time $t$. The two noise terms, $r_{t}$ and $s_{t}$, are Gaussian white noises with zero mean and covariance $\langle r_{t}r_{t^{\prime}}\rangle=\sigma^{2}\delta(t-t^{\prime})$ and $\langle s_{t}s_{t^{\prime}}\rangle=2D\delta(t-t^{\prime})$, respectively. The term $R$ denotes the reflecting boundary condition at $\pm\delta$. The parameter $v_{0}$ characterizes the strength of the self-propulsion666Putting $1/k_{\text{B}}T$ in the force term of eq. (44b) is motivated by the aim to compare with the kinetic MC results..

Similar to the kinetic MC dynamics, in the scaled coordinates $t\to td/\delta$, $\epsilon\to v\delta$, $x\to xd$ the dynamics eq. (44) is described by the Fokker–Planck equation

[TABLE]

Its multi-particle generalization is straightforward. For the inter-particle force eq. (30), the steady-state probability distribution is controlled by three parameters

[TABLE]

The passive limit corresponds to $v_{0}=0$, where the $x$ and $v$ coordinates decouple and the probability $P(x,v)=P(x)P(v)$. In such case, we see that the $P(x)$ follows the usual equilibrium Fokker–Planck equation and the steady-state $P(x)$ is determined by a single control parameter $C/A$, which coincides with $\Gamma_{0}$ for $D=1$. Then, importantly, $P(x)\propto\exp(-Cu(x)/A)$, where $u(x)$ is the rescaled potential and therefore, setting $D=0$ makes the probability distribution uniform and independent of the inter-particle interaction, which corresponds to an infinite temperature.

In the active limit, there are the three control parameters $A$, $B$, and $C$. The only way to reduce the number of parameters, while keeping the particles active and interacting, is by setting $D=0$ and thereby $A=0$. However, this would make the finite-temperature passive limit inaccessible, as we discussed above. This shows that the singular behavior of the MC dynamics, when trying to reduce the number of control parameters to two also happens in this active random acceleration process. A similar singular behavior is present for other models. We have checkedKlamser (2018) that the same statement applies for the active Ornstein–Uhlenbeck process.

From this analysis, we can conclude that in these classes of active dynamics, at least three independent relevant parameters are required to describe the full phase diagram including the passive regime.

V Conclusions

In this work, we presented a kinetic-Monte Carlo perspective on two-dimensional active matter. Within this approach, we established (in extension of our earlier work Klamser, Kapfer, and Krauth (2018)) the presence of the liquid, hexatic, and solid active-matter phases from their constituent decay laws of positional and orientational order (see Table 1). We also ascertained continuity of the active-matter phases in the passive limit and recovered the phases of the corresponding equilibrium system. We have not tested very soft repulsive inverse-power-law potentials ($n<6$ in eq. (4)), but expect on the basis of our findings that for sufficiently steep repulsive inverse-power-law potentials the two-step melting behavior of the equilibrium system is maintained up to high activities, and possibly up to infinite persistence lengths. The stability of the intermediate hexatic phase in the high-activity regime—far above the linear-response regime—is intriguing. It may be due to an underlying symmetry that is yet to be understood. We have not addressed here the nature of the melting transitions in the active region, that is, the question of whether the liquid–hexatic phase transition is of first order or of Kosterlitz–Thouless typeKosterlitz and Thouless (1972, 1973) (as in the equilibrium system Kapfer and Krauth (2015)), and whether this theory continues to apply at all to the hexatic–solid transition in active systems.

Besides the melting phase transitions, we have established the existence of the MIPS phase and identified it as a coexistence between a liquid and a gas: two phases with exponential decays of both correlation functions, but with typically very different densities. The coexistence region was identified for a wide range of interaction potentials including the hard-disk limit. It remains disjoint from the melting lines for all values of $n$ (possibly excluding very soft potentials). An analogy of the MIPS with the liquid–gas coexistence in equilibrium indicates the existenceSiebert et al. (2018) of a critical point, although we have not studied it in detail.

We also discussed the dimensionality of the phase diagram. We showed that the qualitative phase behavior is robust against changes of the maximum step size $\delta$. The steady state reached by the kinetic MC dynamics is fully described by the density, the persistence length, and the maximum step size $\delta$, although a two-dimensional scaling describes the MIPS phase and the melting transitions at high density $\phi$ for small $\delta$. We argued in this direction using the stochastic description of the kinetic MC dynamics, and its formulation in terms of a Langevin dynamics. The common scaling with $\delta$ of the melting transitions and of the MIPS breaks down in the vicinity of the equilibrium phase transition point. This is needed to obtain the one-parameter scaling (for an inverse-power-law potential) in the passive limit.

The detailed phase diagrams found in our present work once more illustrate the rich collective properties in non-equilibrium physics, and the current limits of our understanding of these models and, more generally, of the physics of active matter.

Acknowledgements.

We thank H. Löwen and L. Berthier for helpful discussions. W.K. acknowledges support from the Alexander von Humboldt Foundation.

Appendix A Passive Fokker–Planck equation under diffusive

scaling

In the passive limit, it follows from (16) that $g(\epsilon,\epsilon^{\prime})=1/(2\delta)$. In this case, the Master equation (of eq. (19)) can be recast by using the definition $P_{k}(x)=\int_{-\delta}^{\delta}d\epsilon P_{k}(x,\epsilon)$ and integrating over $\epsilon^{\prime}$ leads to

[TABLE]

Applying the diffusive scaling in eq. (25), it follows

[TABLE]

with $\tilde{P}_{t}(x)=P_{k}(x)$. Expanding in powers of small $\delta$ leads to

[TABLE]

with the definitions

[TABLE]

The expression for $f(x,\epsilon)$ in eq. (6b) can be rewritten as

[TABLE]

with $\Theta(x)$ the Heaviside function. Using this expression in $A(x)$ in eq. (51) and expanding in terms of small $\delta$, it follows

[TABLE]

with the force $F(x)=-\tfrac{\partial}{\partial x}U(x)$. The two cases $F(x)>0$ and $F(x)<0$ lead to the same result for the integral (namely $1/3$), thus giving

[TABLE]

In the same way it follows for $B(x)$ in eq. (51) (the term with $\Theta(..)$ contributes to order $\delta$ and is thus neglected)

[TABLE]

Using this expressions of eqs (52) and (53) in eq. (50) results in the well-known Fokker–Planck equation for a passive particle in a potential

[TABLE]

This corresponds to the standard Langevin description of a passive particle

[TABLE]

The steady-state probability $P(x)\sim\exp\left(-\beta U(x)\right)$, thus the single relevant parameter is $\beta U(x)$.

If the ballistic scaling of eq. (8) was used for the passive case, it would result in

[TABLE]

which clearly does not capture the correct physics.