Systematic extension of the Cahn-Hilliard model for motility-induced phase separation

TL;DR

This paper derives a higher-order extension of the Cahn-Hilliard model to better describe motility-induced phase separation in active particles, establishing a mathematical link and validating the model through numerical comparisons.

Contribution

It introduces a systematic higher-order extension of the Cahn-Hilliard model for active phase separation, connecting it to the mean-field MIPS model and previous phenomenological models.

Findings

Leading order CH model matches full continuum model near MIPS onset

Higher order corrections become relevant beyond certain control parameters

Extended CH model aligns well with recent phenomenological models

Abstract

We consider a continuum model for motility-induced phase separation (MIPS) of active Brownian particles [J. Chem. Phys. 142, 224149 (2015)]. Using a recently introduced perturbative analysis [Phys. Rev. E 98, 020604(R) (2018)], we show that this continuum model reduces to the classic Cahn-Hilliard (CH) model near the onset of MIPS. This makes MIPS another example of the so-called active phase separation. We further introduce a generalization of the perturbative analysis to the next higher order. This results in a generic higher order extension of the CH model for active phase separation. Our analysis establishes the mathematical link between the basic mean-field MIPS model on the one hand, and the leading order and extended CH models on the other hand. Comparing numerical simulations of the three models, we find that the leading order CH model agrees nearly perfectly with the full continuum model near the onset of MIPS. We also give estimates of the control parameter beyond which the higher order corrections become relevant and compare the extended CH model to recent phenomenological models.

Equations

$${\partial }_t\tilde{\rho }=-\nabla \cdot \left[v(\tilde{\rho })\mathbf{p}-D_e\nabla \tilde{\rho }\right],$$

$$v(\tilde{\rho })=v_0-\tilde{\rho }\zeta +{\lambda }^2{\nabla }^2\tilde{\rho }.$$

$${\partial }_t\mathbf{p}=-\nabla P(\tilde{\rho })+D_e{\nabla }^2\mathbf{p}-\mathbf{p},$$

$$P(\tilde{\rho })=\frac{1}{2}v(\tilde{\rho })\tilde{\rho }.$$

$$\tilde{\rho }=\bar{\rho }+\rho .$$

$${\partial }_t\rho$$

$${\partial }_t\mathbf{p}$$

$$+D_e\Delta \mathbf{p}-\mathbf{p},$$

$$\alpha =v_0-R,\beta =\frac{1}{2}(v_0-2R),$$

$$R=\zeta \bar{\rho }.$$

$$\sigma (q)$$

$$=D_2{q}^2-D_4{q}^4+\mathcal{O}(q^6),$$

$$D_2$$

$$D_4$$

$$D_e+\alpha \beta =0$$

$$\frac{1}{2}{v}_0^2-\frac{3}{2}R{v}_0+D_e+R^2=0.$$

$$v_{\ast }$$

$$R_{\pm }=\frac{1}{4}\left[3v_0\pm \sqrt{v_0^2-16D_e}\right].$$

$$R_{\ast }$$

$${\lambda }^2<\zeta v_{\ast }^2/24.$$

$$\rho =\sqrt{\epsilon }{\rho }_1$$

$$v_0=v_{\ast }(1+\epsilon ).$$

$$R=R_{\ast }(1+r_1),\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{mbox}with{r}_1=\sqrt{\epsilon }{\tilde{r}}_1.$$

$${\partial }_x\to \sqrt{\epsilon }{\tilde{\partial }}_x.$$

$${\partial }_t\to {\epsilon }^2{\partial }_{T_1}+{\epsilon }^{5/2}{\partial }_{T_2}.$$

$$\mathbf{p}$$

$${\mathbf{p}}_0$$

$${\mathbf{p}}_1$$

$${\mathbf{p}}_2$$

$${\mathbf{p}}_3$$

$${\mathbf{p}}_4$$

$$-\frac{{\lambda }^2}{2}{\tilde{\partial }}_x\left({\tilde{r}}_1{\rho }_{\ast }+{\rho }_1\right){\tilde{\partial }}_x^2{\rho }_1.$$

$$0=\left({\alpha }_{\ast }{\beta }_{\ast }+D_e\right){\tilde{\partial }}_x^2{\rho }_1.$$

$$0$$

$${\partial }_{T_1}{\rho }_1=-{\tilde{\partial }}_x^2$$

Full text

1 Introduction

