An Infinite Set of Integral Formulae for Polar, Nematic, and Higher Order Structures at the Interface of Motility-Induced Phase Separation

The study of active matter is crucial to our understanding of diverse living matter and driven synthetic systems [1, 2]. In addition to being paramount to our quantitative description of various natural and artificial phenomena, much novel universal behaviour has also been uncovered in active matter systems in the hydrodynamic limits, which ranges from novel phases [3, 4, 5, 6, 7, 8, 9, 10] to critical phenomena [11, 12]. Besides the hydrodynamic limits, interesting emergent phenomena also occur in the microscopic and mesoscopic scales. In the case of motility-induced phase separation (MIPS) [13, 14, 15, 16, 17], one of these emergent phenomena is the polar-nematic ordering behaviour at the liquid-gas interface (Fig. 1) [18, 19, 20, 21, 22]. Understanding this phenomenon will be integral to modelling quantitatively diverse interface-related properties of MIPS that include nucleation [23], wetting [24, 25], negative interfacial tension [26, 27, 28], and reverse Ostwald ripening [29, 30].

Although the interfacial ordering ultimately emerges from the many-body interactions of the particles, similar pattern in fact already occurs in a system of non-interacting self-propelled particles (i.e., an ideal active gas) around an impenetrable wall – particles stuck at the wall tend to point towards it (hence the system is polar), and particles right outside the wall are predominantly moving parallel to the wall (hence nematic), since they constitute particles whose orientations have just rotated enough to be pointing away from the wall [31, 32]. This revelation suggests that the polar-nematic pattern observed in MIPS can be studied using an effective model in which all particle-particle interactions are captured by a velocity field that depends purely on the system’s local configurations, such as the local particle density. Implementing this task is the goal of this paper. Specifically, I first derive the equation of motion (EOM) of the many-particle distribution function under the assumptions that the effective speed and the effective pressure depend solely on the local properties of the particle density. I then obtained the steady state equations that describe the polar, nematic, and higher order tensor fields of the system, and from these an infinite set of integral formulae (IF), one for each tensor field. Moreover, I test these IF on published numerical data obtained from direct particle dynamics simulation [18], and demonstrate that the IF remain valid to a remarkable extent. Finally, I showed in the appendix that half of the IF hold exactly for a wide class of active Brownian particle systems that exhibit MIPS.

In two dimensions, an archetype of active systems exhibiting MIPS consists of a collection of $N$ polar self-propelled particles in a finite volume $V$, such that the particles interact purely sterically via a short-ranged repulsive potential energy $U$, whose exact form is unimportant. Each particle also exerts a constant active force, $f_{a}$, that points to a particular orientation denoted by $\hat{\mathbf{n}}\equiv\cos\Theta\hat{\mathbf{x}}+\sin\Theta\hat{\mathbf{y}}$, and the angle $\Theta$ itself undergoes diffusion with the rotational diffusion coefficient $D_{\rm R}$. Furthermore, the particles can potentially perform Brownian motions characterized by the diffusion coefficient $D_{\rm T}$. Overall, the microscopic equations of motion (EOM) of the system are therefore

where the indices $i,j$ enumerate the particles ($i,j=1,2,\ldots,N$), $\mathbf{R}_{i}(t)$ is the $i$-th particle’s position at time $t$, $\eta$ is the damping coefficient in this overdamped system, and $\mathbf{g}^{\rm T}$’s and $g^{\rm R}$’s are independent Gaussian noise terms with zero means and unit variances, e.g., $\langle g^{\rm R}_{i}(t)\rangle=0$ and $\langle g^{\rm R}_{i}(t)g^{\rm R}_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$.

Since the particles’ identities are irrelevant, we can focus instead on the temporal evolution of the $N$-particle distribution function:

where the angular brackets denote averaging over the noises.

As mentioned before, the goal here is to account for all particle-particle interactions and fluctuations effectively through a speed function and a pressure function that depend solely on the local density: $\rho(\mathbf{r},t)=(2\pi)^{-1}\int{\rm d}\theta\psi(\mathbf{r},\theta,t)$. Specifically, the model EOM of $\psi$ is assumed to be of the form:

where $\mathbf{v}$ is the ‘velocity field’ taken to be

and $u$ and $P$ are the effective density-dependent speed and pressure functions, respectively. Note that the ‘effective pressure’ $P$ does not necessarily correspond to the mechanical pressure in a thermal system, or the Irving-Kirkwood pressure defined for passive systems [33]. Instead, $P$ should be viewed as a density-dependent scalar function whose spatial gradient contributes to the velocity field as described in Eq. (4). I further note that allowing the effective functions $u$ and $P$ to depend also on the spatial derivatives of $\rho$ (e.g., $\nabla^{2}\rho$) will not alter the key results here. In the Appendix, I will relate the approximation adopted here to a formally exact set of hierarchical EOM.

Given the model equations (3,4), we are now ready to systematically construct a reduced set of EOM in which the fields of interest depend on $\mathbf{r}$ and $t$, but not $\theta$. Specifically, these fields will be $m$-th rank tensors of the form:

where the Greek letters enumerate the spatial coordinates and the subtraction of $\rho$ in the integrand ensures that these tensor fields vanish if $\psi$ is isotropic in $\theta$. For instance, the unnormalized polar field ${\bf M}(\mathbf{r},t)$ and nematic field $\mathbf{Q}$ are the first-rank and second-rank tensors, respectively:

These fields are unnormalized since they are not normalized by the density, as opposed to, e.g., the normalized polar field ${\bf M}^{\rm N}=(2\pi\rho)^{-1}\int{\rm d}\theta\hat{\mathbf{n}}\psi$.

