Interparticle torques suppress motility-induced phase separation for rodlike particles

TL;DR

This study investigates how interparticle torques influence motility-induced phase separation (MIPS) in self-propelled particles, revealing that elongated shapes suppress MIPS through collision dynamics, with implications for understanding active matter behavior.

Contribution

The paper introduces a comprehensive simulation and analysis of MIPS suppression in rodlike particles, highlighting the role of particle shape and collision duration in phase behavior.

Findings

Elongated particles suppress MIPS compared to isotropic particles.

A collision duration-based mechanism explains MIPS suppression in rods.

The suppression mechanism also accounts for MIPS enhancement in Vicsek-like particles.

Abstract

To study the role of torque in motility-induced phase separation (MIPS), we simulate a system of self-propelled particles whose shape varies smoothly from isotropic (disks/spheres) to weakly elongated (rods). We construct the phase diagrams of 2D active disks, 3D active spheres and 2D/3D active rods of aspect ratio $l/\sigma=2$. A stability analysis of the homogeneous isotropic phase allows us to predict the onset of MIPS based on the effective swimming speed and rotational diffusion of the particles. Both methods find suppression of MIPS as the particle shape is elongated. We propose a suppression mechanism based on the duration of collisions, and argue that this mechanism can explain both the suppression of MIPS found here for rodlike particles and the enhancement of MIPS found for particles with Vicsek interactions.

Equations

$${\partial }_t{\mathbf{r}}_i$$

$${\partial }_t{\theta }_i$$

$$⟨{\Lambda }_i⟩$$

$$⟨{\Lambda }_i^{\alpha }(t){\Lambda }_j^{\beta }(t^{\prime })⟩$$

$$u_{WCA}(r_{s,ij})=$$

