Phase coexistence of active Brownian particles

We investigate motility-induced phase separation of active Brownian particles, which are modeled as purely repulsive spheres that move due to a constant swim force with freely diffusing orientation. We develop on the basis of power functional concepts an analytical theory for nonequilibrium phase coexistence and interfacial structure. Theoretical predictions are validated against Brownian dynamics computer simulations. We show that the internal one-body force field has four nonequilibrium contributions: (i) isotropic drag and (ii) interfacial drag forces against the forward motion, (iii) a superadiabatic spherical pressure gradient, and (iv) the quiet life gradient force. The intrinsic spherical pressure is balanced by the swim pressure, which arises from the polarization of the free interface. The quiet life force opposes the adiabatic force, which is due to the inhomogeneous density distribution. The balance of quiet life and adiabatic forces determines bulk coexistence via equality of two bulk state functions, which are independent of interfacial contributions. The internal force fields are kinematic functionals which depend on density and current but are independent of external and swim forces, consistent with power functional theory. The phase transition originates from nonequilibrium repulsion, with the agile gas being more repulsive than the quiet liquid.

Figure 1

Overview of phase-separated active Brownian particles. (a) Snapshot from BD simulations. The particles (blue) separate into regions with high and with low density. Due to periodic boundary conditions in $x$ and in $y$, two interfaces form along the short dimension of the simulation box with size $120\sigma$ in $x$ and $24\sigma$ in $y$. (b) Particles with size $\sigma$ at position $\mathbf{r}=(x,y)$ are driven in the (unit vector) direction $\mathit{\omega }$ by the swim force with strength $\gamma s$. The particles interact with the (WCA) repulsion $\varphi \left(r\right)$ with energy scale $\epsilon$. (c) Illustration of the different one-body force density contributions that add up to the total internal force density field ${\mathbf{F}}_{\mathrm{int}}$: the drag force density ${\mathbf{F}}_{\sup ,0}$ acts against the local current direction $\mathbf{J}$, and its magnitude is small (large) in the dilute (dense) phase. The nonspherical drag correction ${\mathbf{F}}_{\sup ,1}$ occurs at the interface, and it acts in direction ${\mathit{\omega }}^{*}$ ($\mathit{\omega }$ mirrored at the $x$ axis). The swim pressure (orange arrow) is due to the polarization of the interface, and it is balanced by the superadiabatic pressure (red arrow), which is low (high) in the dilute (dense) phase. The arrows indicate the direction of the respective negative pressure gradient. The quiet life force field ${\mathbf{F}}_{\sup ,3}/\rho$ compresses the liquid and acts against the adiabatic force field ${\mathbf{F}}_{\mathrm{ad}}/\rho$, which is (solely) due to the density gradient and tends to expand the liquid. (d) Mean scaled forward swimming speed $v_{\mathrm{f}}/s$ as a function of the scaled position $x/\sigma$ across the interface in a phase-separated system of active Brownian particles interacting with the WCA pair potential, with particle size $\sigma$ and energy scale $\epsilon$. Simulation data are shown for $k_{B}T/\epsilon =0.5$ and $s\tau /\sigma =60$, where the timescale is $\tau ={\sigma }^2\gamma /\epsilon$, which corresponds to $\mathrm{Pe}=120$. The aspect ratio of the simulation box is 5, and the number $N$ of particles per system volume $V$ is $N/V=0.7{\sigma }^{-2}$ with $N=2000$; the time step is $\mathrm{\Delta }t/\tau ={10}^{-5}$. Sampling was performed over ${10}^8$ time steps; see Ref. [38] for further simulation details that also apply to the present study. The inset shows the theoretical result (38). The inset axis labels have been omitted for clarity; the scale is identical to that of the main plot.

Figure 2

Nonequilibrium phase diagram as a function of scaled density ${\rho }_{\mathrm{b}}{\sigma }^2$ and Péclet number $\mathrm{Pe}$. Shown are the nonequilibrium binodal (solid line) and spinodal (dashed line) obtained from the present theory, compared to results for the binodal (orange squares) taken from Ref. [22], and for the gas side of the spinodal (blue circles) taken from Ref. [40]. Also shown is the theoretical result for the critical point (open circle).

