Re-entrant phase separation in nematically aligning active polar particles

TL;DR

This paper numerically investigates the complex phase behavior of active polar particles with nematic alignment, revealing re-entrant phase separation and transitions driven by activity and noise levels.

Contribution

It introduces a detailed numerical analysis of phase transitions, including re-entrant behaviors, in active polar particles with nematic alignment, highlighting novel phase phenomena.

Findings

Re-entrant fluid-phase separation observed at varying activity levels.

Distinct mechanisms for phase separation and remelting identified.

Flocking and jamming behaviors influence phase transitions.

Abstract

We present a numerical study of the phase behavior of repulsively interacting active polar particles that align their active velocities nematically. The amplitude of the active velocity, and the noise in its orientational alignment control the active nature of the system. At high values of orientational noise, the structural fluid undergoes a continuous nematic-isotropic transition in active orientation. This transition is well separated from an active phase separation, characterized by the formation of high density hexatic clusters, observed at lower noise strengths. With increasing activity, the system undergoes a re-entrant fluid-phase separation-fluid transition. The phase coexistence at low activity can be understood in terms of motility induced phase separation. In contrast, the re-melting of hexatic clusters, and the collective motion at low orientational noise are dominated by flocking behavior. At high activity, sliding and jamming of polar sub-clusters, formation of grain boundaries, lane formation, and subsequent fragmentation of the polar patches mediate remelting.

Figures

Equations

$${\dot{\mathbf{r}}}_i=v_0{\text{n}}_i(t)+\mu \sum _j{\mathbf{F}}_{ij}+{\xi }_i(t),$$

