Phase transitions and phase coexistence: equilibrium systems versus externally driven or active systems - Some perspectives

ABSTRACT

A tutorial introduction to the statistical mechanics of phase transitions and phase coexistence is presented, starting out from equilibrium systems and nonequilibrium steady-state situations in externally driven systems, such as unmixing of sheared binary fluid mixtures, the driven lattice gas model, and the onset of Rayleigh-Bénard convection. Then, some models for phase separation in models for active systems, where particles possess internal motility, are discussed, emphasizing what one can learn by extending analysis methods to study phase transitions in equilibrium systems by computer simulations to active systems. Specific examples will include colloid-polymer mixtures where the colloids are assumed to be active particles, and active Brownian particles. The extent to which concepts familiar from the study of equilibrium systems are still useful will be critically discussed.

KEYWORDS:

• Active Brownian Particles
• critical exponents
• finite size scaling
• Asakura-Oosawa model

Introduction: Phase Transition in Equilibrium

Active matter has become a hot topic of research in recent years, ^[1] and many phenomena have been found which are reminiscent of phase transitions that are familiar from systems in thermal equilibrium (e.g. ^[2–18]). Both experiment (e.g. ^[6]) and simulation studies of various models (e.g. ^{[4,7,8,10,12,13,15]}) have revealed cases of coexistence between different states of active matter, resembling coexistence between different phases in thermodynamic equilibrium. Is this an accidental analogy, or does there exist a deep connection? Similarly, cases occur where the difference in properties between coexisting states of active matter gradually vanishes, when some control parameter of the active system is varied. Such cases typically involve the observation of large (and slow) fluctuations (e.g. ^[2,10,12]) reminiscent of critical fluctuations associated with continuous-phase transitions in thermal equilibrium. This analogy is particularly intriguing because it is known that such critical behavior is universal for broad classes of systems: this means that the ”critical exponents” ^[19–25] and certain dimensionless combinations of prefactors in the power laws that describe the critical singularities (the “critical amplitudes”) are universal constants, for all systems belonging to such a “universality class.”^[23–26] E.g., the critical behavior of anisotropic magnets (and corresponding theoretical models, such as the Ising model), the vapor-liquid critical behavior and unmixing of binary fluids all belong (with respect to their static properties near criticality) to the same class, and details about the interactions between the Ising spins or the fluid molecules are irrelevant for the universal properties.^[19–26] E.g., experimental evidence for fluids near the vapor-liquid critical point and binary (AB) liquids near the critical mixing point is reviewed in .^[27]

Within the framework of equilibrium statistical mechanics, the renormalization group theory ^[20,22,23] has provided general concepts to understand why universal behavior occurs, and to explain which universality classes exist ^[20,22,23] and identify the parameters that distinguish them: dimensionality (d) of space the system is in; dimensionality (n) of the order parameter that characterizes spontaneous breaking of symmetry associated with the transition from the disordered phase (at temperatures $T$ higher than the critical temperature $T_c$, in case of a thermally driven phase transition) to the ordered phase (at $T<T_c$). Finally, also, the range of the interactions between the elementary degrees of freedom are also involved in the ordering matters: all systems with finite range define one class, while for systems with power-law decay of interactions, different classes may arise.^[24] (Effects due to quenched disorder ^[28,29,30] are disregarded here and also quantum transitions are out of our focus).

This general theoretical framework does not provide information on the nonuniversal properties (such as the location of the critical temperature in relation to the detailed interactions between the elementary degrees of freedom, the values of the critical amplitudes, and the behavior away from the immediate vicinity of the critical point). However, equilibrium statistical mechanics also provides a sound basis for methodologies applied in the context of Monte Carlo and Molecular Dynamics simulations ^[31–34]of model systems, which are suitable to obtain accurate information on all these nonuniversal properties associated with critical behavior and yield in favorable cases also very precise estimates of the universal properties as well.^[24,25] Moreover, such computer simulations are also valuable tools to study coexistence between the different phases of a system, associated with first-order phase transitions, ^[35] including the properties of interfaces separating such coexisting phases.^[36,37]

While, so far, only static properties of systems in thermal equilibrium have been addressed here, it is important to recall that also the dynamics of the (small) thermal fluctuations that occur spontaneously in equilibrium can be addressed similarly. However, there is one important distinction: while with respect to equilibrium properties, in the thermodynamic limit when the number $N$ of the elementary degrees of freedom tends to infinity, the different ensembles of statistical mechanics are strictly equivalent, this is not the case with respect to the dynamics of long-wavelength fluctuations. For instance, in the microcanonical ensemble of a system with a considered volume $V$ both internal energy $E$ and particle number $N$ are strictly constant, i.e. “conserved.” Such conservation laws imply continuity equations for the corresponding densities and gradients. In the canonical NVT ensemble, however, the internal energy is not conserved, and in the grand-canonical $\mu VT$ ensemble where the chemical potential $\mu$ is prescribed and the particle number $N$ may fluctuate, no conservation laws apply. As we shall see later, for the study of static properties of equilibrium systems the freedom in the choice of statistical ensemble is very convenient and useful for the application of computer simulation tools to study phase transitions.^[32] No such freedom occurs for systems far from thermal equilibrium, since no analog of concepts such as thermodynamic potentials and Legendre transformations between them has yet been established.

An interesting aspect of the dynamics of fluctuations in equilibrium states near critical points is also that one static universality class splits into several universality classes with respect to dynamic critical behavior.^[38,39] This is not only due to the conservation laws but also derives from the need to consider the appropriate couplings between the coarse-grained variables needed to describe the dynamics of fluctuations, namely the so-called “Poisson-bracket relations” between the order parameter and the conserved densities. For example, a binary (AB) fluid that undergoes an unmixing transition is considered. There are six conserved densities: the energy density, the mass densities for both components, and the three components of the momentum density vector. Small long-wavelength deviations of these quantities from equilibrium are described by linearized hydrodynamics.^[40] The solution of the corresponding hydrodynamic equations yields five modes (four diffusive modes with relaxation rates proportional to $k^2$ where $k$ is the wavenumber of the mode, and a propagating sound wave (with complex frequency due to sound damping)). The Poisson bracket relation between the total momentum and other densities in the system then implies coupling of modes, which is then crucial for a correct description of the singular behavior of transport coefficients in the system.^[38,39,41] An important feature of fluctuations in equilibrium states is that the correlation function between a fluctuation $\delta Q_1(t)=Q_1(t)-⟨Q_1⟩$ of a quantity $Q_1$ and another quantity $Q_2$, $\delta Q_2(t{ }^{\mathrm{\prime }})$, does not depend on the two times $t$, $t{ }^{\mathrm{\prime }}$ separately, but only on their difference $\mathrm{\Delta }t=t{ }^{\mathrm{\prime }}-t$ (this is called “time translation invariance”).

Kinetics of Phase Changes

One obvious extension of the consideration of the dynamics of fluctuations in equilibrium is the kinetics of phase transitions, where one changes a control parameter of the system (e.g. temperature $T$, or magnetic field $H$, in the case of a ferromagnet) by a finite amount, such that in the phase diagram of the system, a phase boundary is crossed. For example, performing a (instantaneous) quench of temperature in a one-component fluid (or a binary fluid mixture, respectively) such that the final state is inside the two-phase coexistence region, while the initial state was in thermal equilibrium in the (macroscopically homogeneous) one-phase region, a transient non-equilibrium relaxation process occurs, where at intermediate stages of the process, the system is far from equilibrium.^[39,42] When the final state point is only slightly inside the coexistence curve that describes the final macroscopic equilibrium phases, that coexist at this temperature in the thermodynamic limit, the transition starts by nucleation and growth of mesoscopic domains of the minority phase.^[35,43–46] However, both a precise understanding of the lifetime of metastable states (that is limited by such nucleation processes) and of the nature of the saddle point in phase space that is crossed in a nucleation event are still incompletely understood. Further problems concern deeper quenches, when the nucleation barrier due to this saddle point is only of the order of a few times the thermal energy, and a gradual transition to spinodal decomposition (i.e., spontaneous growth of long wave density or concentration fluctuations^{[39,42,47,48]}) takes place. While for systems with long-range forces that are described by mean field theory this transition occurs at a well-defined spinodal curve, ^[49] for systems with short range forces the concept of a spinodal lacks precise meaning.^[35,48,50] However, it is worth remarking that often for the description of the late stages of such phase transition processes the concept of local equilibrium is useful. For example, when after the quench into the two-phase region of a mixture, many nuclei of the minority phase have been formed and grown, so that the supersaturation of the binary mixture has strongly decreased and no more new nuclei can be formed at an observable rate, there is a transient equilibrium between the population of (spherical) domains and the supersaturated background phase: the radius of a domain of the average size is related to the supersaturation of the background with a Kelvin-like relation. This is not a stable equilibrium, since domains larger than the average size will grow while domains smaller than the average size will shrink and ultimately disappear, and thus also the average droplet size grows and the associated supersaturation diminishes. As the time t passed after the quench tends to infinity, also the domain size according to this Lifshitz-Slyozov ^[51] mechanism becomes macroscopic, and true thermal equilibrium is ultimately reached. Irrespective of whether phase separation starts via nucleation processes or spinodal decomposition, characteristic linear dimensions at late stages grow according to power laws with time.^[39,42] Of course, depending on the system that is considered (and the dynamical laws that govern the time evolution of the degrees of freedom of interest), many different kinetic processes associated with phase changes from one equilibrium state to another can be considered;^[39,42,47] while analogous processes occur also in active systems (e.g. ^[52]), we shall not consider such processes here further. Such transient processes are particularly difficult to describe exhaustively, since correlations between dynamical variables at time t and $t+\mathrm{\Delta }t$ do not only depend on the time difference $\mathrm{\Delta }t$ (as in full thermal equilibrium) but also on the time t passed after the instantaneous change of the considered control variable.

