Periodic patterns displace active phase separation

Patterns are ubiquitous in nature. They emerge spontaneously in a plethora of living or inanimate driven systems [1–10], and their esthetic appeal is immediately apparent to all observers [1]. Prototypical patterns, such as stripes, hexagons and traveling waves can be found in various fluid-dynamical, chemical and biological systems and they often play a key role in various processes in nature. To name just a few examples: convection patterns in fluid dynamical systems enhance the thermal transport [2, 6], protein dynamics in E. coli plays a key role in finding the cell center during the cell division process [11–14], and the formation of patterns are the basis of successful survival strategies for vegetation in water-limited regions [8, 9].

The mechanisms that drive nonequilibrium patterns are as diverse as the systems in which they occur. Instabilities leading to patterns are usually divided into types I, II and III, which differ in their preferred wavenumber q₀ and/or in their frequency ω₀ and the wavenumber dependence of the dispersion relation [2]. Despite their very different origins, nonequilibrium patterns of the type I and type III share common universal properties, which are covered by universal Ginzburg–Landau equations [2–5, 15–22]. They describe the dynamics of the envelope (resp. the order parameter field) of stationary and traveling stripe or oscillating patterns.

A number of nonequilibrium demixing phenomena with conserved particle dynamics attract increasing attention recently [23–32]. Ultimately, these are type II instabilities, for which the theory is currently undergoing a lively development. Such nonequilibrium demixing phenomena include cell polarization [23, 24], chemotactically communicating cells [25, 26], motility induced phase separation (MIPS) in self-propelled Brownian particle systems [27–31] and also a very early model for the ion-channel density in a membrane [32, 33].

Demixing phenomena in thermal equilibrium systems are long known and they are described by the Cahn–Hilliard (CH) mean-field model [34–37]. It is very surprising that the CH equation describes also the generic order-parameter field for type II instabilities with particle conservation, i.e. also for non-equilibrium demixing phenomena [38–41]. This was shown by generalizing the perturbation theory used to derive equations for unconserved order parameter fields [2, 3, 17–20, 22] to cases with conserved order parameter fields [38]. Its application to the dynamical equations of chemotactic systems [38, 41], to models of cell polarization [38, 40] and to MIPS [39] results always directly in the generic CH equation. Accordingly, the coefficients of the CH equation now depend on the kinetic parameters of the respective system and it describes nonequilibrium rather than equilibrium transitions.

In the case of MIPS the CH model was extended to the active model B and the active model B+ by higher order nonlinear terms to describe interesting coarse-grained phenomena further away from the onset [42–44] with particular emphasis on the effects of noise. The perturbation theory in reference [38] extended to the next higher nonlinear order and applied to a mean-field theory for MIPS [29, 45] gives the generic structure of AMB+ [39] as well. Moreover, this perturbation-theoretic approach unambiguously links the parameters of AMB+ to the kinetic parameters of the mean field model in references [29, 45]. In terms of this perturbation theory the deterministic part of AMB+ is a generic nonlinear extension of the CH model rather than an extension into nonequilibrium.

The perturbation theory from reference [38] systematically provides, in addition to the generic next higher nonlinear contributions, also a sixth-order spatial derivative of the order parameter field. This higher order derivative is not part of active model B+, but becomes indispensable beyond a finite distance from the threshold of MIPS and at larger amplitudes of the order parameter field. In this range, solutions without the sixth-order derivative then show unbounded growth. The conserved Swift–Hohenberg model+ (CoSH+) introduced in section 2 for a class of systems showing nonequilibrium phase separation includes this higher derivative and shows a novel secondary bifurcation in the nonlinear regime of active phase separation, which is characterized in section 3. There are two complementary scenarios of the secondary bifurcation, and they are described in sections 4 and 5. In the section 6 we summarize our results, classify them and highlight further perspectives of this work.

We introduce and investigate a nonlinear model for a real-valued conserved order parameter field ψ(x, t), which is simultaneously an extension of the CH model [35–37] by the generic next higher order contributions and an extension of the conserved Swift–Hohenberg (CSH) model [46] by the leading order nonlinear gradient terms,

$$\begin{equation}{\partial }_{t}\psi =-{\partial }_{x}^{2}\left[\varepsilon \psi +\left({D}_{4}+{\beta }_{2}\psi \right){\partial }_{x}^{2}\psi -{D}_{6}{\partial }_{x}^{4}\psi +{\beta }_{1}{\psi }^{2}-g{\psi }^{3}+{\beta }_{3}{\left({\partial }_{x}\psi \right)}^{2}\right],\end{equation} \tag{ 1 }$$