Active matter systems are nonequilibrium systems which consume fuel and disspative energy locally. These systems are full of fascinating phenomena and have recently attracted increasing attention in the scientific community Vicsek:2012.1 ; Giardina:2014.1 ; HymanA:2014.1 ; Cates:2015.1 ; Aranson:2015.1 ; Bechinger:2016.1 ; Brangwynne:2017.1 ; Juelicher:2018.1 . Examples range from active molecular processes which are driven by chemical free energy provided by metabolic processes Alberts:2001 up to flocks of birds and schools of fish Vicsek:2012.1 ; Giardina:2014.1 . Various active matter systems also show collective non-equilibrium transitions. On the time scale of these transitions, the involved entities such as proteins, cells or even birds are conserved. Examples include cell polarization EdelsteinKeshet:2007.1 ; Kuroda:2007.1 ; Grill:2013.1 ; Grill:2014.1 ; Alonso:2010.1 ; GoehringNW:2011.1 ; Bergmann:2018.2 , chemotactically communicating cells Tindall:2008.1 ; HillenT:2009.1 ; Romanczuk:2014.1 ; Liebchen:2015.1 , self-propelled colloidal particles Theurkauff:2012.1 ; Palacci:2013.1 ; SpeckT:2013.1 ; Marchetti:2012.1 ; Cates:2013.1 ; Baskaran:2013.2 ; SpeckT:2014.1 , as well as mussels in ecology vandeKoppel:2013.1 .

Self-propelling colloidal particles undergo a non-equilibrium phase transition into two distinct phases - a denser liquid-like phase and a dilute gas-like phase Theurkauff:2012.1 ; Palacci:2013.1 ; SpeckT:2013.1 - if their swimming speed decreases with increasing local density. This is known as motility-induced phase separation (MIPS) Marchetti:2012.1 ; Baskaran:2013.2 ; Cates:2015.1 . It strikingly resembles well-known phase separation processes at thermal equilibrium such as the demixing of a binary fluid. We recently introduced a class of such non-equilibrium demixing phenomena we call active phase separation Bergmann:2018.2 . Among the phenomena identified as members of this class are cell polarization or chemotactically communicating cells. For this class we have shown that the similarities between equilibrium and nonequilibrium demixing phenomena are in fact not coincidental. We have generalized a classical weakly nonlinear analysis near a supercritical bifurcation with unconserved order parameter fields CrossHo to the case of active phase separation with a conserved order parameter field Bergmann:2018.2 . The generic equation describing active phase separation systems turned out to be the classic Cahn-Hilliard (CH) model

• the same generic model that also describes equilibrium phase separation.

In this work, we raise the question whether the recently introduced nonlinear perturbation approach in Ref. Bergmann:2018.2 is also directly applicable to MIPS. We employ this reduction approach to a mean-field description of MIPS provided by Speck et al SpeckT:2014.1 ; SpeckT:2015.1 and show how the MIPS model reduces to the CH model at leading order.

Recently, several phenomenological extensions of the CH model have also been considered as continuum models of MIPS Cates:2014.1 ; Cates:2018.1 . These are extensions of the CH model to the next higher order of nonlinear contributions. In this work, we therefore also introduce an extension of our perturbative scheme that allows us to systematically derive higher order nonlinearities directly from the continuum model for MIPS. Due to our systematic approach, the extended CH model we derive is not a phenomenological model. Instead, we directly map the continuum model for MIPS to the extended CH model. Note that we concentrate on the example of MIPS in this work. However, the extension introduced here can be applied to any system in the class of active phase separation. We thus show in general how both the leading order CH model and its extension describe active phase separation as a non-equilibrium phenomenon.

This work is organized as follows: We first present the mean-field MIPS model and calculate the onset of phase separation in the system. We then introduce the perturbative scheme we use to reduce the MIPS model to the classic CH equation near the onset of phase separation. In the next step, we extend the previous approach to include nonlinearities at the next higher order. Section 5 is an in-depth discussion of the derived leading order and extended CH models including their connection to the mean-field MIPS model and other phenomenological descriptions of MIPS. Finally, in section 6, we present numerical simulations comparing leading order and extended CH to the full mean-field MIPS model to assess validity and accuracy of the reduced models.

2 Model

On a mean-field level, phase separation of active Brownian particles can be described by two coupled equations for the particle density $\tilde{\rho}({\bf r},t)$ and a polarization ${\bf p}({\bf r},t)$ SpeckT:2013.1 ; SpeckT:2015.1 . The evolution of the particle density $\tilde{\rho}$ is determined by

[TABLE]

where $D_{e}$ is the effective diffusion coefficient of the active Brownian particles. $v(\tilde{\rho})$ is the density-dependent particle speed given by

[TABLE]

$v_{0}$ is the speed of a single self-propelled particle. With increasing particle density, the velocity is reduced by $\zeta\tilde{\rho}$ due to interactions with other particles. $\zeta$ is related to the pair distribution function of the individual particles and assumed to be spatially homogeneous SpeckT:2013.1 . The nonlocal contribution in eq. (2) was introduced in Refs. Cates:2013.1 ; Cates:2014.1 . It incorporates the effect that active Brownian particles sample the neighboring particle density on a length scale $\lambda$ larger than the particle spacing. Equation (2) is coupled to a dynamical equation for the polarization SpeckT:2013.1 ; SpeckT:2015.1 ,

[TABLE]

with the “pressure”

[TABLE]

3 Onset of phase separation

The stationary solution of eq. (1) and eq. (3) is any constant density $\bar{\rho}$ and $\mathbf{p}=0$. Therefore, we decompose the particle density into its homogeneous part $\bar{\rho}$ and the inhomogeneous density variation $\rho$:

[TABLE]