$$\{\begin{array}{ll}4ϵ\left[{\left(\frac{\sigma }{r_{s,ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{s,ij}}\right)}^6\right]+ϵ& \mathrm{if}{r}_{s,ij}<2^{1/6}\sigma ;\\ 0& \mathrm{if}{r}_{s,ij}\ge 2^{1/6}\sigma .\end{array}$$

$${\partial }_t{\mathbf{r}}_i$$

$${\partial }_t{\hat{\mathbf{e}}}_i$$

$${\partial }_t\delta \rho (\mathbf{r},t)\approx {\mathcal{D}}_{\delta \rho }(\bar{\rho },v_0){\nabla }^2\delta \rho (\mathbf{r},t).$$

$${\mathcal{D}}_{\delta \rho }(\bar{\rho },v_0)=D_t^{\text{eff}}+\frac{v^{\text{eff}}(2v^{\text{eff}}-v_0)}{d(d-1)D_r^{\text{eff}}},$$

$$⟨{\dot{\mathbf{r}}}_i(t)\cdot {\hat{\mathbf{e}}}_i(t)⟩$$

$$⟨{\hat{\mathbf{e}}}_i(t)\cdot {\hat{\mathbf{e}}}_i(0)⟩$$

$${\partial }_t\psi =$$

$$-\mathcal{R}\cdot \left(\beta D_r\mathcal{T}-D_r\mathcal{R}\psi \right),$$

$$\mathcal{F}({\mathbf{r}}_1,{\hat{\mathbf{e}}}_1,t)\equiv$$

$$\int \text{d}{\mathbf{r}}_2\text{d}{\hat{\mathbf{e}}}_2\left(-{\nabla }_1{V}_{{\hat{\mathbf{e}}}_1,{\hat{\mathbf{e}}}_2}({\mathbf{r}}_1,{\mathbf{r}}_2)\right){\psi }_{{\hat{\mathbf{e}}}_1,{\hat{\mathbf{e}}}_2}^{(2)}({\mathbf{r}}_1,{\mathbf{r}}_2,t);$$

$$\mathcal{T}({\mathbf{r}}_1,{\hat{\mathbf{e}}}_1,t)\equiv$$

$$\int \text{d}{\mathbf{r}}_2\text{d}{\hat{\mathbf{e}}}_2\left(-{\mathcal{R}}_1{V}_{{\hat{\mathbf{e}}}_1,{\hat{\mathbf{e}}}_2}({\mathbf{r}}_1,{\mathbf{r}}_2)\right){\psi }_{{\hat{\mathbf{e}}}_1,{\hat{\mathbf{e}}}_2}^{(2)}({\mathbf{r}}_1,{\mathbf{r}}_2,t).$$

$${\psi }_{{\hat{\mathbf{e}}}_1{\hat{\mathbf{e}}}_2}^{(2)}({\mathbf{r}}_1,{\mathbf{r}}_2,t)=\psi ({\mathbf{r}}_1,{\hat{\mathbf{e}}}_1,t)\psi ({\mathbf{r}}_2,{\hat{\mathbf{e}}}_2,t)g_{{\hat{\mathbf{e}}}_1{\hat{\mathbf{e}}}_2}({\mathbf{r}}_1,{\mathbf{r}}_2,t)$$

$$\tilde{\mathcal{F}}({\mathbf{r}}_1,{\hat{\mathbf{e}}}_1,t)\equiv$$

$$\int \text{d}{\mathbf{r}}_2\text{d}{\hat{\mathbf{e}}}_2\left(-{\nabla }_1{V}_{{\hat{\mathbf{e}}}_1,{\hat{\mathbf{e}}}_2}({\mathbf{r}}_1,{\mathbf{r}}_2)\right)\psi ({\mathbf{r}}_2,{\hat{\mathbf{e}}}_2,t)g_{{\hat{\mathbf{e}}}_1,{\hat{\mathbf{e}}}_2}({\mathbf{r}}_1,{\mathbf{r}}_2,t).$$

$$\tilde{\mathcal{F}}(\mathbf{r},\hat{\mathbf{e}},t)\cdot \hat{\mathbf{e}}=-\rho (\mathbf{r},t)\zeta (\bar{\rho },v_0).$$

$$v^{\text{eff}}=v_0-\beta D_t\rho (\mathbf{r},t)\zeta (\bar{\rho },v_0).$$

$$\mathcal{T}(\mathbf{r},\hat{\mathbf{e}},t)\approx b{\mathcal{R}}_1\psi (\mathbf{r},\hat{\mathbf{e}},t).$$

$${\partial }_t\psi =-\nabla \cdot \left(v^{\text{eff}}\hat{\mathbf{e}}\psi -D_t\nabla \psi \right)+D_r^{\text{eff}}\mathcal{R}\cdot \mathcal{R}\psi ,$$

$$\rho (\mathbf{r},t)$$

$$m_{\alpha }(\mathbf{r},t)$$

$$S_{\alpha \beta }(\mathbf{r},t)$$

$${\partial }_t\rho =$$

$${\partial }_t{m}_{\alpha }=$$

$$-(d-1)D_r^{\text{eff}}{m}_{\alpha };$$

$${\partial }_t{S}_{\alpha \beta }=$$

$$-d(d-1)D_r^{\text{eff}}{S}_{\alpha \beta }.$$

$$\rho (\mathbf{r},t)$$

$$\mathbf{m}(\mathbf{r},t)$$

$$\mathbf{S}(\mathbf{r},t)$$

$$\mathbf{\delta m}(\mathbf{r},t)=-\frac{1}{d(d-1)D_r^{\text{eff}}}\nabla \left(v^{\text{eff}}(\mathbf{r},t)\rho (\mathbf{r},t)\right).$$

$$\mathbf{\delta m}(\mathbf{r},t)=-\frac{1}{d(d-1)D_r^{\text{eff}}}\left(v_0-2\beta D_t\rho (\mathbf{r},t)\zeta \right)\nabla \rho (\mathbf{r},t).$$

$${\partial }_t\delta \rho (\mathbf{r},t)={\mathcal{D}}_{\delta \rho }(\bar{\rho },v_0){\nabla }^2\delta \rho (\mathbf{r},t),$$

Full text

I Introduction

Recently, there has been an increased interest in the thermodynamics of what has been coined active matter: systems formed by agents that can convert ambient or external energy into kinetic energy. These systems are diverse. Many are biological in nature: birds, fish, cells and bacteria all convert some form of ambient chemical energy into kinetic energy. Synthetic systems also exist in the form of colloidal particles that self-propel, typically by diffusiophoresis Theurkauff et al. (2012); Palacci et al. (2013); Buttinoni et al. (2013); Ginot et al. (2018). All these active-matter systems are driven so far out of equilibrium that current theories of statistical thermodynamics fail to describe many of their properties. The dynamics of active-matter systems can be modeled quite easily. However, it has unfortunately proven to be very difficult to relate the dynamics to their steady states or to any kind of probability distribution. If we could apply the concepts of equilibrium statistical thermodynamics to active systems, it would greatly improve our ability to describe and predict the structure and behaviour of these systems. Some progress has been made in this regard. For instance, for self-propelled disks, the pressure has been defined Solon et al. (2015a); Patch et al. (2017), the glass transition has been investigated Fily et al. (2014a), and equations of state have been constructed Levis et al. (2017). Furthermore, the work of Cates and Tailleur contributes significantly in constructing effective free energies for these nonequilibrium systems Cates and Tailleur (2015).

For active 2D disks or 3D spheres, two well-studied model systems, there is often a parameter regime in which the system demixes into a dense and a dilute region. This phase separation closely resembles the well-known gas-liquid coexistence found in, for instance, water or Lennard-Jones systems. Unlike the gas-liquid phase separation, however, the clustering in active matter occurs because of the motility rather than the particle-particle attractions. Hence, the phenomenon has been coined motility-induced phase separation (MIPS). Recently, MIPS has been studied extensively: it has been identified for both active Brownian Fily and Marchetti (2012); Redner et al. (2013a); Fily et al. (2014b); Marchetti et al. (2016) and run-and-tumble particles Tailleur and Cates (2008), its nucleation has been studied Redner et al. (2016), its interface has been shown to allow for a negative surface tension Bialké et al. (2015), it has been derived from equations of state Patch et al. (2017); Levis et al. (2017), and from nonequilibrium thermodynamics theories Cates and Tailleur (2015).

All of the above studies logically constrained themselves to the simplest possible model systems, in which particles interact either through hard-particle excluded-volume interactions or through short-range repulsions. Importantly, such models contain no torques. Studies that do include torques typically fall into two categories. The first uses particles with Vicsek-like alignment interactions Barré et al. (2015); Sese-Sansa et al. (2018), which mimic a visual alignment mechanism, such as for birds or fish. The second uses particles with an anisotropic, typically rodlike shape Peruani et al. (2006); Ginelli et al. (2010); Peruani (2016); Weitz et al. (2015); Abaurrea Velasco et al. (2018). This most closely mimics bacteria, whose alignment arises simply from bumping into one another. While studies of active rods reveal a zoo of nonequilibrium phases, they do not exhibit MIPS; there seems to be no parameter regime for which there is a separation into dense and dilute regions without strong alignment. Naturally, this raises questions such as: why does MIPS occur for 2D disks and 3D spheres, but not for 2D and 3D rods? How anisotropic or rod-like must a particle be for MIPS to disappear? In this paper, we will address these questions by both simulations and theory.

To address these questions numerically, we need a model system which exhibits MIPS, and a means to identify MIPS when it occurs. Section II describes both the Active Brownian Particle model we use, and the modified cluster algorithm we apply to identify MIPS. In Section III we present an analytical criterion for the onset of MIPS, on the basis of a stability analysis of density fluctuations in the homogeneous isotropic phase, with the full derivation presented in the Appendix. In Section IV.1 we discuss the phase diagrams for the 2D disks, 3D spheres and 2D and 3D rods, showing unambiguously that MIPS indeed disappears for increasing aspect ratio. Subsequently, we discuss the mechanism behind this suppression in Section IV.2. Section V then concludes this paper by discussing the influence of torque on MIPS in a more general context.

II Computational methods

II.1 Active Brownian Particles

Using Brownian Dynamics (BD) simulations, we study a system of $N$ spherocylinder-shaped active Brownian particles (ABP) of head-to-tail length $l$ and diameter $\sigma\leq l$ in a periodic area $A$, self-propelling with a velocity $v_{0}$ along their long axis $\bm{\hat{e}}$. The particles are subject to rotational and translational noise, with rotational diffusion constant $D_{r}$ and translational diffusion tensor $\bm{D}_{t}=D_{\|}\bm{\hat{e}}\bm{\hat{e}}+D_{\bot}(\bm{\mathcal{I}}-\bm{\hat{e}}\bm{\hat{e}})$, with parallel and perpendicular components $D_{\|}$ and $D_{\bot}$, respectively. For such a 2D system, shown schematically in Fig. 1, the overdamped Langevin equations are given by

[TABLE]

where $i=1,...,N$ is the particle label, $\bm{r}_{i}$ is the position of particle $i$, $\bm{\hat{e}}_{i}=(\cos\theta_{i},\sin\theta_{i})$ the particle orientation, and $\beta=1/k_{B}T$. The force $\bm{\mathcal{F}}_{ij}$ and torque $\mathcal{T}_{ij}$ are due to particle-particle interactions. We assume fluctuation-dissipation to hold on the scale of individual particles, such that the translational and rotational noise terms $\Lambda^{t,\alpha}_{i}$ and $\Lambda^{r}_{i}$, respectively, are Gaussian distributed random numbers with zero mean and unit variance, i.e.:

[TABLE]

To describe excluded-volume interaction between particles $i$ and $j$, we let the forces $\bm{\mathcal{F}}_{ij}=(\partial u_{WCA}(r_{s,ij})/\partial r_{s,ij})\bm{\hat{r}}_{s,ij}$ be the result of a short-range pairwise repulsive WCA potential $u_{WCA}(r_{s,ij})$ acting on the shortest distance $r_{s,ij}$ between particle cores:

[TABLE]

For disks ($l/\sigma=1$), the distance $r_{s,ij}$ is simply the distance between their centers. For $l/\sigma>1$ the cores of the particles are no longer points, but lines. The distance $r_{s,ij}$ then corresponds to the shortest distance between these two line segments. The torques $\mathcal{T}$ are calculated from the forces by $\bm{\mathcal{T}}_{ij}=\bm{a}_{ij}\times\bm{\mathcal{F}}_{ij}$, where $\bm{a}_{ij}$ is the lever arm for the applied force $\bm{\mathcal{F}}_{ij}$ on rod $i$ by rod $j$. For each pair of particles, both the shortest distance $r_{s,ij}$ and the lever arms $\bm{a}_{ij}$ are calculated using the algorithm described in Ref. Vega and Lago (1994). In 2D this torque always points out of plane, so we only need to consider its scalar magnitude $\mathcal{T}$ in the equations of motion.

This 2D model easily generalizes to 3D: aside from vectorial quantities now being three- rather than two-dimensional, we must now also consider the direction of the torque. For convenience, we also switch to vector notation in the orientational equation of motion. The equations of motion in 3D are thus:

[TABLE]

We nondimensionalize the 2D and 3D system by expressing all distances in units of the particle diameter $\sigma$, all energies in terms of the thermal energy $k_{B}T$, and all units of time in terms of $\tau=1/D_{r}$.

II.2 Choice of model parameters and additional assumptions

For our investigation, we will study the influence of four parameters: the dimensionality $d=2$ and $d=3$, the aspect ratio $p=l/\sigma$, the packing fraction $\phi=N((\pi/4)\sigma^{2}+(l-\sigma)\sigma)/A$ ($\phi=N((\pi/6)\sigma^{3}+(\pi/4)(l-\sigma)\sigma^{2})/V$ in 3D) and the Péclet number $\text{Pe}=v_{0}/\sigma D_{r}$. Note that literature sometimes defines the Péclet number in terms of the translational diffusion instead. The diffusion constants $D_{\|}$ and $D_{\bot}$ for rodlike particles can be calculated from simulations including hydrodynamics as in e.g. Ref. Bet et al. (2017) or, for short spherocylinders, approximated by the exact results for ellipsoids (Perrin (1934), see also SI). We found that the influence of this change of diffusion constants is negligible for the aspect ratio range we look at, so for simplicity we will set $D_{\|}=D_{\bot}=D_{t}=\sigma^{2}D_{r}/3$ from now on. This choice corresponds to the correct ratio between translational and rotational diffusion for 3D spheres.

Some care is required in the way we vary the Péclet number. The most straightforward way is to simply vary it by changing the self-propulsion velocity $v_{0}$. However, if we do this and keep the pair interaction strength fixed, the ratio between active and interaction forces will depend on the Péclet number. The result of changing this ratio is that the particle interaction effectively becomes softer as the Péclet number increases. In the extreme case, MIPS may even disappear for high enough Péclet numbers. Earlier work has remarked on this subtlety of varying the Péclet number Stenhammar et al. (2014); Bialké et al. (2013). As our aim is not to provide quantitative but only qualitative data on the phase behaviour, we nevertheless use the straightforward approach by fixing $\epsilon=24k_{B}T$ and changing the Péclet number by varying $v_{0}$.

II.3 Identifying motility-induced phase separation by clustering regions of similar density

MIPS is a separation of a system of self-propelled particles into a dense and a dilute region. While it can be identified quite readily from visual inspection of particle configurations, it is also useful to have a more quantitative method. Two of these methods are common. The first is to measure the distribution of the local density: for a homogeneous system, such a distribution is unimodal, while for a phase-separated system it is bimodal Suma et al. (2014); Speck (2016); Bruss and Glotzer (2018). However, such a distribution can not tell us whether the system has separated into one or into multiple domains, which means it cannot distinguish between micro- and macrophase separation. This distinction becomes important for rods.

The other method is to group particles together into clusters based on a distance cutoff and to determine the fraction $f_{\text{cl}}$ of particles in the largest cluster Fily and Marchetti (2012); Buttinoni et al. (2013); Prymidis et al. (2016). Since MIPS eventually forms one large, dense cluster in a very dilute background gas, $f_{\text{cl}}\rightarrow 1$ for MIPS for large enough systems, while for a homogeneous fluid $f_{\text{cl}}\rightarrow 0$. This latter method requires a cutoff distance that specifies whether particles are close enough to belong to the same cluster. In practice, we found that there is no single cutoff distance that yielded reasonable results for the resulting cluster fraction across all shapes and densities we wish to study.

To solve this problem, we developed a slightly different clustering method, shown schematically in Fig. 2. From the particle positions (Fig. 2a), we make a Voronoi construction. This provides us not only with a parameter-free way to define neighbouring particles, but also with a means of measuring the local packing fraction: $\phi_{l}=((\pi/4)\sigma^{2}+(l-\sigma)\sigma)/A_{v}$ or $((\pi/6)\sigma^{3}+(\pi/4)(l-\sigma)\sigma^{2})/V_{v}$, with $A_{v}$ and $V_{v}$ the area (2D) or volume (3D) of the Voronoi cell (colors in Fig. 2b). Our requirements for two particles to belong to the same cluster are then that (a) their Voronoi cells are connected and (b) they both have a local packing fraction that is either lower or higher than the mean packing fraction $\phi$ by a certain cutoff $\Delta\phi$. Using this method, we create clusters of similar local density (Fig. 2c). We choose $\Delta\phi=0.025$, as we found through trial and error that this cutoff allows us to meaningfully distinguish between homogeneous states with $f_{\text{cl}}<0.5$ and phase-separated states with $f_{\text{cl}}\geq 0.5$ for all aspect ratios and Péclet numbers of interest and for nearly all densities, both in $d=2$ and $d=3$. Note that $f_{\text{cl}}$ is not guaranteed to go to zero in the homogeneous phase when using this definition of clusters due to density fluctuations, but $f_{\text{cl}}=0.5$ still offers a reasonable threshold.

III An analytical criterion for the onset of MIPS

Having described the means to obtain and identify MIPS numerically, we now describe an analytical criterion for the onset of MIPS. We are aware of three ways to obtain such a criterion: by considering the particle flux balance between a dense cluster and a dilute gas phase Redner et al. (2013a, b), by constructing an effective free energy and proceeding as in equilibrium Cates and Tailleur (2013, 2015), and by a stability analysis of density fluctuations of the homogeneous isotropic phase Bialké et al. (2013); Speck et al. (2015). All three methods have previously been used for torque-free systems. We extend the mean-field-like third method laid out in Speck et al. (2015) to 3D systems with torque. The derivation of this extension is given in the Appendix. In short, we map our system to an active ideal gas, where the effect of the many-body forces and torques is subsumed into a modified, effective swim speed $v^{\text{eff}}$, rotational diffusion $D_{r}^{\text{eff}}$, and translation diffusion $D_{t}^{\text{eff}}$. These effective constants then depend on the mean density $\bar{\rho}$ and input swimming speed $v_{0}$. By doing this mapping, we effectively make two approximations: the only effect of the interparticle forces $\bm{\mathcal{F}}$ is to slow particles down, and the only effect of the torques $\bm{\mathcal{T}}$ is to change the rate at which particles change their orientation. The former is a good approximation in the absence of structural order, the latter in the absence of orientational order. Both approximations become poorer at higher densities, where structure and alignment become important.

As will be derived in Appendix 2, the evolution of long-range density perturbations $\delta\rho(\bm{r},t)$ for this active ideal gas is given by

[TABLE]

This is a diffusion equation with a diffusion coefficient $\mathcal{D}_{\delta\rho}$ given by:

[TABLE]

where $d$ now indicates the dimensionality. In agreement with literature and as detailed in the SI, we confirm from our simulations that $D_{t}^{\text{eff}}\approx D_{t}$ for spheres to a reasonable approximation, similar to what was reported in Ref. Speck et al. (2015) for disks. This result extends to 3D rods as well. Thus, we set $D_{t}^{\text{eff}}=D_{t}$ from here on.

The effective constants $v^{\text{eff}}$ and $D_{r}^{\text{eff}}$ can now be found in two ways: we can either formulate closed-form equations for these effective constants, or we could measure them in some way. We choose the latter method, and determine their value from the following correlation functions:

[TABLE]

which measure the effective velocity in the direction of self-propulsion and how quickly a particle loses its orientation, respectively. In other words, we can measure $v^{\text{eff}}$ and $D_{r}^{\text{eff}}$ by simulating a (small) system in the homogeneous isotropic phase.

IV Results & Discussion

To explore the MIPS-related phase behaviour, we performed Brownian Dynamics simulations with $N=10^{4}$ particles in the packing-fraction range $0.1\leq\phi\leq 0.7$, and the Péclet-number range $1\leq\text{Pe}\leq 100$ (2D) and $1\leq\text{Pe}\leq 150$ (3D). This spans the entire density range from the fluid regime to just below the hexatic/solid regimes Bates and Frenkel (2000); Digregorio et al. (2018). The Péclet range spans from below the MIPS critical point to high enough Péclet that the MIPS region attains a near-constant width in density Speck et al. (2014); Stenhammar et al. (2014). We also performed smaller simulations to measure the effective constants $v^{\text{eff}}$ and $D_{r}^{\text{eff}}$. The initial state for all simulations was one with random positions and orientations. Using only a limited number of particles ($N=100$) ensures that the system remains in the homogeneous isotropic phase regardless of density or activity. Of course, these smaller simulations suffer from finite-size effects. The supplementary information contains a more detailed analysis of these finite-size effects.

We now turn to answer the questions posed in the introduction. How elongated do particles have to be to not display MIPS? And what is the mechanism that suppresses MIPS for rods?

IV.1 Phase diagrams of self-propelled disks, spheres and rods

Before we can appreciate how the phase behaviour changes with aspect ratio, we must first establish the relevant features of MIPS for isotropic particles. Let us start in 2D.

While the phase boundaries of MIPS for disks have been studied by a number of authors Redner et al. (2013a); Speck et al. (2014); Solon et al. (2015a); Speck et al. (2015); Levis et al. (2017), a comprehensive study that also includes the high-density hexatic and solid phases has only appeared quite recently Digregorio et al. (2018). In this study, the authors report not only the commonly-reported U-shaped MIPS region (in the density-activity plane), but also that there a narrow density regime wherein a hexatic phase can be found. This regime spans from the passive system ($\text{Pe}=0$) to connect to the MIPS region ($\text{Pe}\sim 100$). Note that Ref. Digregorio et al. (2018) defines the Péclet number in terms of the active force, while we express it in terms of $v_{0}$ and $D_{r}$—this shifts the scale by a constant factor of $D_{r}\sigma^{2}/D_{t}=3$. We can expect to find qualitatively similar features here. Quantitatively, the phase boundaries will be shifted somewhat because of differences in the repulsive pair potential: Ref. Digregorio et al. (2018) uses $U\propto(\sigma/r)^{64}-(\sigma/r)^{32}$, while we use $U\propto(\sigma/r)^{12}-(\sigma/r)^{6}$. Our softer potential decreases the size of the liquid-hexatic coexistence region Kapfer and Krauth (2015). Given the similar temperatures ($k_{B}T/\epsilon=1/24$, versus $k_{B}T/\epsilon=1/20$ in Ref. Digregorio et al. (2018)), the slightly longer range of the WCA potential (cutoff radius $r/\sigma=2^{1/6}\sim 1.12$, as opposed to $r/\sigma=2^{1/32}\sim 1.02$ in Ref. Digregorio et al. (2018)) will shift the solid phases to slightly lower packing fractions. However, an in-depth comparison of the high-density phase boundaries falls outside the scope of this paper. Instead, we will merely identify the solid-like phases by looking at where the effective velocity $v^{\text{eff}}$ becomes vanishingly small ($v^{\text{eff}}<0.1\sigma D_{r}$). Although this is not a very accurate measure, but it serves to crudely distinguish the solid or hexatic phase from the fluid and MIPS phases, at least at low self-propulsion. We use this criterion for all phase diagrams throughout this paper. With this information in mind, let us now consider the phase diagrams in Figure 3.

Figure 3 shows phase diagrams in the Péclet number Pe - packing fraction $\phi$ representation for 2D disks and rods that show both the MIPS region predicted on the basis of the stability analysis (blue-tinted region, $\mathcal{D}<0$), and the MIPS region found in the simulations using $N=10^{4}$ (black points). Both methods seem to indicate MIPS in roughly the same region, but there are a few notable differences. On the low density side, we also find MIPS outside of the predicted spinodal region. Making the analogy with the gas/liquid phase separation, we would expect MIPS in this region to then occur through nucleation and growth. Is this also the case?

A simple way to see if MIPS forms through a nucleation process is to look at domain growth, which we can track using a time series of cluster fraction $f_{\text{cl}}$ defined in Section II.3. If the system immediately decays from an isotropic to a MIPS state, this fraction will likewise increase immediately. If, on the other hand, the system stays in the fluid state for a prolonged period of time, only to later transition into MIPS through a nucleation process, $f_{\text{cl}}$ will retain the value corresponding to the fluid for a finite time.

Figure 4 compares the time evolution of the fraction $f_{\text{cl}}$ of particles in the largest cluster for a number of Péclet numbers at two different densities: one on the low density side of the MIPS regime at $\phi=0.25$ and one on the high density side at $\phi=0.7$. On the low density side and outside of the predicted spinodal region, the cluster fraction can stay constant for a significant amount of time ($t>30\tau$) before transitioning to a MIPS state. On the high density side of the MIPS region, such a delay is absent. The stability analysis predicts spinodal decomposition in this regime, and the cluster growth agrees. This asymmetry is consistent with the findings of Speck et al. Speck et al. (2015), who report that the MIPS transition is discontinuous at low densities, but continuous at high densities.

There is also a discrepancy between the stability analysis and the large-scale simulation at low Péclet numbers. This is to be expected: in this region the fluid-MIPS transition is continuous, and the difference in density between the coexisting phases is small when we close to a critical point. Consequently, distinguishing between clusters of particles is difficult, and the exact choice of cluster fraction threshold $f_{\text{cl}}$ can shift the boundary quite a bit in this region.

Having identified the most important features of the phase diagram for active disks, let us now turn to rods and see how these features change. Figure 3b shows the phase diagram in the density-activity representation for rods with an aspect ratio of $l/\sigma=2$, using the same density and activity ranges as for the disks. The most obvious difference with the rods is that the MIPS region is now both shifted to higher densities and much narrower. The predictions of the stability analysis are worse for the rods: the predicted spinodal now lies in the middle of the simulated MIPS region. We find that the transition from fluid to MIPS now appears to be completely continuous—the system always starts clustering immediately, without any nucleation-like transient period. As can be seen from Figure 5a, the suppression is continuous with increasing aspect ratio, and it eventually pushes the fluid-MIPS transition into the regime where solid phases typically emerge.

Let us now see whether the 3D case is similar. Figure 6 displays phase diagrams in the $(\phi,\text{Pe})$ representation, in Fig. 6a for 3D spheres and in Fig. 6b for 3D rods with $l/\sigma=2$. Somewhat unsurprisingly, they are similar to their 2D counterparts. The most important feature is retained: MIPS disappears when the aspect ratio is increased. The fluid gap we found in between the solid and MIPS phases is also present for the active spheres. However, there are notable differences between the 2D and 3D cases.

In contrast to the 2D case, we observe no region for the active spheres where the MIPS transition is discontinuous. All simulations that form MIPS appear to undergo immediate spinodal decomposition. This does not necessarily mean that there is no binodal region: it may simply be quite small or have low nucleation barriers. The density regime of the metastable region for 3D active spheres is not well understood. We are only aware of one comparable simulation study by Stenhammar et al. Stenhammar et al. (2014), who looked at 2D and 3D active disks/spheres to study the influence of dimensionality. However, their binodal lines were defined as the density at which a high-Péclet system phase separated, which is not directly comparable to the metastable region we define here. Hence, further studies are needed to explain the difference in the width of the metastable region between $d=2$ and $d=3$.

Another difference is at high Péclet number, where the predicted MIPS region for the spheres continues to shift towards higher density, instead of moving towards a constant one. We believe this to be the behaviour that we discussed in Section II.2: for higher Péclet numbers the particles can approach each other closer due to the active forces, causing the effective diameter of the particles to decrease. This effect appears to be stronger in 3D than in 2D, presumably due to the increased coordination of each particle.

The final difference between the 2D and 3D cases is perhaps the most notable one: for the rods, MIPS has disappeared completely. Whatever mechanism suppresses MIPS, it appears to be stronger in 3D than in 2D. Curiously enough, the stability analysis still predicts MIPS in a significant portion of the phase diagram. This discrepancy, combined with its 2D counterpart, suggests that our theoretical approach breaks down for longer aspect ratios. We will see why this is the case in the next section, where we discuss the suppression mechanism.

Armed with the knowledge of these phase diagrams, can we now answer the first question posed in our introduction: “How rodlike must a particle be for MIPS to disappear?” Only partially, unfortunately. Determining the exact aspect ratio where MIPS disappears turns out to be quite difficult. We now know that the nature of the suppression stems from the fluid-MIPS transition shifting to higher densities, but unfortunately our methods to identify MIPS are less reliable at higher densities. More importantly, however, when the particle interactions are not isotropic, MIPS is no longer defined unambiguously and multiple types of clustered phases are possible which all fit the present criteria. When we identify MIPS according to a) the system phase-separating into a single dense cluster in a background gas and b) this cluster having no net orientational order, there are still multiple realizations of such a system (see SI, figure S5), such as a dense cluster with large domains of oppositely oriented particles ($l/\sigma=1.1$) or a percolating cluster with low orientational order and many voids ($l/\sigma=1.3,\ 2.0$). Therefore, establishing the boundaries of MIPS at these higher densities requires a more careful consideration of both hexatic Digregorio et al. (2018) and orientational order Shi and Chaté (2018). We leave this investigation to future work and instead, having established that MIPS is suppressed when particles become elongated, we now turn to finding out why.