Phase Transitions in Externally Driven Systems

Here, we shall rather focus on systems that are far from equilibrium but in a steady state that does not change. Considering again fluid binary mixtures, a problem of great importance in materials sciences is phase separation under steady state shear flow.^[39] Is there an extension of the concepts from equilibrium statistical mechanics to such systems? The simplest case to consider is a flow profile of the form $\dot{\gamma }y{\vec{e}}_x$, i.e. the flow direction is taken to be along the x-axis, ${\vec{e}}_x$ denotes the unit vector along the x-axis. The velocity field varies in the y-direction (i.e., the direction of the shear gradient), while the z-axis denotes the so-called vorticity direction. The shear rate $\dot{\gamma }$ characterizes the strength of the resulting deformation of the system. The concentration fluctuations, which in a system without shear would on average be isotropic in space, become now very anisotropic: this deformation is controlled by the so-called Deborah number $D=\dot{\gamma }{\tau }_{\xi },$ with ${\tau }_{\xi }$ being the characteristic time of the critical fluctuations in equilibrium, when the correlation length is $\xi$, and ${\tau }_{\xi }\propto {\xi }^3$.^[38,39] While in the absence of shear, the starting point for the description of critical dynamics is a Langevin equation known as “model $H$”, ^[38]

(1) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}\mathrm{\Psi }=-\mathrm{\nabla }\cdot (\mathrm{\Psi }\vec{v})+L_o{\mathrm{\nabla }}^2\frac{\delta }{\delta \mathrm{\Psi }}(\mathcal{H}/k_{B}T)+\mathrm{\Theta }$(1)

where $\mathrm{\Psi }$ is the order parameter field, $\vec{v}$ the velocity field, $L_o$ a kinetic coefficient, and $\mathrm{\Theta }$ a random force (satisfying the fluctuation-dissipation theorem). In the presence of shear, a term $-\dot{\gamma }y\mathrm{\partial }\mathrm{\Psi }/\mathrm{\partial }x$ needs to be added. While in mean field theory (without shear) the static structure factor is ($q$ is the wavenumber, $r_o>0$ is the distance from criticality in the phase diagram, for a state point in the one-phase region)

(2) $S(\vec{q})=[r_o+q^2{]}^{-1},$(2)

in the presence of shear, the isotropy in reciprocal space, which applies to this well-known Ornstein-Zernike result, gets replaced by an anisotropic and singular expression ^[53]

(3) $S(\vec{q})\approx [r_o+c[(\dot{\gamma }/L_o)|q_x|{]}^{2/5}+q^2{]}^{-1},$(3)

where $c$ is a constant. A consequence of this result is that the “upper critical dimensionality” $d_u$ (for $d>d_u$ mean field critical exponents apply) gets lowered from the standard result ^[20–23] $d_u=4$ to a value $d_u=12/5$. In addition, a nontrivial shift of the critical point due to the shear occurs.

Equations (1(1) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}\mathrm{\Psi }=-\mathrm{\nabla }\cdot (\mathrm{\Psi }\vec{v})+L_o{\mathrm{\nabla }}^2\frac{\delta }{\delta \mathrm{\Psi }}(\mathcal{H}/k_{B}T)+\mathrm{\Theta }$(1) –3(3) $S(\vec{q})\approx [r_o+c[(\dot{\gamma }/L_o)|q_x|{]}^{2/5}+q^2{]}^{-1},$(3) ) already illustrate the generic features common to many nonequilibrium phase transitions: one must base the treatment on a description of the dynamics of the corresponding equilibrium system (Equation 1(1) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}\mathrm{\Psi }=-\mathrm{\nabla }\cdot (\mathrm{\Psi }\vec{v})+L_o{\mathrm{\nabla }}^2\frac{\delta }{\delta \mathrm{\Psi }}(\mathcal{H}/k_{B}T)+\mathrm{\Theta }$(1) ) even if one wishes just to calculate structural properties, such as $S(\vec{q})$, for the nonequilibrium system, and typically the nonequilibrium perturbation (the term $-\dot{\gamma }y\mathrm{\partial }\mathrm{\Psi }/\mathrm{\partial }x$ in the present case) creates an essential anisotropy of the system. Note that this extra term cannot be incorporated into the functional derivative of the effective free energy of the system – there is no counterpart to “thermodynamic potentials” far from equilibrium, at least not in general. We mention, however, that for shear induced phase separation of polymer solutions and polymer blends one can obtain a quasi-thermodynamic theory, by adding a stored elastic energy contained in the sheared system to the standard (Flory-Huggins^[54]) free energy of the system, describing thus the shift of the critical point due to shear.^[55] However, the validity of this approach is doubtful.

There is a rich variety of nonequilibrium phase transitions occurring in systems that are externally driven. One of the simplest systems is the lattice gas model driven by an electric field (the particles carry an electrical charge, but as in the standard lattice gas model, that is equivalent to the Ising ferromagnet, the only interparticle interaction is the nearest-neighbor attraction between the particles).^[56,57] The density of this lattice gas is conserved, and the dynamics is associated by allowing hopping of the particles from their sites to adjacent empty sites on the lattice. Due to the periodic boundary conditions, a steady-state situation (with a constant current of particles along the field direction) is established. While for the equilibrium lattice gas model, the rate for attempted jumps (in the case of infinite temperature, $T\to \mathrm{\infty }$, when the attraction does not matter) is the same for all neighboring sites, say $\mathrm{\Gamma }$, in the presence of the field the rates are largest in the direction of the field $({\mathrm{\Gamma }}_{+})$, smallest against it $({\mathrm{\Gamma }}_{-})$, and in between in transverse direction $(\mathrm{\Gamma })$. In the absence of the field, this model is the Kawasaki spin exchange kinetic Ising model ^[58] suitable to describe (qualitatively) the dynamics of interstitial alloys (in $d=3$ dimensions) or diffusion in physisorbed monolayers at crystal surfaces (in $d=2$ dimensions). The extension including the electric field is sometimes used to model fast ionic conductors (also called superionic conductors).^[59] Often, this model is considered in the limit of a very strong electric field, where no hopping against the field is possible. The model is clearly formulated as appropriate for studies using Monte Carlo simulation methods;^[31,32,33] but while the latter have yielded very precise and reliable results for the lattice gas/Ising model in thermal equilibrium,^[24,25] despite longstanding attempts ^[57,60] the description of the critical behavior of this model has met many difficulties. Note that due to the anisotropy created by the driving field one needs to distinguish two correlation lengths ${\xi }_{||},{\xi }_{\mathrm{\perp }}$ parallel and perpendicular to the field direction (which also diverge with different exponents ${\nu }_{||},{\nu }_{\mathrm{\perp }}$). While for standard critical phenomena in isotropic systems finite size scaling methods ^[61] are well developed and rather straightforward to apply ^[31,32] in systems with anisotropic critical behavior, the linear dimensions $L_{||},L_{\mathrm{\perp }}$ have to be deliberately chosen consistent with this anisotropic critical behavior.^[62,63] While in equilibrium systems, the principle that, in the thermodynamic limit, the different ensembles of statistical thermodynamics (that are related to each other, as is well known, by Legendre transformations) yield equivalent results allows the choice of the ensemble which for finite size scaling methods is most convenient, no such principle exists for systems far from equilibrium! In the equilibrium lattice gas, the choice of the grand canonical ensemble allows to study the phase transition from the disordered phase to the ordered phase without being hampered by interfacial effects associated with the two-phase coexistence: the latter cannot be avoided in nonequilibrium problems, such as the driven lattice gas.