Accordingly, we investigate the following dynamical equations for $\rho$ and ${\bf p}$:

[TABLE]

where

[TABLE]

with the density parameter

[TABLE]

We assume $\zeta$ and $D_{e}$ to be constant SpeckT:2015.1 .

The homogeneous basic solution $\rho=0,\leavevmode\nobreak\ {\bf p}=0$ is unstable if the perturbations $\rho,{\bf p}=\hat{\rho},\hat{\bf p}\exp(\sigma t+iqx)$ grow, i.e. if the growth rate $\sigma$ is positive. Solving the linear parts of eqs. (6) with this perturbation ansatz, we find the dispersion relation

[TABLE]

where

[TABLE]

$D_{2}$ changes its sign as a function of $v_{0}$. Assuming $D_{4}>0$, the growth rate $\sigma$ becomes positive in a finite range of $q=[0,q_{max}]$, when $D_{2}>0$. Note that the range of wavenumbers $q$ with positive growth rate extends down to $q=0$. The related instability condition

[TABLE]

provides a quadratic polynomial for the critical mean density $\bar{\rho}$ (represented by the density parameter $R$) and the respective particle speed $v_{0}(R)$:

[TABLE]

For particle speeds $v_{0}>v_{\ast}$, where

[TABLE]

this polynomial has two real solutions

[TABLE]

This corresponds to a critical value $R_{*}$ of the density parameter:

[TABLE]

Note that the assumption $D_{4}>0$ is fulfilled if $\lambda^{2}<2\zeta\alpha\beta^{2}/R$, i.e. for sufficiently small $\lambda$. At the critical point, $v_{0}=v_{*}$ and $R=R_{*}$, this condition simplifies to

[TABLE]

For particle velocities below $v_{\ast}$, the homogeneous solution is stable for any value of the density parameter $R=\zeta\bar{\rho}$. For $v>v_{\ast}$ and $R_{-}<R<R_{+}$ (shaded region in fig. 1) the homogeneous particle density becomes unstable with respect to perturbations.

4 Derivation of Cahn-Hilliard models

In this chapter, we will apply the systematic pertubative scheme introduced recently in Ref. Bergmann:2018.2 to the mean field model, eqs. (6), and reduce them near onset to the well-known Cahn-Hilliard (CH) model. In a second step, we will then expand the pertubative scheme to include higher order contributions.

The transition from the homogenous state of eqs. (1) and (3) to MIPS is either supercritical or slightly subcritical. In both cases, cubic nonlinearities limit the growth of density modulations - as we also confirm in this work a posteriori. Therefore, the amplitudes of the density modulations near MIPS are small and we write

[TABLE]

with a small parameter $\varepsilon$ and $\rho_{1}\sim{\it O}(1)$. Thereby $\varepsilon$ measures the distance from the critical velocity $v_{\ast}$:

[TABLE]

This also allows an expansion of $R_{\pm}(v_{0})$ in eq. (15) near $R_{\ast}$. At leading order, we find $R_{\pm}\simeq R_{*}(1\pm\gamma_{1}\sqrt{\varepsilon})$ with $\gamma_{1}=\sqrt{2}/3$. This suggests the following parameterization of $R$ in the ranges $v_{0}>v_{\ast}$ and $R_{-}<R<R_{+}$ near $R^{\ast}$:

[TABLE]

According to the dispersion relation in eq. (3), the fastest growing mode is given by $q_{e}^{2}=D_{2}/(2D_{4})$. The largest growing wavenumber $q_{max}$ (calculated from $\sigma=0$) is $q_{max}^{2}=D_{2}/D_{4}$. Thus, both $q_{e}^{2}$ and $q_{max}^{2}$ scale with the factor $D_{2}/D_{4}$. Using the previously introduced definitions and expanding for small values of the control parameter $\varepsilon$, we find $D_{2}/D_{4}\propto\varepsilon$ at leading order. Thus, both $q_{e}$ and $q_{max}$ are of the order $\sqrt{\varepsilon}$, i.e. perturbations of the homogeneous basic state vary on a large length scale. Accordingly, we introduce the new scaling $\tilde{x}=\sqrt{\varepsilon}x$, resulting in the following replacement of the differential operator:

[TABLE]

From $q^{2}$ of order $\mathcal{O}(\varepsilon)$ and $D_{2}\propto\varepsilon$ follows that $\sigma\propto\varepsilon^{2}$ according to eq. (3). Thus, the growth of these long wavelength perturbations is very slow. Accordingly, we introduce the slow time scale $T_{1}=\varepsilon^{2}t$. In order to capture the dynamics at the next higher order, we also introduce a second slow time scale $T_{2}=\varepsilon^{5/2}t$. This suggests the following replacement of the time derivatives:

[TABLE]

Since we expressed the density $\rho$ as a multiple of $\sqrt{\varepsilon}$, see eq. (18), we also expand the polarization field ${\bf p}$ in orders of $\sqrt{\varepsilon}$:

[TABLE]

We insert these scalings into the dynamic equations (6) and collect terms of the same order $\sqrt{\varepsilon}^{n}$. The polarization follows the density field adiabatically. Thus, the contributions to the polarization in increasing orders up to $\varepsilon^{5/2}$ are:

[TABLE]

With these solutions, we can systematically solve the equations for the density $\rho_{1}$ in the successive orders of $\sqrt{\varepsilon}$. In the lowest order $\mathcal{O}(\varepsilon^{3/2})$, we find

[TABLE]

This equation, however, is trivially satisfied due to the instability condition ${\alpha_{*}\beta_{*}+D_{e}=0}$.

At order $\mathcal{O}(\varepsilon^{2})$, we get

[TABLE]

With the definition of $R_{*}$ in eq. (14) it follows that $\alpha_{*}+\beta_{*}=0$. Thus, eq. (30) is again trivially fulfilled.

At order $\mathcal{O}(\varepsilon^{5/2})$, we finally get a dynamic equation for $\rho_{1}$:

[TABLE]

Note that we used the expressions in eq. (14) to eliminate $R_{*}$ and $D_{e}$. Equation (31) has the form of the well-known Cahn-Hilliard (CH) equation Bray:1994.1 ; Desai:2009 . This shows that MIPS is a further example of the non-equilibrium demixing phenomenon which shares the universal CH model with classic phase separation. Recently, the notion active phase separation was coined for these types of non-equilibrium phenomena Bergmann:2018.2 . Other recently discussed examples of active phase separation are cell polarization or chemotactically communicating cell colonies Bergmann:2018.2 . All of these very different systems can be reduced to the same universal equation near the onset of phase separation. They thus share generic features as expressed in their common representation via the CH equation.

In the next step, we extend the reduction scheme introduced in Ref. Bergmann:2018.2 to include higher order nonlinearities. Continuing the expansion above to the next order $\mathcal{O}(\varepsilon^{3})$, we obtain:

[TABLE]

We will discuss these new contributions in detail in Chapter 5.2 below.

Equations (31) and (32) can be combined into a single equation by reconstituting the original time scale via $\partial_{t}\rho_{1}=\varepsilon^{2}\partial_{T_{1}}\rho_{1}+\varepsilon^{5/2}\partial_{T_{2}}\rho_{1}$. In addition, we go back to the original spatial scaling by setting $\tilde{\partial}_{x}=\partial_{x}/\sqrt{\varepsilon}$, to the original density $\rho$ via eq. (18), and $r_{1}$ as defined in eq. (20). The complete extended amplitude equation for the density variations $\rho$ then reads:

[TABLE]

In this equation, contributions with the coefficients $\alpha_{i}$ originate from the leading order and are given by

[TABLE]

In other words, eq. (33) with $\beta_{i}=0$ is the rescaled version of eq. (31). The coefficients $\beta_{i}$ signal the new contributions from the next higher order. They are given by

[TABLE]

5 Discussion of the derived Cahn-Hilliard models

In this section, we will discuss the results obtained in the previous section 4. At first we consider the classic CH equation that resulted at leading order of our perturbative analysis. We then take a closer look at the higher order corrections $\propto\beta_{i}$ in eq. (33). We also focus on the relation of the higher order coefficients $\beta_{i}$ to the parameters of recently introduced phenomenological extensions of the CH model for MIPS Cates:2014.1 ; Solon:2018.1 ; Cates:2018.1 .

5.1 Classic CH equation at leading order

For $\beta_{i}=0$, the leading order of eq. (33),

[TABLE]

corresponds to the asymmetric version of the Cahn-Hilliard (CH) equation, see e.g. Refs. Bray:1994.1 ; Desai:2009 , The coefficients $\alpha_{i}$ are given in eqs. (34). Note that the quadratic nonlinearity implies a broken $\pm\rho$-symmetry. This is usually not included in the classic representation of the CH equation since it can be removed by adding a constant to the density: $\rho\rightarrow\rho+\rho_{h}$. In any case, the quadratic nonlinearity vanishes for $\alpha_{3}=0$. For MIPS, this is fulfilled for $r_{1}=0$, or $\bar{\rho}=\rho_{*}$ accordingly. This special case has also been considered in SpeckT:2015.1 where they found a CH equation with coefficients consistent with $\alpha_{i}$ above.

Equation (36) can be derived from the energy functional

[TABLE]

via

[TABLE]

On first glance this is a surprising result since the two initial dynamical equations for the density, eq. (1), and the polarization, eq. (3), do not follow potential dynamics and therefore cannot be derived from a functional. Nevertheless, this specific property has been seen for other non-equilibrium systems: The evolution equation for the envelope of spatially periodic patterns also follows potential dynamics while the dissipative starting equations do not CrossHo ; Cross:2009 .

5.2 Extended CH model

We now take a closer look at the CH model extended to the next higher order, eq. (33) with coefficients $\beta_{i}$ given in eqs. (35). The contributions $\beta_{1}$, $\beta_{2}$ and $\beta_{3}$ are corrections to the coefficients $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ of the leading order CH equation. Note, however, that according to eqs. (35a) and (35c), $\beta_{1}$ and $\beta_{3}$ are functions of $\varepsilon$ and thus both increase with the distance $\varepsilon$ from phase separation onset. Notably, $\beta_{3}$ - the correction to the quadratic nonlinearity - is not a function of the relative deviation $r_{1}$ from the critical density parameter $R_{*}$. Thus, while for $r_{1}=0$ the CH model at leading order is $\pm\rho$-symmetric, the symmetry is always broken at higher order.