with D₄, D₆, g > 0. Equation (1) is symmetric with respect to the transformation (ψ, β₁, β₂, β₃) → −(ψ, β₁, β₂, β₃). It can be written in the form of a conservation law ∂_t ψ = −∂_x j with j = j(ψ, ∂_x ψ) and contains a number of models as special cases. For D₆ = β₁ = β₂ = β₃ = 0 equation (1) reduces to the classic CH model in one spatial dimension [34–37]. Using an ɛ-scaling the contributions including D₆, β₂, β₃ are shown to vanish in the limit of small control parameter values ɛ in appendix A.

For D₆ = 0 the model corresponds to the nonlinearly extended CH model for active phase separation [39], which was recently derived by a systematic perturbation calculation [38] from a dissipative mean-field model for MIPS in active colloids [45]. Equation (1) reduces for D₆ = β₁ = 0 to the one-dimensional version of the so-called active model B+ with K₁ ≠ 0 in reference [44]. In the case of vanishing coefficients D₆ = g = β₁ = β₂ = 0 the equation (1) reduces to the conserved Kuramoto–Sivashinsky model [47]. For the parameter set β_i = 0 (i = 1, 2, 3), D₆ = 1, D₄ = −2 and ɛ = r − 1 equation (1) describes the CSH model for spatially periodic patterns [46], also known as the phase-field crystal model [48].

The 6th order derivative ${D}_{6}{\partial }_{x}^{6}\psi$ in the linear operator in equation (1) has been added in contrast to the models in references [39, 44]. Without this higher order contribution model (1) will become structurally unstable for larger values of the control parameter. Above the onset of phase separation one has $\psi \propto \sqrt{\varepsilon }$ and ψ takes positive and negative values equally likely. Therefore, the prefactor D₄ + β₂ ψ may become small ${D}_{4}+{\beta }_{2}\psi \sim \sqrt{\varepsilon }$ for negative values of β₂ ψ < 0 or even negative for D₄ > 0. Moreover, in this parameter range close the threshold one has slow variations of the fields, i.e. ${\partial }_{x}^{2}\psi \sim O\left(\varepsilon \right)$, ∂_x ∼ O(ɛ^1/4) and all contributions to equation (1) are of the same order ∝ɛ², including ${\partial }_{x}^{6}\psi$. This illustrates why equation (1) would be incomplete without the contribution ${\partial }_{x}^{6}\psi$. In this sense equation (1) is a consistent generic model being structurally stable also for small or negative values ${D}_{4}+{\beta }_{2}\psi \sim \sqrt{\varepsilon }$. Equation (1) is a systematic extension of the CSH model by leading order nonlinear gradient terms and therefore, we call it CoSH+ model.

Equation (1) has the additional interesting property that its solutions obey a gradient dynamics

$$\begin{equation}{\partial }_{t}\psi ={\partial }_{x}^{2}\left(\frac{\delta \mathcal{F}}{\delta \psi }\right),\end{equation} \tag{ 2 }$$

for the parameter combination β₂ = 2β₃ with the functional

$$\begin{equation}\mathcal{F}\left(\psi \right)=\int \mathrm{d}x\left[\frac{{\psi }^{4}}{4}-\frac{{\beta }_{1}{\psi }^{3}}{3}-\frac{\varepsilon }{2}{\psi }^{2}+\left(\frac{{D}_{4}}{2}+{\beta }_{3}\psi \right){\left({\partial }_{x}\psi \right)}^{2}+\frac{{D}_{6}}{2}{\left({\partial }_{x}^{2}\psi \right)}^{2}\right].\end{equation} \tag{ 3 }$$

In this work we focus either on the case without a gradient dynamics, especially on the parameter range β₂ ∝ O(1) with β₃ = 0 or small values of β₃. In a second case we focus on the range with a gradient dynamics and its neighborhood β₂ ∼ 2β₃. In both ranges we find two different scenarios of a secondary transition taking place above the onset of phase separation.

Numerically we solve equation (1) with periodic boundary conditions on a spatial domain (0, L) (see appendix B for specifics). Unless stated otherwise, the following set of parameters is used throughout: ${\beta }_{2}=\sqrt{2}$, D₄ = 0.15, D₆ = 1/6 and g = 1.

In the range of small ɛ the CoSH+ model (1) reduces to the classic CH-model as demonstrated in appendix A. This is also true for its solutions. It can be seen in figure 1(a) for a small ɛ = 0.005, where a stationary solution is shown that has the same form as typical monotonic domain-wall solutions of the CH-model connecting the plateau values ${\psi }_{0}={\pm}\sqrt{\varepsilon }$ [35, 36].