IV.2 Torque-induced suppression of motility-induced phase separation

Since the main difference between the disk and rod systems is the presence of torque, it is likely that the suppression of MIPS must arise there. In our stability analysis, the only effect of torques is to modify the rotational diffusion. Looking at Eq. (33), it might be possible to suppress MIPS if $D_{r}^{\text{eff}}$ is increased enough to make the second term smaller than $D_{t}^{\text{eff}}$. Is this the case? Is the rotational diffusion perhaps enhanced so much that we effectively end up with a thermal system again?

Not so. Looking at Figs. 7e and 7f, we can see that $D_{r}^{\text{eff}}$ is indeed increased significantly where MIPS vanishes. However, when we insert the actual values of $v^{\text{eff}}$ and $D_{r}^{\text{eff}}$, we see that this is not the case: the typical values of $v^{\text{eff}}$ are simply too large. So if it is not $D_{r}^{\text{eff}}$, it must be $v^{\text{eff}}$ that contains the key information that allows us to predict MIPS or its suppression. After all, the stability analysis does correctly predict that MIPS is suppressed for high aspect ratios. Comparing the effective swimming speeds of different aspect ratios (Figs. 7a and 7c, or 7b and 7d), we see that the rods slow down less with increasing density than the disks. In other words, the rods hinder each other’s movement less than the disks do. Why is this? Again we must look to the main difference between the two systems: torque.