Figure 3

Representative Fourier coefficients obtained from BD simulation (left panels) and theory (right panels). Shown are Fourier components of order $n=0,1,2$ (as indicated) as a function of $x/\sigma$ for (a) the density ${\rho }_n{\sigma }^2$, (b) the $x$ component of the current, $J_n^x\sigma \tau$, and (c) the $y$ component of the current, $J_n^y\sigma \tau$, where the “molecular” timescale is $\tau =\gamma {\sigma }^2/\epsilon$. The simulation parameters are identical to those of Fig. 1.

Figure 4

Spherical drag force density ${\mathbf{F}}_{\sup ,0}$ and nonspherical drag force density ${\mathbf{F}}_{\sup ,1}$, obtained from (a) simulations via correlators, (b) via kinematic functionals using simulation data as input, and (c) stand-alone theory. Shown are the $x$ and $y$ components of the respective force density field: $F_{\sup ,0}^{\left(x\right)}$ (first column) $F_{\sup ,0}^{\left(y\right)}$ (second column), $F_{\sup ,1}^{\left(x\right)}$ (third column), and $F_{\sup ,1}^{\left(y\right)}$ (fourth column) in units of $\epsilon /{\sigma }^3$ and as a function of distance $x/\sigma$ across the interface and angle $\phi$ with respect to the interface normal (pointing towards the liquid). The plus and minus signs indicate the sign of the force density fields.

Figure 5

(a) Superadiabatic spherical pressure profile ${\mathrm{\Pi }}_2\left(x\right)$, swim pressure profile $P_{\mathrm{swim}}\left(x\right)$, and the sum ${\mathrm{\Pi }}_2\left(x\right)+P_{\mathrm{swim}}$, in units of $\epsilon /{\sigma }^2$ as a function of distance $x/\sigma$ across the interface. Results are obtained from the correlator expressions using simulation data as input (blue dotted lines), from the kinematic functionals with simulation data input (black dashed lines) and from the stand-alone theory (red solid lines). (b) Superadiabatic spherical chemical potential profile ${\nu }_3\left(x\right)$, the sum of adiabatic excess and ideal chemical potential profile ${\mu }_{\mathrm{id}}\left(x\right)+{\mu }_{\mathrm{ad}}\left(x\right)$, and the total chemical potential, i.e., the sum ${\nu }_3+{\mu }_{\mathrm{ad}}+{\mu }_{\mathrm{id}}$ in units of $\epsilon$ and as a function of $x/\sigma$. Here the ideal chemical potential is ${\mu }_{\mathrm{id}}=k_{B}Tln\left({\eta }^{\prime }\right)$ with the local rescaled packing fraction ${\eta }^{\prime }=0.8{\rho }_{\mathrm{b}}/{\rho }_{\mathrm{jam}}$. For clarity the results for the sum have been shifted upwards by two units. The superadiabatic results are obtained using the kinematic functional with either the simulated (dashed black lines) or stand-alone (solid red lines) Fourier coefficients.

Figure 6

Isotropic Fourier component of the scaled density profile, ${\rho }_0{\sigma }^2$, as a function of $x/\sigma$ obtained from BD simulations for $s\tau /\sigma =60$ and $k_{B}T/\epsilon =0.5$. (a) For different values of the average density $N{\sigma }^2/V=0.5,0.7$, and 0.9 (as indicated) with simulation box aspect ratio $\mathcal{A}=5$. (b) For different aspect ratios $\mathcal{A}=2.5,5,$ and 10, with average density $N{\sigma }^2/V=0.7$.

Figure 7

Fourier coefficients ${\rho }_n{\sigma }^2$ of the density distribution as a function of $x/\sigma$ across the interface. Shown are the theoretical results for $n=0,1,2$ (dashed lines, as indicated) for different values of the parameter $\xi =0$ (a) and 2100 (b). As a reference the results for $\xi =700$ are also shown (solid black lines); these results are identical to the data in the right panel of Fig. 3.