Zoom In Zoom Out Reset image size

By increasing ɛ as in figure 2(a) up to ɛ = 0.06, a domain wall solution as in figure 1 is deformed for finite β₂ > 0 (β₂ < 0) as damped oscillations develop at the lower (upper) plateau near the domain wall. In addition, the lower (upper) plateau widens and the amplitudes of both plateaus are shifted upwards (downwards).

To investigate whether these changes will eventually lead to an instability of domain wall solutions, as depicted in figure 2, we use the following ansatz

$$\begin{equation}\psi \left(x,t\right)={\psi }_{0}\left(x\right)+{\psi }_{1}\left(x,t\right),\end{equation} \tag{ 4 }$$

and separate the primary state, ψ₀(x), from contributions ψ₁(x) of a possible secondary state to ψ(x). This gives the nonlinear dynamical equation for ψ₁:

$$\begin{align}\hfill {\partial }_{t}{\psi }_{1}& =-{\partial }_{x}^{2}\left[\left(\varepsilon +{D}_{4}{\partial }_{x}^{2}-{D}_{6}{\partial }_{x}^{4}-3{\psi }_{0}^{2}+2{\beta }_{1}{\psi }_{0}+{\beta }_{2}\left({\partial }_{x}^{2}{\psi }_{0}+{\psi }_{0}{\partial }_{x}^{2}\right)+2{\beta }_{3}\left({\partial }_{x}{\psi }_{0}\right){\partial }_{x}\right){\psi }_{1}\right.\hfill \\ \hfill & \quad \left.-3{\psi }_{0}{\psi }_{1}^{2}-{\psi }_{1}^{3}+{\beta }_{1}{\psi }_{1}^{2}+{\beta }_{2}{\psi }_{1}{\partial }_{x}^{2}{\psi }_{1}\enspace +{\beta }_{3}{\left({\partial }_{x}{\psi }_{1}\right)}^{2}\right].\hfill \end{align} \tag{ 5 }$$

Further away from the domain walls and at small values of ɛ the constant amplitude at the plateaus is approximately ψ_* $={\pm}\sqrt{\varepsilon }$. In a first step we neglect the spatial variation of ψ₀(x) near the domain walls, assume at the end of the plateaus periodic boundary conditions and the constant amplitude ψ = ψ_* at the whole plateau. In this case equation (5) simplifies with β₁ = β₃ = 0 to an equation with constant coefficients and its linear part can be solved by the ansatz ψ₁ = F exp(σt + iqx). The resulting growth rate has the following q-dependence

$$\begin{equation}\sigma \left(q\right)=\left(\varepsilon -3{\psi }_{{\ast}}^{2}\right){q}^{2}-\left({D}_{4}+{\beta }_{2}{\psi }_{{\ast}}\right){q}^{4}-{D}_{6}{q}^{6}.\end{equation} \tag{ 6 }$$

The prefactor in front of ∝q⁴ may change its sign in the range ψ_* β₂ < −D₄ even for D₄ > 0. In the range of a vanishing or even negative coefficient of q⁴ the stabilizing higher-order derivative ${D}_{6}{\partial }_{x}^{6}\psi$ in equation (1) becomes crucial and ensures with D₆ > 0 the damping of short-scale perturbations in equation (6).

Zoom In Zoom Out Reset image size

Beyond the onset of phase separation one has ɛ > 0 and therefore a negative damping coefficient of q² in equation (6). In the case of a sign change of the prefactor of q⁴ the growth rate σ(q) develops besides q = 0 a second extremum as indicated in figure 1(b). The wavenumber at the finite-q extremum is

$$\begin{equation}{q}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{2}=\frac{1}{3{D}_{6}}\left[-{\beta }_{2}{\psi }_{{\ast}}-{D}_{4}+\sqrt{{\left({\beta }_{2}{\psi }_{{\ast}}+{D}_{4}\right)}^{2}+3{D}_{6}\left(\varepsilon -3{\psi }_{{\ast}}^{2}\right)}\right],\end{equation} \tag{ 7 }$$

if ψ_* β₂ > (3D₆ ɛ − D₄). For a sufficiently negative D₄ + β₂ ψ_* the growth rate becomes positive in a finite q-range around its maximum σ(q_max), as indicated by the solid line in figure 1(b). In the CSH model the q-dependence of the perturbation growth rate with respect to the homogeneous basic state shows a similar pattern forming behavior [46].

For the approximation ${\psi }_{{\ast}}=-\sqrt{\varepsilon }$ and β₂ > 0 the critical ɛ-value, above which σ(q_max) becomes positive, is given by