For disks, one can derive the linear decrease of the velocity with increasing density $v^{\text{eff}}=v_{0}(1-\rho/\rho^{*})$ from mean-field theory and kinetic arguments Stenhammar et al. (2013); Fily et al. (2014a); Bruss and Glotzer (2018). This is done by assuming that particles slow down at low density due to time spent in binary collisions, which leads to $v^{\text{eff}}(\rho)\simeq v_{0}(1-\tau_{c}/\tau_{f})$, where $\tau_{c}$ is the time spent in a collision and $\tau_{f}=1/(\sigma v_{0}\rho)$ the mean free time between collisions. At low density, we expect the mean free time $\tau_{f}$ to be mostly unaffected by the presence of torques as long as there are no significant short-range orientational correlations. The duration of collisions $\tau_{c}$, however, can change significantly when torques are involved. For disks, the duration of their collision—of their hindrance—is determined by how long it takes for them to slide around each other. Rods, however, will rotate to reorient their swimming directions away from the combined center of mass of the collision. This will decrease the collision duration. Since collisions are now shorter, the rods spend more time moving freely: less hindered. Furthermore, this reorientation leads to an enhanced rotational diffusion—exactly what we find.

Interestingly, this suggests that an inverse mechanism might also exist. If the torques between two colliding particles cause the particles to rotate towards their center of mass, collisions would be prolonged and MIPS would be enhanced. Precisely this inverse effect was reported earlier in Refs. Barré et al. (2015); Sese-Sansa et al. (2018): MIPS is enhanced for self-propelled particles that align through Vicsek interactions. In binary collisions, the Vicsek torques always rotate particles towards the combined center of mass, increasing the duration of collisions, increasing hindrance and thus enhancing MIPS.