The coefficients $\beta_{5}$ and $\beta_{6}$ are the prefactors of higher order nonlinearities. These new contributions $\propto\partial_{x}^{2}(\partial_{x}\rho)^{2}$ and $\propto\partial_{x}^{4}\rho^{2}$ are structurally different compared to the terms in the leading order CH model. In general, an additional nonlinearity $\propto\partial_{x}^{2}\rho^{4}$ is of the same order as these two contributions. However, in the exemplary case of MIPS we analyze here this term does not appear. Note, however, that the higher order extension of the CH model presented here can also be applied to other active phase separation systems. We expect the additional nonlinearity of the form $\propto\partial_{x}^{2}\rho^{4}$ to be relevant in other examples such as cell polarization or chemotaxis.

In the context of MIPS, a contribution $\propto\partial_{x}^{2}(\partial_{x}\rho)^{2}$ has been introduced via a phenomenological approach in Ref. Cates:2014.1 . The CH model extended by this term has been called Active Model B. It was considered as a non-equilibrium extension of the CH model and minimal model for MIPS. We would like to reiterate that the CH model as given by eq. (36) (without any additional nonlinear terms) is the leading order description of the non-equilibrium phenomenon of active phase seperation Bergmann:2018.2 . As we have shown here, this also includes MIPS. All higher order nonlinearities vanish for $\varepsilon\rightarrow 0$ (see also the discussion in sec. 5.4). In that respect Active Model B is a nonlinear extension of the CH model - not an extension of the CH model to non-equilibrium systems. Our systematic approach reveals the existence of the additional higher nonlinearity $\propto\partial_{x}^{4}\rho^{2}=2\partial_{x}^{2}[(\partial_{x}\rho)^{2}+\rho\partial_{x}^{2}\rho]$. It includes the nonlinear correction to the CH model, $\propto\partial_{x}^{2}(\partial_{x}\rho)^{2}$, that leads to the Active Model B Cates:2014.1 ; Cates:2015.1 . The second part of the new nonlinear correction term, $\propto\partial_{x}^{2}(\rho\partial_{x}^{2}\rho)$, has recently been included in a further CH extension for MIPS called Active Model B+ Solon:2018.1 ; Cates:2018.1 . Note that the contribution $\propto\beta_{6}$ in eq. (33) vanishes for $\lambda=0$. Active Model B and Active Model B+ also do not include the quadratic nonlinearity $\propto\beta_{3}\rho^{2}$. Our analysis shows, however, that the coefficients $\beta_{i}$ in general are not independent from each other and $\beta_{2}$ in fact always appears simultaneously with the nonlinearity $\propto\beta_{5}$. The broken $\pm$-symmetry and the resulting asymmetric phase separation profiles depend on the distance $\varepsilon$ from threshold (see $\beta_{3}$ in eq. (35c). It is an important qualitative feature of the system behavior above threshold.

As discussed in sec. 5.1, the leading order CH model can be derived from an energy potential. For the extended Ch model, eq. (33), integrability depends on the coefficients of the additional higher order contributions: For arbitrary values of $\beta_{5}$ and $\beta_{6}$, the extended CH model is non-integrable. In the special case $\beta_{6}=-\beta_{5}$, however, eq. (32) can be derived from the energy functional

[TABLE]

For the MIPS model, eqs. (6), this condition is fulfilled for

[TABLE]

Note, however, that the linear stability analysis in sec. 3 introduced a condition for $\lambda$: $\lambda^{2}<\zeta v_{*}^{2}/24$ in eq. (17). This condition and eq. (40) cannot be fulfilled simultaneously. Thus, the integrability of the extended CH model depends on the exact parameter choices. For the MIPS continuum model we investigate here, there do not seem to be suitable parameter choices that enable integrability. But note again that our approach can be applied to other systems showing active phase separation. For these other models, the coefficients of the extended CH model could allow for integrability.

5.3 Comparison of linear stability

As a first step to assess the quality of our derived reduced equations, we analyze the linear stability of the homogeneous basic state $\rho=0$, and compare to the stability of the full MIPS model. As discussed in sec. 3, the instability condition for the full MIPS system is given by eq. (12). Using $v_{0}=v_{*}(1+\varepsilon)$, $R=R*(1+r_{1})$ and the definitions of $D_{e}$ and $R_{*}$ as given by eqs. (14), we find

[TABLE]

for the onset of phase separation. Thus, in the symmetric case $r_{1}=0$ the threshold is $\varepsilon_{c}=0$. For $r_{1}\neq 0$ the onset of phase separation is shifted to larger values of $\varepsilon$. Larger particle velocities $v_{0}$ are thus required to trigger the demixing process.

Similarly, we can analyze the linear stability of both the leading order CH equation, eq. (36), and its higher order extension, eq. (33). The threshold calculated from the linear parts of eq. (36) is given by

[TABLE]

Comparing this to $\varepsilon_{c}$ in eq. (41), we find that the shifting of the threshold due to finite $r_{1}$ is represented up to leading order of $r_{1}$. Assuming $r_{1}>0$, $\varepsilon_{c,\text{lead}}$ significantly overestimates the real threshold $\varepsilon_{c}$. For the extended CH equation, eq. (33), we find the threshold

[TABLE]

This is in agreement with the threshold for the full model, eq. (41), up to the order $\mathcal{O}(r_{1}^{3})$. The threshold is therefore only slightly underestimated compared to the full model. Keeping these different threshold values in mind is particularly important for the numerical comparison of the MIPS model, eqs. (1) and (3), to its two reductions, eqs. (36) and (33) in sec. 6. All three equations only provide the exact same threshold, namely $\varepsilon_{c}=0$, in the special case $r_{1}=0$.

The linear stability analysis also provides the dispersion relation for the perturbation growth rate $\sigma$. For the full model, it is given by eq. (3). Expanding for small perturbation wavenumbers $q$, the general form of the growth rate is

[TABLE]

The coefficients $D_{2}$ and $D_{4}$ are given in eqs. (10) and (11), respectively. Using the definitions introduced in the course of the perturbative expansion, $D_{2}$ can be rewritten to

[TABLE]

Good agreement between the full MIPS model and its reduction to eq. (33) can only be expected if the reduced equations are able to reproduce the basic form of this growth rate. The linear part of eq. (33) leads to a growth rate of the form

[TABLE]

where

[TABLE]

$G_{2}$ is in agreement with $D_{2}$ of the full model equations up to linear order in $\varepsilon$. $D_{2}$ only includes an additional term of order $\mathcal{O}(\varepsilon^{2})$: $D_{2}=G_{2}-v_{*}^{2}\varepsilon^{2}/2$. $G_{4}$ exactly reduces to $D_{4}$ in the case $\varepsilon=r_{1}=0$. In the limit $\varepsilon\rightarrow 0$ but $r_{1}\neq 0$, the two terms agree up to linear order in $r_{1}$. As discussed in sec. 3, the coefficient $D_{4}$ has to be positive for the instability condition to hold and to ensure damping of short wavelength perturbations. The same applies to the coefficient $G_{4}$. The condition $G_{4}>0$ is fulfilled if

[TABLE]

Note the similarity to the previously derived condition in eq. (17).

5.4 Significance of nonlinear corrections

In this section, we discuss the importance of the higher order nonlinearities compared to the leading order terms of the classic Cahn-Hilliard model in eq. (36). For this comparison we focus on the case with $\pm$-symmetry at leading order, i.e. $r_{1}=0$. We rescale time, space and amplitude in eq. (33) via $t^{\prime}=\tau_{0}\varepsilon^{2}t$, $x^{\prime}=\xi_{0}\sqrt{\varepsilon}x$ and $\rho^{\prime}=\rho_{0}\rho/\sqrt{\varepsilon}$, respectively, where

[TABLE]

This allows us to rewrite eq. (33) in the following form:

[TABLE]

where

[TABLE]

The first line in eq. (51) is the parameter-free, $\pm\rho$-symmetric version of the Cahn-Hilliard model as described, e.g., in Refs. Bray:1994.1 ; Desai:2009 . The additional three contributions are the first higher order corrections as gained above via a systematic reduction of the continuum model for MIPS. These three corrections are proportional to $\sqrt{\varepsilon}$ and thus vanish when approaching the onset of active phase separation ($\varepsilon\to 0$). In the limit $\varepsilon\to 0$ the classic CH model thus fully describes the non-equilibrium mean-field dynamics of MIPS. With increasing $\varepsilon$, the higher order contributions become more and more important.

Note that eq. (51) was derived under the assumption $r_{1}=0$. As discussed in sec. 5.1, the CH model at leading order is $\pm\rho$-symmetric in this case. The three higher order contributions in eq. (51), however, break the $\pm\rho$-symmetry with increasing $\varepsilon$. Moreover, in the case of the MIPS model we analyze here, the coefficient $\gamma_{1}$ does not depend on any of the system parameters at all. Thus, there is in fact no special case in which this contribution can be neglected.

The coefficients of the other two higher order nonlinearities, $\gamma_{2}$ and $\gamma_{3}$, are functions of the system parameters, especially of $\lambda$. Typical parameter choices for the continuum model in eq. (6) are such that $v_{*}$ and $\zeta$ are of order $\mathcal{O}(1)$. Accordingly, $\lambda$ has to be small to fulfill the condition in eq. (17). Therefore, an expansion of $\gamma_{2}$ and $\gamma_{3}$ in terms of small $\lambda$ is appropriate:

[TABLE]

In the limit $\lambda=0$ the coefficient $\gamma_{2}$ vanishes, i.e. $\gamma_{2}=0$, and $\gamma_{3}$ simplifies to $\gamma_{3}=\gamma_{1}$. For finite $\lambda$, $\gamma_{2}$ also becomes finite. But since according to eq. (53) $\gamma_{2}$ is proportional to $\lambda^{2}$, it will be much smaller than $\gamma_{3}$ for small $\lambda$. For MIPS as described by the mean-field model in eqs. (6), the impact of the nonlinearity $\propto\partial_{x}^{2}(\partial_{x}\rho)^{2}$ thus seems to overshadow the term $\propto\partial_{x}^{4}\rho^{2}$. This predominance of $\gamma_{3}$, however, is specific to MIPS. For other examples of active phase separation such as cell polarization or chemotactically communicating cells, we expect that the nonlinearities described by $\gamma_{1}$ or $\gamma_{2}$ can be of similar order as $\gamma_{3}$. As mentioned earlier, for both examples of active phase separation we also expect an additional higher order correction $\propto\partial_{x}^{2}\rho^{4}$ which is completely absent in the MIPS model.

6 Numerical comparison

In this section, we compare numerical simulations of the full MIPS model, eqs. (6), to both the leading order CH equation, eq. (36), as well as the extended version including higher nonlinearities, eq. (33). On the one hand, this allows us to assess the quality and validity range of our reduction scheme in general. On the other hand, comparing the leader order and the extended CH model also gives us information about the importance of higher order nonlinearities in MIPS.

All simulations were performed using a spectral method with a semi-implicit Euler time step. The system size was $L=100$ with periodic boundary conditions and $N=256$ Fourier modes were used.

We first analyze the special case $r_{1}=0$, i.e. $\bar{\rho}=\rho_{*}$. This is the case in which the $\pm$-symmetry-breaking quadratic nonlinearity vanishes at leading order. We choose $v_{*}=1$ and $\zeta=1$ throughout all of the following simulation results. As discussed in sec. 5.4, $\lambda$ has to be small and is thus not expected to significantly influence the results. We thus set $\lambda=0$.

Figure 2 shows the steady state profiles for the three models (full MIPS model, leading order CH and extended CH) at $\varepsilon=0.01$. The profiles are typical for phase separation solutions: We find two distinct regions where the mean density is either increased ($\rho>0$) or decreased ($\rho<0$). In each of the regions $\rho$ is essentially spatially constant, creating two distinct density plateaus $\rho_{\text{min}}$ and $\rho_{\text{max}}$. The two plateaus are smoothly connected at their boundary, resembling a hyperbolic tangent function. Note that the mean density in the system is conserved. Thus, the areas under the positive and negative parts of $\rho(x)$ are equal.

The solution for the full system is represented as the outline of the grey shaded area. We first compare this to the leading order CH equation (blue dashed line). As predicted, the leading order CH equation results in a symmetric phase separation profile, i.e. the two plateaus have the same absolute value: $\rho_{\text{max}}=\lvert\rho_{\text{min}}\rvert$. This does not accurately represent the solution for the full system, which is already slightly asymmetric. However, the leading order CH equation gives a good approximation of the plateau values with a deviation of less than $7\%$ from the real value. Extending the CH equation to the next higher order (red dotted line in fig. 2), we can almost perfectly reproduce the profile for the full MIPS model. It accurately represents the asymmetry of the phase separation profile. The deviation in the plateau values shrinks to less than $2\%$.

Figure 3 shows the absolute plateau values $\lvert\rho_{\text{min}}\rvert$ and $\rho_{\text{max}}$ as a function of $\varepsilon$ - the distance from phase separation onset. The bifurcation to active phase separation is supercritical in this case: starting at $\varepsilon_{c}=0$, the plateau values increase monotonically. Considering only the leading order approximation (blue dashed line), we again find the system to be symmetric for all values of $\varepsilon$. In reality, the full system (black solid lines) becomes more and more asymmetric for increasing $\varepsilon$. This is very accurately represented by the higher order approximation (red dotted lines). It only starts to deviate from the full model further from threshold. Importantly though, close to the onset of mobility induced phase separation, as $\varepsilon$ becomes smaller, the full model becomes more and more symmetric. All three models then are in increasingly good agreement. This again underlines the fact that the classic CH model is the simplest generic model for active phase separation. All active phase separation phenomena of this type can be reduced to the CH model close to onset. Higher order nonlinearities only come into play further from threshold.

If we allow $r_{1}\neq 0$, phase separation is asymmetric even at leading order. This can be seen in fig. 4 which shows the steady state profiles for the full MIPS model, leading order CH and extended CH at $\varepsilon=0.02$. Here, the leading order CH equation (dashed blue line) results in an asymmetric solution. However, the predicted plateau values deviate about $20\%$ from the full system (outlines of shaded grey region). The extended CH model, meanwhile, is still able to accurately predict the full system solution with a deviation of less than $6\%$.

Looking at the plateau values as a function of $\varepsilon$ (see fig. 5) solidifies this impression: The leading order CH model gives a good qualitative representation of the full system. Going to the extended CH model provides very good quantitative agreement with the full model even for larger values of $\varepsilon$. As discussed earlier in sec. 5.3, the onset of phase separation (i.e. the $\varepsilon$-value at which the homogeneous solution $\lvert\rho_{\text{min}}\rvert=\rho_{\text{max}}=0$ becomes unstable) is shifted to finite values of $\varepsilon$ in the case $r_{1}\neq 0$. For the given system parameters, the threshold for the full system is shifted to $\varepsilon_{c}\approx 0.013$. The leading order CH model significantly overestimates this threshold, shifting to $\varepsilon_{c}\approx 0.02$. The extended CH model only very slightly underestimates the real threshold. Note that above this threshold, the plateau values immediately jump to finite values. Thus, the transition from the homogeneous to the phase-separated state is no longer smooth. On the other hand, fig. 5) also shows that the branches of finite density plateau values extend below the thresholds noted above. This creates a range of bistability - a range of control parameter values in which both the homogeneous and the phase-separated state are stable simultaneously. All of these characteristics indicate that bifurcation from the homogeneous state to active phase separation is now subcritical.