$$\begin{equation}{\varepsilon }_{qc}=\frac{{D}_{4}^{2}}{{\left(2\sqrt{2{D}_{6}}-{\beta }_{2}\right)}^{2}},\end{equation} \tag{ 8 }$$

with ɛ = ɛ_qc at q_c = q_max. The formula shows that the threshold for an instability of the lower plateau value increases with ${D}_{4}^{2}$. If we account for the numerically observed upward shift in the lower plateau value, as in figure 2(a), by an increase of approximately 25% compared to ${\psi }_{{\ast}}\simeq -\sqrt{\varepsilon }$, then the threshold of the instability is reduced by about 50% according to equation (6) compared to ɛ_qc .

The spatial variation of ψ₀(x) near the domain wall causes spatially varying coefficients in the nonlinear equation (5). Such spatially varying coefficients can be represented by a Fourier series with locally varying Fourier amplitudes as demonstrated in reference [49]. It is known, that spatial modulations of parameters with a wavenumber twice that of the wavenumber of the pattern at onset lead to a reduction of the pattern threshold proportional to the amplitude of the modulation due to a so-called 1:2 resonance [50, 51]. Here the spatially varying coefficients in equation (5) include Fourier modes with a wavenumber twice that of the wavenumber q of ψ₁, but restricted to the range near the domain wall. These modulations are locally in 1:2 resonance with ψ₁ in conjunction with the mentioned shift of ψ_*.

In figure 2(a) the two-domain solution at the control-parameter value ɛ = 0.06 is stable. This ɛ-value is far below the threshold ɛ_qc ≃ 0.334 of an instability of a constant plateau solution ψ_* $=-\sqrt{\varepsilon }$ (for parameters as in figure 2). However, the shift of the plateau values ψ_* and the local resonance lead to a shift of the onset of a secondary pattern such that a small ɛ-jump to ɛ = 0.065 already destabilizes the two-domain solution in figure 2(a) leading to the nonlinear evolution as shown in figure 2(c).

According to the arguments above, the domain wall solution becomes unstable far below the predicted linear stability threshold ɛ_qc ≃ 0.334 for a constant plateau in the absence of spatial variations near the domain wall. As can be seen in figure 2(c) the periodic pattern indeed nucleates near the domain walls and then evolves to a system filling extended periodic pattern as shown in figure 2(b). It is quite surprising that the periodic solution spreads throughout the entire system instead of only over the lower plateaus. Note that for β₂ < 0 the upper plateau becomes unstable. The scenario of a periodic instability restricted to the lower plateau for β₂ > 0 (to the upper plateau for β₂ < 0) is described in section 5.

The extended spatially periodic patterns which displace active phase separation beyond the secondary instability as in figure 2(c), are multistable in a finite wavenumber band, as we show in section 4.1. In addition, we demonstrate in section 4.2 that the transitions between stable periodic patterns of different wavenumbers exhibit the same universal behavior as those in classical pattern forming systems [2, 52–60]. Furthermore, the transition from phase separation to periodic patterns and back are shown to be hysteretic.

4.1. The Eckhaus stability band and transitions between periodic patterns

To address multistability of periodic patterns of different wavenumbers we first determine a stationary periodic, nonlinear solution ψ_s by using an N-mode truncated Fourier ansatz

$$\begin{equation}{\psi }_{\mathrm{s}}\left(x\right)=\sum\limits _{\vert \ell \vert {\leqslant}N}{c}_{\ell }{\mathrm{e}}^{\mathrm{i}\ell qx},\end{equation} \tag{ 9 }$$

with the conjugate symmetry ${c}_{\ell }={c}_{-\ell }^{{\ast}}$ for the real field ψ_s and c₀ = 0 as a consequence of conservation. Using the ansatz (9) and projecting equation (1) onto the modes exp(inqx) for all |n| ⩽ N yields coupled nonlinear equations for the Fourier amplitudes c_n . The nth equation reads

$$\begin{equation}0=L\left(n,0\right){c}_{n}-\sum\limits _{\vert n-m-j\vert {\leqslant}N}{c}_{n-m-j}{c}_{m}{c}_{j}+\sum\limits _{\vert n-m\vert {\leqslant}N}g\left(n,m\right){c}_{n-m}{c}_{m},\end{equation} \tag{ 10 }$$

where for convenience we have introduced:

$$\begin{align}\hfill L\left(n,k\right)& =\varepsilon -{D}_{4}{\left(k+nq\right)}^{2}-{D}_{6}{\left(k+nq\right)}^{4},\hfill \\ \hfill g\left(n,m\right)& ={\beta }_{1}-{\beta }_{2}{\left(mq\right)}^{2}-{\beta }_{3}\left(n-m\right)m{q}^{2}.\hfill \end{align} \tag{ 11 }$$