While the driven lattice gas (as well as phase separation in a sheared binary mixture) starts with a system that already has a phase transition in equilibrium, and one then asks the question how is the transition changed when the external driving is present, there are also many cases of genuine non-equilibrium phenomena, not related to any phase transition in equilibrium. Examples include hydrodynamic instabilities such as the onset of the convective Rayleigh-Bénard instability,^[64–66] or “Absorbing Phase Transitions” in models for reaction-diffusion processes.^[67] Here we shall briefly describe the Rayleigh-Bénard problem only. Here, one considers a fluid held between parallel plates, where the upper plate is held at temperature $T$ while the lower plate is held at temperature $T+\mathrm{\Delta }T$, the distance between the (infinitely large) plates is $h$. Due to this temperature gradient $\mathrm{\Delta }T/h$, heat is transported from the lower plate to the upper one, maintained by a power density $\dot{Q}$ of heating, and one is interested in the effective thermal conductivity ${\lambda }_{eff}=\dot{Q}h/\mathrm{\Delta }T$. If no convection is present, only standard heat conduction occurs, ${\lambda }_{eff}=\lambda$, the standard thermal heat conductivity. The instability arises because, in this setup, the colder liquid on top is slightly heavier than the warmer liquid on the bottom, against the action of gravity. The characteristic control variable for this problem is the “Rayleigh number” $R$, defined as

(4) $R=g\alpha h^3\mathrm{\Delta }T/(k{\nu }_{kin})$(4)

with gravitational acceleration $g$, isobaric thermal expansion coefficient $\alpha$, thermal diffusivity $k$ , and kinematic viscosity ${\nu }_{kin}$. It turns out that convection in the form of a roll pattern sets in ^[64–66] when $R=R_c$, while for $R>R_c$ the “Nusselt number” $\mathcal{N}\equiv {\lambda }_{eff}/\lambda$ exceeds unity. Thus, the region $R<R_c$ plays the role of the “disordered phase” of this nonequilibrium transition, while $\mathcal{N}-1$ is a kind of order parameter.

To estimate $R_c$, one uses the standard hydrodynamic equations, but assumes that the temperature dependence of $\alpha$, k, ${\nu }_{kin}$ can be neglected. Denoting by $\vec{U},\tilde{T},$ and $p$ the deviations of the velocity, temperature and pressure profiles across the layer from the profiles in pure heat conduction, these equations are

(5) $\frac{\mathrm{\partial }\vec{U}}{\mathrm{\partial }t}+(\vec{U}\cdot \mathrm{\nabla })\vec{U}=-\mathrm{\nabla }p+{\mathrm{\nabla }}^2\vec{U}+\sqrt{R}\tilde{T}\hat{Z},\mathrm{\nabla }\cdot \vec{U}=0,$(5)

(6) $\frac{{\nu }_{kin}}{k}(\frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }t}+\vec{U}\cdot \mathrm{\nabla }\tilde{T})={\mathrm{\nabla }}^2\tilde{T}+\sqrt{R}(\vec{U}\cdot \hat{Z}),$(6)

$\hat{Z}$ being a unit vector in the direction of the gradient. In Equations (5)(5) $\frac{\mathrm{\partial }\vec{U}}{\mathrm{\partial }t}+(\vec{U}\cdot \mathrm{\nabla })\vec{U}=-\mathrm{\nabla }p+{\mathrm{\nabla }}^2\vec{U}+\sqrt{R}\tilde{T}\hat{Z},\mathrm{\nabla }\cdot \vec{U}=0,$(5) ,(6(6) $\frac{{\nu }_{kin}}{k}(\frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }t}+\vec{U}\cdot \mathrm{\nabla }\tilde{T})={\mathrm{\nabla }}^2\tilde{T}+\sqrt{R}(\vec{U}\cdot \hat{Z}),$(6) ), lengths are measured in units of the layer thickness $h$, times in units of $h^2/{\nu }_{kin}$, temperature in units of $\mathrm{\Delta }T{\nu }_{kin}/(k\sqrt{R})$, so we deal here with macroscopic variables, not a description on molecular scales. One can add random force terms ^[66] to describe statistical fluctuations, but there are no fluctuation–dissipation relations linking them to transport coefficients of the fluid. The critical value $R_c$ is then estimated applying the well-known concept of linear stability analysis ^[68] of Equations 5(5) $\frac{\mathrm{\partial }\vec{U}}{\mathrm{\partial }t}+(\vec{U}\cdot \mathrm{\nabla })\vec{U}=-\mathrm{\nabla }p+{\mathrm{\nabla }}^2\vec{U}+\sqrt{R}\tilde{T}\hat{Z},\mathrm{\nabla }\cdot \vec{U}=0,$(5) ,6(6) $\frac{{\nu }_{kin}}{k}(\frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }t}+\vec{U}\cdot \mathrm{\nabla }\tilde{T})={\mathrm{\nabla }}^2\tilde{T}+\sqrt{R}(\vec{U}\cdot \hat{Z}),$(6) , using boundary conditions $U_Z=0$, $\tilde{T}=0$, $\mathrm{\partial }{U}_x/\mathrm{\partial }z=\mathrm{\partial }{U}_y/\mathrm{\partial }z=0$ at both the lower and upper plates.

By studying the decay rate of modes for various wavevectors $\vec{k}=(k_x,k_y)$, one finds that an instability sets in for $R=R_c=27{\pi }^4/4$, and the convective roll pattern is characterized by a wavenumber $k_c=\pi /\sqrt{2}$. The amplitude and shape of the pattern can only be obtained from a study of the full nonlinear problem.^[66] While the driven lattice gas model concerns a nonequilibrium phase transition where the elementary degrees of freedom (the hopping particles) represent certain atoms of a superionic conductor, for the Rayleigh-Bénard problem, one considers pattern formation on a macroscopic scale. The active systems, which we will consider next, study degrees of freedom also of macroscopic objects (e.g. flocks of birds) or at least mesoscale objects, such as active Brownian particles.

Selected Models for Active Matter Systems

While the constituent particles of standard condensed matter range in size from atoms or molecules to mesoscopic size of a few $\mu m$ (e.g. the colloidal particles in colloid-polymer mixtures where liquid-gas type phase separation in a phase rich in polymers and a phase rich in colloids can be observed ^[69–72]), the “particles” constituting “active matter” range in size from mesoscopic objects (colloids, ^[9,73–75] living bacteria, ^[6,76] etc.) to macroscopic living animals (flocks of birds ^[77] etc.). The distinctive feature of these “particles” in active matter is that they are self-propelled, i.e. they convert free energy into directed motion. In all these systems, energy has to be supplied continuously to maintain a steady state: this can be due to chemical processes, e.g. when a spherical colloidal particle has two different hemispheres (“Janus particle”) so that the chemical reaction occurs on one hemisphere only.^[78] However, this energy can also be supplied by external fields (light, magnetic fields, etc.) or by mechanical shaking of granular particles.^[79] In all such active systems, there is no detailed balance, unlike systems in thermal equilibrium. In the present article, we shall not at all dwell on the mechanisms that create such a motility of particles, but shall only give a qualitative discussion of the consequences of this motility for the behavior of a few popular models.

The most well known of these models is the Vicsek model, ^[2–5] which was originally intended to understand the collective behavior of flocks of birds or other animals. No realistic description at all was aimed at, and thus, particles moving in a two-dimensional plane with periodic boundary conditions were assumed. Each particle is characterized by its position ${\vec{r}}_i$ and velocity ${\vec{v}}_i$, pointing in a direction characterized by an angle ${\vartheta }_i$ relative to the x-axis, at time $t$. $|{\vec{v}}_i|=v_0=const$, and at each time step $\mathrm{\Delta }t$ the direction of the motion is updated according to ^[2,4]

(7) ${\vartheta }_i(t+\mathrm{\Delta }t)=⟨\vartheta (t){⟩}_{s(i)}+\zeta ,$(7)

where the noise $\zeta$ is a random variable with a uniform distribution in the interval $[-\eta /2,\eta /2]$ and $⟨\cdots {⟩}_{s(i)}$ denotes an average over all particles in the local neighborhood of particle $i$ (i.e., all particles within a circle of radius $R$ centered at the position of the i’th particle). The locations of the particles are updated as ^[2,4]

(8) ${\vec{r}}_i(t+\mathrm{\Delta }t)={\vec{r}}_i(t)+{\vec{v}}_i(t)\mathrm{\Delta }t.$(8)