Is the changing density dependence of $v^{\text{eff}}$ with increasing anisotropy enough to completely describe the suppression of MIPS? If we would have a system of self-propelled particles with some arbitrary shape and we would know how the effective swim speed depends on density, could we then predict whether and where it will undergo MIPS? Unfortunately, no. As we can see from the rod phase diagrams in Figs. 3b and 6b, our stability analysis predicts the right qualitative trend, but its quantitative prediction is poor. This is probably due to neglecting alignment effects in the stability analysis. As the rod length increases, nematic and polar alignment of the particles start playing a more significant role in their phase behaviour, which is not captured by our theory. For instance consider Fig. 8, where we show a snapshot of rods at $\phi=0.5$, $\text{Pe}=100$, just outside the MIPS region, and colour particles according to nematic orientation. The clusters formed by the rods have significant short-range nematic order. Incorporating the dynamics of the polarization and nematic fields using theory developed for active nematics Baskaran and Cristina Marchetti (2010); Bertin et al. (2015) might allow for more accurate predictions for the onset of MIPS for longer rods.

V Conclusions & outlook

In this paper, we showed that motility-induced phase separation does not occur for rodlike particles when they become sufficiently anisotropic. This disappearance is observed both from many-particle simulations and from a stability analysis of the homogeneous isotropic phase. The latter provides a simple criterion for the onset of MIPS by considering the effective swimming speed of the particles and their effective rotational diffusion. Both methods agree qualitatively in that MIPS is pushed to higher densities for increasing rod aspect ratio, and they agree quantitatively for short rods that deviate only slightly from disks or spheres. For longer aspect ratios the quantitative agreement is lost, presumable due to alignment interactions that are present, but not taken into account in the stability analysis.