7 Conclusion

Starting from the mean-field model for active Brownian particles in Refs. SpeckT:2013.1 ; SpeckT:2015.1 , we applied a perturbative approach introduced in Ref. Bergmann:2018.2 . We showed that the non-equilibrium phenomenon motility-induced phase separation (MIPS) is described near its onset at leading order by the Cahn-Hilliard (CH) model Cahn:58.1 ; Cahn:1961.1 ; Bray:1994.1 ; Desai:2009 . This is in agreement with a recent observation that the CH model describes the system-spanning behavior of a number of very different demixing phenomena in active and living systems far from thermal equilibrium Bergmann:2018.2 . The results in this work show that MIPS also belongs to this class of active phase separation. Thus, even though the CH model was originally introduced to describe phase separation of binary mixtures in thermal equilibrium, our analysis shows that it is also the generic leading order description of active phase separation

• a non-equilibrium phenomenon.

We also extended the perturbative scheme introduced in Ref. Bergmann:2018.2 beyond the CH model to next higher order nonlinearities. In this work, we used the continuum model for MIPS as a framework to establish this concept. The extension of our nonlinear expansion, however, can also be applied to other systems showing active phase separation (with a conserved order parameter field) such as cell polarization and clustering of chemotactically communicating cells. Having a $\pm$-symmetric CH model at onset of active phase separation, we find that in general four nonlinear terms come into play at the next higher order. Two of them have the same form as contributions suggested in previous phenomenological extensions of the CH model for MIPS Cates:2014.1 ; Cates:2015.1 ; Solon:2018.1 ; Cates:2018.1 . These phenomenological models are thus related to the extended CH model that our perturbative scheme provides. Our approach, however, is non-phenomenological: It establishes a direct mathematical link between the coefficents of the extended CH model and the full mean-field description of MIPS (or any other basic model of active phase separation in general). It shows in addition, that the coefficients of the additional contributions in the extended CH model are in general not independent from each other, as often assumed in phenomenological approaches. Furthermore, these coefficents are system-specific and cannot be removed by rescaling as in the case of the leading order CH model. It is also important to reiterate that these nonlinear extensions become negligible when approaching the onset of MIPS or other examples of active phase separation. Therefore, the leading order CH model already covers the universal behavior of MIPS (as a non-equilibrium phenomenon) near its onset. Higher order nonlinearities mainly improve accuracy and become relevant further from threshold. They should thus not be seen as the key to expand the CH model to non-equilibrium systems.