The motivation for the choice of interactions among particles implicit in Eq. 7(7) ${\vartheta }_i(t+\mathrm{\Delta }t)=⟨\vartheta (t){⟩}_{s(i)}+\zeta ,$(7) is the analogy with isotropic ferromagnets, where exchange interactions align the directions of neighboring spins, while here the directions of the neighboring particles are preferentially aligned, depending on the strength $\eta$ of the noise and the density $\rho$ of the particles. The “fluid” described by Equations 7(7) ${\vartheta }_i(t+\mathrm{\Delta }t)=⟨\vartheta (t){⟩}_{s(i)}+\zeta ,$(7) , 8(8) ${\vec{r}}_i(t+\mathrm{\Delta }t)={\vec{r}}_i(t)+{\vec{v}}_i(t)\mathrm{\Delta }t.$(8) is far from equilibrium also on small timescales, since only the density is conserved, while there is no momentum conservation nor Galilean invariance.^[80] One may introduce the analog of an order parameter by considering the magnitude of the average momentum of the N particles, $\varphi \equiv (1/N)|\sum {\vec{v}}_j|$. Depending on the strength of the noise $\eta$, one finds (for $N\to \mathrm{\infty }$) that a steady state develops, with a disordered phase for $\eta >{\eta }_c$, and an ordered phase with $\varphi >0$ for $\eta <{\eta }_c$, even though one starts the system with an initially disordered configuration, $\varphi (t=0)\approx 0$. Note that a system at initially uniform density and $\varphi (t=0)$ develops into a state of non-uniform density, i.e. phase separation into an ordered cluster with large $\varphi$ and a low-density disordered phase occurs (Fig. 1).^[4] While for a XY ferromagnet in $d=2$ dimensions, there is no long-range order,^[81] in this nonequilibrium phase transition, obviously there does occur a spontaneous breaking of a continuous symmetry in $d=2$. This destruction of the spontaneous magnetization in two-dimensional isotropic magnets is due to long wavelength spin wave excitations, which cost little energy ${ϵ}_k$ (${ϵ}_k\propto k^2$ at wavevector $\vec{k}$), and occur with a large statistical weight, proportional to the Boltzmann factor $(\propto exp(-{ϵ}_k/k_{B}T)$, where $T$ is absolute temperature and $k_B$ Boltzmann’s constant). However, for nonequilibrium systems such as described by the Vicsek model, there is no quantity that would correspond to temperature, and the Boltzmann factor there has no meaning. While in equilibrium, critical phenomena are universal, and hence details on small scales are irrelevant and can be eliminated by a suitable coarse-graining, ^[20–23] for active matter systems it is by no means a priori clear which features of a system are irrelevant details and which features are generic. Thus, one should not be misled by superficial analogies!

Figure 1. Snapshot pictures of the time development of a system with $N=4000$ particles in $L\times L$ box with $L=40$ and $v_0=0.01$ at 50 (a) and 3000 (b) steps. Adapted from Czirok et al.^[4] © IOP publishing. Reproduced with permission. All rights reserved

Our second example are active Brownian particles (ABP).^[8,82–86] Recall that passive Brownian particles are ubiquitous in colloidal suspensions, where we may have spherical silica particles with radius R much larger than the size of the fluid molecules in the suspension, and also their mass is orders of magnitude larger than the masses of the fluid molecules.^[71,87] Due to this disparity of scales, the motions of the colloidal particles are orders of magnitude slower than the motions of the molecules: on the time scale of colloidal motion, interactions with molecules can simply be modeled as uncorrelated random collisions, and so one is led to the standard Langevin equation description for the diffusive motion of colloidal particles.^[87,88] Both translational diffusion and rotational diffusion can be considered. Again, the statistical thermodynamics of equilibrium enters on a fundamental level, linking the noise and the friction constant via the fluctuation dissipation theorem.^[88] In the absence of any (effective) attractive interactions between such passive colloidal particles, no phase transitions whatsoever are expected in such a colloidal suspension in equilibrium, apart from entropy-driven crystallization ^[87] at high densities of the colloids (when the typical distance between the particles is of the order of the colloid diameter) and potentially an intermediate hexatic phase in two-dimensional systems.^[89–91]

However, the situation changes fundamentally when we give the Brownian particles a velocity ${\vec{v}}_i$ of fixed absolute value $v_0$, like in the Vicsek model. The direction of this speed decorrelates gradually via rotational diffusion, as in thermal equilibrium, with angular diffusion constant $D_r$. The equation of motion for an active Brownian particle is then postulated as a Langevin equation for the position ${\vec{r}}_i(t)$ of particle $i$ ^[86]

(9) $\frac{d{\vec{r}}_i}{dt}=v_o{\vec{e}}_i(t)+\frac{1}{\tilde{\gamma }}[{\vec{F}}_i(t)+{\vec{\mathrm{\Gamma }}}_i(t)]$(9)

where the first term on the right-hand side describes that the particle is propelled with constant velocity $v_0$ along the direction of its orientation vector ${\vec{e}}_i(t)$. The translational friction coefficient here is denoted by $\tilde{\gamma }$, ${\vec{F}}_i(t)=-\mathrm{\nabla }{U}_i$ is the force acting on the particle due to the total potential $U_i$ acting on the particle, and ${\vec{\mathrm{\Gamma }}}_i(t)$ is a random force as in thermal equilibrium, with zero average $(⟨{\vec{r}}_i(t)⟩\equiv 0$) and satisfying the usual fluctuation-dissipation theorem ^[88]

(10) $⟨{\mathrm{\Gamma }}_{\alpha ,i}(t){\mathrm{\Gamma }}_{\beta ,j}(t{ }^{\mathrm{\prime }})⟩=2\tilde{\gamma }{k}_{B}T{\delta }_{\alpha \beta }{\delta }_{ij}\delta (t-t{ }^{\mathrm{\prime }})$(10)

with $\alpha ,\beta \in (x,y,z)$ in $d=3$ dimensions or $\alpha ,\beta \in (x,y)$ in $d=2$ dimensions. Note that the translational diffusion constant (in equilibrium) then simply is $D_t^0=k_{B}T/\tilde{\gamma }$; henceforth we may use units such that $k_{B}T=1,$ $\tilde{\gamma }=1$ and then $D_t^0$ is unity as well. A crucial point now is that the direction ${\vec{e}}_i(t)$ undergoes rotational diffusion, with a diffusion coefficient $D_r$. The importance of noise for the motion of the active particle is measured by the Péclet number, $P_e\equiv 3v_0/(\sigma D_r)$, $\sigma$ being the diameter of the (spherical) particles. In the context of simulation studies, it is often convenient to work with a strongly repulsive Weeks-Chandler-Andersen pair potential ^[92] modeling the interaction between the particles and then $\sigma$ in $P_e$ gets renormalized by an effective hard sphere or hard disk parameter, respectively.^[16]

Clearly, in the absence of self-propulsion, this model is essentially a fluid of hard disks $(d=2)$ or hard spheres $(d=3)$, and no vapor-liquid-type phase transition can occur. Thus, the occurrence of a vapor-liquid-type critical point and liquid-vapor phase separation in this model for strong self-propulsion has been a surprise (Fig. 2) .^[16]