Equation (10) are solved numerically for a given set of system parameters and fixed values of ɛ and q. In equation (9) modes with even and uneven values of ℓ are coupled in accordance with the broken ± symmetry of the equation as seen in figure 2(b). The solutions of interest are dominated mostly by the amplitude c₁. While still significant, higher modes decay quickly allowing us to use N = 10 for the calculation of ψ_s.

Similar as for ψ₀(x) above we also investigate the linear stability of the stationary solution ψ_s(x) with respect to a small perturbation ψ_p(x, t). In this case one replaces in equation (5) ψ₀ by ψ_s. With ψ_s given by equation (9) the part of equation being linear in ψ₁ (5) can be solved with a Floquet ansatz for ψ₁ = ψ_p,

$$\begin{equation}{\psi }_{\mathrm{p}}\left(x,t\right)={\mathrm{e}}^{\sigma t}\enspace {\mathrm{e}}^{\mathrm{i}kx}\sum\limits _{\vert m\vert {\leqslant}N}{a}_{m}\enspace {\mathrm{e}}^{\mathrm{i}mqx},\end{equation} \tag{ 12 }$$

and real Floquet exponent k. Using this ansatz for the linear part of equation (5) and projecting this equation onto exp[i(k + nq)x], we obtain the following linear eigenvalue problem for the coefficients a_n :

$$\begin{align}\hfill \sigma {a}_{n}& ={\left(k+nq\right)}^{2}L\left(n,k\right){a}_{n}-3{\left(k+nq\right)}^{2}\sum\limits _{\vert n-m-j\vert {\leqslant}N}{c}_{n-m-j}{c}_{m}{a}_{j}\hfill \\ \hfill & \quad +\sum\limits _{\vert n-m\vert {\leqslant}N}h\left(n,m,k\right){c}_{n-m}{a}_{m},\hfill \end{align} \tag{ 13 }$$

where we have introduced in addition:

$$\begin{align}\hfill h\left(n,m,k\right)& =2{\beta }_{1}-{\beta }_{2}{\left(k+nq\right)}^{2}\left[{\left(k+mq\right)}^{2}+{\left(n-m\right)}^{2}{q}^{2}\right]\hfill \\ \hfill & \quad -2{\beta }_{3}{\left(k+nq\right)}^{2}\left(k+mq\right)\left(n-m\right)q.\hfill \end{align} \tag{ 14 }$$

By determining the eigenvalues via equation (13), we identify the parameter ranges wherein one finds stable periodic solutions ψ_s(x), i.e. wherein the eigenvalue with the largest real part σ_m (k) is still negative:

$$\begin{equation}{\mathrm{max}}_{k}\enspace \mathrm{R}\mathrm{e}\left\{{\sigma }_{m}\left(k\right)\right\}{< }0.\end{equation} \tag{ 15 }$$

The so-called Eckhaus stability boundary is determined by ${\mathrm{max}}_{k}\enspace \mathrm{R}\mathrm{e}\left\{{\sigma }_{m}\left(k\right)\right\}=0$. It is displayed as the solid line in figure 3 for the parameter set given in figure 1. Between the Eckhaus stability boundary and the dashed lines we still find numerically stationary periodic solutions, but they are unstable.

Zoom In Zoom Out Reset image size

Spatially periodic patterns ψ_s with wavenumbers chosen from within the Eckhaus band, i.e. the area enclosed by the solid line in figure 3, are all linearly stable. This type of multistability is known to occur in a number of other pattern forming systems beyond a primary instability [2, 52–61]. The insight that stationary, spatially periodic solutions predicted to occur beyond a secondary instability in systems with conserved order parameter fields share this generic property of multistability with a great number of other systems [2, 52–61] is central to this work.

In contrast to other systems with spatially periodic patterns the shape of the stable Eckhaus band at its bottom in figure 3 is not parabolic. Most of the commonly known spatially periodic patterns emerge above a primary instability out of a homogeneous basic state that is symmetric with respect to spatial reflections and translations. In these cases the Eckhaus stability band has a parabolic shape near its minimum. Here in this work the secondary bifurcation to a spatially periodic patterns takes place after a primary bifurcation to phase separation and in the state of phase separation both symmetries are absent.

4.2. Transition scenario to and between extended periodic patterns

The similarities to other nonequilibrium stripe pattern are highlighted further in the way the solutions of the CoSH+ model in equation (1) respond to changes in the control parameter or the wavenumber beyond the secondary instability. This is demonstrated in figure 4 for a generic property of periodic patterns.

