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

Motility-induced phase separation (MIPS) is a purely non-equilibrium phenomenon in which self-propelled particles phase separate without any attractive interactions. One surprising feature of MIPS is the emergence of polar, nematic, and higher order structures at the interfacial region, whose underlying physics remains poorly understood. Starting with a model of MIPS in which all many-body interactions are captured by an effective speed function and an effective pressure function that depend solely on the local particle density, I derive analytically an infinite set of integral formulae (IF) for the ordering structures at the interface. I then demonstrate that half of these IF are in fact exact for a wide class of active Brownian particle systems. Finally, I test the IF by applying them to numerical data from direct particle dynamics simulation and find that all the IF remain valid to a great extent.


Introduction
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 noninteracting self-propelled particles (i.e., an ideal active gas) around an impenetrable wall Figure 1. Schematic of the generic polar-nematic ordering observed at the interface of MIPS. The shade of blue depicts the particle density of the system. The particles' orientations are isotropic in the bulks of both liquid and gas phases. As one approaches the interface from the liquid (condensed) phase, the average orientation of the particles at the interfacial region becomes predominantly polar (pointing towards the liquid phase as indicated by the black arrow), and then predominantly nematic (i.e., particles' orientations tend be approximately vertical, as indicated by the red double-arrow.) Higher order tensorial structures beyond the polar and nematic fields also emerge at the interfacial region. Here, using an model in which the effective speed of the particles and the effective local pressure in the system depend solely on the local density, an infinite set of integral formulae is obtained for all these tensorial structures.
-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 polarnematic 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.

Effective model
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 byn ≡ cos Θx + sin Θŷ, and the angle Θ itself undergoes diffusion with the rotational diffusion coefficient D R . Furthermore, the particles can potentially perform Brownian motions characterized by the diffusion coefficient D 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, . . . , N ), R i (t) is the i-th particle's position at time t, η is the damping coefficient in this overdamped system, and g T 's and g R 's are independent Gaussian noise terms with zero means and unit variances, e.g., g R i (t) = 0 and g R i (t)g R j (t ) = δ ij δ(t − t ). 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: ρ(r, t) = (2π) −1 dθψ(r, θ, t). Specifically, the model EOM of ψ is assumed to be of the form: where v is the 'velocity field' taken to be v(r, θ, t) = u(ρ(r, t))n(θ) − ρ(r, t) −1 ∇P (ρ(r, t)) , 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 ρ (e.g., ∇ 2 ρ) 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 r and t, but not θ. Specifically, these fields will be m-th rank tensors of the form: where the Greek letters enumerate the spatial coordinates and the subtraction of ρ in the integrand ensures that these tensor fields vanish if ψ is isotropic in θ. For instance, the unnormalized polar field M(r, t) and nematic field 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 M N = (2πρ) −1 dθnψ. From (3,4), the EOM of the density and polar fields can be obtained: where ∂ α ≡ ∂/∂r α and repeated indices are summed over. The EOM of higher order tensor fields can be derived similarly, as illustrated next.
3. Integral formulae at the steady state 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···x . For later convenience, I will further denote T 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: dθ cos m θ∂ 2 θ ψ = dθ (m − 1)m cos m−2 θ − m 2 cos m θ ψ , 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 ψ is isotropic in both x and θ deep in the liquid and gas phase, M x , Q xx , O m , E m , and dP/dx all vanish when |x| 0. Therefore, by integrating both sides of (10) over the whole x domain, we arrive at the followings: , and the even-rank tensor fields E m (c), with respect to x at the steady state of a system of self-propelled particles with repulsive interactions that exhibits MIPS. The simulation data used for these curves are from a previous study [18].
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).

Testing the integral formulae on a microscopic model of MIPS
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(π/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: η = 1, f a = 100, D T = 0, and D 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 ∼ 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.

Summary & Outlook
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].
Acknowledgments I thank an anonymous referee for their suggestion to look into the exactness of the IF, which led to the realisation that half of the IF can be shown to be exact from first principles without resorting to the approximation adopted in the main text.
Appendix A. Relating the approximation to an exact hierarchical EOM In this appendix, I will relate the key approximation made in the main text to a formally exact set of EOM of the tensor fields. To do so, I will start with the fluctuating Nparticle distribution function: where the superscript (f ) specifies that ψ (f ) is a fluctuating quantity (hence different from ψ (2)) due to the lack of noise averaging (hence without the angular brackets). Following standard procedures [15,50,51,52,53,20], the model EOM of ψ (f ) is of the form: u 0 ≡ f a /η is the constant 'active speed', and with ρ (f ) being the 'fluctuating' density: , the steady state of the average tensor fields can be obtained similarly as before: and for odd m > 1: (A.7) and for even m > 2: (A.8) The above steady state equations are similar to those in the main text (10), except for two crucial differences: (1) the active speed u 0 is a constant here instead of being dependent on the density ρ, and (2) the previous effective pressure gradient terms are now replaced by the correlation functions between w Deep in the bulks of the liquid and gas phases, we expect that that the correlation ψ (f ) (r, θ)ρ (f ) (r ) cannot distinguish left from right. Hence where H ± are constants at the x → ±∞ limits. While the form in (A.14) expectedly satisfies the symmetry in (A.12), I cannot show that the integrals on the L.H.S. are exactly proportional to cos θ, and their demonstrations from first principles will be a very interesting future direction.