Figure 2. Simulation results for the phase diagram of active Brownian particles in $d=2$ dimensions, plotting the state of the system in the plane of variables packing fraction $\varphi$ and Péclet number $Pe$. Within the region confined by the binodal (dashed and solid line), two-phase coexistence occurs between a gas-like phase and a liquid like-phase, as indicated. Open circles show the coexisting packing fractions of the gas $({\varphi }_g)$ and liquid $({\varphi }_{\mathrm{\ell }})$, respectively, and also the packing fractions at the rectilinear diameter ${\varphi }_d=({\varphi }_g+{\varphi }_{\mathrm{\ell }})/2$ are indicated (red circles). The location of the critical point (diamond) was estimated by a subsystem finite size scaling analysis. Reprinted figure with permission from J.T. Siebert et al., physical review E 98, 030601 (2018) ^[16] (https://doi.org/10.1103/PhysRevE.98.030601). Copyright 2018 by the American Physical Society

At first sight, Fig 2 looks like an equilibrium-phase diagram, $P_e$ taking the role that inverse temperature plays in equilibrium. However, this analogy is extremely superficial: in a thermal equilibrium system, in the thermodynamic limit, the properties of the system do not depend on the boundary conditions that are chosen, and all (large) subsystems of the system are at rest, there are no large-scale relative motions of subsystems, so choosing (in a theoretical or simulational study) periodic boundary conditions does not affect the findings on the “bulk” properties of the system. This is completely different here: When we choose periodic boundary conditions here, the system develops at high densities an alignment of the particle velocities $\{{\vec{v}}_i(t)\}$ which differ from the orientations of the active forces, ${\vec{v}}_i(t)\ne v_0{\vec{e}}_i(t)$. This phenomenon is shown in Fig 3, and has no analog whatsoever in equilibrium.^[18] It is also noteworthy that the internal structure of the dense phase is rather complex consisting of patches with hexagonal order and bubbles with rather low densities.^[93] The special character of active fluids emerges already in the dilute limit, where pair correlations between active particles display a special anisotropic structure (“depletion wings” ^[94]), unlike equilibrium fluids formed from point like particles, where pair correlations are isotropic, of course. Due to the active motion, a new characteristic length, called “persistence length” $l_p=v_0/D_r$, matters. Experiments on Janus particles have verified some of the theoretical consequences of this description.

Figure 3. (a) Snapshot configuration of the two-dimensional active Brownian particle model, for a large volume fraction $\eta =0.64$, in a square $L\times L$ box with periodic boundary conditions and $D_r=0.2$, $v_0=5.0$. (b) Enlarged portion of the high density cluster. Colors encode the direction of the self-propulsion. (c,d) Windows from the bulk where the velocities of the particles are shown by blue arrows, showing aligned and vortex domains, respectively. So one can see that there is little correlation between these directions for neighboring particles, even if strong correlations between their velocities (shown by arrows) occur. Reprinted figure with permission from L. Caprini et al., physical review letters 124, 078001 (2020) ^[18] (https://doi.org/10.1103/PhysRevLett.124.078001). Copyright 2020 by the American Physical Society

When, for a system in thermal equilibrium, one has a phase diagram of the type in Fig 2, a state within the miscibility gap means two-phase coexistence, with both phases having then the same temperature $T$ and pressure $p$. It is already a nontrivial issue to give these quantities a precise and unique meaning (recall that also for the Rayleigh-Bénard problem), ^[66] we had to deal with a space-and-time-dependent temperature field, but due to the macroscopic character of the problem, the formulation of Equations 5(5) $\frac{\mathrm{\partial }\vec{U}}{\mathrm{\partial }t}+(\vec{U}\cdot \mathrm{\nabla })\vec{U}=-\mathrm{\nabla }p+{\mathrm{\nabla }}^2\vec{U}+\sqrt{R}\tilde{T}\hat{Z},\mathrm{\nabla }\cdot \vec{U}=0,$(5) ,6(6) $\frac{{\nu }_{kin}}{k}(\frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }t}+\vec{U}\cdot \mathrm{\nabla }\tilde{T})={\mathrm{\nabla }}^2\tilde{T}+\sqrt{R}(\vec{U}\cdot \hat{Z}),$(6) rested on a local equilibrium assumption, so there was no problem to define what the temperature field $\tilde{T}(\vec{r},t)$ meant). For the mesoscopic active particles of interest in the present section, we have a well-defined temperature $T$ for the background fluid of the colloid suspension, responsible for the noise terms in the Brownian motion (cf. Equation 10(10) $⟨{\mathrm{\Gamma }}_{\alpha ,i}(t){\mathrm{\Gamma }}_{\beta ,j}(t{ }^{\mathrm{\prime }})⟩=2\tilde{\gamma }{k}_{B}T{\delta }_{\alpha \beta }{\delta }_{ij}\delta (t-t{ }^{\mathrm{\prime }})$(10) ), but this temperature has nothing to do with the effective temperature related to the average translational kinetic energy of the active particles, $T_{eff}=\frac{1}{2}m⟨v_{\alpha }^2⟩$ ($m$ denotes the mass of the active particles). Now it turns out that for the model defined by Equations 9(9) $\frac{d{\vec{r}}_i}{dt}=v_o{\vec{e}}_i(t)+\frac{1}{\tilde{\gamma }}[{\vec{F}}_i(t)+{\vec{\mathrm{\Gamma }}}_i(t)]$(9) 10 $T_{eff}$ indeed is constant in the steady states of the system, but this is not the case for a slightly different model, ^[17] where one considers underdamped rather than overdamped Brownian motion, where an acceleration term $md{\vec{v}}_i(t)/dt$ is added to the friction term $\tilde{\gamma }{\vec{v}}_i(t)$ in Equation 9(9) $\frac{d{\vec{r}}_i}{dt}=v_o{\vec{e}}_i(t)+\frac{1}{\tilde{\gamma }}[{\vec{F}}_i(t)+{\vec{\mathrm{\Gamma }}}_i(t)]$(9) , and also the moment of inertia I needs to be accounted for orientational relaxation.^[17] In this case, particles in the gas phase are typically much faster than in the dense phase. It has been shown that the coexistence of phases is possible which have different kinetic temperatures!

Also, the meaning of the pressure $P$ needs to be carefully discussed.^[86] In equilibrium, for a system enclosed by hard walls, there is an obvious mechanical interpretation of the pressure in terms of the forces exerted by the particles on the walls. Being interested in the simulation of systems with periodic boundary conditions, the standard approach is to apply the Virial theorem, which can be extended to active Brownian particles as $({\rho }_N=N/V=$ density of the particles) ^[86]

(11) $p=\rho k_{B}T+\frac{1}{2dV}\sum _{i=1}^N\sum _{\genfrac{}{}{0ex}{}{j=1}{(j\ne i)}}^N⟨{\vec{F}}_{ij}\cdot ({\vec{r}}_i-{\vec{r}}_j)⟩+p_{swim}.$(11)

Here, the first term represents the ideal gas contribution, and the (standard) second term is due to the forces ${\vec{F}}_{ij}$ acting between particles $i$ and $j$. The so-called “swim pressure” is due to the activity of the particles,

(12) $p_{swim}=\frac{\tilde{\gamma }{v}_0\rho }{2d{D}_r}v,v=\frac{1}{N}\sum _{i=1}^N⟨{\vec{v}}_i\cdot {\vec{e}}_i⟩,$(12)

$v$ being the average velocity of the particles. As an example, the swim pressure is shown to be a function of volume fraction for various Péclet numbers in Fig 4.^[86] A detailed analysis of the total pressure shows that the phase transition from vapor to liquid is accompanied by a pressure discontinuity, ^[86] unlike what happens in equilibrium. Recall that in thermal equilibrium two bulk coexisting phases have the same pressure, and an isotherm of the density versus pressure has a horizontal part (a state in this horizontal part means that both phases are present in the system, their amount being given by the lever rule). Mean field theories in equilibrium rather yield a “van der Waals loop,” but the correct description of phase coexistence is then found by a Maxwell construction. The latter also applies to simulation results for finite systems (which are numerically exact, apart from statistical errors^[32,33]), where a loop occurs as a finite size effect due to interfacial contributions (see, e.g., ^[50]). However, there is no Maxwell construction for active systems, not even within a mean field framework (see, e.g., ^[95] and references therein).

Figure 4. The swim pressure of a three-dimensional system of active Brownian particles plotted vs packing fraction $\varphi$, for various Péclet numbers: Pe = 9.8 (purple), 29.5 (olive), 44.3 (green), 59.0 (blue), and 295.0 (red). The values are normalized by the ideal (bulk) pressure of individual active Brownian particles as at the same Péclet number. From Winkler et al. ^[86] – Published by the royal society of chemistry

While neither the Vicsek model nor the model of active Brownian particles exhibits a liquid-gas phase separation when the motility is shut off, it is also possible to modify models, which show a transition already in equilibrium, by adding motility to them, in the desire to study how the properties of the transition then change. An example of this approach is colloid-polymer mixtures, where already in the passive limit, phase separation may occur in a polymer-rich phase and a colloid-rich phase.^[69–72] In thermal equilibrium, the generic model for such systems is the Asakura-Oosawa model,^[69,71,72] where the colloids are hard spheres of diameter ${\sigma }_c$ and the polymers are soft spheres of diameter ${\sigma }_p$. While polymers can overlap without energy cost, overlap of colloids and polymers is strictly forbidden. This choice of interactions creates the entropic “depletion attraction” ^[71] among the colloids. One wants to understand for which choices of the polymer and colloid packing fractions ${\eta }_p=(\pi {\sigma }_p^3/6){\rho }_p$ (with polymer density ${\rho }_p=N_p/V),{\eta }_c=(\pi {\sigma }_c^3/6){\rho }_c$ (with ${\rho }_c=N_c/V)$ this depletion attraction leads to phase separation into a polymer-rich and a colloid-rich phase.