Zoom In Zoom Out Reset image size

Figure 4(a) depicts two examples of transition scenarios: in the first one, A → C, a linearly stable periodic solution inside the Eckhaus band is quenched to the unstable range between the band and existence estimate. The system reacts to this change by reducing the wavenumber of the pattern via a series of transient states B to arrive at C back inside the Eckhaus band. Figure 4(b) shows the corresponding initial stable state A, an example of a transient B and the final state C. In the second scenario A' → C', we jump upwards into the unstable range and the system reacts by increasing the pattern wavenumber via transients B' to arrive at C' within the band. Figure 4(c) again shows the corresponding initial stable, the final stable periodic solution and an example of a transient B'. Note that in both scenarios the shift in ɛ is accompanied by a shift in the amplitude, which increases and decreases with the control parameter.

These transition scenarios of periodic states after parameter changes as shown in figure 4 are also observed for spatially periodic patterns in models with unconserved order parameters [53], in electrohydrodynamic convection in nematic liquid crystals [52] and in experiments with periodic patterns in Taylor-vortex flows [55]. Note that in these examples the transition scenarios have been observed for patterns that emerge out of a homogeneous basic state. In contrast, periodic patterns in our example emerge from a phase-separated state with broken translational and the reflection symmetry. Nevertheless, the transition scenarios for periodic patterns in a conserved model as shown in figure 4 are generic as in all the previous examples.

Due to the quadratic nonlinearities including β₂ and β₃ in equation (1), the maximal amplitudes ψ_max and |ψ_min| differ in both the phase separated state and the periodic state, cf figures 2(a) and (b). Plotting them as functions of increasing ɛ one finds at the transition from phase separation (circles) to the periodic state (triangles) a jump in both amplitudes as shown in figure 5. Repeating the same by starting with a periodic pattern and decreasing ɛ one finds the hysteresis of the transition as shown in figure 5.

Zoom In Zoom Out Reset image size

By increasing ɛ, the transition scenario from phase separation to periodic states is similar as shown in figure 2(c). The wavenumber of the final state is always within the Eckhaus-stable band in figure 3 and may depend on the spatial structure of phase separation, i.e. on the distribution of the domain walls just before the transition takes place.

Starting with a periodic pattern, a sufficient reduction of ɛ leads latest at the lower boundary of the Eckhaus stability band to a transition to phase separation. Since the spatial structure of the state of phase separation is multifaceted in long systems, the transition dynamics from periodic pattern to phase separation is in general less universal and rather complex.

The hysteretic transition scenario shown in figure 5 for β₁ = β₃ = 0 will not change very much for finite but not to large parameter values β₃ and β₁. However, in the range β₂ ∼ 2β₃ we find a different transition scenario as described in the next section 5.

For the parameter combination β₂ = 2β₃ model (1) has the surprising property of a gradient dynamics as described by equations (2) and (3). The stability analysis of a constant plateau solution in section 3 is independent of β₃ and so the threshold ɛ_qc remains unaltered. The dynamics and stationary solutions above ɛ_qc , however, differ greatly from those described in sections 3 and 4 for β₃ = 0. Note that for β₂ < 0 the upper instead of the lower plateau becomes unstable, but the scenario is otherwise unchanged.

We find in a range around β₂ ∼ 2β₃ a further novel transition from the state of phase separation, with many plateau areas of different length, to a hybrid state, where constant plateaus alternate with spatially periodic patterns, as shown for example in a shorter system with only two domain walls in figure 6. Similar areas are occupied by the plateaus and the periodic states according to the conservation of the field ψ(x, t), but by increasing the distance ɛ − ɛ_qc , the region occupied by a periodic state increases slightly as indicated in figure 6 from part (b) to (c).

Zoom In Zoom Out Reset image size

In the range β₂ ≃ 2β₃ the secondary periodic pattern occurs again at first locally near the domain walls via wavy overshoots caused by the higher order spatial derivative. Therefore, localized finite amplitude patterns occur already below the threshold ɛ_qc , as indicated in figure 6(a) at ɛ = 0.32 < ɛ_qc = 0.334 (for the parameters used), but do not invade the lower plateau range. Patterns localized near the domain walls in the range ɛ ≲ɛ_qc , as in figure 6(a), decay exponentially with the distance from the domain walls and the decay length increases with decreasing values of D₆.