We also propose a more intuitive explanation for the suppression mechanism. MIPS relies on particles slowing down sufficiently with increasing density Cates and Tailleur (2015). This hindrance is closely linked to the duration of collisions between particles Stenhammar et al. (2013); Fily et al. (2014a); Bruss and Glotzer (2018). Excluding torques, the duration of collisions is determined by how long it takes for them to slide along one another. Including torques can dramatically decrease the duration of collisions by rotating the forward axes of the self-propelled particles away from each other. Formulated in this way, we can also explain the results of Refs. Barré et al. (2015); Sese-Sansa et al. (2018), where MIPS is enhanced for particles with Vicsek interactions. Simply put, Vicsek torques prolong particle collisions, while rodlike excluded volume torques shorten them. Intriguingly, this provides us with a particle design tool to enhance or suppress MIPS. MIPS is enhanced for Vicsek-like interactions Barré et al. (2015); Sese-Sansa et al. (2018), for faceted, concave and/or rough particles Ilse et al. (2016); Prymidis et al. (2016); Cugliandolo et al. (2017), while it is suppressed for smooth particles and rodlike shapes Shi and Chaté (2018). In addition to steric interactions, hydrodynamic interactions between active particles also play an important role in whether or not MIPS can form. While hydrodynamics seems to usually suppress MIPS Zöttl and Stark (2014); Matas-Navarro et al. (2014); Blaschke et al. (2016), the details depend on whether particles are “pushers” or “pullers” and on the dimensionality Oyama et al. (2016, 2017).