As is well known, equilibrium phase behavior is discussed most elegantly in terms of intensive rather than extensive thermodynamic variables (or the densities of the latter). Thus, for the standard liquid–vapor transition for a simple fluid, we use temperature and pressure rather than entropy (or enthalpy) and volume. For the colloid-polymer mixture, the intensive variables of interest are the chemical potentials of polymers ${\mu }_p$ and colloids ${\mu }_c$, respectively. It has become common practice ^[72] to use the so-called polymer reservoir packing fraction ${\eta }_p^r=(\pi {\sigma }_p^3/6){\mathrm{\Lambda }}_p^{-3}exp({\mu }_p/k_{B}T)$ instead of ${\mu }_p$ as a control variable. The phase transition then occurs at a transition line in the $({\mu }_c,{\eta }_p^r)$ plane, ending in a critical point. Of course, in equilibrium we can use Legendre transformations to map out the phase diagram also in the $({\eta }_p,{\eta }_c)$-plane or $({\eta }_p^r,{\eta }_c)$ plane, respectively (Fig. 5 ^[96]). The latter choice has the advantage that the tie lines between coexisting vapor-liquid and liquid-like phases must be strictly horizontal, and hence the critical point must be at the minimum of the miscibility gap in Fig. 5 (insert). In contrast, the tie lines in the ${\eta }_p-{\eta }_c$ representation have nontrivial slopes and finding the location of the critical point is a difficult task.

Figure 5. Phase diagram of the AO model for ${\sigma }_p/{\sigma }_c=0.8$ in the plane of variables ${\eta }_c$ and ${\eta }_p$, as well as using ${\eta }_c$ and ${\eta }_p^r$ (insert). Continuous curves show the so-called free volume approximation (which implies a mean-field-type critical behavior) while symbols show Monte Carlo data (which imply a critical behavior compatible with the Ising model universality class). The square indicates the location of the critical point. Reprinted from R.L.C. Vink and J. Horbach, journal of chemical physics, 121, 3253, (2004), ^[96] with the permission of AIP Publishing

In order to extend the AO model to the case of active colloids, a variant of the original model was used, where the hard-core potentials were mimicked by smooth potentials of the WCA type, ^[97] and instead of $U_{pp}(r)\equiv 0$ a soft repulsion was chosen.^[97] The change of the phase diagram of this model in comparison with Fig. 5 is almost negligibly small. For this model, it is then possible to use a Langevin-type equation of motion similar to Equation 9(9) $\frac{d{\vec{r}}_i}{dt}=v_o{\vec{e}}_i(t)+\frac{1}{\tilde{\gamma }}[{\vec{F}}_i(t)+{\vec{\mathrm{\Gamma }}}_i(t)]$(9) ^[10,12]

(13) $m\frac{d^2{\vec{r}}_i(t)}{d{t}^2}=-\tilde{\gamma }{\vec{v}}_i(t)+{\vec{F}}_i(t)+\sqrt{2\tilde{\gamma }{k}_{B}T}{\vec{\mathrm{\Gamma }}}_i(t),$(13)

but including an inertial term. Unlike the active Brownian particles, a force ${\vec{F}}_i(t)$ is chosen that aligns the velocity of the considered particle with the velocities of its neighbors, as in the Vicsek model ^[10,12]

(14) ${\vec{F}}_i=f_A⟨{\vec{v}}_j{⟩}_R/⟨|{\vec{v}}_j|{⟩}_R$(14)

where $f_A$ is an amplitude factor, which is a control variable (typically $f_A=10$), and $⟨\dots {⟩}_R$ denotes an average over all colloids within a sphere of radius $R$ (typically twice the colloid diameter). While a configuration snapshot of phase coexistence in the active system (Fig. 6a) looks like its counterpart of a passive system, and also the average density profiles in the direction perpendicular to the interfaces (Fig. 6b) resemble equilibrium systems, other properties are very different. Indeed, it was found that using Equations 13(13) $m\frac{d^2{\vec{r}}_i(t)}{d{t}^2}=-\tilde{\gamma }{\vec{v}}_i(t)+{\vec{F}}_i(t)+\sqrt{2\tilde{\gamma }{k}_{B}T}{\vec{\mathrm{\Gamma }}}_i(t),$(13) , 14(14) ${\vec{F}}_i=f_A⟨{\vec{v}}_j{⟩}_R/⟨|{\vec{v}}_j|{⟩}_R$(14) , the miscibility gap became enlarged, in comparison with the passive model.^[10,12] In addition, the model exhibits very pronounced nonequilibrium effects, since the velocities of the particles exhibit huge fluctuations both in sign and in magnitude. In the regime where phase coexistence occurs, the clustered active particles travel mostly in the same direction, and the kinetic temperature in the colloid-rich phase deviates from the temperature chosen for the thermostat by more than an order of magnitude (Fig. 6) .^[10]

Figure 6. (a) Snapshot picture of the active polymer-colloid mixture at the state point ${\eta }_c=0.15$, ${\eta }_p=0.20$ with $f_A=10$, for a box with linear dimensions $L_x\times L_y=12$, $L_z=48$. Colloids are shown in yellow, polymers in black. (b) Profiles of ${\eta }_p$ and ${\eta }_c$ in z-direction. (c) Corresponding profiles for the various components of T(z) as obtained from the equipartition theorem. Reprinted from B. Trefz et al., journal of chemical physics, 144, 144902, (2016), ^[10] with the permission of AIP publishing

The AO model ^[69] already in equilibrium is very useful to understand another important aspect of phase behavior, namely, the competition between the liquid–vapor transition and the liquid-solid transition of fluids.^[70–72] If the attractive forces (in equilibrium) between the particles are completely absent, we still have (for spherical particles of diameter $\sigma$) a liquid – solid transition. This transition in $d=3$ is of first order; we have a two-phase coexistence region from a packing fraction $\eta =(\frac{\pi }{6}{\sigma }^3)\cdot \rho =0.494$ to $0.545$, where a face-centered cubic solid takes over. Attractive interactions (due to flexible polymer coils) are characterized by two basic parameters, their strength and their range ${\sigma }_p$. Only when ${\sigma }_p$ is large enough can we observe a vapor-liquid-like transition: for small ${\sigma }_p$ one rather observes that the boundary between the fluid and the two-phase region with increasing interaction strength bends over to rather small $\eta$. For intermediate attraction strength, a metastable vapor, liquid-like transition is found (in mean-field treatments) within the fluid-solid two-phase region. Interestingly, there is also evidence that the $d=3$ ABP model belongs to this scenario, with a metastable phase boundary separating the vapor-like phase from the vapor-liquid two-phase region INSIDE of the vapor-crystal two-phase coexistence region .^[98] In $d=2$, however, for small Peclet numbers rather a coexistence region between an active fluid and a hexatic phase, followed by a pure hexatic phase, is found.^[99] Note, however, that the existence of hexatic phases is due primarily to entropic effects: there seems to exist evidence ^[100] that in the presence of strong enough attractive forces the fluid-solid transition in $d=2$ is of first order, as in $d=3$ dimensions.

How Can We Learn from Simulation Studies about the Critical Behavior of Active Systems?

When strong phase separation occurs, as shown in Fig. 6a,b extracting the densities ${\eta }_p,{\eta }_c$ in either of the coexisting phases can be done rather straightforwardly, directly from the average profiles (6b). However, near the critical point, this naive analysis breaks down, since the width of the interfaces increases and strong fluctuations in the remaining bulk phases make quantitatively reliable direct estimations almost impossible.

This problem would arise also if we study the equilibrium phase transitions by simulations and analyze the data in this straightforward way. However, for equilibrium phase transitions, it is well established that the most efficient approach to analyze simulations is finite size scaling.^{[23,31,32,61,101]} The simplest case is the Ising model of ferromagnets, where at zero magnetic field $H=0$ upon lowering the temperature $T$ a phase transition occurs, in the thermodynamic limit, from the paramagnetic state at $T>T_c$ to the ferromagnetic state at $T<T_c$, where states of nonzero spontaneous magnetization $(\pm m_s$; we normalize the magnetization per spin) occur. Of course, in a finite system there cannot be any spontaneous symmetry breaking: in equilibrium the distribution $P_L(m)$ of the magnetization per spin, for a d-dimensional lattice of linear dimension $L$, in thermal equilibrium is always symmetric with respect to the sign of $m$. To a good approximation, $P_L(m)$ for $T\gg T_c$ is simply a Gaussian centered at $m=0$, as one can justify simply from the law of large numbers; for $T\ll T_c$ one finds an appearance of two peaks, centered near $\pm m_s$. As one approaches $T_c$ from above, the width of the distribution broadens and its form becomes distinctly nongaussian, when the correlation length of spin correlations $\xi$ (which diverges at $T_c$ for $L\to \mathrm{\infty }$) becomes comparable to the chosen $L$. Near $T_c$ a gradual change of the distribution from single peak form to double peak form occurs. Thus, as a function of $T$ both $P_L(m)$ and all its moments $⟨m^2{⟩}_L,⟨m^4{⟩}_L\dots$ are nonsingular functions of $T$. Now, finite size scaling ^[61,101] states that near $T_c$ and for large enough $L$, $P_L(m)$ depends on the three variables $L,m$ and $T$ not separately but only on two scaled variables, $L/\xi$ and $m/m_s$. Since $\xi \propto (T/T_c-1{)}^{-\nu }$, $m_s\propto (1-T/T_c{)}^{\beta }$, where $\nu$ and $\beta$ are the critical exponents of the correlation length and order parameter, one can choose also $m{\xi }^{\beta /\nu }$ instead of $m/m_s$ as a scaling variable. Normalization of the probability then requires that $P_L(m)$ can be written as $(ϵ\equiv 1-T/T_c$ here) ^{[31,32,61,101]}