A numerical solution of equation (1) shows that the amplitude of the periodic pattern in the middle between the two domain walls increases or decreases continuously $\propto \sqrt{\varepsilon -{\varepsilon }_{qc}}$ by increasing or decreasing ɛ. Therefore, the transition from the state of phase separation to the hybrid state is continuous and not hysteretic in spite of its local onset near the domain walls.

We have confirmed this numerical observation by a two-mode approximation ψ_hyb = ψ_* + A cos(q_max x). By choosing this ansatz we assume a long lower plateau at ψ_* and periodic boundary conditions at its ends, i.e. we neglect the effects of the domain walls at the ends of a long plateau. In this case we obtain in agreement with our numerical observation that the amplitude A is real with A² > 0 only in the range of ɛ > ɛ_qc . One finds ψ_hyb with A ≠ 0 only in the range ɛ > ɛ_qc . The constant plateaus with A = 0 corresponds to a local maximum, i.e. the functional has for finite A an lower value than for A = 0.

We discovered two new secondary bifurcation scenarios in systems described by conserved order parameter fields. Both scenarios occur by increasing the control parameter further after the primary bifurcation into active phase separation. In the first case, phase separation is stopped and completely replaced by an extended, system-filling, spatially periodic pattern. It is stable for different values of the wavenumber that belong to the so-called Eckhaus stability band, which we have determined for our model as shown in figure 3. The multistability of patterns with different wavenumbers is universal [2, 52–60]. It has been found experimentally and theoretically in a considerable number of nonequilibrium pattern forming systems with unconserved order parameter fields [2, 52–60]. As we have shown here, this multistability also applies for periodic patterns beyond a secondary bifurcation in systems with a conserved order parameter field.

The second bifurcation scenario describes a transition from active phase separation to a novel and stable nonlinear hybrid state. These hybrid solutions are spatial alternations between sections of spatially constant and spatially periodic states. Surprisingly, this spatially alternating coexistence of different spatial structures is stable. The subranges occupied by the homogeneous and periodic states may vary in space, but due to the conservation law both occupy similar areas immediately above the transition to the hybrid state. With increasing distance of the control parameter from the secondary threshold, the regions occupied by the periodic states increases slightly.

The basis of our novel results is the CoSH+ model in equation (1), which we introduced with this work. This model is a further step within a recent systematic and dynamic development for nonequilibrium systems described by conserved order parameter fields that we mention shortly. It has been found recently, that the weakly nonlinear behavior of several demixing phenomena in systems far from equilibrium are also described at leading order by the classic CH equation [38–41]. This was demonstrated exemplarily for collective demixing phenomena in chemotactic systems [38, 41], for models of cell polarization [38, 40] and MIPS in active Brownian particle systems [39]. Previously, for modeling MIPS further away from its onset an extension of the CH model to the next higher nonlinear order, i.e. ${\partial }_{x}^{2}\left(\psi {\partial }_{x}^{2}\psi \right)$ and ${\partial }_{x}^{2}{\left({\partial }_{x}\psi \right)}^{2}$, was suggested with the introduction of active model B+ [42, 44]. One of the two nonlinear extensions contributes to the prefactor D₄ + β₂ ψ of ${\partial }_{x}^{4}\psi$ in the active model B+ and in the CoSH+ model (1). This prefactor may change its sign at larger amplitudes of the order parameter field. In this range the contribution ${D}_{6}{\partial }_{x}^{6}\psi$ is according to our perturbation theory of the same order as the nonlinear extensions. It is indispensable in this range to avoid dangerous large wavenumber instabilities after sign change of D₄ + β₂ ψ. Our model in (1) includes the sixth order derivative, similar as the CSH model in reference [46], and it is in this sense an extension of the conserved SH model by generic leading order nonlinear gradient terms. Accordingly, we call it in section 2 the CoSH+ model.

The spatially varying sign change of the prefactor D₄ + β₂ ψ in equation (1) is the origin of the two novel bifurcation scenarios. This sign change is qualitatively similar to the Lifshitz point in pattern formation, found near the so-called oblique-roll instability in electroconvection in nematic liquid crystals, where the coefficient of a higher order spatial derivative undergoes a similar change in sign, but spatially homogeneously [62–65].

Secondary bifurcations of, for example, stripe patterns are known in the field of classical pattern formation with unconserved order parameter fields [2, 66–72]. However, to the best of our knowledge there is so far no further example of a secondary bifurcation that evolves in a conserved nonequilibrium system out of active phase separation into stable nonlinear periodic pattern.