Despite recent advancements, the role of torque in active systems is still not well understood. Much of the developed theory has been restricted to the torque-free regime, but recent numerical studies suggest that torque, either from boundaries Solon et al. (2015b) or from particle interactions Barré et al. (2015); Sese-Sansa et al. (2018), can have a significant effect on the structure and dynamics of active matter systems. In order to understand active matter beyond torque-free model systems, more theoretical work is needed to elucidate the influence of torques in active systems.

Supplementary Material

See supplementary material for an analysis of the significance of the effective translational diffusion $D_{t}^{\text{eff}}$, the hydrodynamic friction coefficients for short spherocylinders, a finite-size analysis of the effective self-propulsion velocity $v^{\text{eff}}$ and effective rotational diffusion $D_{r}^{\text{eff}}$, and additional simulation snapshots.

Acknowledgements

This work is part of the D-ITP consortium, a program of the Netherlands Organisation for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW). We acknowledge financial support from an NWO-VICI grant. We would like to acknowledge Bram Bet for useful discussions and information on the friction coefficients of short spherocylinders.

Appendix A Stability analysis of the homogeneous isotropic phase, including torque

There is currently no general theory for out-of-equilibrium statistical physics from which we can derive the onset of MIPS. In equilibrium, one can derive spinodal and binodal lines by considering (derivatives of) the free energy. Out of equilibrium, effective free energies can only be constructed under specific circumstances Cates and Tailleur (2015). Nevertheless, there are other ways to derive criteria for the onset of MIPS. We are aware of three ways to obtain such criteria: by constructing an effective free energy and proceeding as in equilibrium Cates and Tailleur (2013, 2015), by looking at the particle flux balance between a dense cluster and a dilute gas phase Redner et al. (2013a, b), and by a stability analysis of density fluctuations of the homogeneous isotropic phase Bialké et al. (2013); Speck et al. (2015). The first method cannot be applied directly to our system, as one of its underlying assumptions is that no torques act between the particles. The second method is also likely to fail for rods, as it relies on the assumption that the orientations of particles in the boundary of the dense cluster evolve diffusively. Thus, we opt for the third method: deriving a criterion for the (in)stability of the homogeneous isotropic phase to density fluctuations, by extending the mean-field-like method from Ref. Speck et al. (2015) to 3D systems with torques. This is done by capturing the effect of the torques into a modified, effective rotational diffusion $D_{r}^{\text{eff}}$. Here we describe this extended derivation.

A.1 Effective Smoluchowski equation

To render the problem analytically tractable, our first goal is to simplify the effect of the pairwise forces and torques. We will do this using a mean-field-like approximation. Following the same procedure as Refs. Bialké et al. (2013); Speck et al. (2015), we start from the Smoluchowski equation for the one-particle probability density $\psi(\bm{r},\bm{\hat{e}},t)$, given by

[TABLE]

where $\bm{\nabla}$ are the 2D and 3D gradient operators and $\bm{\mathcal{R}}$ is the rotation operator, defined as $\mathcal{R}=\partial_{\theta}$ in 2D and $\bm{\mathcal{R}}=\bm{\hat{e}}\times\bm{\nabla}_{\bm{\hat{e}}}$ in 3D. Note that similar to our numerical model we neglect the influence of particle shape anisotropy on the translational diffusion and simply set $\bm{D}_{t}=D_{t}\bm{\mathcal{I}}$. The pairwise force density $\bm{\mathcal{F}}$ and torque density $\bm{\mathcal{T}}$, which arise due to the particle-particle interactions of a pair potential $V_{\bm{\hat{e}}_{1},\bm{\hat{e}}_{2}}(\bm{r}_{1},\bm{r}_{2})$, can then be written in terms of the two-body probability density $\psi^{(2)}_{\bm{\hat{e}}_{1},\bm{\hat{e}}_{2}}(\bm{r}_{1},\bm{r}_{2},t)$ as

[TABLE]

[TABLE]

In order to close this hierarchy, the force and torque densities $\bm{\mathcal{F}}$ and $\bm{\mathcal{T}}$ need to be expressed in terms of the one-body PDF. To do so, we first use the identity