(15) $P_L(m)=L^{\beta /\nu }\hspace{0.17em}\tilde{P}(ϵL^{1/\nu },m{L}^{\beta /\nu }),⟨m^k{⟩}_L=L^{-k\beta /\nu }{\tilde{M}}_k(ϵL^{1/\nu })$(15)

where $k=2,4,\dots$. Choosing the arguments of the scaling functions $\tilde{P}$ and ${\tilde{M}}_k$ in this particular way is advantageous because it brings out that $\tilde{P}$ and ${\tilde{M}}_k$ are nonsingular functions of their variables, since in a finite system $P_L(m)$ and $⟨m^k{⟩}_L$ are not singular as function of $T$ and hence singularities only arise when the limit $L\to \mathrm{\infty }$ is taken before $ϵ\to 0$.

From Equation 15(15) $P_L(m)=L^{\beta /\nu }\hspace{0.17em}\tilde{P}(ϵL^{1/\nu },m{L}^{\beta /\nu }),⟨m^k{⟩}_L=L^{-k\beta /\nu }{\tilde{M}}_k(ϵL^{1/\nu })$(15) the well-known cumulant intersection method ^[101] emerges, i.e., $U_L$ is a function of $ϵL^{1/\nu }$ only and hence for $ϵ=0$ all cumulants for different choices of $L$ must intersect at a common point $U^{\ast }$,

(16) $U_L=1-⟨m^4{⟩}_L/(3⟨m^2{⟩}_L^2)=\tilde{U}(ϵL^{1/\nu }),U^{\ast }=\tilde{U}(0).$(16)

The location of this common intersection point hence yields $T_c$, with no need to know anything else in beforehand. At $T_c$, the second moment varies as $⟨m^2{⟩}_L\propto L^{-2\beta /\nu }$, hence allows to estimate the ratio $\beta /\nu$ of the critical exponents. Since $\tilde{U}$ is a continuous function of $T$, we must have $\tilde{U}(ϵL^{1/\nu })\approx U^{\ast }+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}\cdot ϵL^{1/\nu }+\dots$ near $ϵ=0$, studying the L-dependence of the slope of $U_L$ at $U^{\ast }$ hence yields $1/\nu$. Since finite size scaling implies ^[101] also the validity of the “hyperscaling” relation ^[21–23] $d\nu =\gamma +2\beta$, the susceptibility exponent $\gamma$ follows (as well as the other exponents via scaling laws^[19–21]).

This account of finite size scaling in a nutshell is not the whole story, of course: (i) Equations 15(15) $P_L(m)=L^{\beta /\nu }\hspace{0.17em}\tilde{P}(ϵL^{1/\nu },m{L}^{\beta /\nu }),⟨m^k{⟩}_L=L^{-k\beta /\nu }{\tilde{M}}_k(ϵL^{1/\nu })$(15) ,16(16) $U_L=1-⟨m^4{⟩}_L/(3⟨m^2{⟩}_L^2)=\tilde{U}(ϵL^{1/\nu }),U^{\ast }=\tilde{U}(0).$(16) hold only strictly in the asymptotic limit $L\to \mathrm{\infty }$, for practically useful values of $L$ corrections to finite size scaling need careful consideration.^[24,25,32] (ii) If one would apply a local Monte Carlo algorithm ^[31,32] e.g. single spin flips, the accuracy of the sampling is strongly hampered by critical slowing down,^[38,39] i.e. the sequence of generated system configuration exhibits strong “dynamic” correlations with a relaxation time $\tau \propto {\xi }^z$, where (in the absence of any conservation laws) the dynamic exponent $z$ is ^[32,38] $z\approx 2$. Accurate results for large $L$ could be obtained ^[24,25] thanks to the invention of cluster algorithms ^[32] which strongly reduce critical slowing down.

However, obtaining accurate information on the critical behavior of models for fluids exhibiting a vapor-liquid transition has been much more difficult and possible only at the beginning of the present century, ^[102] due to several reasons: (i) for off-lattice systems cluster algorithms are much less powerful than for the Ising model (ii) While in the Ising model there is the special symmetry that the internal energy is invariant against spin reversal while the order parameter changes its sign, these two quantities are in a sense orthogonal to each other, unlike fluids where no special symmetry exists, and energy density and particle number density are “coupled operators” (“field mixing” effect ^[103]). Therefore, it is a nontrivial task not only to locate the critical temperature, but also to locate the critical density (which in the lattice gas fluid, that corresponds to the Ising model, by symmetry is ${\rho }_c=1/2$) and the chemical potential ${\mu }_c$ and pressure $p_c$ at criticality. A consistent analysis ^[102] needs to consider the coupling between the three “scaling fields” $ϵ,\mu$ and $p$ in order to correctly extract critical properties.

Still, the task is manageable since in equilibrium one can perform the simulation in the grand canonical ensemble where $\mu$ is used as a control variable, and moves involving the insertion and deletion of particles are carried out, and thus the density distribution $P_L(\rho )$, which is the analog of Equation 15(15) $P_L(m)=L^{\beta /\nu }\hspace{0.17em}\tilde{P}(ϵL^{1/\nu },m{L}^{\beta /\nu }),⟨m^k{⟩}_L=L^{-k\beta /\nu }{\tilde{M}}_k(ϵL^{1/\nu })$(15) , can be sampled. However, for systems of active particles, no analog of the grand-canonical ensemble exists; rather one is forced to work in the canonical ensemble where the average density of active particles in the system is strictly constant; and it is just a heuristic speculation to assume that the finite size scaling concept can be transferred to phase transitions in active matter!

Of course, also in equilibrium systems one can try an extension of the finite size scaling method to the canonical ensemble where $\rho =N/V=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}$. This can be done by considering sub-systems of the considered total system: exchange of particles through the subsystem boundary into the neighborhood is then possible, and for each subsystem, the rest of the system acts like a reservoir of particles.

This approach was tried by Siebert et al., ^[16] using $\mathrm{\ell }\times \mathrm{\ell }$ subsystems of a simulation box which has extensions $6\mathrm{\ell }\times 2\mathrm{\ell }$. When the system is in the two-phase coexistence region, the interfaces between the coexisting phases of low density and high density will be oriented perpendicular to the long direction (the x-direction). In each configuration that is analyzed one can then locate these interfaces (in the region where phase separation occurs) and record the properties of subsystems that are located in the “pure” phases, to avoid artifacts due to “contamination” of the bulk properties with interfacial properties (Fig. 7a). If one includes subsystems that contain the interfaces, one does not obtain useful results as one can verify already for the Ising model.^[16,104] With the present variant of the method, the critical temperature of the Ising model can be located with reasonable accuracy (Fig. 7b), and also good estimates of the critical exponents could be extracted .^[16] Unfortunately, for the ABP system, the observed cumulant crossing is less impressive (Fig. 7c). Of course, for the Ising system using a total magnetization of zero trivially ensures that at $T_c$ one reaches critical conditions; the choice of density for the ABP model requires already some estimate for the critical density, which is unknown beforehand, and an inaccurate choice clearly may be the source of an uncontrolled error. Neither in this Ising model with conserved magnetization nor in the ABP model could the region very close to criticality be reliably probed (unlike studies of the Ising model in the grand canonical ensemble). This fact is due to the difficulty that in the Ising model with conserved order parameter the dynamic exponent $z=2+\gamma /\nu$ is much larger than in the non-conserved case,^[38] critical slowing down hence is a much more serious problem for simulations. So while Siebert et al.^[16] found that the exponents differ significantly from those of the Ising model, one must make the caveat that the range of $\mathrm{\ell }$ that could be explored is only $\mathrm{\ell }=10$ to $\mathrm{\ell }=17.5$, i. e. it spans less than a factor 2, while for accurate estimates of equilibrium systems, one needs about a decade.^[24,25,32] Also, a direct study of the order parameter as a function of the distance $\xi$ from criticality $(ϵ=P_{ec}/P_e-1)$ here could only probe the range $0.3\le ϵ\le 0.7$, while in equilibrium the proper exponents are seen for $ϵ\ll 0.3$ only!