The existence of stability bands of periodic patterns, which we found for our first secondary bifurcation scenario, is a universal and robust phenomenon in nonlinear physics. To add or to remove a periodic unit (wavelength) of a pattern requires strong local deformations of the pattern, which cannot be driven by small inherent fluctuations. Therefore, in the absence of strong external deformations, periodic patterns remain stable at a wavenumber within the Eckhaus band. It is an important result of our work that this robust behavior of periodic patterns applies also to patterns beyond the secondary bifurcation in conserved systems. This robust stability of periodic patterns in conserved systems is also highly relevant for modeling cell division, where it possibly helps to find the center of a cell with great certainty, before it divides into two equally sized daughter cells carrying the same genetic information [73, 74].

In some former studies active phase separation was changed into a bifurcation to spatially periodic patterns by the introduction of additional kinetic effects that violate conservation laws and therefore change the bifurcation type [75–77]. For such a kinetically changed primary bifurcation from phase separation to periodic patterns, the notion arresting phase separation was coined [77]. Note that this is a crucial difference between these examples and our work. The transition to periodic patterns occurs in our work via a secondary bifurcation under the retention of the conservation conditions.

An extension of the model in equation (1) to two spatial dimensions is a further promising perspective. It is very likely that our results obtained for one spatial dimension lead in two dimensions to further interesting secondary bifurcation scenarios, for example, due to the broken up–down symmetry one expects at even lower values of the control parameter a secondary bifurcation to hexagonal patterns.

For pattern forming systems with unconserved order parameters a number of interesting finite size, inhomogeneity, noise and boundary effects on periodic patterns are known [49–51, 69, 73, 78–81]. In contrast, such effects on periodic patterns in conserved systems, including those beyond secondary bifurcations, are unexplored so far.

Candidates of specific systems where the described phenomena are expected as well, are, for example, generalized chemotactic systems (see e.g. [41]), Brownian particle systems showing MIPS in the presence of quorum sensing [82] and cell biology, where periodic patterns play a key role in finding the cell center so that cells divide into equal parts with the same amount of genetic information [74].

Support by the Elite Study Program Biological Physics and discussions with Andre Förtsch, Samuel Grimm, Markus Hilt, Mirko Ruppert and Winfried Schmidt are gratefully acknowledged.

The data that support the findings of this study are available upon reasonable request from the authors.

In this appendix, we show as in reference [39] that contributions to our model (1), which go beyond the classic Cahn–Hilliard model, vanish in the limit of small values of the distance ɛ from the onset of phase separation. For this we rescale time, space and amplitude in equation (1) such that the part of the classic CH in equation (1) becomes dimensionless: t' = ɛ² t, ${x}^{\prime }=\sqrt{\varepsilon }x$ and ${\psi }^{\prime }=\psi /\sqrt{\varepsilon }$. In addition we choose β₁ = ɛβ₁', because the quadratic term at leading order can be removed by adding a constant. This allows us to rewrite equation (1) with a dimensionless classic CH part in the following form:

$$\begin{align}\hfill {\partial }_{{t}^{\prime }}{\psi }^{\prime }& =-{\partial }_{{x}^{\prime }}^{2}\left[{\psi }^{\prime }+{\partial }_{{x}^{\prime }}^{2}{\psi }^{\prime }-{\psi }^{\prime 3}\right]\hfill \\ \hfill & \quad +\sqrt{\varepsilon }{\partial }_{{x}^{\prime }}^{2}\left[{\beta }_{2}{\psi }^{\prime }{\partial }_{{x}^{\prime }}^{2}{\psi }^{\prime }+{\beta }_{3}{\left({\partial }_{{x}^{\prime }}{\psi }^{\prime }\right)}^{2}+{\beta }_{1}^{\prime }{\psi }^{\prime 2}\right]+\varepsilon {\partial }_{{x}^{\prime }}^{6}{\psi }^{\prime }.\hfill \end{align} \tag{ A.1 }$$

Therefore, all contributions which do not belong to the leading order classic Cahn–Hilliard model vanish in the limit of small ɛ. This behavior explains, why in figure 1(a) the solution of model (1) has in the limit of small ɛ the same form as domain wall solutions of the classic CH model.

We perform spatio-temporal simulations of equation (1) on a spatial domain (0, L) with periodic boundary conditions. We focus on stationary solutions and use a pseudo-spectral method with semi-implicit Euler time-stepping, treating the linear differential operator implicitly and the nonlinear contributions explicitly. The interval is discretized uniformly into an N-point grid with resolution Δx = L/N, where N (i.e. the number of Fourier modes) is adapted according to the number of periods per system length. For figures 2, 4 and 5 we used L = 100, N = 256 with time-step size Δt = 1 × 10⁻³. To resolve the more rapid spatial variations in the functional case displayed in figure 6 we used N = 1024 with Δt = 1 × 10⁻⁴.