[TABLE]

to rewrite Eq. (13) as $\bm{\mathcal{F}}=\bm{\tilde{\mathcal{F}}}\psi$, where

[TABLE]

To obtain a closure, we make the following assumptions. First, we assume that the force $\bm{\mathcal{F}}$ acts along the direction of self-propulsion, i.e. $\bm{\mathcal{F}}=(\bm{\mathcal{F}}\cdot\bm{\hat{e}})\bm{\hat{e}}$. Whereas this is exact in a homogeneous, isotropic bulk as dictated by symmetry, in general we neglect a possible second component that is perpendicular to $\bm{\hat{e}}$. In Ref. Speck et al. (2015), Speck et al. consider this second component to act along the gradient of the one-particle PDF i.e. $\bm{\mathcal{F}}=(\bm{\mathcal{F}}\cdot\bm{\hat{e}})\bm{\hat{e}}+a\nabla\psi$. This additional assumption leads to a modified translational diffusion $D_{t}^{\text{eff}}=(1-\beta a)D_{t}$. We measured the magnitude of this modification for 3D spheres and rods, and found that the modification provided by $\beta a$ is of negligible influence on the location of the phase boundaries. Thus, we do not consider this additional component here and simply set $a=0$, i.e. $D_{t}^{\text{eff}}=D_{t}$, from now on. We did not explicitly check the validity of this assumption in the 2D case, but see no reason to assume a difference.

To continue our derivation, we make the second assumption that $\bm{\tilde{\mathcal{F}}}\cdot\bm{\hat{e}}$ is linear in the local density $\rho(\bm{r},t)=\int\text{d}\bm{\hat{e}}\psi(\bm{r},\bm{\hat{e}},t)$ and has no further dependence on $(\bm{r},\bm{\hat{e}},t)$:

[TABLE]

Here the constant $\zeta$ is independent of $(\bm{r},\bm{\hat{e}},t)$, but can still depend on the mean density $\bar{\rho}=N/A$ (or $N/V$ in 3D) and the self-propulsion strength $v_{0}$. In this way, using Eq. (12), the effect of the interaction forces can be absorbed into a modified self-propulsion velocity $v^{\text{eff}}$, which is given by

[TABLE]

For the torques, we make the approximation that its only influence is to modify the rotational diffusion i.e.

[TABLE]

Proceeding on, we assume the corresponding constant to be homogeneous and isotropic, depending only on the mean density and self-propulsion: $b=b(\bar{\rho},v_{0})$. With these approximations we can simplify Eq. (12) as the Smoluchowski equation for an active ideal gas:

[TABLE]

where $v^{\text{eff}}$ and $D_{r}^{\text{eff}}$ are now the effective self-propulsion and rotational diffusion constant, respectively.

A.2 Stability analysis of the homogeneous isotropic phase

Now that we have reduced the full Smoluchowski Eq. (12) into the ideal-gas form of Eq. (20), we can perform a linear stability analysis on the homogeneous isotropic phase. We start by defining the relevant moments of the one-particle PDF $\psi(\bm{r},\bm{\hat{e}},t)$,

[TABLE]

Here, the Greek indices label the Cartesian vector- or tensor components, and in the following we shall employ the Einstein summation convention. Considering the same moments of the ideal gas Smoluchowski equation (20) yields the following evolution equations:

[TABLE]

Here $B$ is the next (third) order moment. The structure of this hierarchy of time-evolution equations (24)-(26) is such that the time-derivative of each moment depends linearly on itself and lower order ones, and on the next one. However, as we shall see, moments beyond $m$ are irrelevant for the instability we wish to consider.

A steady-state solution to Eq. (20) is the homogeneous isotropic phase: $\psi(\bm{r},\bm{\hat{e}},t)\propto\bar{\rho}$. Expressed in terms of the moment equations (24)-(26), this gives $\rho(\bm{r},t)=\bar{\rho}$ and $\bm{m}(\bm{r},t)=\bm{S}(\bm{r},t)=0$. To obtain a criterion for the stability of this solution, we investigate the behaviour of small perturbations to the homogeneous state:

[TABLE]

Since MIPS is a macroscopic phase separation, we should study the instability with respect to long-range perturbations i.e. perturbations with small spatial gradients. In this limit, the dynamics are dominated by the terms in Eqs. (24)-(26) with the fewest gradients. Of the three moments, it is $\rho$ whose time evolution is slowest. Its timescale is of order $\nabla^{-1}$, while $\bm{m}$ and $\bm{S}$ evolve as $(D_{r}^{\text{eff}})^{-1}\sim\nabla^{0}$. As we are interested in the evolution of the density perturbations, i.e. of the slow variable, we can assume that at any given time, the higher moments $\bm{m}$ and $\bm{S}$ are given by their steady-state solutions that correspond to the density profile $\rho(\bm{r},t)$ at that instant. Solving Eq. (26) for its steady-state solution $\delta S_{\alpha\beta}$ reveals that it scales as $\mathcal{O}(\nabla^{1})$. Therefore, its contribution to the evolution of polarization perturbations (Eq. (25)) is of higher order. To leading order, the evolution of polarization perturbations is then given by

[TABLE]

Recalling that $v^{\text{eff}}=v_{0}-\beta D_{t}\rho(\bm{r},t)\zeta(\bar{\rho},v_{0})$, we can take this gradient explicitly and obtain

[TABLE]

Using this result, the equation for the time evolution of density perturbations becomes

[TABLE]

which is a diffusion equation with diffusion constant

[TABLE]

Whenever the diffusion constant $\mathcal{D}_{\delta\rho}$ is negative, density perturbations $\delta\rho(\bm{r},t)$ will grow. Therefore, the region in $(\bar{\rho},v_{0})$-space where the homogeneous isotropic phase becomes unstable is then given by the condition $\mathcal{D}_{\delta\rho}(\bar{\rho},v_{0})<0$. This can only occur for self-propulsion velocities $v_{0}$ above the critical threshold $v^{*}=2\sqrt{2}\sqrt{d(d-1)}\sqrt{D_{t}D_{r}^{\text{eff}}}$.