$$U(r_{ij})=\{\begin{array}{ll}4ϵ\left\{(\sigma /r_{ij}{)}^{12}-(\sigma /r_{ij}{)}^6\right\},& \text{if }{r}_{ij}\le 2^{1/6}\sigma \\ 0,& \text{otherwise}\end{array}$$

$${\theta }_i^{t+\delta t}=arg\left[\sum _{j\in R_i}\mathrm{sign}[\cos ({\theta }_j^t-{\theta }_i^t)]e^{i{\theta }_j^t}\right]+{\zeta }_i(\eta ).$$

Full text

I Introduction

Active particles transduce stored or ambient energy into motion, generating non-equilibrium drive from the smallest scales of a system Marchetti2013 ; Romanczuk2012 ; Vicsek2012 . Examples in biology range from molecular motors, cytoskeleton, cells, tissues Schliwa2003 ; Alberts2007 ; Howarda ; Surrey2001 ; Schaller2010 , bacterial suspension Cates2012 ; Sokolov2007 to ant trails, bird flocks, animal groups and human crowd Toner2005 ; Ballerini2008 ; Couzin2007 ; Couzin2003 ; Katz2011 ; Attanasi2014 ; Helbing1997 ; Helbing2000 . In nonliving matter, examples include colloidal Janus particles Zottl2016 , and vibrated asymmetric granular matter Weber2013 . The seminal model by Vicsek, and its subsequent extensions Vicsek1995 ; Chate2006 ; Chate2008 ; Vicsek2012 , considered active point particles that orient the direction of their active velocity ferromagnetically, i.e., parallel to that of their neighbours. The collective properties of active particles depend only on a few microscopic details. The self propulsion can be polar Vicsek1995 ; Gregoire2004 or apolar Chate2006 ; Gruler1999 . The direction of propulsion may align parallel to or in a nematic manner with the neighboring particles Vicsek1995 ; Ginelli2010 . The alignment may arise from hydrodynamic coupling Sokolov2007 ; Katz2011 , visual cues Nagy2010 ; Ballerini2008 , or direct collision between extended objects Narayan2007 ; Peruani2006 ; Weber2013 ; Weitz2015 ; Shi2018 ; Grossman2019 . With increasing tendency of alignment and activity, the coupling of velocity fluctuations with density lead to flocking, collective motion, and giant number fluctuations Toner1995 ; Toner1998 ; Narayan2007 ; Dey2012 ; Gregoire2004 ; Peruani2011b ; Solon2013 ; Peruani2013 .

Another important aspect in the emergent phase behavior of active particles is the impact of volume exclusion. Its interplay with activity is most extensively studied in active Brownian particles (ABP) Fily2012 ; Redner2013 ; Solon2015 ; Stenhammar2013 ; Cates2015 . In ABPs the orientation of self propulsion of each particle undergoes Brownian rotation, independent of others. This persistent motion suppresses the change in velocity direction after collision between ABPs, unlike in elastic scattering. The resultant slow-down generates a positive feedback effect allowing accumulation of more particles to form clusters, eventually leading to a motility induce phase separation (MIPS) Cates2015 . The high density clusters could be in either a hexatic or solid phase, determined by the amount of activity and mean density Digregorio2018 . It has been shown recently that at low activity, the equilibrium two stage defect mediated melting scenario of two-dimensional solid Kosterlitz1973 ; Halperin1978 ; Young1979 ; Kapfer2015 ; Bernard2011 ; Gasser2010 ; Thorneywork2017 persists, with the melting points shifting to higher densities with increasing activity Klamser2018 ; Digregorio2018 .

The two aspects of active alignment and exclusion have been combined in a number of recent studies analyzing the impact of volume exclusion on ferromagnetically aligning active polar particles Peruani2011b ; Weber2014 ; Martin-Gomez2018 ; Sese-Sansa2018 . The role of topological defects in the melting of such active solids was demonstrated Weber2014 . At lower densities, the system shows formation of micro or macro clusters, lanes and bands Martin-Gomez2018 . A recent experiment observed interrupted phase separation in active Janus colloids where impact of alignment and excluded volume are simultaneously present Linden2019 . Active granular disks with built-in asymmetry showed polar activity under external shaking. These disks repel each other in a circularly symmetric fashion, and their active velocities showed approximate nematic alignment Weber2013 . Active point particles with nematic alignment were known to undergo nematic to isotropic phase transition and associated clustering Ginelli2010 .

In this paper, we investigate the effect of volume exclusion on the collective properties of nematically aligning active polar particles. Using numerical simulations, we present a comprehensive study of the active orientational phase transition and spatial phase separation with changing values of the active velocity and orientational noise. We fix the mean density at a value that is slightly lesser than the equilibrium melting point. With reduction of alignment noise, the active fluid shows a continuous isotropic to nematic phase transition. However, across this transition, structurally, the system remains in a homogeneous fluid phase, unlike the point particles. With further reduction of the noise, the active particles undergo a phase separation showing coexistence of high density hexatic clusters and low density fluid. This emerges due to the interplay of activity and excluded volume interaction. At an intermediate noise, the system demonstrates a remarkable re-entrant transition from fluid - phase coexistence - fluid with increasing active velocity. We present the detailed phase diagram, supported by a comprehensive characterization of the respective phases and phase-transitions in this paper. As we show, the onset of clustering at low active velocities is like MIPS, mediated by an effective slowdown of particles with density. The remelting of these active hexatic clusters at higher active velocities takes place by formation of grain boundaries and fragmentation of local polar patches. The cluster size distribution changes qualitatively at the onset of the phase separation. The hexatic clusters display complex dynamics – from sliding stripes, to jamming, to formation of counter-propagating channels. A finite size scaling analysis shows the robustness of the clustering transition, the mean cluster size remains a finite fraction of the overall system size even in the largest systems.

In Section II, we describe the model and simulation details. We present the complete phase diagram and the details of the phase transitions in Section III. Finally we conclude in Section IV presenting a summary and outlook.

II Model

We consider a two-dimensional system of $N$ active polar particles in a volume $A=L_{x}L_{y}$, performing active Brownian motion with a self-propulsion speed $v_{0}$, in direction ${\textbf{n}_{i}}\,=(\cos\theta_{i},\sin\theta_{i})$ for the $i$-th particle, where $\theta_{i}$ denotes the angle ${\textbf{n}_{i}}$ subtends with the $x$-axis. The positional evolution of the $i$-th particle is described by an overdamped Langevin equation

[TABLE]

where, $\mu$ is the mobility, $\xi_{i}(t)$ is a Gaussian white noise with $\langle\xi_{i}(t)\rangle=0$ and $\langle\xi_{i}(t)\xi_{i}(0)\rangle=2D\delta(t)$ so that $D=\mu k_{B}T$, ${\bf F}_{ij}=-\nabla_{i}U(r_{ij})$, is the repulsive force between particles given by the WCA potential Weeks1971

[TABLE]

where, $r_{ij}=\mid{\bf r}_{i}-{\bf r}_{j}\mid$ is the inter-particle separation. We assume a nematic alignment of the polar orientation ${\textbf{n}_{i}}$ of the $i$-th particle with its neighbors. It aligns (anti-) parallel to the local mean active orientation, if the relative angle between the two orientations is (obtuse) acute. We use the evolution Ginelli2010 ,

[TABLE]

The alignment gets randomized by the orientational noise $\zeta_{i}(\eta)$ that follows uniform distribution between $(-\eta\pi,\eta\pi)$ with the noise strength $\eta\in(0,1)$.

The unit of length, energy, and time are set by $\sigma$, $\epsilon=k_{B}T$, and $\tau=\sigma^{2}/D$, respectively. The activity of the particles is controlled by a dimensionless Péclet number $Pe=v_{0}\sigma/D$, and the noise strength $\eta$. The numerical simulation is performed using Euler integration. We use integration time step $\delta t=10^{-5}\tau$ for $Pe\leq 50$, and a smaller time-step $\delta t=5\times 10^{-6}\tau$ for higher $Pe$, to avoid inter-particle overlap.

In the absence of activity, the two-dimensional solid is expected to display a solid- hexatic- liquid equilibrium phase-transition Bernard2011 ; Kapfer2015 ; Kosterlitz1973 ; Halperin1978 ; Young1979 , with the solid melting point at the dimensionless density $\rho=(N/A)(\pi/4)\sigma^{2}\approx 0.71$. We present results for $\rho=0.6$, a density slightly below the equilibrium melting point. At this density, the excluded volume interaction plays an important role, while not dominating over activity.

III Phase transition

In this section, we present a detailed phase diagram, and extensive discussion of the corresponding phase transitions.

III.1 Phase diagram

Fig.1($a$) shows the phase diagram in the $Pe$-$\eta$ plane, identifying the isotropic, and the nematic fluid, and the region of active phase separation. With decreasing $\eta$, the active orientations of isotropic fluid gets aligned in the nematic fashion, captured by the time-averaged scalar order parameter $S=\frac{1}{\tau_{0}}\int_{0}^{\tau_{0}}dt\,|{\bf s(t)}|$, where the mean nematic order of a given configuration is ${\bf s}=(1/N)\sum_{k=1}^{N}\exp(2i\theta_{k})$. The increasing convection of orientational information with increasing $Pe$, requires a larger $\eta$ to destroy the nematic order. However, due to the volume exclusion, the range over which the nematic alignment may get extended due to convection is restricted at a given density. Thus the isotropic- nematic phase boundary shows an initial increase with $Pe$ before eventual saturation (Fig.1($a$) ). The transition points denoted by the filled blue circles in Fig.1($a$) are evaluated using $N=1024$ particles. They are identified by the maximum of fluctuations in the order parameter $S$. For comparison, we performed simulations for a larger system of $N=16384$ particles at $Pe=10$ and $20$. The red filled circles show the isotropic-nematic transition points corresponding to this system size. They compare well with $N=1024$ particle results, within error bars. The determination of the transition points and the corresponding error bars are discussed in detail in Sec. III.2.

At further lower $\eta~{}(<0.06)$, the active system undergoes phase separation between a nematic fluid and denser hexatic clusters. The phase coexistence is captured by the bimodal distribution function of the local density $\rho({\bf r})$, and the amplitude of the local hexatic bond-orientational order $\vec{\psi_{l}}({\bf r})=\langle\vec{\psi}_{6}^{k}\rangle_{k\in{n_{l}}}$, where $n_{l}$ denotes the number of particles in a coarse-grained volume $l_{x}\times l_{y}$ around ${\bf r}$. The bond-orientational order associated with the $k$-th particle is defined as $\vec{\psi}_{6}^{k}=\sum_{j=1}^{n_{v}}(l_{j}/l)\,\exp(6i\phi_{kj})$ Mickel2013 ; Steinhardt1983 , where, $\phi_{kj}(t)$ is the angle subtended by the bond between the $k$-th particle and its $j$-th Voronoi neighbor on the positive $x$-axis; $l_{j}$ is the length of the Voronoi edge corresponding to the $j$-th neighbor, and $l=\sum_{j=1}^{n_{v}}l_{j}$ where $n_{v}$ is the total number of Voronoi topological neighbors. The solid line denotes the onset of phase coexistence of the nematic fluid and hexatic clusters. The non-monotonic variation of this boundary as a function of $Pe$ indicates a re-entrant transition. For example, keeping the orientational noise strength fixed at $\eta=0.05$, change in $Pe$ from $5$ to $30$ leads a nematic fluid to the hexatic- nematic fluid phase coexistence. The hexatic clusters remelt into a homogeneous nematic fluid by increasing $Pe$ further to $60$. The dashed line $\eta\sim Pe^{1/2}$ shown at the phase boundary is an estimate obtained using the argument of MIPS, as outlined later in Sec.III.3. The clusters in the rest of the phase coexistence regime are dominated by alignment and flocking. Before going into the quantitative analysis, we illustrate the local phase behavior in terms of system configurations and local hexatic order.

Fig.1($b$)-($f$) show typical instantaneous configurations of particles shaded by colors denoting the projection $h_{k}=\vec{\psi}_{6}^{k}\cdot\vec{\psi}_{6}$, of their individual hexatic orientations $\vec{\psi}_{6}^{k}$ to the instantaneous system averaged bond-orientational order $\vec{\psi}_{6}$. $h_{k}$ is bounded between $-1$ and $1$. Fig.1($b$) shows a typical configuration in the two phase regime. The high density clusters of particles are observable as contiguous colored patches. The large red patches identify regions with local bond orientational order aligned along the mean instantaneous hexatic order. On the other hand, the blue patches denote the regions where local bond orientational orders are anti-parallel to that of the mean hexatic order. A magnified portion of one such high density cluster is shown in the upper right inset in Fig.1($a$). It shows local formation of solid-like triangular lattice structure. This region of high six-fold bond orientational order is a jammed state of two sets of active particles with opposite polarity, denoted by red and blue colors in this figure. In the upper left inset of Fig.1($a$), we show the same configuration in terms of Voronoi tessellation, color coding the disclinations where the particles do not have six Voronoi neighbours, denoting particles with five neighbors in blue and that with seven neighbors in red. In a pure triangular lattice, all particles have six neighbors. The presence of bound $5-7$ disclination pairs, well-separated from each other, identify free dislocations, a characteristic of the hexatic phase. Fig.1($c$) shows a typical configuration in the isotropic fluid phase. The two phase points at which configurations Fig.1($b$),($c$) are plotted are denoted by pink circles in the phase diagram Fig.1($a$). The black dashed line in the phase diagram Fig.1($a$) follows a remelting transition. Typical configurations corresponding to the three points denoted by $\blacktriangle$ in Fig.1($a$) are presented in Fig.1($d$)-($f$), showing significant hexatic clustering in Fig.1($e$). The details of these phase transitions are quantified in the following.

III.2 Isotropic Nematic Transition

The nematic-isotropic transition is quantified by following the reduction of the nematic scalar order parameter $S$ with the increase in the orientational noise strength $\eta$. This is shown for seven different $Pe$ values in Fig.2, using a system of $N=1024$ particles. The fluctuation of the order parameter $(\Delta S)^{2}=[\langle S^{2}\rangle-\langle S\rangle^{2}]$ shows a pronounced maximum at the transition point $\eta_{c}$. A representative behavior is shown in the inset of Fig.2 for $Pe=10$. This maximum is identified as the transition point. To determine the exact location of the maximum $\eta_{c}$ and its standard deviation $\Delta\eta_{c}$, we use a Gaussian fitting around the maximum (red solid line in the inset of Fig.2). These constitute the nematic-isotropic phase transition points and their error-bars plotted in Fig.1($a$). With increasing system size, the transition becomes sharper, however, the transition point remains within the error bar. One data set for a larger system of $N=16384$ particles at $Pe=10$ is shown using the filled blue circles in Fig.2. The observed transition is between a quasi-long ranged ordered (QLRO) nematic phase to a disordered isotropic phase. The QLRO within the nematic phase is illustrated in Appendix-A in terms of an algebraic decay of the scalar order parameter with system size. A further discussion on the nature of $S(\eta)$ curves, and their data-collapse in an intermediate $Pe$ regime are presented in Appendix-B.

The nature of the nematic-isotropic phase transition is further qualified using the probability distribution of the scalar order parameter, $P(S)$, across the transition. At a given activity, with increasing $\eta$ the probablity distribution remains unimodal, only the peak position shifts from $S\approx 1$ to values close to $S=0$ (Fig. 3). This unimodal nature of the distribution functions signifies the absence of any metastable phase on the other side of the transition, a characteristic of continuous transitions. The distributions broaden near the phase transition point capturing the enhanced fluctuations.

The continuous nature of the transition is further characterized using the two-point correlation of the normalized nematic order parameter $C_{ss}(r)=\langle\sum_{j,k}\cos[2(\theta_{j}-\theta_{k})]\delta(r-r_{jk})/\sum_{j,k}\delta(r-r_{jk})\rangle$. The correlations are plotted for different $\eta$ at $Pe=40$ in the inset in Fig. 3. Within the QLRO nematic phase, the correlations decay as power law $r^{-\nu}$, with the value of the exponent $\nu$ increasing with $\eta$. At the transition point the exponent $\nu\approx 1/4$, in approximate agreement with the prediction of equilibrium Kosterlitz-Thouless transition Kosterlitz1973 .

Across the nematic-isotropic transition, the volume exclusion maintains a homogeneous density distribution. This behavior is quite unlike that observed for point particles, where increased orientational disorder leads to density inhomogeneity resulting in clustering, and band formation Ginelli2010 . This remains the zero density or vanishing exclusion interaction limit of our model. At further smaller $\eta$, enhanced alignment reduces orientational diffusion, thereby, increasing persistence. As a result, as we show in the following section, the interplay of activity and volume exclusion leads to phase separation.

III.3 Phase separation

The structural order in the coexisting phases are analyzed using a system of $N=16384$ particles in a volume of $L_{x}\times L_{y}$ with $L_{x}=157.37\,\sigma$ and $L_{y}=136.28\,\sigma$ fixing the dimensionless density to $\rho=0.6$. The phase behavior is studied varying $Pe$ and $\eta$ values. At $\eta\lesssim 0.06$, the system undergoes phase separation with increasing $Pe$. We characterize the phase coexistence using the density distribution function $P(\rho)$ and the distribution of local hexatic order $P(\psi_{l})$, where $\psi_{l}({\bf r})=|\vec{\psi}_{l}({\bf r})|$ are calculated over local coarse-grained volumes of $\ell_{x}\times\ell_{y}=5.25\,\sigma\times 4.5\,\sigma$ around ${\bf r}$. The system is relaxed up to $10^{7}\delta t$, before analyses are performed using $4000$ configurations saved in intervals of $10^{4}\delta t$ obtained from simulations over further $4\times 10^{7}\delta t$.

III.3.1 Re-entrant phase separation with $Pe$

In the regime of small orientational noise $\eta\lesssim 0.05$, a re-entrant phase-separation is observed with increasing $Pe$. A detailed characterization is presented in Fig.4 for a fixed $\eta=0.04$. At low activities, e.g., $Pe=5$, $P(\rho)$ exhibits a single maximum at the mean density $\rho=0.6$ corresponding to the homogeneous fluid (see Fig.4($a$) ). With increasing $Pe$, initially the distribution broadens, capturing the increased fluctuations in density. As the activity is further increased to $Pe=10$, clusters form, and the system phase separates into high and low density regions.

The phase separation is most pronounced at $Pe=20$. The phase separation persists up to $Pe=50$. As activity is increased further to $Pe=60$, the high density clusters get destabilized and the bimodality disappears identifying remelting. The new unimodal distribution is significantly broad with respect to the fluid before the onset of phase separation at small $Pe$. It shows a shallow peak near $\rho=0.48$ (Fig.4($a$) ). This is due to the larger density fluctuations associated with the stronger activity. The peak position(s) of the local density distributions for $\eta=0.04$ are plotted as a function of $Pe$ in Fig.4($b$). This clearly demonstrates phase coexistence over an intermediate regime of $Pe$, between $10$ to $50$.

Associated with the clustering, a local hexatic order $\psi_{l}({\bf r})$ emerges. Three representative probability distributions $P(\psi_{l})$ at $Pe$ values corresponding to the two homogeneous phases and the phase coexistence region, are shown in Fig.4($c$). At $Pe=5,\,60$ the distributions show a single maximum at a low hexatic order $\psi_{l}\approx 0.15$. At intermediate $Pe$ values, e.g., $Pe=20$, a bimodality in $P(\psi_{l})$ is observed, showing emergence of a new large $\psi_{l}$ peak near $\psi_{l}=0.94\pm 0.02$, associated with the emergence of the high density hexatic clusters. However, even within the phase coexistence region, $Pe=20$, the larger peaks of the distributions remain at a small $\psi_{l}$ value, reflecting the fact that the homogeneous fluid covers the larger spatial fraction of the system. The peak positions of the probability distributions $P(\psi_{l})$ are plotted in Fig.4($d$). This again captures the re-entrant phase separation, as in Fig.4($c$). With increasing $Pe$, the system shows a re-entrant transition from a single fluid to a fluid-hexatic coexistence that finally returns to a single fluid.

III.3.2 Phase separation with $\eta$

To demonstrate the dependence of the phase separation on the strength of orientational fluctuations $\eta$, we consider the system at a sufficiently high activity $Pe=40$ and vary $\eta$. As we have seen above, at this activity the system remains phase separated at $\eta=0.04$. In Fig.5($a$) we show the variation of $P(\rho)$ with $\eta$. At low $\eta$ ($=0.02$, blue up triangles), the system remains phase separated, captured by the bimodality of $P(\rho)$. With increasing $\eta$, the peaks in the distribution slowly come together and merge at $\eta=0.05$ (red circles). For larger $\eta$, the density distribution remains unimodal (e.g., $\eta=0.12$, black squares) with the peak at the mean density $\rho=0.6$, signifying a homogeneous fluid phase. The peak positions of $P(\rho)$ are plotted against $\eta$ in Fig.5($b$). This shows how the coexistence lines merge together as a function of $\eta$ at a critical point. This behavior is similar to the merger of liquid-gas coexistence lines as the critical temperature is approached from below. The detailed properties of such a critical point has been explored for the ABPs recently Siebert2018 , and remains to be studied in the current context in future. The distributions of the local hexatic order $P(\psi_{l})$ show pronounced multi-modality at the phase coexistence. The peak near $\psi_{l}=0.96$ captures the significantly large hexatic order of the active clusters. The other peak near $\psi_{l}=0.17$ signifies the coexisting homogenoeus fluid. In addition, we observe a peak at $\psi_{l}=0$, which corresponds to the voids generated at the cost of clustering. The non-zero peaks of $P(\psi_{l})$ are plotted against $\eta$ in Fig.5($d$), recapturing the loss of phase coexistence with increasing $\eta$, as in Fig.5($b$).

Using the above analysis it is possible to construct a combined phase- coexistence diagram in the $Pe-\eta$ plane as is shown in Fig.6. The $\rho_{m}$ values denote the binodals extracted from the peaks of the probability distribution of the local density. The apex of the plot at large $\eta=0.056$ denotes the critical point below which phase coexistence becomes possible. The lines indicate the boundary of the onset of phase coexistence. As is clear from the plot, the coexisting densities separate out more strongly at lower $\eta$. At a given $\eta$, after the onset of active phase separation appears at a low $Pe$, coexisting densities first separate out to finally merge as $Pe$ is increased. At lower $\eta$, large activity supports formation of flocks. As a result this merger happens at a larger $Pe$. Moreover, lowering of $\eta$ reduces the orientational diffusion supporting MIPS at small $Pe$, thereby reducing the $Pe$ value at which phase separation appears. The flocking and MIPS mechanisms are discussed further in the following section.

III.3.3 Cluster formation

In the absence of alignment, volume exclusion induces slow down of active motion with increasing local density generating a positive feedback mechanism that mediates MIPS. However, the clustering we observe is deep within an orientationally ordered nematic phase. This raises the question whether the MIPS- like mechanism still plays any role in clustering, or the phase coexistence is due to the onset of flocking related to orientational fluctuations. It was argued before that alignment for polar disks reduces rotational diffusion, and thereby may favor MIPS Martin-Gomez2018 ; Sese-Sansa2018 . On the other hand, steric alignment of self propelled rods beyond an asymmetry parameter produces destabilizing torques to break MIPS clusters Grossman2019 ; Shi2018 ; Weitz2015 . We note that, for our self propelled disks in the nematic phase with antiparallel active orientations, crowding effects may slow down the propulsion speed. In fact, as we show below, this is the dominant mechanism of clustering in the small $Pe$ regime.

Here we compute the self-propulsion speed $v_{a}(\rho)$ reduced by the inter-particle repulsion at a local density $\rho$. We proceed by calculating the coarse- grained quantity $(1/n_{r})\mid\sum_{i\in r}^{n_{r}}[v_{0}{\bf n}_{i}-\mu\nabla_{i}U(r_{ij})\,]\cdot{\bf n}_{i}\mid$ in cells of volume $\ell_{x}\times\ell_{y}$ having $n_{r}$ particles Sese-Sansa2018 , such that the instantaneous density is $(\pi/4)n_{r}\sigma^{2}/\ell_{x}\ell_{y}$. Within the MIPS paradigm, it is expected that $v_{a}(\rho)=v_{0}(1-a\rho)$ where $a$ depends on $v_{0}$, a stall time at collision, and a scattering cross-section Stenhammar2013 . Fig.7 shows the variation of $v_{a}(\rho)$ along the phase separation boundary. In the small $Pe$ regime, between $Pe=2$ to $20$, $v_{a}(\rho)$ decreases linearly with $\rho$ with a single slope $a\approx 1$, in agreement with the MIPS expectation. However, at higher $Pe$, the behavior starts to get dominated by flocking, shown by the slow down of the decay of $v_{a}(\rho)$ at higher $\rho$. Note the cross-over to a second slope $a\approx 1/2$ at large $\rho$ for $Pe>20$. This happens as a cluster starts to get an effective polar alignment such that the active velocities of particles do not completely block each other. Similar cross-overs have been observed earlier for ferromagnetically aligning ABPs Sese-Sansa2018 .

At small $Pe$, colliding particles slow down as in MIPS. The size of a hexatic cluster is maintained by the absorption of particles from the surrounding fluid having density $\rho_{l}$, with a rate $r_{\rm on}=\rho_{l}v_{0}/\pi$ and evaporation from the cluster surface due to reorientation of the direction of self-propulsion with a rate $r_{\rm off}=\kappa D_{\theta}(\eta)/\sigma$, where $\kappa$ denotes the total number of particles leaving the cluster in each escape  Redner2013 . The orientational relaxation is controlled by the noise strength $\eta$, such that the related rotational diffusion constant $D_{\theta}\sim\eta^{2}$ (see Appendix-C). The mean density of the system in the phase coexistence regime can be expressed as $\rho=f\rho_{h}+(1-f)\rho_{l}$, where $\rho_{h}$ ($\rho_{l}$) denotes the mean density of the hexatic (liquid) fraction, and $0\leq f<1$ denotes the amount of hexatic fraction. At the onset of phase coexistence, $f=0$ and $\rho=\rho_{l}$. Thus, equating $k_{\rm on}=k_{\rm off}$ and setting $\rho_{l}=\rho$, we get $v_{0}=\pi\kappa D_{\theta}(\eta)/\sigma\rho$. This leads to the condition $Pe=v_{0}\sigma/D=\pi\kappa D_{\theta}(\eta)/\rho D$. Using $D_{\theta}\sim\eta^{2}$ we obtain the prediction for MIPS- like phase boundary $\eta\sim Pe^{1/2}$, plotted in the phase diagram Fig.1($a$). This captures the behavior of the low $Pe$ boundary for the onset of active phase separation.

To illustrate the melting process of hexatic clusters at high $Pe$, in Fig.8, we show the evolution of a steady state configuration after quench from within the phase coexistence regime at $\eta=0.04$ and $Pe=10$ to $Pe=70$, an activity beyond the remelting boundary. The plots show velocity orientation (color code) and speed (arrow) of particles, e.g., right moving particles (orientation [math]) are denoted by green and the left moving ones by red (orientation $\pm\pi$). If counter moving polar patches do not block each other, they slide past. In other places, e.g., in the largest cluster in the upper- left corner at the initial configuration $t=0\,\tau$ (Fig.8), a left moving patch creates channel and squeezes through the majority right moving particles. The polar patches eventually detach, destabilizing the big clusters to smaller clusters. Similar ejection of polar clusters were observed before in self propelled rods Weitz2015 . The ejection process continues to finally melt the hexatic clusters, e.g., see the last configuration after a time lapse of $5\tau$ in Fig.8. The detailed dynamics can be observed in the movie added in the supplemental material 111The video shows the time evolution of the local density and velocity fields over $5\,\tau$ after a quench from $Pe=10$ to $Pe=70$ at $\eta=0.04$. The color code in the video follows that in Fig.8..

III.4 Topological defects

In this section we investigate the role of the topological defects in the melting of the active hexatic clusters. We define the coordination number of each particle $n_{v}$ as the number of Voronoi neighbours it has. The $n_{v}\neq 6$ particles constitute the topological defects. Within hexatic clusters, particles with $n_{v}=5$ and $n_{v}=7$ stay as bound pairs. Such $5-7$ bound pairs, well separated from each other, identify the presence of dislocations in hexatic. For example, see the upper left inset in Fig.1($a$) displaying several $5-7$ bound pairs in a portion of a hexatic cluster of the active system. Unbinding of such $5-7$ pairs, generate free disclinations. As we show, interconnected chains of $5-7-5-7..$ defects form grain boundaries, held by local polar patches of the hexatic clusters.

We analyze the variation of the fraction of particles belonging to different defect types within the hexatic clusters, with increasing $\eta$ at $Pe=40$ Qi2014 . A typical configuration at $\eta=0.02$ is shown in Fig.9. In Fig.9($a$) we display the hexatic order associated with all the $k$ particles $\mid\vec{\psi}_{6}^{k}\mid^{2}$. The large clusters have an overall large hexatic order $\mid\vec{\psi}_{6}^{k}\mid^{2}\approx 1$. The local drops in the values of $\mid\vec{\psi}_{6}^{k}\mid^{2}$ within the large clusters are associated with formation of topological defects. In Fig.9($b$) we show the three different defect types separately – dislocations (red points), disclinations (blue points) and grain boundaries (string of black points) within the three largest clusters colored with cyan (largest cluster), grey (2nd largest cluster) and pink (3rd largest cluster). To understand the localization of the grain boundaries, we display the projection of the local polar orientation $\hat{n}_{k}=(\cos\theta_{k},\sin\theta_{k})$ of the particles onto the orientation of instantaneous overall nematic order ${\hat{s}}={\bf s}/\mid{\bf s}\mid$ in Fig.9($c$). A comparison between Fig.9($c$) with ($b$) shows that the grain boundaries appear at the interface of anti-parallel polar patches of activity. At high activity, such as $Pe=40$, the opposing active push from these neighboring polar regions leads to the formation of the grain boundaries within the clusters.

To quantify the behaviour of the topological defects across the melting transition, we compute the fraction $\psi_{d}$ of particles associated with the dislocations, grain boundaries and disclinations within the largest cluster of a given configuration. The averaging is performed over more than $3000$ well separated configurations in the steady state. In Fig.10 we show the variation of these three average quantities with increasing $\eta$, at $Pe=40$. In the active clusters, all the disclinations are thrown out to the hexatic- fluid interfaces. The fraction of the disclination within the hexatic cluster, thus, remains largely unchanged. The fraction of dislocations increases with $\eta$, as well as the fraction of grain boundaries. Note that the dislocations characterize hexatic phase, and do not destabilize it. Thus the melting of the active hexatic at large $Pe$ is mediated by grain boundary formation along the interface of oppositely aligned polar patches within a cluster.

III.5 Cluster size and dynamics

As has been mentioned before, the coexistence at the active phase separation is characterized by the presence of large correlated clusters. They are identified using a clustering algorithm, defining particles within a separation $\sigma$ to be part of a single cluster Allen1989 . The number of particles belonging to a cluster determines the cluster size $c$. Three typical steady state cluster size distributions $P(c)$ in the homogeneous fluid phase and at the active phase coexistence are shown in Fig.11.

The distribution of small clusters are approximately independent of the phase, and the weight decays as $c^{-2.5}$. The decay of the distribution for large clusters, on the other hand, changes qualitatively with the change in phase. In the homogeneous fluid phase, it decays exponentially. In contrast, in the regime of phase separation, the tail of the distribution shows a power law decay $P(c)\sim c^{-1.5}$. This appears in addition to a small but approximately uniform weight associated with the largest clusters. The change of the tail of the cluster size distribution from exponential decay to a power law is qualitatively similar to that observed in the conserved mass aggregation model (CMAM) Majumdar1998 , however, the exponents in our case differ significantly. Moreover, for the largest clusters structural ordering competes with ejection of polar fractions, not allowing stabilization of infinite-aggregates, unlike in the CMAM.

A measure of the spatial extension of the clusters is the correlation length obtained from the two point density correlation function $g(r)=\langle\rho(0)\rho(r)\rangle/\langle\rho\rangle^{2}$. In Fig.12 the variation of $g(r)$ is shown in the fluid and the phase coexistence region. The typical correlation length, $l_{\rho}$, is estimated as the value of the separation $r$ beyond which the correlation decays to $1\pm 0.01$. As is shown in Fig.12, in the regime of phase separation ($Pe=10,\,\eta=0.02$) one gets a large $l_{\rho}$ ($\approx 35\,\sigma$) corresponding to the large extension of the clusters. However, in the homogeneous fluid phase the correlation length becomes small ($l_{\rho}\approx 5\sigma$). If the appearance of phase separation is not a finite size effect, the cluster size $l_{\rho}$ should remain a finite fraction of the system size $L_{0}=(L_{x}\times L_{y})^{1/2}$, in the thermodynamic limit of infinite $N$. To test this we perform a system size scaling. As we show in the inset of Fig.12, indeed, the ratio $l_{\rho}/L_{0}$ does not decrease with system size $N$.

The collective motion is dominated by flock- like behavior at large $Pe$ as well as small $\eta$, at which alignment dominates. This consists of (i) structurally ordered stripes sliding past each other, (ii) jamming of polar clusters, (iii) which eventually break via formation of counter-propagating channels. The behaviors are shown in Fig.13($a$)-($c$) with the help of local vector fields identifying the orientations and amplitudes of the locally coarse grained velocity of particles. The coarse- graining is performed over cells of size $\ell_{x}\times\ell_{y}$, and the velocity fields are calculated from displacements obtained over $0.1\,\tau$. The orientations of velocity are further highlighted by color codes explained in the figure. The snapshots are shown at three different well separated time instances in the steady state, deep inside the phase coexistence regime. For the ease of illustration, we have chosen $Pe=10$ and $\eta=2\times 10^{-4}$, as the clusters at such a low $\eta$ are system spanning in extension and have long life-times.

In Fig.13($a$) we plot the first steady-state snapshot of local velocity fields. It shows two large counter propagating clusters that glide past each other, in the up and down directions ($\approx\pm\pi/2$). The velocity field at the interface of the two clusters vanishes, before changing direction. With time, the velocity field from one domain may penetrate the other due to dynamical fluctuations along the interface. In Fig.13($a$), such events are identified by a bulge of vectors marked within a circle. This generates instability towards the formation of lanes. At a later time, $t=41\tau$, as is shown in Fig.13($b$), particles in this region starts to flow leftward (direction $\approx\pm\pi$). They create path for other particles to trickle along. As time progresses, this local trickling takes the form of a stream of particles forming counter propagating lanes, as is shown in Fig.13($c$) at $t=216.6\,\tau$. The crossovers between these different regimes happen faster at higher $Pe$ and $\eta$.

IV Discussion

We have studied the combined effect of volume exclusion and nematic alignment on collective behavior of active polar particles. Using numerical simulations, we have presented a phase diagram encapsulating the structural and orientational transitions with changing activity $Pe$ and noise $\eta$. In the equilibrium limit of vanishing activity our system behaves as a high density fluid. On the other hand, the limit of vanishing interaction reduces it to the nematically aligning point particles Ginelli2010 . The nematic-isotropic transition in activity obtained by increasing $\eta$ does not change the homogeneous nature of the fluid phase of the excluded volume particles. This behavior is in striking contrast to the properties observed for point particles Ginelli2010 , where the nematic-isotropic transition was accompanied by phase separation. The volume exclusion in the current model, in addition to a significant rotational diffusivity, $D_{\theta}\sim\eta^{2}$, even at the nematic-isotropic transition, suppresses such fluctuations.

However, we do find clustering transition in our system, but at a lower strength of $\eta$ associated with larger persistence, $D_{\theta}^{-1}$. As we have shown, the competition of persistence and volume exclusion led to a kinetic slowdown with density. This provided a MIPS mechanism towards phase separation at small $Pe$. The nematic fluid undergoes a transition from the homogeneous phase to an active phase separation characterized by the formation of high density hexatic clusters. We found a re-entrant fluid- phase coexistence- fluid transition with increasing $Pe$. At the small $Pe$ regime, volume exclusion dominates – the formation of clusters is mediated by the slow-down of colliding particles providing a positive feedback, as in MIPS. The phase boundary in this regime follows a form $\eta\sim Pe^{1/2}$, that could be motivated using the MIPS phenomenology.

With increasing $Pe$, the opposing push from the counter propagating polar sub-regions of a nematic cluster leads to formation of grain boundaries. Once the active push gets large enough, opposing polar regions squeeze past each other forming lanes. At higher $Pe$, the clusters eject polar patches that form co-moving flocks. At high enough $Pe$, this breaking and ejection of polar patches continue to dissolve larger clusters into smaller ones, finally leading to remelting. As in the high $Pe$ regime, at the lowest orientational noise strengths, alignment dominates leading to flock- like collective motion that consists of sliding clusters, jamming, and lane formation.

An experiment on active colloidal disks, recently, examined the combined impact of volume exclusion and alignment in the phase separation kinetics at high activity Linden2019 . Similar experiments, may potentially verify our predictions. A coarse grained hydrodynamic description combining active alignment and density dependent slowdown remains to be developed, and may provide further insight into the phase transitions of such systems.

Acknowledgements

The simulations were performed on SAMKHYA, the high performance computing facility at IOP, Bhubaneswar. DC thanks ICTS-TIFR, Bangalore for an associateship, and SERB, India, for financial support through grant number EMR/2016/001454.

Appendix A Quasi long range order

Here we demonstrate the nature of the nemtaic order before the nematic-isotropic transition. To this end, we perform a system size analysis by calculating the scalar order parameter $S$ for system sizes $N=128$ to $9216$, at a fixed $Pe=40$, and varying $\eta=0.02-0.10$ (Fig.14). With increasing $N$, we find an algebraic decay $S=S_{0}\,N^{-\alpha}$, where both $S_{0}$ and the exponent $\alpha$ depend on $\eta$. Even at the smallest $\eta$, the algebraic decay of the order parameter with system size is evident (see inset of Fig.14). This is a signature of the quasi- long ranged order (QLRO) of the active nematic phase mentioned in Sec.III.2. The exponent $\alpha$ increases with $\eta$.

Appendix B Scaling of the $S$-$\eta$ curves

We identify three different regimes, depending on how the function $S(\eta)$ behaves at different $Pe$ (Fig.2). For $Pe\leq Pe_{\rm min}(=2.0)$, the curves fall on top of each other implying that a minimum activity $Pe_{\rm min}$ is required for the system to show effects of $Pe$ in the macroscopic scale. This $Pe_{\rm min}$ vanishes with increase in system size(data not shown). Increase in $Pe$, allows a single particle to cover larger space, nematically aligning with larger number of particles over a span of time. Thus a larger orientational noise $\eta$ is required to destabilize the nematic phase, increasing the transition noise $\eta_{c}$ with $Pe$. However, for a given density, there is a limit up to which $Pe$ can suppress the nematic to isotropic transition. This limit is decided by the available free space, determined by the mean density. For $Pe>Pe_{\rm max}(=20)$, again, the $S(\eta)$ curves fall on top of each other, making $\eta_{c}$ independent of $Pe$. In between $Pe=Pe_{\rm min}$ and $Pe_{\rm max}$, the variations of $S$ with $\eta$ show data collapse when plotted against $(\eta-\eta_{c})/\eta_{c}$ (Fig.15). Within this restricted window of activity, we find a dependence of transition point $\eta_{c}$ with $Pe$ that shows an approximate scaling form $\eta_{c}(Pe)\sim Pe^{0.2}$ (inset in Fig.15).

Appendix C Orientational relaxation

In the presence of the nematic alignment, the active orientation of the particles undergo orientational relaxation controlled by the noise of strength $\eta$. In Fig.16 we show that the orientational diffusion constant increases as $D_{\theta}\sim\eta^{2}$, and remains essentially unchanged with variation of $Pe$.