Within the systematics of pattern formation theory, the work we introduced in Ref. Bergmann:2018.2 and extended here is a weakly nonlinear analysis and reduction method for active phase separation described by conserved order parameter fields. It can be seen as a yet unexplored counterpart to the weakly nonlinear analysis of (non-oscillatory) spatially periodic patterns with unconserved order parameter fields and its numerous applications CrossHo ; Pismen:2006 ; Cross:2009 ; Meron:2015 ; Misbah:2016 .

Our generic approach for active phase separation opens up several pathways for further system-spanning investigations. Coarsening dynamics in large systems, and especially the role of higher nonlinearities in this context, have already been of particular interest to the scientific community (see, e.g., Ref. Cates:2018.1 for MIPS). Other active phase separation phenomena such as cell polarization, on the other hand, take place in very small systems where coarsening plays a less important role EdelsteinKeshet:2013.1 . For these systems, spatial constraints may significantly influence the behavior instead. Studies on spatially periodic patterns have already shown that confinement may trigger various interesting generic effects (see e.g. Rapp:2016.1 ) and even induce patterns in small systems which are unstable in larger systems (Bergmann:2018.1 and references therein). On the basis of our results, it will be interesting to investigate finite size effects on non-equilibrium phase transitions with conservation constraints.

Acknowledgements.

Support by the Elite Study Program Biological Physics is gratefully acknowledged.

Author Contribution Statement

All authors contributed to the design of the research, to calculations, the interpretation of results and the writing of the manuscript. LR performed numerical simulations.

Data vailability statement

The simulation datasets used in this article are available from the corresponding author on request.