From (3,4), the EOM of the density and polar fields can be obtained:

where $\partial_{\alpha}\equiv\partial/\partial r_{\alpha}$ and repeated indices are summed over. The EOM of higher order tensor fields can be derived similarly, as illustrated next.

I will now focus exclusively on the steady state in the phase separated regime, where there is a single liquid-gas interface located around $x=0$. By construction, spatial variations occur only along the $x$ direction, and we thus need only to consider tensor fields with all subscripts being $x$, i.e., tensor fields of the form $T^{(m)}_{x\cdots x}$. For later convenience, I will further denote $T^{(m)}_{x\cdots x}$ by $O_{m}$ when $m$ is odd and by $E_{m}$ when $m$ is even. Specifically, from the definitions of the tensor fields (5), we have

where

At the steady state, the equations satisfied by the tensor fields are, from (3,4):

The L.H.S. of (10) follow from the identities below:

and for the L.H.S. of (10d), I have additionally used the identity:

which can be readily derived from the definition of $K_{m}$ (9).

Since $\psi$ is isotropic in both $x$ and $\theta$ deep in the liquid and gas phase, $M_{x},Q_{xx},O_{m}$, $E_{m}$, and ${\rm d}P/{\rm d}x$ all vanish when $|x|\gg 0$. Therefore, by integrating both sides of (10) over the whole $x$ domain, we arrive at the followings:

where

In fact, one more simplification can be made: starting with (13a) and using (13b), one can readily prove by mathematical induction that

The integral formulae (IF) expressed in (13c) and (15) are the key results of this paper. In the appendix, I will further show that for even $m$, the IF (13c) are in fact exact in a wide class of active Brownian particle systems.

Introducing the following notation:

the plot of $I_{m}/I_{1}$ vs. $m$ is shown in Fig. 2 (red circles). Note in particular the slow decay of $I_{m}/I_{1}$ for odd $m$, which goes asymptotically to 0 as $m^{-1/2}$ (9).

As a test of the validity of the IF to realistic active systems exhibiting MIPS, I now re-analyse the numerical results obtained from simulating a microscopic model of MIPS [18]. Specifically, the simulation is on a system of 300 particles confined in a 2D channel of height $10\sin(\pi/3)$ and length $50$. A periodic boundary condition is used for the vertical direction, a hard wall is placed on the left of the channel, and particles exiting the right wall will stay at the right wall but with its orientation randomized. These particles interact via short-ranged repulsive interactions in the form of the Weeks–Chandler–Andersen potential [34]:

Other parameters are: $\eta=1,f_{a}=100,D_{\rm T}=0$, and $D_{\rm R}=3$.

At the steady-state, the density and various tensor field profiles are shown in Fig. 3. Given that the liquid phase is on the left (Fig. 3a), the odd-rank tensor fields $O_{m}$ are all negative around the interface (Fig. 3b), as the particles’ orientations are mostly pointing towards the liquid phase (black arrow in Fig. 1). For the even-rank tensor fields $E_{m}$ (Fig. 3c), they first become positive when approaching the interface from deep in the liquid phase, indicating that their orientations are predominantly horizontal, consistent with the orientation of the polar field. These fields then subsequently become negative as one exits the interface, indicating that the particles’ orientations become predominantly vertical (red double arrow in Fig. 1). The “oscillatory” feature of $E_{m}$ around the interface is of course necessary in order for the overall integrals of $E_{m}$ to be zero (13c). Note also that while all tensor fields go to zero as $m$ goes to infinity, they apparently do so rather slowly (Fig. 3(b)-(c)). This is again consistent with the analytical result that $I_{m}\sim m^{-1/2}$ for odd $m$. Overall, this slow decays indicates that a quantitative theory focusing on the interfacial region will need to incorporate a large number of higher order tensor fields. Finally, performing the integrals to this set of data as prescribed in (16), the numerical result (black crosses and pluses in Fig. 2) demonstrates that the IF are indeed valid to a great extent.

Starting with a generic microscopic model of active particle system with only repulsive interactions, and adopting the simplifying assumptions that all many-body physics and fluctuations can be captured by an effective speed function and an effective pressure that depend solely on the local particle density, I have derived a set of EOM that describes the polar, nematic, and higher order tensor fields of the system. Focusing then on the steady state with a single liquid-gas interface, I obtained an infinite set of IF, one for each tensor field. I then tested these IF on data obtained from direct particle dynamics simulation of a microscopic model of MIPS, and found that the IF remain valid to a high accuracy. Furthermore, I showed in the appendix that half of these IF (even $m$) are in fact exact for a wide class of active Brownian particle systems. The overall validity of all the IF suggests that the model assumptions are appropriate when studying interfacial properties of MIPS.

This work opens up a number of interesting future directions: (1) For finite systems, the IF can be made more precise by incorporating the exact boundary conditions at the two limits of the spatial integrals. (2) The IF can be readily modified to sum rules for lattice models of MIPS [35, 36, 37, 38]. (3) It would be interesting to consider how these IF can be generalized to a system of active particles with alignment interactions. Indeed, the various types of travelling bands observed in polar active matter [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49] are reminiscent of the soliton solutions in the Korteweg–De Vries equation, which also admits an infinite number of integrals of motion. (4) Finally, the analytical treatment here provides a basis for developing a quantitative theory that elucidates how the interfacial ordering impacts upon other interface-related phenomena in MIPS, such as the Gibbs-Thomson relation [18], wetting [24, 25], and reverse Ostwald ripening [29, 30].