Figure 7. (a) Schematic representation of the simulation box used for the cumulant analysis. Simulations are done at medium packing fraction of particles with an edge length ratio of 3:1 and periodic boundary conditions. This results in a two-phase coexistence with a slab geometry of the liquid domain. The slab is always aligned with the short axis. Two $\mathrm{\ell }\times \mathrm{\ell }$ subboxes are then placed in the liquid domain, such that the x-coordinates of their center of mass coincides with the x-coordinate of the liquid domain. The other two subboxes are shifted by $3\mathrm{\ell }$ along the x-axis, so they are centered in the vapor region. (b) Cumulant ratio $Q_{\mathrm{\ell }}=⟨m^2{⟩}_{\mathrm{\ell }}^2/⟨m^4{⟩}_{\mathrm{\ell }}$ for several choices of $\mathrm{\ell }$, namely $\mathrm{\ell }=8,10,12$ and 15 lattice spacings, for the Ising/lattice gas model in $d=2$ dimensions, plotted vs. temperature $T$ (in units of the exchange constant $J/k_B$). (c) Cumulants for the two-dimensional active Brownian particle model plotted vs. the Péclet numbers for four choices of $\mathrm{\ell }$ (in units of $\sigma$) as indicated. A crossing of $Q_{\mathrm{\ell }}(Pe)$ can be seen for roughly $Pe=40\pm 2$. The dashed lines are guides to the eye only. Reprinted figure with permission from J.T. Siebert et al., physical review E 98, 030601 (2018) ^[16] (https://doi.org/10.1103/PhysRevE.98.030601). Copyright 2018 by the American Physical Society

Thus, it is no surprise that doubts on the validity of the results of Siebert et al.^[16] have been raised.^[105] In particular, Partridge and Lee ^[105] proposed a variant of APB on a triangular lattice, where particles with an orientation pointing in the lattice direction can either undergo rotational diffusion or jump (with preference in the direction of their orientation) to a neighboring site, provided this site is empty. Applying the same type of finite size scaling analysis as Siebert et al. ^[16] for $18\le L\le 36$, evidence for Ising-like behavior was found, while indications are somewhat less clear for another active lattice model ^[106] as investigated, recently, in Ref.^[107] 2d-Ising behavior has also been reported of late for off-lattice active Ornstein-Uhlenbeck particles.^[108]

In our view, these discrepancies between ^[16] and ^[105] do not imply that the described finite size scaling analysis fails in principle; rather, the difficulties are reminiscent of the early days of Monte Carlo studies of equilibrium critical behavior, ^[31,32,101] where due to the use of too small system sizes and too short simulation runs the desired accuracy of estimates could not be reached. Thus, it is hoped that future work can resolve these problems.

Also, in the active version of the AO model (polymers + active colloids) ^[10,12] the situation allows to draw rather tentative conclusions only: while again a subsystem finite size scaling approach is the method of choice, one needs to analyze the distribution of two densities, $P_{\mathrm{\ell }}({\eta }_p,{\eta }_c)$. Using subsystems in the disordered phase near criticality gave results ^[12] compatible with the Ising universality class, but (in $d=3$ dimensions) only rather small values of $\mathrm{\ell }$ were accessible $(3\le \mathrm{\ell }\le 6)$, and the location of the critical point could only be very roughly estimated^[12]). Thus, while in principle simulations should be able to clarify the critical behavior of phase transitions in active matter systems, clearly more work is still necessary.

Conclusions

Active matter systems undergoing phase transitions exhibit a wealth of interesting phenomena, with similarities to systems undergoing phase transitions in thermal equilibrium, or in systems under external drive creating a steady flow of heat or particles through the system, but there also occur very interesting differences. We have shown that phenomena such as coherent spontaneous large-scale spontaneous flow patterns of particles can emerge (e.g. Fig. 3), the effective temperatures can differ between coexisting phases (Fig. 6), and so on. An interesting feature that we have not yet addressed is the fact that the interface tension between coexisting phases, when defined from the anisotropy of the pressure tensor, can become negative.^[14,37,109] In equilibrium fluid systems, the interfacial tension equals the interfacial stiffness, which controls the capillary wave broadening of the interface (and must be positive, while the interface would be unstable against roughening). Bialké et al.^[109] did find a positive interfacial stiffness together with a negative interfacial tension: this result again shows that many relations one knows to hold in equilibrium do not carry over to nonequilibrium systems. Note that in thermal equilibrium the interfacial free energy can be defined mechanically (interfacial stress is defined as an excess term computed from the stress tensor^[110]) or thermodynamically (e.g. considering the variation of the free energy with interfacial area^[110]), and both definitions must yield the same result. This agreement is lost in active systems (see, e.g. ^[95]), where the mechanical route (which is still well defined) yields a negative interfacial tension. Hence, alternative definitions of interfacial tension are still a matter of debate (e.g. ^{[93, 111, 112]}), since a uniquely defined free energy does not exist for active systems.

A description related to these issues is the so-called “active model B+ “ and the “cluster phases” following from it.^[113] In the classification of dynamic critical phenomena, ^[40] “model B” describes critical dynamics of binary mixtures where only the concentration is conserved. Ref.^[113] suggests that an active version of this model is obtained by adding suitable terms to the current that break the local time reversal symmetry. The strength of these terms is described by two coupling parameters $\lambda$ and $\zeta$, and in the space of these two parameters, regions can be observed where “reverse Ostwald ripening” for droplets or bubbles (but never for both of them) occurs, in addition to the standard (Lifshitz-Slyozov^[51]) Ostwald ripening (i.e., droplets coarsen with time), that is observed near equilibrium. Tjhung et al.^[113] suggest that in this way one can obtain a generic phase diagram for all variants of “motility induced phase separation (MIPS)” as a generalization of the standard theory of phase separation kinetics leading to thermal equilibrium.

There are very interesting attempts to develop a theoretical understanding of phase transition in active matter systems by formulating suitable analytical theories, extending (e.g. ^[9,14]) theories on the dynamics of critical phenomena ^[38,39] or of the kinetics of transitions caused by quenching experiments.^[35,42,48] In this tutorial perspective article, we have not attempted to describe them for several reasons: (i) lack of expertise of the authors; (ii) basic assumptions are sometimes not consistent with features of the specific models that we describe in Sec. II: e.g., in ^[14] it is assured that the vectorial degrees of freedom corresponding to the orientations of the active particles are fast degrees of freedom and hence there is the dynamics sufficiently described by the continuity equation for the density $d\rho /dt=-\mathrm{\nabla }\cdot \vec{J}$, and the current vanishes in homogeneous phases.^[14] We feel, however, that in such a description interesting phenomena (such as seen in Fig. 3) are disregarded from the outset. (iii) While the extensions of Cahn-Hilliard theory ^[35,48] to transitions in active matter clearly render a mathematically very elegant description, ^[14] we recall that for equilibrium systems many predictions of Cahn-Hilliard theory (initial exponential growth of unstable fluctuations, critical-like singularities at the spinodal, etc.) are misleading ^[35,50] apart from fluids with extremely weak and extremely long range attractive forces.^[47,49] We have also not described attempts to obtain an approximate mapping of active systems to passive ones ^[10] where integral equation theories are used to derive effective forces between the active particles from the observed pair correlations; we again expect that such approaches are of limited usefulness. A very interesting subject that we could not cover here are phase transitions of active nematics ^[114] and active polymer networks, ^[115] which may play a role in biological contexts.

Thus, many open questions exist in the field; we hope that the present article helps to formulate such questions more precisely, and stimulates to develop further models and concepts for their analysis.

Acknowledgments

The authors had the privilege of a stimulating collaboration with S. K. Das, F. Dittrich, S. A. Egorov, F. Schmid, J. T. Siebert, T. Speck and B. Trefz in some preceding original work [10, 12, 16]. It is a pleasure to thank them for many insightful discussions, and for permitting them to reproduce some figures. PV gratefully acknowledges the financial support by the Deutsche Forschungsgemeinschaft within priority program SPP 1726 (Grant VI 237/5-2). The authors gratefully acknowledge the computing time granted on the supercomputer Mogon at Johannes Gutenberg University Mainz (hpc.uni-mainz.de).

Additional information

Funding

This work was supported by the Deutsche Forschungsgemeinschaft [VI 237/5-2].