††thanks: Present address: Kavli Institute for Theoretical Physics, University
of California Santa Barbara, Santa Barbara, CA 93106, USA

#  Coarsening and wavelength selection far from equilibrium: a unifying
framework based on singular perturbation theory

Henrik Weyer  Fridtjof Brauns<EMAIL_ADDRESS>Arnold Sommerfeld Center for
Theoretical Physics and Center for NanoScience, Department of Physics, Ludwig-
Maximilians-Universität München, Theresienstraße 37, D-80333 München, Germany
Erwin Frey<EMAIL_ADDRESS>Arnold Sommerfeld Center for Theoretical Physics and
Center for NanoScience, Department of Physics, Ludwig-Maximilians-Universität
München, Theresienstraße 37, D-80333 München, Germany  Max Planck School
Matter to Life, Hofgartenstraße 8, D-80539 Munich, Germany

(August 27, 2024)

###### Abstract

Intracellular protein patterns are described by (nearly) mass-conserving
reaction–diffusion systems. While these patterns initially form out of a
homogeneous steady state due to the well-understood Turing instability, no
general theory exists for the dynamics of fully nonlinear patterns. We develop
a unifying theory for nonlinear wavelength-selection dynamics in (nearly)
mass-conserving two-component reaction–diffusion systems independent of the
specific mathematical model chosen. Previous work has shown that these systems
support an extremely broad band of stable wavelengths, but the mechanism by
which a specific wavelength is selected has remained unclear. We show that an
interrupted coarsening process selects the wavelength at the threshold to
stability. Based on the physical intuition that coarsening is driven by
competition for mass and interrupted by weak source terms that break strict
mass conservation, we develop a singular perturbation theory for the stability
of stationary patterns. The resulting closed-form analytical expressions
enable us to quantitatively predict the coarsening dynamics and the final
pattern wavelength. We find excellent agreement with numerical results
throughout the diffusion- and reaction-limited regimes of the dynamics,
including the crossover region. Further, we show how, in these limits, the
two-component reaction–diffusion systems map to generalized Cahn–Hilliard and
conserved Allen–Cahn dynamics, therefore providing a link to these two
fundamental scalar field theories. The systematic understanding of the length-
scale dynamics of fully nonlinear patterns in two-component systems provided
here builds the basis to reveal the mechanisms underlying wavelength selection
in multi-component systems with potentially several conservation laws.

##  I Introduction

Across many non-equilibrium systems, small constituents self-organize into
macroscopic patterns on much larger length scales. One of the most intriguing
aspects of such patterns far from equilibrium is the emergence of intrinsic
length scales independent of the system size or other spatial cues but
entirely determined by the local interaction of the constituents. Examples
include chemical systems far from equilibrium [1, 2, 3, 4], especially
intracellular protein patterns [5, 6, 7, 8], collective states in active
matter [9, 10, 11, 12, 13, 14, 15, 16], and phase separation of chemically
active species [17, 18, 19] or species undergoing population dynamics [20, 21,
22]. In contrast, phase-separation processes approaching thermodynamic
equilibrium usually develop toward full separation via a continuous growth of
the average domain size, a process termed “coarsening”. The coarsening process
can be interrupted, and a wavelength selected in close-to-equilibrium
processes if the system allows for long-range interactions [23, 24, 25].
Recently, the mechanisms underlying length-scale selection were explained for
non-equilibrium extensions of classical phase separation, advancing the
understanding of active phase separation important for intracellular
condensates and active matter systems [12, 18, 21, 19].

Reaction–diffusion systems are governed by (chemical) reactions that induce
transitions between particle states with different diffusivities, instead of
physical interactions that induce the demixing of phases. Originally, Alan
Turing proposed such systems to explain the patterning during morphogenesis,
i.e., the development of organisms [1]. By now, reaction–diffusion systems are
found to describe diverse biological processes, for example, intracellular
protein pattern formation [26, 27, 8, 7, 28] or signaling via trigger waves
[29]. In the intracellular context, the timescale of protein production and
degradation is long compared to the timescale of pattern formation, and the
total number of proteins can be assumed to be (approximately) constant. The
proteins just switch between different states, for example, a membrane-bound
and a cytosolic state. These dynamics are captured by mass-conserving
reaction–diffusion (McRD) systems which describe the concentrations of the
different protein states [30, 31, 32, 33, 34, 5, 35, 36, 7, 37, 38]; for an
introduction to the theory of McRD systems, see the lecture notes Ref. [39].
Importantly, McRD systems do not conserve each individual component but the
total number density of each protein summing over the concentrations of its
conformational states. McRD systems also describe, for example, granular media
[40], precipitation patterns [41], and braided polymers [42].

Two-component reaction–diffusion (2cRD) systems _lacking_ mass conservation
generally exhibit patterns with an intrinsic length scale [43, 2, 3, 4, 44].
In contrast, the fully nonlinear patterns in two-component _mass-conserving_
reaction–diffusion (2cMcRD) systems have been observed to undergo coarsening
until completion, that is, the pattern length scale grows until the system
size is reached [30, 45, 36, 46]. This distinct behavior of two-component
reaction–diffusion systems without and with mass conservation calls for a
detailed analysis of the length-scale dynamics in nearly mass-conserving 2cRD
systems. We started to investigate this question in Ref. [46] by proposing
that wavelength selection in 2cRD systems can be understood as an interruption
of the coarsening process due to (weak) source terms that break strict mass
conservation. Together with the splitting of pattern domains at larger
wavelengths, this gives rise to a broad range of stable wavelengths. Here, we
systematically analyze under what conditions this reasoning applies in general
2cMcRD systems and how the final pattern wavelength is selected dynamically.

The coarsening process observed in mass-conserving 2cRD systems is reminiscent
of the typical dynamics found in phase-ordering and phase-separation kinetics.
In these systems, which are close to thermal equilibrium, the dynamics is
governed by gradient flows in a free energy landscape toward the minimum of
the respective free energy functional. The free-energy cost due to interfaces
between the different phases leads to dynamics that continuously minimizes the
surface area of the interfaces (Model A/B dynamics [47]) [48]. Thus, small
domains with a high surface-to-bulk ratio collapse in favor of larger domains,
and the characteristic pattern length scale grows uninterruptedly until the
fully phase-ordered or phase-separated state is reached [49, 50, 51, 52, 53].
Intriguingly, the uninterrupted coarsening in (inherently far-from-
equilibrium) 2cMcRD systems exhibits regimes where the coarsening resembles
either Cahn–Hilliard dynamics (bulk-diffusion-controlled phase separation
[54]) [46, 55] or conserved Allen–Cahn dynamics (interface-controlled kinetics
[53, 50, 56]) [57, 58, 59], two classical models of phase separation.
Moreover, stationary states of 2cMcRD systems can be analyzed similarly to a
Maxwell construction [38], and close to the supercritical onset of pattern
formation 2cMcRD systems reduce to an amplitude equation which agrees with the
Cahn–Hilliard equation [60]. While some 2cMcRD systems with a specific
mathematical form of the reaction term allow for an abstract mapping onto a
gradient flow for an effective free-energy functional [61, 62], a general and
unifying understanding of the similarities between reaction–diffusion and
phase-separation dynamics is lacking.

Similarly, the changed phenomenology giving rise to wavelength selection in
the presence of source terms is reminiscent of interrupted coarsening observed
if non-equilibrium extensions are included in classical phase-separation
dynamics. In binary phase separation, for example, the coarsening process can
be arrested if chemical reactions are introduced that convert particles from
one of the phase-separating species into the other [17]. This mechanism has
been of increasing interest in recent years to describe intracellular
condensates [63, 18]. Moreover, in the fields of active matter and
nonreciprocal systems, wavelength selection in non-equilibrium settings is of
increasing interest and interrupted coarsening has been frequently studied as
well [12, 64, 65, 21, 66]. It would therefore be telling to find out whether
common principles underlie coarsening and wavelength selection in these
systems taking the form of active phase separation and 2cRD models.

Several approaches have been developed to analyze the length-scale dynamics in
systems of (active) phase separation and reaction–diffusion systems. Close to
a _supercritical_ onset of pattern formation, the _amplitude equation_
formalism provides a powerful tool to study the pattern properties and
dynamics, including the question of wavelength selection [43, 4, 24]. As this
approach critically depends on a small pattern amplitude, it is not applicable
to fully developed patterns of large amplitude. If the patterns only evolve on
length scales much larger than the typical pattern wavelength, an approach
based on _phase equations_ , which are conceptually closely related to
amplitude equations, can be applied to obtain the long-time dynamics in the
highly nonlinear regime [67, 68, 69]. In both the amplitude and phase equation
approach, conservation laws play a critical role and must be accounted for
explicitly [43, 70].

(Nearly) mass-conserving two-component reaction–diffusion systems as well as
thermodynamic systems exhibiting phase separation generally show a
_subcritical_ onset of pattern formation such that the regime of (spontaneous)
lateral instability (termed as _spinodal_ regime in the language of phase
separation) is surrounded by a multistable regime which shows stimulus-induced
pattern formation (_nucleation-and-growth_ regime) [38, 24]. In these
reaction–diffusion systems, the coarsening process proceeds by mass
redistribution, which is fastest on the shortest distances, that is, between
neighboring pattern domains [46]. The dominant dynamic process thus acts on
the length scale of the patterns. As a result, amplitude and phase equations
cannot capture the dynamics of the highly nonlinear patterns in nearly mass-
conserving reaction–diffusion systems. The same holds for systems describing
(close-to-equilibrium) phase separation if the dynamics is not analyzed close
to the critical point. Therefore, different methods have been developed and
have been first used to describe the long-time dynamics of close-to-
equilibrium phase separation. In the mathematical literature, _slow-manifold
theory_ is applied to phase-separating systems [71, 72, 73, 74, 75]. Moreover,
Lifshitz, Slyozov, and Wagner [49, 50] developed their classical theory of
Ostwald ripening (LSW theory), building on the physical properties of single
(quasi-)stationary droplets [quasi-steady-state (QSS) approximation]. The
resulting theory describes the particle exchange between droplets of different
sizes of an immiscible minority phase sparsely distributed in the majority
phase. LSW theory was extended to analyze wavelength selection in systems
introducing chemical reactions between the species undergoing phase separation
[63, 18, 21]. In contrast, this physical reasoning has not been employed to
describe reaction–diffusion dynamics. Instead, in 2cRD systems _singular
perturbation theory_ [57] was adopted to analyze wavelength selection [76, 77,
59, 58, 44]. Also, numerical bifurcation analysis was applied to discuss the
wavelength of stable stationary patterns when tuning parameters [78, 66]. In
Ref. [46], we have employed a QSS approximation inspired by LSW theory in
(nearly) mass-conserving reaction–diffusion systems, gaining an in-depth
understanding of the coarsening dynamics and the stabilization of finite
pattern wavelengths in these systems. However, a systematic justification for
the QSS approximation in reaction–diffusion systems is still missing.
Providing such a justification in terms of singular perturbation theory
supported by physical arguments is the central goal of the present work.

Our theoretical analysis in Ref. [46] has shown that the length-scale dynamics
in reaction–diffusion systems can be understood—analogous to active phase
separation—as a coarsening process which is arrested by counteracting
processes above a particular length scale. Coarsening is explained on the
basis of a mass-competition instability between neighboring pattern domains
(the droplets in LSW theory) and wavelength selection by suppression of this
instability. Here, we complement the physical reasoning by a singular
perturbation analysis for (weakly) mass-conserving 2cRD systems, which
provides explicit expressions for the growth rates of the different processes
involved in the mass-competition instability. Previous works have focused on
specific mathematical forms of the reaction term [77, 44, 79, 80] and the
reaction-limited regime [59, 58]. Our results are independent of the specific
mathematical form of the reaction terms as they build on geometric reasoning
in phase space [38]. Moreover, our results apply in both diffusion- and
reaction-limited regimes, including the crossover between them. These
mathematical results are explained by a detailed analysis of how mass is
exchanged between pattern domains, revealing a systematic link of the
reaction–diffusion dynamics to diffusion-limited (bulk-diffusion-controlled)
and reaction-limited (interface-controlled) phase-separation kinetics. We then
derive the wavelength dynamics of patterns in large systems from the rate
expressions for the competition between neighboring peaks. This dynamics
reveals that the _threshold_ of interrupted coarsening selects the final
pattern wavelength. Taken together, our results underline that wavelength
selection in both the reaction–diffusion and phase-separating systems is
driven by the mass exchange between domains and the mass exchange with a
reservoir described by the source terms. Consequently, we expect that the
analysis provides the basis to analyze wavelength-selection dynamics in multi-
component reaction–diffusion and active matter systems governed by one, or
possibly multiple, (approximate) conservation laws.

The remainder of this paper is organized as follows. We first describe in
Section II the phenomenology of the pattern dynamics we set out to explain,
that is, the phenomenology of coarsening and its arrest. In Section III, we
then focus on the approximate conservation law governing 2cRD system with weak
source terms. We also introduce two standard models exhibiting phase
separation—(generalized) Cahn–Hilliard (CH) and conserved Allen–Cahn (cAC)
models—which account for mass conservation by either a local or a global
constraint. These classical systems also serve to compare the
reaction–diffusion with (active) phase-separation dynamics. Section IV
discusses the stability properties of the homogeneous steady state which
highlights the connections found between these three models in their initial
pattern-forming instabilities. Afterward, we construct the stationary patterns
of all three systems in a unifying picture (Sec. V). On the basis of this
classification of stationary patterns, we then focus in section VI on strictly
mass-conserving systems. The mass-competition instability is described, which
underlies the self-amplifying mass transport between neighboring pattern
domains, and which causes the growth of larger droplets at the expense of
smaller droplets. This instability is the elementary motive underlying the
uninterrupted coarsening process observed in the mass-conserving systems, and
we use the derived growth rate of the instability in section VII to obtain a
scaling law for the time evolution of the coarsening process. This completes
the description of uninterrupted coarsening in mass-conserving systems, and we
turn in section VIII to discuss the influence of weak source terms that break
the strict conservation law. From the growth rate of the mass-competition
instability under the influence of weak source terms, we explain in this
section the central criterion for the wavelength selected by interrupted
coarsening. In particular, we discuss how the suppression of the mass-
competition instability by weak source terms above a threshold pattern
wavelength determines the wavelength-selection dynamics. After explaining the
underlying processes, the comparison with several examples analyzed
numerically verifies the found relations. The discussion of our findings and
future applications are found in Sec. IX.

Throughout this work, the main text focuses on the discussion of the results
from the singular perturbation analysis, i.e., their implications and physical
interpretations while the formal mathematical derivations are deferred to the
appendices. Moreover, we restrict the analysis to one-dimensional systems to
avoid mathematical complications resulting from the system’s geometry that
would hamper the analysis of the underlying mechanisms. Additional effects
appearing in higher dimensions are mentioned in the discussion.

Figure 1:  The phenomenology of coarsening. (a) The kymograph shows the
density profile $\rho$ (grayscale, black corresponds to large densities) for
the simulation of the second model of Ref. [30], a peak-forming 2cMcRD system
defined by ${f=\rho^{2}\eta-\rho}$, where $\rho=u+v$ is the total density and
$\eta=v+du$ is the mass-redistribution potential (see Sec. III.1). The
parameters are ${D_{u}=1}$, ${D_{v}=10}$ and ${\bar{\rho}=10}$. (b) A single
peak participating in the coarsening process [(red dashed line in panel (a);
colored (dark to light gray) profiles in (b) correspond to the times indicated
in (a)] is well described by the stationary peak profile (black dashed,
analytical form was derived in Ref. [30]). (c,e) A mesa-forming model shows
much slower coarsening [cubic model ${f=(\eta-\rho^{3}+\rho)/(1-d)}$; see Sec.
VIII.4.1 and Appendix H.1]. At low average total density $\bar{\rho}$ [panel
(c)], coarsening proceeds mainly via competition [blue (dark gray) arrows]
while coalescence [red (light gray) arrows] dominates at large average density
[panel (e)]. The parameters are $D_{u}=1$, $D_{v}=10$ and $\bar{\rho}$ as
indicated. (d) Again, the single interface profiles are well approximated by
the stationary profile. All simulations employ a system length $L=1000$ and
periodic boundary conditions.

##  II Phenomenology of the coarsening process and its arrest

Let us start with an overview of the typical phenomenology of (nearly mass-
conserving) 2cRD systems before getting into their detailed mathematical
analysis. Consider the general reaction–diffusion dynamics of two species $u$
and $v$. We decompose the governing equations into a mass-conserving core
system and source terms of strength $\varepsilon$ (cf. [81, 82, 46])

$\displaystyle\partial_{t}u(x,t)$
$\displaystyle=D_{u}\,\nabla^{2}u+f(u,v)+\varepsilon s_{1}(u,v)\,,$ (1a)
$\displaystyle\partial_{t}v(x,t)$
$\displaystyle=D_{v}\,\nabla^{2}v-f(u,v)+\varepsilon s_{2}(u,v)\,,$ (1b)

defined on a $D$-dimensional spatial domain $\Omega$ with no-flux boundary
conditions for both $u$ and $v$. Although we here introduce the general
system, our analysis in the following chapters focuses on the one-dimensional
case ${D=1}$. For specificity, we choose the relative diffusion constant as
${d:=D_{u}/D_{v}<1}$.

Thus, $u(x,t)$ describes the density of a slowly diffusing species, which we
will interpret in the context of intracellular pattern formation as proteins
attached onto a cell membrane. In contrast, $v(x,t)$ describes the density of
the fast-diffusing species, which may be a different conformational state of
the proteins that has detached from the membrane and undergoes fast diffusion
in the cytosol of the cell (${D_{v}\gg D_{u}}$ for intracellular pattern
formation). The reaction term $f(u,v)$ accounts for the mass-conserving
conversion between these two species, that is, between the two conformational
states. In this intracellular context, $f$ is typically given as attachment-
detachment kinetics ${f(u,v)=a(u)v-b(u)u}$, where the attachment rate $a(u)$
and detachment rate $b(u)$ depend on the membrane concentration $u$. These
terms model, for instance, nonlinear recruitment and enzyme-driven detachment
[34, 39, 83].

As the chemical reactions represented by $f(u,v)$ correspond to conversions
between the two protein states, they locally conserve the total density
${\rho(x,t)=u+v}$ of proteins at each point $x$ in the domain $\Omega$.
Consequently, in the case of mass conservation (${\varepsilon=0}$), the
density $\rho(x,t)$ changes only due to the redistribution of proteins within
the system. We stress that—different from conserving systems of Cahn–Hilliard
type—mass conservation in the 2cMcRD system [Eq. (1) with ${\varepsilon=0}$]
only holds for the total density $\rho(x,t)$ while the dynamics of the single
protein states $u$ and $v$ is not conserving. In contrast, $s_{1}$ and $s_{2}$
are source terms that break strict mass conservation of the total density
$\rho$; their strength is given by the dimensionless parameter $\varepsilon$.
These source terms, for instance, account for the production or degradation of
proteins, and they change the total number of proteins.

###  II.1 Long-time dynamics of 2cMcRD systems

We study two-component mass-conserving reaction–diffusion (2cMcRD) systems for
their ability to form spatially heterogeneous protein concentration patterns.
Such patterns form if homogeneous steady states (HSS) with uniform protein
concentrations within the whole domain exhibit a lateral instability. This
instability drives the exponential growth of small fluctuations in the protein
concentration fields around the homogeneous steady state. In 2cMcRD systems,
this process results in the formation of peak or mesa patterns (Fig. 1). Both
types of patterns subsequently undergo a coarsening process, but typically on
a much larger timescale, during which smaller peaks or mesas vanish while
larger domains grow.

Figure 1(a) exemplifies this dynamics for peak-forming patterns. Small
perturbations around the HSS grow exponentially and rapidly form a series of
density peaks. This is followed by a slow coarsening process which continues
until just one peak is left. During the coarsening process, the dynamics of
each individual peak is well described by a (quasi-)stationary profile that
slowly evolves over time [Fig. 1(b)]. The elementary motif of the coarsening
process is the competition for mass between peaks: We observe that between two
peaks of almost equal size, mass is redistributed from the smaller peak toward
the larger peak. This process destabilizes the symmetric state of equally
sized peaks and leads to the collapse of the smaller peak. Because diffusive
mass redistribution from one to the other peak is slower over longer length
scales, the competition is strongest between _neighboring_ peaks. We call this
elementary instability between neighboring peaks the _mass-competition
instability_ [46].

Figure 2:  The coarsening scenarios. (a) The competition scenario is
illustrated for mesa (top) and peak (bottom) patterns: One of the two
mesas/peaks is growing, and the other is shrinking [indicated by (red) arrows]
due to competition for mass. (c) Coarsening can also proceed by coalescence of
mesas (top) or peaks (bottom): Neighboring high-density domains move toward
each other until they merge. This coalescence scenario can also be understood
as competition for mass between low-density regions (‘troughs’) of the
mesa/peak patterns.

Similar dynamics as for peak patterns occurs in mesa-forming systems [Fig.
1(c,d)]. However, the coarsening process is much slower, and at high average
protein densities, it is dominated by the _coalescence_ , i.e. relative
motion, of mesas rather than their competition for mass [examples marked by
red (light gray) arrows in Fig. 1(e)]. Competition for mass, in contrast, is
(mainly) observed in peak-forming systems as well as at lower densities in
mesa-forming systems [blue (dark gray) arrows in Fig. 1(e)]. At coalescence
events, the mesa number decreases not because a mesa loses all its mass and
vanishes but because mesas shift their positions and merge. The interaction of
fronts resulting in coarsening dynamics has been analyzed in diverse mesa-
forming systems [51, 71, 72, 73, 74, 75, 84]. We will use a linear stability
analysis of periodic patterns to discuss the mass-competition instability and
the resulting dynamics in 2cMcRD systems.

As in the case of competition for mass between the mesas/peaks, the
coalescence process is also based on an instability of the periodic patterns:
A high-density mesa which is, say, closer to its left neighbor than to its
right neighbor moves even farther to the left until it merges with its left
neighbor. Figure 2 compares the different coarsening scenarios where
coarsening is driven by competition for mass [panel (a)] and coalescence
[panel (b)], respectively. This comparison illustrates that the dynamics in
the coalescence scenario is actually driven by (inverted) mass competition of
the low-density domains, that is, competition for (negative) mass between the
‘troughs’ of the pattern. For example, if a peak moves toward a neighboring
peak [Fig. 2(b), bottom], the trough on one side grows, and the trough on the
other side collapses similarly as the larger high-density mesa grows while the
smaller one collapses during the mesa-competition process [Fig. 2(a), top].
Consequently, both competition and coalescence are driven by destabilizing
mass redistribution between pattern domains, i.e., a mass-competition
instability between domains of high or low density.

In the next sections, we analyze mass competition in an isolated compartment
containing only two ‘half’ peaks or mesas, that is, one period of the
stationary pattern (gray vertical lines in Fig. 2). This corresponds to the
elementary motif of nearest-neighbor competition and it allows us to isolate
each of the two coarsening scenarios by placing the no-flux boundaries of the
compartment such that they reflect the symmetry of the respective perturbation
mode of the pattern [cf. (red) arrows in Fig. 2.

Below, we will use linear stability analysis to describe the evolution of
small perturbations from the periodic patterns.111Note that the linear
stability analysis is performed for the dynamics linearized around fully
nonlinear periodic patterns. This is different from the typical Turing
analysis where the dynamics is linearized around the HSS. Mass competition and
coalescence then correspond to distinct unstable eigenmodes of the linearized
dynamics close to the periodic patterns. In this linear regime, deviations
from the periodic pattern grow exponentially ${\sim\operatorname{e}^{\sigma
t}}$ in time, where $\sigma$ is the growth rate of the corresponding
eigenmode. The sign of $\sigma$ determines whether the associated coarsening
mode is stable (${\sigma<0}$) or unstable (${\sigma>0}$), i.e., whether the
perturbations decay or grow.222In general, $\sigma$ is complex, and the sign
of the real part of the growth rate $\Re[\sigma]$ determines the stability.
For the competition and coalescence modes, it will turn out that $\sigma$ is
real. We will learn from this analysis why 2cMcRD systems always exhibit
uninterrupted coarsening. Moreover, the magnitudes of the growth rates
$\sigma$ for the different coarsening scenarios illustrated in Fig. 2 will
explain that (almost) no coalescence events occur for peak patterns because
peak competition is much faster than peak coalescence. Moreover, for mesa
patterns, the magnitude of both growth rates for mesa competition and mesa
coalescence will turn out to be strongly reduced compared to the rate of peak
competition, revealing why the coarsening process is much slower for mesa than
for peak patterns. In addition, we will show that the relative strength of
mesa competition and mesa coalescence depends on the average density
$\bar{\rho}$, explaining why the coarsening process for mesa patterns is
dominated by competition at low average densities $\bar{\rho}$ and coalescence
at high average densities [cf. Fig. 1(c,e)]. To describe the coarsening
dynamics in a large system containing many peaks or mesas (cf. Fig. 1), we use
the growth rate $\sigma$ to determine the temporal law of the coarsening
dynamics, i.e., the time evolution of the average peak (mesa) separation or of
the peak (mesa) number (see Ref. [46] and Sec. VII).

Figure 3: Weak source terms interrupt the coarsening process by suppression of
the mass-competition instability at large wavelengths. (a) The kymograph
obtained from numerical simulations shows mesa patterns (grayscale, $-1$ to
$1$) that collapse and coalesce for an intermediate period of time [blue (dark
gray) and red (light gray) arrow]. At later times ${t\gtrsim 10^{6}}$, the
domains merely rearrange into the final periodic stationary state [green (two
right-most) arrows]. As an example, the cubic model
${f=(\eta-\rho^{3}+\rho)/(1-d)}$ ($\rho$ and $\eta$ defined as in the caption
of Fig. 1) was simulated with the source terms ${(s_{1},s_{2})=(0,-\rho)}$ on
a domain of length ${L=20\,000}$ (see Sec. VIII.4.1 and Appendix H.1). The
parameters are ${D_{u}=1}$, ${D_{v}=10}$ and ${\varepsilon=10^{-6}}$. (b) A
given source strength $\varepsilon$ is strong enough to stabilize patterns at
length scales ${\Lambda>\Lambda_{\mathrm{stop}}(\varepsilon)}$ (top edge of
the shaded region, blue line). Patterns with a shorter characteristic length
undergo coarsening [(blue-) shaded region].

###  II.2 Weak source terms interrupt coarsening

As we have just demonstrated phenomenologically, strictly mass-conserving two-
component reaction–diffusion systems exhibit uninterrupted coarsening: The
pattern length scale grows until it reaches the system size, and no intrinsic
length scale is selected. In contrast, classical non-mass-conserving 2cRD
systems, including production and degradation terms, result in patterns with a
fixed length scale [3, 4, 44]. To bridge the phenomenological gap between
mass-conserving 2cRD systems and those without mass conservation, we analyze
systems with weak source terms. Such weak coupling to a reservoir is also a
first step to generalize the analysis toward reaction–diffusion systems with
more than two components, which are important as models of pattern formation
in complex biochemical reaction networks [5, 7, 6, 85].

Figure 3(a) shows the time evolution of a mesa-forming model with weak source
terms. Initially, the mesa pattern develops out of perturbations around the
homogeneous steady state. These mesas then undergo coarsening [blue (second
from the left) and red (left-most) arrow in Fig. 3(a)], which halts after some
time. After that, no more mesas (and troughs) collapse, and the remaining
mesas rearrange themselves into a periodic pattern. During this rearrangement,
smaller mesas grow, and larger mesas shrink so that their masses balance out
[green (two right-most) arrows in Fig. 3(a)]. These observations suggest that
the mass-competition instability must be suppressed by the weak source terms,
and the mass competition between two neighboring peak/mesas no longer results
in the collapse of the smaller peak/mesa. Rather, the mass-competition process
is reversed and stabilizes the symmetrical configuration of two equally sized
domains.

The slowdown of the coarsening process as it progresses [see Fig. 1] indicates
that the mass-competition instability weakens as the typical distance
$\Lambda$ between peaks/mesas increases. Moreover, we will show below that
source terms of strength $\varepsilon$ [cf. Eqs. (1)] lead to stabilizing
effects ${\sim\varepsilon}$. Thus, a critical pattern length scale
${\Lambda_{\mathrm{stop}}(\varepsilon)}$ exists where the source terms are
sufficiently strong to suppress the mass-competition instability [see Fig.
3(b)]. For ${\Lambda<\Lambda_{\mathrm{stop}}(\varepsilon)}$ the growth rate
${\sigma^{\varepsilon}(\Lambda)}$ of at least one coarsening mode under the
influence of weak source terms stays positive, giving rise to a mass-
competition instability, and driving coarsening. For patterns of larger length
scale ${\Lambda>\Lambda_{\mathrm{stop}}(\varepsilon)}$ the growth rates for
all coarsening scenarios will fulfill ${\sigma^{\varepsilon}<0}$ which
signifies that deviations from the symmetric pattern of equally sized peaks or
mesas are not amplified any longer but relax back toward the symmetric state.
Below, we determine this threshold $\Lambda_{\mathrm{stop}}(\varepsilon)$
where $\sigma^{\varepsilon}=0$.

##  III Models

In the last section, we introduced the general form of (nearly) mass-
conserving two-component reaction–diffusion dynamics. Here, we will rewrite
this dynamics to make the role of mass conservation explicit. Then, we discuss
the Cahn–Hilliard and conserved Allen–Cahn models for phase separation and
recast these models into a form that underlines their structural similarity to
mathematical descriptions of reaction–diffusion dynamics. Importantly, by
pointing out how the three models are akin at a mathematical level, we are not
claiming that these models describe the same physical processes, but rather we
are trying to explain in what sense the phenomena they exhibit are related.

###  III.1 Two-component reaction–diffusion system

We consider the general 2cRD dynamics Eqs. (1). If the source terms are
switched off by setting ${\varepsilon=0}$, the system locally conserves the
total density ${\rho(x,t)=u+v}$ of proteins. Including the source terms one
has in a well-mixed reaction compartment

$\partial_{t}\rho(t)=\varepsilon s_{\mathrm{tot}}(u,v)\,,$ (2)

with the total source term ${s_{\mathrm{tot}}:=s_{1}+s_{2}}$.

In contrast, in a spatially extended system the total density $\rho(x,t)$
additionally changes by redistribution of proteins, and its time evolution is
governed by the (modified) continuity equation [add up Eqs. (1), and, for ease
of notation denote $s_{\mathrm{tot}}(u(\rho,\eta),v(\rho,\eta))$ by
$s_{\mathrm{tot}}(\rho,\eta)$]

$\partial_{t}\rho=D_{v}\nabla^{2}\eta+\varepsilon
s_{\mathrm{tot}}(\rho,\eta)\,.$ (3)

Here we defined ${\eta(x,t):=v(x,t)+d\,u(x,t)}$, termed the _mass-
redistribution potential_ [30, 45, 42, 38]. In the absence of source terms
(${\varepsilon=0}$) the spatially averaged density
${\bar{\rho}=1/|\Omega|\int_{\Omega}\mathrm{d}x\,\rho}$ is conserved. Thus,
the mass-conserving chemical reactions entail that the _total_ protein density
$\rho$ follows a locally conserved dynamics given by a continuity equation,
akin to the Cahn–Hilliard equation and gradient dynamics of scalar field
theories for conserved order parameters (“Model B”) [54, 47, 24]. The mass-
redistribution potential here plays a similar role as the chemical potential
in the Cahn–Hilliard equation.

Unlike the chemical potential in thermal equilibrium systems, however, the
mass redistribution potential is _not_ given by the gradient of a free energy
functional but follows its own time evolution [using Eqs. (1) and the
definition $\eta=v+du$]

$\partial_{t}\eta=\big{(}D_{v}{+}D_{u}\big{)}\nabla^{2}\eta-
D_{u}\nabla^{2}\rho-\tilde{f}(\rho,\eta)+\varepsilon\big{(}ds_{1}{+}s_{2}\big{)}\,,$
(4)

with ${\tilde{f}(\rho,\eta):=(1-d)\,f(u(\rho,\eta),v(\rho,\eta))}$. While the
continuity equation, Eq. (3), describes the redistribution of the total
density and its production or degradation, Eq. (4) describes the local
reactions between $u$ and $v$ that adjust their ratio given a prescribed
total-density profile $\rho$. These reactions induce the relaxation of the
mass-redistribution potential $\eta=v+du$ toward the reactive equilibrium
$\tilde{f}(\eta^{*},\rho)=0$.333Consequences of bistability of the reaction
kinetics are discussed in [38].  Here, the nullcline $\eta^{*}(\rho)$ gives
the family of reactive equilibria for different total densities $\rho$. If the
density profile is not uniform, the densities $u$ and $v$ show gradients that
lead to particle diffusion. The (differential) diffusion of $u$ and $v$ is
accounted for by the gradient terms in Eq. (4).

###  III.2 Generalized Cahn–Hilliard equation

The classical Cahn–Hilliard (CH) equation [54] describes the equilibrium
dynamics of an incompressible binary mixture of two types of particles A and B
which undergo phase separation due to a stronger affinity between particles of
the same type than between particles of distinct types. The thermodynamics of
such phase-separating systems are determined by the free-energy functional
[54, 24]

$\mathcal{F}[\phi]=\int\mathrm{d}x\bigg{[}\frac{\kappa}{2}(\nabla\phi)^{2}+g(\phi)\bigg{]}\,,$
(5)

where the scalar field $\phi(x)$ corresponds to the local composition of the
binary mixture and is linearly related to the density of the A particles. The
square-gradient term in $\mathcal{F}[\phi]$ accounts for the free energy costs
of gradients in the composition. We denote the local free-energy density by
$g(\phi)$ to avoid confusion with the term $f(u,v)$ describing chemical
reactions in the 2cRD model. This free-energy density $g(\phi)$ can be
obtained from symmetry arguments [86, 24] or, alternatively, derived, for
example, from a lattice model [87, 88, 89]. Below the critical point at
sufficiently low temperatures, the free-energy density is a double-well
potential that models both entropic effects and the effectively repulsive
interactions between particles of different types. For symmetric binary
mixtures, one often takes the simple quartic form
${g(\phi)=-\frac{r}{2}\phi^{2}+\frac{u}{4}\phi^{4}}$ as obtained from a
Ginzburg-Landau expansion close to the critical point at ${r=0}$ and
${\phi=0}$. This simple form also captures the qualitative shape of the free-
energy density further away from the onset of phase separation.

To describe the dynamics of phase separation, we need an evolution equation
associated with the free-energy functional $\mathcal{F}[\phi]$. Since the
total number of A and the total number of B particles are each conserved, the
local composition $\phi(x,t)$ must obey a continuity equation
${\partial_{t}\phi(x,t)=-\nabla J(x,t)}$, and for the system to relax into
thermal equilibrium, the particle current $J(x,t)$ must be proportional to the
gradient of a corresponding chemical potential $\mu(x,t)$ [90, 91, 24]. These
general concepts of non-equilibrium thermodynamics lead to the following
equations444We consider only the deterministic dynamics here. Effects of noise
during the pattern-formation process (phase-separation process) are discussed
shortly in the discussion Sec. IX.

$\displaystyle\partial_{t}\phi(x,t)$
$\displaystyle=\nabla\big{[}\Gamma(\phi)\nabla\mu\big{]}\,,$ (6a)
$\displaystyle\mu(x,t)$
$\displaystyle=\frac{\delta\mathcal{F}[\phi]}{\delta\phi}=-\kappa\nabla^{2}\phi+\partial_{\phi}g(\phi)\,,$
(6b)

where $\Gamma(\phi)$ denotes the mobility, i.e., an Onsager coefficient that
may depend on the composition $\phi$. These equations correspond to Model B in
the classification scheme of Hohenberg and Halperin [47]. Beyond phase
separation, equations of this form are also used to describe, for example, the
dynamics of fluid thin films [92].

The classical CH equation follows for constant mobility
${\Gamma(\phi)=\Gamma}$ and the quartic free-energy density $g(\phi)$. It
reads

$\displaystyle\partial_{t}\phi$ $\displaystyle=\Gamma\nabla^{2}\mu\,,$ (7a)
$\displaystyle\mu$ $\displaystyle=-\kappa\nabla^{2}\phi+\phi^{3}-\phi\,,$ (7b)

where we eliminated the parameters in the free-energy density by rescaling.
This equation yields phase separation into A- and B-rich domains with
${\phi_{\pm}\approx\pm 1}$ which undergo coarsening and grow until the fully
phase-separated state is reached [52].

As an extension of the equilibrium phase-separation dynamics, reactions
$R(\phi)$ between the particles involved in the phase separation have been
considered [17, 18, 21]. The resulting model couples Model A (reactions) and
Model B (phase separation) dynamics

$\partial_{t}\phi=\nabla\big{[}\Gamma(\phi)\nabla\mu\big{]}+R(\phi)\,,$ (8)

with the chemical potential $\mu$ again given by Eq. (6b). For example, a
linear reaction that allows for conversion between A and B particles yields
${R(\phi)=k(\phi_{0}-\phi)}$ where $k$ is the reaction rate, and $\phi_{0}$
denotes the chemical equilibrium composition balancing the forward and
backward reactions. It was shown that such reactive dynamics may lead to an
interruption of the coarsening process found in pure Model B dynamics and thus
a selection of a finite domain size [17]. This motivates the study of such
models to understand the size control of intracellular condensates which are
thought to compartmentalize biochemical reactions in living cells [18]. This
system also captures models for active matter considering birth and death [20]
and the Oono–Shiwa equation [93]. Overall, Eq. (8) defines a paradigmatic non-
equilibrium field theory if the free energies for the Model A and Model B
dynamics are chosen independently [21]. Allowing for the interaction with
different diffusing species through the reaction term, a related model was
considered to analyze phase separation in cell membranes coupled to protein-
pattern formation [94].

To work out the similarities with the 2cRD system, we rewrite the generalized
CH equation, Eq. (8), in the form

$\displaystyle\partial_{t}\rho$
$\displaystyle=\nabla\big{[}D_{v}(\rho)\nabla\eta\big{]}+\varepsilon
s_{\mathrm{tot}}(\rho,\eta)\,,$ (9a) $\displaystyle 0$
$\displaystyle=-D_{u}\nabla^{2}\rho-\big{[}\eta-\eta^{*}(\rho)\big{]}\,.$ (9b)

To underline the close similarities, we are using the same nomenclature as for
the reaction–diffusion system [Eqs. (3), (4)]: ${\phi\leftrightarrow\rho}$,
${\mu\leftrightarrow\eta}$, ${\Gamma(\phi)\leftrightarrow D_{v}(\rho)}$,
${R(\phi)\leftrightarrow\varepsilon s_{\mathrm{tot}}(\rho,\eta)}$,
${\kappa\leftrightarrow D_{u}}$, and
${\partial_{\phi}g(\phi)\leftrightarrow\eta^{*}(\rho)}$.

Comparing the 2cRD dynamics, Eqs. (3), (4), with the rewritten CH equation
including chemical reactions, Eqs. (9), one notices the following similarities
and differences. First, in the 2cRD system, $\rho$ describes the total
molecule density which may change only due to redistribution of molecules
within the system (current $-D_{v}\nabla\eta$) or production and degradation
via $\varepsilon s_{\mathrm{tot}}$ [cf. Eq. (3)]. Similarly, $\phi$ in the
generalized CH equation is a measure for the local fraction of A particles,
and the local amount of A particles can also only change by redistribution of
A particles [and ensuing redistribution of B particles due to the
incompressibility constraint; driven by the current $-(\Gamma\nabla\mu)$] or
by conversion of particles between A and B type [described by $R(\phi)$].
Therefore, $\rho$ and $\phi$ follow analogous dynamics with the mass-
redistribution potential $\eta$ corresponding to the chemical potential $\mu$
and the total source term $\varepsilon s_{\mathrm{tot}}$ taking the place of
the conversion reactions $R(\phi)$. The timescale of redistribution is set by
the, possibly density-dependent, mobility $\Gamma$ in the extended CH system.
In the 2cRD system, it is set by $D_{v}$, the larger of the two diffusion
coefficients, which is typically density-independent.

Rearranging the defining equation for the chemical potential, Eq. (6b), into
Eq. (9b) suggests that the value of the chemical potential for the CH equation
is fixed by the steady-state balance of a diffusion [$D_{u}\nabla^{2}\rho$]
and a reaction term. One can choose ${\tilde{f}=\eta-\eta^{*}(\rho)}$ as
reaction term in the 2cRD system [Eq. (4); for the classical CH equation we
have ${\eta-\eta^{*}(\rho)=\eta+r\rho-u\rho^{3}}$]. Then, Eq. (9b) corresponds
to the (quasi-)stationary equation for Eq. (4), that is, the equation
following from ${\partial_{t}\eta=0}$, up to the additional diffusion term
${{\sim}\nabla^{2}\eta}$. The continuity equation, Eq. (3), shows that the
strength of the diffusion term agrees with the strength of mass transport.
Consequently, the diffusion term is negligible if mass transport is slow in
comparison to the local reactions. This indicates that a connection between
the reaction–diffusion and CH dynamics can be found under this condition.
Moreover, the identification ${\tilde{f}=\eta-\eta^{*}(\rho)}$ shows that the
steady-state mass-redistribution potential of the homogeneous (well-mixed)
2cMcRD system determined by the reactive equilibria ${\tilde{f}=0}$
corresponds to the chemical potential ${\mu=\partial_{\phi}g=\eta^{*}}$ of a
uniform composition.

Finally, this comparison shows that the rigidity $\kappa$ plays the role of an
effective diffusion constant, a correspondence well known for non-conserved
dynamics where the Allen–Cahn system (Model A) is compared to the Schlögl
reaction–diffusion model [95, 24]. In phase separation, the rigidity $\kappa$
determines the width of the interface between the two separated phases.
Similarly, the diffusion constant $D_{u}$ determines the width of the
transition region between high- and low-density domains in the (nearly) mass-
conserving 2cRD system [38].

###  III.3 Conserved Allen–Cahn equation

If the redistribution of material between different droplets is fast during a
phase-separation process, the dynamics can be limited by the growth of
individual droplets and not by the redistribution of matter between the
different phase-separating domains [50]. Examples are sublimation–deposition
processes where diffusion through the gaseous phase is orders of magnitude
faster than the growth of individual domains [56]. Model B, describing only
the mass-redistribution dynamics, cannot account for this process. Instead,
the growth of individual domains has to be described.

Because one again describes phase separation, the model is based on the same
square-gradient free energy functional Eq. (5) as Model B. However, the time
evolution is different. The growth and shrinking of single domains proceeds to
minimize the free-energy functional $\mathcal{F}[\phi]$. During this process,
the domains are coupled to a global pool which represents the fast-diffusing
phase (the gaseous phase). Thereby, mass conservation is enforced as a global
conservation law. Consequently, the conserved Allen–Cahn (cAC) dynamics is
derived using relaxational Model A dynamics restricted to the mass-conserving
(MC) submanifold in phase space [53, 56]

$\displaystyle\partial_{t}\phi(x,t)$
$\displaystyle=-\left[\Gamma(\phi)\frac{\delta\mathcal{F}[\phi]}{\delta\phi}\right]_{\mathrm{MC}}$
(10)
$\displaystyle=-\left[\Gamma(\phi)\frac{\delta\mathcal{F}[\phi]}{\delta\phi}-\frac{1}{|\Omega|}\int_{\Omega}\mathrm{d}x\,\Gamma(\phi)\frac{\delta\mathcal{F}[\phi]}{\delta\phi}\right]$
$\displaystyle=\Gamma(\phi)\big{(}\kappa\nabla^{2}\phi-\partial_{\phi}g\big{)}+\frac{1}{|\Omega|}\int_{\Omega}\mathrm{d}x\,\Gamma(\phi)\partial_{\phi}g\,,$

again implemented on a system domain $\Omega$ with no-flux boundary conditions
for $\phi$. In the second line, the global mass conservation constraint is
written out explicitly as a Lagrange-multiplier term. Originally, the model
was introduced in Ref. [53] to describe phase separation dominated by viscous
effects. The cAC model also finds application in granular media [96] and the
simulation of incompressible two-fluid flow [97, 98].

We are interested in the cAC system because in the 2cRD system, the cytosolic
density ${v(x,t)=v(t)}$ is spatially constant in the limit of fast cytosolic
diffusion ${D_{v}\to\infty}$. The species $v$ then effectively acts as a
global pool that ensures instantaneous redistribution of molecules over the
whole domain (like the gaseous phase in a sublimation/deposition process). For
reaction–diffusion systems this limit is called the shadow limit [57]. In the
limit ${D_{v}\to\infty}$, the time evolution of the global pool
${v(t)=\eta(t)}$ (because ${d\to 0}$) follows from integration of Eq. (4)
(${\varepsilon=0}$) over the whole domain $\Omega$ with no-flux or periodic
boundary conditions and yields

$\partial_{t}v=\partial_{t}\eta=-\frac{1}{|\Omega|}\int_{\Omega}\mathrm{d}x\,\left(\tilde{f}-\varepsilon
s_{2}\right).$ (11)

Inserting the continuity equation, Eq. (3), into Eq. (4), and identifying
$f_{\mathrm{cAC}}=\tilde{f}-\varepsilon s_{2}$, we find the time evolution of
the density profile as

$\displaystyle\partial_{t}\rho$
$\displaystyle=D_{u}\nabla^{2}\rho+f_{\mathrm{cAC}}(\rho,\eta)-\frac{1}{|\Omega|}\int_{\Omega}\mathrm{d}x\,f_{\mathrm{cAC}}(\rho,\eta)$
$\displaystyle\quad+\varepsilon s_{\mathrm{tot}}(\rho,\eta)\,,$ (12)

which has the standard form of the cAC equation, Eq. (10), with an additional
source term $\varepsilon s_{\mathrm{tot}}$ and the density $\eta$ in the
global pool as an additional ‘parameter’. If the 2cRD system is mass-
conserving, the source term drops out. If, additionally, one sets
${\tilde{f}=\eta-\eta^{*}(\rho)}$ (as for the extended CH dynamics), the
dependence on $\eta$ drops out and the shadow limit of the 2cMcRD system, Eq.
(III.3), agrees exactly with the cAC model, Eq. (10).

##  IV Instability of the homogeneous steady-state

To start the comparison of the 2cRD dynamics with the CH and cAC dynamics, we
begin with the analysis of the homogeneous steady states (HSS) and their
stability properties. The instability of the homogeneous steady state against
spatial modulations describes the onset of pattern formation due to the growth
of small fluctuations around the uniform density distribution.555Outside the
regime of lateral instability, a multistable regime allows for pattern
formation through a finite perturbation of a stable homogeneous steady state
in (nearly) mass-conserving 2cRD systems [38] (nucleation-and-growth regime in
phase-separation dynamics [24]). We begin our analysis with mass-conserving
2cRD systems, compare these with the CH and cAC systems, and then discuss the
implications of broken mass conservation.

The homogeneous steady states $(\rho_{\mathrm{HSS}},\eta_{\mathrm{HSS}})$ of
the 2cMcRD system depend on the average density
${\bar{\rho}=\rho_{\mathrm{HSS}}}$ as a control parameter because the total
mass in the system is not determined by the system dynamics but fixed by the
initial condition. With the nullcline $\eta^{*}(\rho)$ defined as the set of
reactive equilibria ${\tilde{f}(\rho,\eta^{*}(\rho))=0}$, the homogeneous
steady state at average total density $\bar{\rho}$ can be written as
$[\bar{\rho},\eta^{*}(\bar{\rho})]$ [see Eqs. (3), (4) with
${\varepsilon=0}$]. By the identification
${\partial_{\phi}g\leftrightarrow\eta^{*}(\rho)}$ introduced above [cf. Eqs.
(9)], the same expression describes the homogeneous steady states in the mass-
conserving CH and cAC systems [see Eqs. (9) and (III.3)]. We assume that these
homogeneous steady states are linearly stable against homogeneous
perturbations, i.e., we demand that
${\partial_{\eta}\tilde{f}(\rho_{\mathrm{HSS}},\eta_{\mathrm{HSS}})>0}$ [cf.
Eqs. (4), (9), (11)]. Otherwise, already the dynamics of the well-mixed system
would leave the considered, unstable steady state under an infinitesimal
perturbation, and the system would evolve into a different stable steady
state.

Performing a linear stability analysis of Eqs. (3), (4) in terms of Fourier
modes for small perturbations around the homogeneous steady state
$(\rho_{\mathrm{HSS}},\eta_{\mathrm{HSS}})$ yields the dispersion relation
$\sigma_{\mathrm{HSS}}(q)$ that determines the growth rates for each mode with
wavenumber $q$ (see Appendix A). One finds that 2cMcRD models show a lateral
instability if and only if [38]

$\partial_{\rho}\eta^{*}(\rho){\big{|}}_{\rho=\rho_{\mathrm{HSS}}}<0\,.$ (13)

Then the dispersion relation shows a band of unstable modes
${0<q<q_{\mathrm{max}}}$ with growth rates ${\sigma_{\mathrm{HSS}}(q)>0}$ and
a fastest-growing mode $q_{\mathrm{c}}$ (Fig. 4).

The condition Eq. (13) can be heuristically understood as follows: If it is
satisfied, the mass-redistribution potential $\eta$ is lowered in regions of
higher total density $\rho$. As gradients in $\eta$ drive mass redistribution,
the decrease of $\eta$ in regions of high total density $\rho$ leads to
additional mass transport toward these regions, thus amplifying the initial
perturbation. The same condition holds for the (mass-conserving) CH model. For
these thermal-equilibrium systems, the ‘nullcline’ $\eta^{*}(\phi)$ is the
derivative $\partial_{\phi}g(\phi)$ of the free-energy density $g$ such that
the above criterion agrees with the well-known curvature criterion
$\partial_{\phi}^{2}g<0$ which defines the spinodal regime [89].

For the cAC model, the dispersion relation agrees with the one obtained for
the standard Allen–Cahn model for all modes with wavenumbers $q>0$ because the
(linearized) integral terms in Eqs. (11), (III.3) vanish for modes with a
finite wavenumber:

$\sigma_{\mathrm{HSS}}(q)=-D_{u}q^{2}+\partial_{\rho}\tilde{f}(\rho_{\mathrm{HSS}},\eta_{\mathrm{HSS}})\,.$
(14)

Since the diffusion term is negative, a lateral instability
(${\sigma_{\mathrm{HSS}}>0}$) can only be induced by the reaction term, and
the instability condition reads
${\partial_{\rho}\tilde{f}(\rho_{\mathrm{HSS}},\eta_{\mathrm{HSS}})>0}$. With
${\partial_{\rho}\eta^{*}=-\partial_{\rho}\tilde{f}/\partial_{\eta}\tilde{f}}$,
which follows from the definition ${\tilde{f}=0}$ of the nullcline, and the
stability of the homogeneous steady state against homogeneous perturbations
[${\partial_{\eta}\tilde{f}(\rho_{\mathrm{HSS}},\eta_{\mathrm{HSS}})>0}$] we
recover the same instability criterion Eq. (13) also for the cAC system.

Figure 4: Dispersion relations obtained from a linear stability analysis of
the homogeneous steady state (HSS). (a) The dispersion relations
$\sigma_{\mathrm{HSS}}(q)$ as a function of wavenumber $q$ for the mass-
conserving 2cMcRD (${\varepsilon=0}$), CH and cAC models show a band of
linearly unstable modes (${\sigma_{\mathrm{HSS}}>0}$ for
${0<q<q_{\mathrm{max}}}$). The 2cMcRD model has two eigenvalue branches of
which only the unstable branch is shown (see Appendix A). For small
wavenumbers $q$, the dispersion relation for the 2cMcRD dynamics is well
approximated by the dispersion relation of the CH equation, while in the
regime of large wavenumbers, the dispersion relation of the cAC model is found
to be a good approximation. For the 2cMcRD system, the fastest-growing mode
$q_{\mathrm{c}}$ is indicated in the graph. (b) Weak source terms
(${\varepsilon\neq 0}$) stabilize the homogeneous mode
[${\sigma_{\mathrm{HSS}}(q=0)<0}$] in the 2cRD, the CH and the cAC model. As a
result, in the 2cRD and CH models only modes with wave numbers larger than
${q_{\mathrm{min}}\sim\sqrt{\varepsilon/D_{v}}}$ are unstable.

Further analysis shows that not only do the instability criteria agree, but
that the mass-conserving CH and cAC models also approximate the dispersion
relation of the 2cMcRD system for small and large wavenumbers, respectively
(see Fig. 4). We refer to these limits as the _diffusion-_ and _reaction-
limited regimes_ of the reaction–diffusion dynamics (cf. [38, 46]). The two
regimes arise because the 2cMcRD dynamics given by Eqs. (3), (4)
(${\varepsilon=0}$) contains two distinct processes. At low wavenumbers $q$,
the wavelength of the Fourier modes is large and mass is transported over long
distances during the dynamics. Over these large distances, mass transport by
gradients in the mass-redistribution potential—the process described by the CH
model—is the rate-limiting process and determines the growth rate. In
contrast, at large wavenumbers, mass transport proceeds on short distances and
is fast compared to the local reaction dynamics that drives the conversion
between the individual species $u$ and $v$. Thus, for large $q$, the dynamics
can be approximated by the shadow limit ${D_{v}\to\infty}$ where mass
redistribution becomes instantaneous, keeping only the dynamics of the local
reactions between $u$ and $v$. This corresponds to the cAC model. The
mathematical analysis of the dispersion relations is provided in Appendix A
(see also Ref. [38]).

In 2cRD models with finite source terms ($\varepsilon>0$), the homogeneous
steady states (HSS) have to fulfill the two conditions:
${0=\tilde{f}+\varepsilon(s_{2}+ds_{1})}$ and ${0=s_{\mathrm{tot}}}$ where the
second condition replaces the mass constraint of mass-conserving systems. The
first condition entails that the HSS under the influence of weak source terms
$0<\varepsilon\ll 1$ is, at the lowest order in the source strength
$\varepsilon$ (exact for the CH model), given by the HSS of the mass-
conserving system $[\rho_{\mathrm{HSS}},\eta^{*}(\rho_{\mathrm{HSS}})]$. The
second condition arises since in the absence of mass conservation the total
density $\rho_{\mathrm{HSS}}$ is not determined by the initial condition.
Instead, the total density evolves until production and degradation balance at
steady state (${s_{\mathrm{tot}}=0}$).

In the dispersion relation of the 2cRD and CH models, small but finite source
terms (${\varepsilon>0}$) shift the dispersion relation at long wavelength
$2\pi/q$ down such that ${\sigma_{\mathrm{HSS}}(q{=}0)<0}$ and modes with
wavenumbers ${q<q_{\mathrm{min}}}$ are stabilized [see Fig. 4(b) and Appendix
A]. Due to the quadratic form of the dispersion relation at small $q$, one
obtains ${\lambda_{\mathrm{min}}=2\pi/q_{\mathrm{min}}\sim\varepsilon^{-1/2}}$
[see Fig. 4(c)], independent of the specific form of the reaction term
$\tilde{f}$. In contrast, the scaling relationship between $\varepsilon$ and
the wavelength selected by interrupted coarsening of fully developed patterns
depends on the specific reaction term (see Sec. VIII.2). Due to the singular
form of the dispersion relation of the cAC system, weak source terms only
shift the growth rate for homogeneous perturbations (${q=0}$) to negative
values but do not introduce a band of stable modes at small, finite
wavenumbers (see Fig. 4).

In summary, we have shown that the linear instabilities in the homogeneous
steady state of the 2cRD, CH, and cAC models are closely related. We
identified a diffusion- and a reaction-limited regime in which the dynamics of
the 2cRD model is well approximated by corresponding CH or cAC models,
respectively. Later, we will also identify parallels between the instability
of the homogeneous steady state that drives pattern formation and the mass-
competition instability that drives coarsening, and the two different regimes
will re-emerge. In the following section, we first construct the stationary
patterns in the 2cRD system to get an overview of the nonlinear patterns
formed.

##  V Stationary patterns

Our focus in this section is a conceptual understanding of the final
stationary patterns employing a geometric construction in phase space based on
local equilibria theory [38]. These stationary density profiles will play an
important role in the analysis of the coarsening process. There, they serve as
quasi-steady states (QSS) to locally approximate the pattern profile during
the dynamics at asymptotically long times (see Sec. II.1). Such an approach is
commonly used in phase-separating systems where the _shape_ of the interface
profile between the two phases can be approximated as stationary [52, 24, 48].
A more detailed discussion and analysis of the stationary states, which will
later also be used for the mathematical analysis of the mass-competition
process, is summarized in Appendix B for mass-conserving systems and in
Appendix F for systems including weak source terms.

The stationary patterns
$[\rho_{\mathrm{stat}}^{\varepsilon}(x),\eta_{\mathrm{stat}}^{\varepsilon}(x)]$
of the 2cRD system satisfy the equations [stationarity of Eqs. (3), (4)]

$\displaystyle 0$
$\displaystyle=D_{v}\partial_{x}^{2}\eta_{\mathrm{stat}}^{\varepsilon}+\varepsilon
s_{\mathrm{tot}}^{\varepsilon}\,,$ (15a) $\displaystyle 0$
$\displaystyle=D_{u}\partial_{x}^{2}\rho_{\mathrm{stat}}^{\varepsilon}+\tilde{f}(\rho_{\mathrm{stat}}^{\varepsilon},\eta_{\mathrm{stat}}^{\varepsilon})+\varepsilon(s_{1}^{\varepsilon}+ds_{2}^{\varepsilon})\,,$
(15b)

where the superscript $()^{\varepsilon}$ indicates that we consider the system
in the presence of source terms. These stationary equations have different
types of solutions. First, they are satisfied by the homogeneous steady state
discussed in the previous section. Second, they allow for spatially periodic
stationary patterns, which we analyze further below because we are interested
in the stability properties of patterns consisting of many equally sized peaks
or mesas (see Sec. II). Third, asymmetric stationary patterns can be
constructed as well. We do not consider these asymmetric states further
because they are unstable [99, 100].

The periodic stationary patterns can be constructed from _elementary
stationary patterns_ which comprise half a period of the periodic pattern.
Because the system is parity symmetric, the spatial profiles of the left and
right interfaces of peaks or mesas are mirror images of each other (cf. Fig.
2). Thus, the elementary stationary pattern for the periodic pattern with
wavelength $\Lambda$ is the monotonic solution to Eqs. (15) on a domain with
no-flux boundary conditions for $\rho$ and $\eta$ and length $\Lambda/2$ [see
Figs. 5(b,f)].666Next to a stable elementary pattern, an unstable elementary
pattern exists in the multistable regime in which both the homogeneous steady
state and an elementary stationary pattern are stable [38]. We do not consider
those unstable patterns here.

For the generalized CH model, the stationary states are determined by the
stationarity of the field equations, Eqs. (9). Thus, they satisfy Eqs. (15)
with ${\tilde{f}=\eta-\eta^{*}(\rho)}$ and ${s_{1}=-ds_{2}}$. The first
identity shows that the deviation of the chemical potential ${\mu\equiv\eta}$
from its associated fixed point value $\eta^{*}(\rho)$ for a uniform density
profile plays the role of a reaction term. The choice ${s_{1}=-ds_{2}}$ is
required to cancel the source terms in Eq. (15b) that do not appear for the
generalized CH model [Eqs. (9)]. Similarly, the limit ${D_{v}\to\infty}$ of
the stationary equations, Eqs. (15), also determines the stationary states of
the cAC model, that is, the shadow limit of the 2cMcRD system, Eqs. (11),
(III.3). Therefore, we now focus on the 2cRD system only.

The comparison shows that the uniformity of the mass-redistribution potential
in steady state for 2cMcRD systems corresponds to the condition of equal
chemical potential in coexisting phases at thermal equilibrium. The condition
of equal pressure will be obtained in a generalized form below by showing that
the stationary patterns of 2cMcRD systems fulfill a geometric condition
similar to a Maxwell construction. Importantly, for 2cMcRD systems the
geometric construction corresponds to a balance of reactive fluxes instead of
following from a common tangent construction for an underlying free-energy
potential [38, 24].

Figure 5:  Phase-space construction of the elementary stationary mesa and peak
patterns. (a) An $\mathsf{N}$-shaped nullcline $\eta^{*}(\rho)$ (NC; black)
intersects the flux-balance subspace ${\eta={\eta_{\mathrm{stat}}}(M)}$ [FBS;
orange (light gray)] three times, corresponding to the low-density plateau
$\rho_{-}$, the inflection point at density $\rho_{\mathrm{infl}}$, and the
high-density plateau $\rho_{+}$ [(blue) dots]. The representation of the
stationary pattern in phase space [blue (dark gray) line] falls onto the FBS
and (approximately) connects the three FBS-NC inflection points. The position
of the FBS is determined by total turnover balance represented approximately
by a balance of the (red-) shaded areas. (b) In real space, the elementary
stationary pattern for an $\mathsf{N}$-shaped nullcline has the form of a
half-mesa. An interface of width $\ell_{\mathrm{int}}$ connects the upper and
lower plateaus with lengths $L_{\pm}$. The interface position is determined by
the mesa mass $M$ [compare the dark (dark gray) and pale blue (light gray)
profiles]. (c,d) In more detail, on an infinite domain, the stationary profile
$[\rho_{\mathrm{stat}}^{\infty}(x),\eta_{\mathrm{stat}}^{\infty}]$ approaches
the plateau densities $\rho_{\pm}$ exponentially [green (thin) line].
Therefore, on a finite domain, the pattern maximum
$\rho_{\mathrm{stat}}(\Lambda/2)$ and minimum $\rho_{\mathrm{stat}}(0)$
deviate from the plateau densities by $\delta\rho_{\pm}$. This changes the
areas representing total turnover balance, shifting
${\eta_{\mathrm{stat}}}(M)$ relative to $\eta_{\mathrm{stat}}^{\infty}$.
${\eta_{\mathrm{stat}}}(M)$ approaches $\eta_{\mathrm{stat}}^{\infty}$
exponentially as the plateau lengths $L_{\pm}$ increase (pale blue and red
(pale gray) construction). The same shown here for the high-density plateau
applies to the low-density plateau. (e) A $\mathsf{\Lambda}$-shaped NC
$\eta^{*}(\rho)$ (black) intersects the FBS ${\eta={\eta_{\mathrm{stat}}}(M)}$
[orange (light gray)] only twice. The FBS-NC intersection points correspond to
the low-density plateau $\rho_{-}$ and the pattern inflection point
$\rho_{\mathrm{infl}}$ [(blue) dots]. Again, the position of the FBS is
determined by total turnover balance [represented by the (red-) shaded areas].
An increase in the pattern mass increases the pattern amplitude
$\rho_{\mathrm{stat}}(\Lambda/2)$, and the FBS shifts downwards (pale blue and
red (pale gray) construction). (f) Systems with $\mathsf{\Lambda}$-shaped
nullclines sustain peak-shaped patterns [blue (dark gray)]. The ‘half-peak’
width is denoted by $\ell_{\mathrm{int}}$. The peak amplitude grows with the
peak mass $M$ [pale blue (light gray) profile]. (g) The stationary mass-
redistribution potential ${\eta_{\mathrm{stat}}}$ is a decreasing function of
the peak mass $M$.

###  V.1 Mass-conserving case

In the mass-conserving case, the average total density $\bar{\rho}$ of the
final stationary pattern is set by the initial conditions which fix the total
amount of molecules in the system. How does $\bar{\rho}$ enter in the
stationarity equations Eqs. (15) of the mass-conserving system
(${\varepsilon=0}$)?

We denote the stationary patterns of the mass-conserving system by
$[\rho_{\mathrm{stat}}(x),{\eta_{\mathrm{stat}}}]$. Balanced mass
redistribution, i.e., the _continuity equation_ , Eq. (15a), ensures that
${\eta_{\mathrm{stat}}}$ is spatially constant and enters only as a parameter
in the _profile equation_ , Eq. (15b). For each value of
${\eta_{\mathrm{stat}}}$ we can calculate the pattern profile and the average
total density $\bar{\rho}$ using Eq. (15b). This establishes the relation
$\bar{\rho}({\eta_{\mathrm{stat}}})$. Reciprocally, the average density
$\bar{\rho}$ fixes the value of the stationary mass-redistribution potential
${\eta_{\mathrm{stat}}}(\bar{\rho})$.777For the family of stable elementary
stationary patterns, the relation $\bar{\rho}({\eta_{\mathrm{stat}}})$ is
bijective and may therefore be inverted (see supplementary material of Ref.
[46]).

For a fixed average density $\bar{\rho}$ (or, equivalently, fixed
${\eta_{\mathrm{stat}}}$) we are now interested in the stationary pattern
profile $\rho_{\mathrm{stat}}(x)$. For this, we consider the
$(\rho,\eta)$-phase space [see Fig. 5(a,e)]. In phase space, the stationary
pattern $[\rho_{\mathrm{stat}}(x),{\eta_{\mathrm{stat}}}]$ is represented by a
curve parametrized by the spatial position $x$. As noted, the mass-
redistribution potential ${\eta_{\mathrm{stat}}}$ must be constant because all
redistribution processes have to balance. Thus, a curve corresponding to a
stationary pattern $[\rho_{\mathrm{stat}}(x),{\eta_{\mathrm{stat}}}]$ must be
restricted to a horizontal subspace, the _flux-balance subspace_ (FBS)
${\eta=\eta_{\mathrm{stat}}}$ [38]. Now consider the nullcline (NC)
$\eta^{*}(\rho)$ of chemical equilibria ${\tilde{f}=0}$. The corresponding
curve in phase space represents the homogeneous steady states of the system,
and the parts of a pattern which are plateaus, i.e., approximately flat in
real space, must lie on (close to) this nullcline. Thus, FBS-NC intersection
points determine the densities of the pattern plateaus [38]. The profile
equation, Eq. (15b), shows that the inflection point
(${\partial_{x}^{2}\rho|_{\mathrm{infl}}=0}$) also corresponds to a FBS-NC
intersection point. For the plateaus to be stable, the corresponding
homogeneous steady state has to be linearly stable, i.e., in geometrical terms
the nullcline slope needs to be positive (see Sec. IV). In contrast, the
inflection point corresponds to a FBS-NC intersection point with a negative
nullcline slope [38].

Depending on the shape of the nullcline we distinguish two qualitatively
different pattern types [see Fig. 5].888We consider only systems for which
$\rho$ stays bounded from below as is the case for example if $\rho$ describes
a concentration. If the nullcline is $\mathsf{N}$-shaped, three FBS-NC
intersection points exist [see Fig. 5(a)]. The resulting pattern is mesa-
shaped with a low- and a high-density plateau connected by an interface [see
Fig. 5(b)]. The densities $\rho_{+}$ and $\rho_{-}$ of the upper and lower
plateaus are given by the outer FBS-NC intersection points, the inflection
point density $\rho_{\mathrm{infl}}$ by the middle FBS-NC intersection point.
The average density $\bar{\rho}$ sets the position of the interface, i.e., the
lengths $L_{\pm}$ of the upper and lower plateaus [see Fig. 5(b)]. To make the
mesa size explicit, we define the (surplus) mesa mass $M$ with respect to the
background density $\rho_{-}$ as

$M=2\int_{0}^{\frac{\Lambda}{2}}\mathrm{d}x\,\big{(}\rho_{\mathrm{stat}}(x)-\rho_{-}\big{)}\,.$
(16)

To analyze the coarsening process, we will exploit that the mesa mass $M$ can
be used instead of the average density $\bar{\rho}$ to parametrize the
stationary patterns of different total masses.

In contrast, if the nullcline is $\mathsf{\Lambda}$-shaped, only two FBS-NC
intersection points exist which correspond to a low-density plateau $\rho_{-}$
and the inflection point $\rho_{\mathrm{infl}}$ [see Fig. 5(e)]. For this
setting, the pattern does not saturate in a high-density mesa but attains a
peak-shaped profile [see Fig. 5(f)]. We define the peak mass $M$ in the same
way as for mesa patterns. While for mesa patterns only the width depends on
the average density $\bar{\rho}$, it determines both the peak height and width
for peak patterns. Again, the peak mass $M$ can be used instead of the average
density $\bar{\rho}$ to parametrize the different stationary peak profiles
(see supplementary material of Ref. [46]). Finally, we note that one also
observes peak patterns in systems with highly asymmetric $\mathsf{N}$-shaped
nullclines if the mass $M$ is low and the maximal pattern density does not
saturate in the high-density plateau given by the third FBS-NC intersection
point $\rho_{+}$.

In the context of intracellular pattern formation, one observes membrane-bound
protein patterns that exhibit narrow interfaces [101, 102] or strongly
localized high-density clusters (concentration peaks) [85]. Sharp interfaces
can form in the concentration profiles of membrane-associated proteins because
protein diffusion on the membrane is very slow. The diffusion constant for
membrane-associated proteins, corresponding to $D_{u}$ in the 2cRD system, is
about two to three orders of magnitude slower than diffusion in the cytosol
[103, 36]. Specifically, in the 2cMcRD system, the interface or ‘half-peak’
width $\ell_{\mathrm{int}}$ [see Fig. 5(b,f)] is set by the diffusion length
${\ell_{\text{diff}}:=\sqrt{D_{u}/r}}$, where the reaction rate $r$ describes
the reactive timescale at the interface [38] [see also Eq. (63)]. Thus, slow
diffusion allows for steep gradients in the protein density. Motivated by this
characteristic feature of membrane-bound protein patterns, we will consider
peak and mesa patterns in the limit999This is called the singular limit in the
mathematical study of reaction–diffusion systems because gradients at the
interface become infinitely steep.

$\frac{\ell_{\mathrm{int}}}{\Lambda}\ll 1\,.$ (17)

This will allow one to make a _sharp-interface approximation_ and apply
asymptotic theory to quantitatively describe the profiles and dynamics of
fully nonlinear patterns (cf. for example [104, 77, 71]).

Recall that the stationary mass-redistribution potential
${\eta_{\mathrm{stat}}}$, that is, the position of the FBS is fixed by the
average density $\bar{\rho}$. We have also introduced the peak- or mesa-mass
$M$, which can alternatively be used as a parameter of the stable elementary
patterns (see supplementary material of Ref. [46]). The key to understanding
the coarsening process will be how the stationary mass-redistribution
potential ${\eta_{\mathrm{stat}}}$ depends on the domain mass $M$. Therefore,
let us work out the relation ${\eta_{\mathrm{stat}}}(M)$.

Integration of the profile equation, Eq. (15b), yields under consideration of
the no-flux boundary conditions

$0=\int_{0}^{\frac{\Lambda}{2}}\mathrm{d}x\,\tilde{f}(\rho_{\mathrm{stat}},{\eta_{\mathrm{stat}}})\,.$
(18)

In the (biologically relevant) limit ${D_{v}\gg D_{u}}$ this condition for the
value ${\eta_{\mathrm{stat}}}$ is easy to interpret since we have
${\eta\approx v}$: Molecules detach from (attach onto) the membrane and
increase (decrease) the cytosolic density ${v\approx\eta}$ until attachment
and detachment balance, that is, the _total reactive turnover_ vanishes.

To evaluate this balance condition in phase space, we can multiply the profile
equation, Eq. (15b), by $\partial_{x}\rho_{\mathrm{stat}}$ before integrating
and find the modified condition [38]

$0=\int_{\rho_{\mathrm{stat}}(0)}^{\rho_{\mathrm{stat}}(\frac{\Lambda}{2})}\mathrm{d}\rho\,\tilde{f}(\rho,{\eta_{\mathrm{stat}}})\,.$
(19)

This total turnover balance is qualitatively represented by an (approximate)
balance of the areas enclosed between NC and FBS [(red-) shaded areas in Fig.
5(a,e)], similarly to a Maxwell construction (see Ref. [38] and Appendix B).
Total turnover balance agrees exactly with the area equality, i.e., the
Maxwell construction, if the reaction term has the simple form
${\tilde{f}\sim\eta-\eta^{*}(\rho)}$ as in the CH system.

Total turnover balance Eq. (19) determines the relation
${\eta_{\mathrm{stat}}}(M)$: For mesa patterns in the sharp-interface limit, a
change of the mesa mass $M$ changes the lengths $L_{\pm}$ of the upper and
lower plateaus [see Fig. 5(b)]. Now, linearization of the profile equation,
Eq. (15b), around the plateau densities $\rho_{\pm}$ shows that the pattern
profile approaches these plateau densities exponentially [see Fig. 5(d)].
Since ${\eta_{\mathrm{stat}}}$ depends on the maximal and minimal pattern
densities $\rho_{\mathrm{stat}}(\Lambda/2)$ and $\rho_{\mathrm{stat}}(0)$,
${\eta_{\mathrm{stat}}}$ changes exponentially with increasing length of the
upper or lower plateau [see Fig. 5(c); detailed analysis in Appendix B and
supplementary material of Ref. [46]]:

$\partial_{L_{\pm}}{\eta_{\mathrm{stat}}}\propto\mp\exp\left(-\frac{2L_{\pm}}{\ell_{\pm}}\right),$
(20)

where $\ell_{\pm}$ describes the length scale of the exponential approach of
the upper or lower plateau density. Because a change ${\delta
M\approx(\rho_{+}-\rho_{-})\delta L}$ in the mesa mass is proportional to a
length change $\delta L$ of the plateaus, ${\eta_{\mathrm{stat}}}$ also
changes only exponentially slowly with the mesa mass $M$.

For peak patterns the peak amplitude increases with the peak mass $M$ [see
Fig. 5(f)]. Thus, the stationary mass-redistribution potential
${\eta_{\mathrm{stat}}}$ has to decrease to balance the total turnover [see
Fig. 5(e,g) and Ref. [46]], and we obtain

$\eta_{\mathrm{stat}}(M)\sim M^{-\alpha}\,,$ (21)

if the reaction term $\tilde{f}$ [and thus the nullcline $\eta^{*}(\rho)$] is
of power-law form at large densities $\rho$ (see the scaling analysis in
Appendix B.3.1).

###  V.2 The effect of weak source terms

The profile of the stationary pattern
$[\rho_{\mathrm{stat}}^{\varepsilon}(x),\eta_{\mathrm{stat}}^{\varepsilon}(x)]$
under the influence of weak source terms (${0<\varepsilon\ll 1}$) shows only
small deviations (of order ${\sim\varepsilon}$) from the stationary profile
$[\rho_{\mathrm{stat}}(x),{\eta_{\mathrm{stat}}}]$ of the mass-conserving
system (${\varepsilon=0}$).

Importantly, as for the homogeneous steady state, the average total density
$\bar{\rho}$ is no longer a control parameter of the system with source terms
as it is not set by the initial condition. Instead, the average density
evolves as [integrate Eq. (3)]

$\partial_{t}\bar{\rho}=\varepsilon\frac{2}{\Lambda}\int_{0}^{\frac{\Lambda}{2}}\mathrm{d}x\,s_{\mathrm{tot}}(\rho(x),\eta(x))\,.$
(22)

Thus, the average density $\bar{\rho}$, that is, the ‘size’ of the stationary
peak or mesa adapts until the overall production and degradation of molecules
in the system are balanced. The mass of the stationary pattern is thus fixed
by the source balance condition ${\partial_{t}\bar{\rho}=0}$.

In summary, we have shown that (nearly) mass-conserving 2cRD systems form peak
or mesa patterns. The type of pattern is determined by the shape of the
nullcline $\eta^{*}(\rho)$. In the mass-conserving system, the stationary
mass-redistribution potential is constant and its value depends on the domain
mass $M$. It decreases exponentially slowly with mass for mesa patterns, and
(typically) like a power law for peak patterns. In both cases, we have
${\partial_{M}{\eta_{\mathrm{stat}}}<0}$. Weak source terms fix the peak or
mesa mass because the stationary pattern size has to balance the overall
production and degradation in the system.

##  VI Growth rate of the mass-competition instability

Having classified the possible patterns, we now study the growth rate $\sigma$
of the mass-competition instability in 2cMcRD systems in terms of the
interaction of two neighboring elementary patterns, each approximated by a
quasi-steady-state density profile. This instability underlies the coarsening
dynamics (see Sec. II). Detailed, systematic derivations of the results
presented in the following section can be found in Appendix E.

###  VI.1 Assumptions and approximations

To determine the growth rates, we build on two assumptions. First, we employ
the sharp-interface approximation (${\ell_{\mathrm{int}}\ll\Lambda}$) which is
justified for well-separated peaks and interfaces (see section V). It allows
the use of singular perturbation theory to determine growth rates of
perturbations around stationary patterns in the fully nonlinear regime.
Second, we assume that mass competition is slow compared to the relaxation of
single peaks and mesas toward their stationary profile. With this separation
of timescales, we account for the observation that—during the coarsening
process—the pattern profiles of single peaks or mesas are well approximated by
the stationary profiles (fast ‘regional’ relaxation) [see Fig. 1(b,d)]. Thus,
we use a QSS approximation for the density profile of single peaks and
interfaces but not for the value of $\eta$ in the peak or interface region (as
used in Ref. [46]). As a result of the QSS approximation, mass competition is
described in terms of a few _collective variables_ : peak masses and positions
or interface positions.

Quantitatively the timescale separation between local relaxation and mass
redistribution can be analyzed using singular perturbation theory.
Specifically, in Appendix D, we derive the relaxation rates of the relaxation
modes of the elementary pattern that redistribute mass from the peak or
interface region to the plateau regions.101010Other relaxation modes exist
that locally deform the peak or interface profile. These modes are not
accessible via singular perturbation methods where one assumes the peak or
interface as infinitely sharp. Within the sharp-interface approximation,
however, we expect these modes to relax rapidly because mass redistribution
within the narrow interface or peak regions will be fast. Also, we find
numerically that those modes are relaxed during the coarsening process [see
Fig. 1(b,d)]. These growth rates describe how fast a single peak or mesa
relaxes toward its steady state profile if mass is added (or removed) from its
plateaus. If this plateau relaxation is fast compared to mass redistribution
between different peaks or mesa interfaces (cf. Fig. 2), we can assume that
both the pattern plateaus and the peak or mesa mass $M$ are fully relaxed
during the competition process (timescale separation, see also Appendix E).
The detailed analysis shows that local relaxation is fast compared to the
mass-competition process if the following condition is satisfied

$\frac{1}{\Lambda\,\partial_{{\eta_{\mathrm{stat}}}}\rho_{\pm}}\gg\left|\partial_{M}{\eta_{\mathrm{stat}}}\right|\,,$
(23)

where $\rho_{\pm}$ denotes the densities in the high- and low-density plateau.
Heuristically this condition can be understood as a consequence of the fact
that the strength of changes in the mass-redistribution potential $\eta$
determines how fast mass redistribution proceeds. Consider adding a small
amount of mass $\delta M$ in a plateau [l.h.s. of Eq. (23)] or a peak/mesa
[r.h.s. of Eq. (23)]. Within the sharp-interface approximation, additional
mass in the plateaus (of length ${\approx\Lambda}$) leads to the change
${\delta{\eta_{\mathrm{stat}}}\approx\partial_{\rho_{\pm}}\eta^{*}\,\delta
M/\Lambda}$ of the (stationary) mass-redistribution potential. Using the
implicit function theorem this can be written as
${\delta{\eta_{\mathrm{stat}}}\approx\delta
M/(\Lambda\partial_{\eta_{\mathrm{stat}}}\rho_{\pm})}$. Similarly, additional
mass in a domain induces a change
${\delta{\eta_{\mathrm{stat}}}\approx(\partial_{M}{\eta_{\mathrm{stat}}})\,\delta
M}$. Thus, the condition Eq. (23) ensures that the $\eta$ gradients induced by
redistribution of mass into the pattern plateaus are large compared to the
gradients arising during mass competition between neighboring domains.

Importantly, the properties of stationary patterns in 2cMcRD systems ensure
that the condition Eq. (23) is generically fulfilled for sufficiently large
peaks and mesas because $\partial_{M}{\eta_{\mathrm{stat}}}$ is strongly
suppressed at large peak and mesa masses
${M=\Lambda\left(\bar{\rho}-\rho_{-}\right)}$ [see Eqs. (20), (21)].

Conveniently, the condition Eq. (23) not only allows us to neglect the time
evolution of the pattern plateaus and assume these as fully relaxed to a
(quasi-)steady state, but we can also neglect the total mass stored in the
background plateau $\rho_{-}$ during the mass-competition process.
Specifically, the total mass of an elementary pattern changes with the peak or
mesa mass $M$, that is, with the stationary mass-redistribution potential
${\eta_{\mathrm{stat}}}$ by [see Eq. (16)]

$\partial_{\eta_{\mathrm{stat}}}\int_{0}^{\frac{\Lambda}{2}}\\!\mathrm{d}x\,\rho_{\mathrm{stat}}=\partial_{\eta_{\mathrm{stat}}}\big{(}\Lambda\rho_{-}\,{+}\,M\big{)}\approx\partial_{\eta_{\mathrm{stat}}}M\,,$
(24)

where the last estimate follows by inverting the above condition Eq. (23).
Therefore, redistribution of mass into the background plateau $\rho_{-}$ does
not play a significant role in the mass-competition dynamics if condition Eq.
(23) is fulfilled. The same simplification arises in LSW theory for coarsening
in phase-separating systems [49].

###  VI.2 Instability growth rates for competition between two peaks/mesas

Figure 6:  The mass-competition instability of two half peaks. (a) The
stationary pattern $[\rho_{\mathrm{stat}}(x),{\eta_{\mathrm{stat}}}]$ to
analyze peak competition is symmetric and composed of two equally sized half
peaks of mass $M/2$ each [cf. Fig. 2(a) bottom]. (b) The antisymmetric mass-
competition eigenmode [blue (dark gray), orange (light gray)] destabilizes
this symmetric stationary state. The density profile $\delta\rho$ [blue (dark
gray)] is localized to the peak regions. In these regions, it is well
approximated by the ‘mass mode’ $\partial_{M}\rho_{\mathrm{stat}}(x)$ (black
dashed). In the sharp interface limit, the mass-redistribution potential
$\delta\eta$ [orange (light gray), strongly magnified] is approximately linear
(red dashed). The linear stability analysis is exemplified in this figure for
the peak-forming model ${\tilde{f}=\eta-10\rho/(1+\rho^{2})}$ with parameters
${D_{u}=1}$, ${D_{v}=10}$, ${\bar{\rho}=4}$, and ${\Lambda=200}$ (see Appendix
H.3).

Applying the approximations discussed in the previous section, we now explain
the processes that underlie the mass-competition instability. We find that the
mass-competition process is composed of three substeps: Particle release at
one peak or mesa, diffusive mass transport through the cytosol, and particle
incorporation into the other domain. The timescale of mass competition is
shown to be the sum of the timescale of reactive conversion (release and
incorporation) and mass transport. Depending on which of the two contributions
is rate-limiting, we find, alike to the instability of the HSS (see Sec. IV),
a diffusion- and a reaction-limited regime which correspond to CH and cAC
dynamics, respectively.

The mass-competition process. — For all coarsening scenarios (see Fig. 2),
mass competition follows the same principles. In essence, a gradient
$\delta\eta(x)$ induced in the mass-redistribution potential between the two
competing domains drives the mass-competition instability: A mass increase in
one peak/mesa reduces the mass-redistribution potential there [cf. Eqs. (20),
(21)], and the resulting gradient leads to mass redistribution that enhances
the initial perturbation [46]. For concreteness, we discuss the competition
between peaks (see Fig. 6). The mathematical analysis of the eigenmode
corresponding to peak competition [see Fig. 6(b)] is given in Appendix E.

Figure 7: The three steps of mass competition between two peaks (similarly for
mesas) in the limit ${D_{v}\gg D_{u}}$. In this limit, particles accumulate
into sharp peaks on the membrane while the cytosolic concentration is nearly
uniform at each peak due to the much faster cytosolic diffusion. Therefore,
most of the peak mass $M$ accumulates in the membrane species $u$, and one has
${\rho\approx u}$. Moreover, it holds ${\eta=v+du\approx v}$. During mass
competition, chemical reactions convert particles from the slow into the fast
species at the smaller peak, where the half-peak mass is reduced by $-\delta
M/2$ [particle release; red (dark gray) arrow]. The resulting change in the
cytosolic density yields a gradient that transports these particles toward the
other peak [orange (light gray) arrow] where they are reactively re-converted
to the slow species and increase the peak mass of the larger peak [particle
incorporation; red (bold, dark gray) arrow]. The rate of release and
incorporation is given by $\sigma_{\mathrm{R}}$ while the rate of mass
transport between the peaks is given by $\sigma_{\mathrm{D}}$.

To understand the competition process in detail, we focus on the (biologically
relevant) limit ${D_{u}\ll D_{v}}$. In this limiting case, the pattern mainly
forms in the slow-diffusing (membrane-bound) species $u$ while, in comparison,
the density profile in the species $v$ is smoothed out by fast diffusion. At
the same time, redistribution of mass between the peaks mainly proceeds
through the fast-diffusing species $v$; note that
${\partial_{x}\eta\approx\partial_{x}v}$ in the plateau between the peaks. As
a consequence of these different roles of $u$ and $v$ in the dynamics,
redistribution of mass from smaller into larger peaks, i.e., the mass
competition process driving coarsening, must involve three subsequent steps
illustrated in Fig. 7: First, particles at the smaller peak detach from the
membrane into the cytosol (“particle release”, left part of Fig. 7). Second,
these additional cytosolic particles create a gradient ${\delta
v(x)\approx\delta\eta(x)}$ in the fast-diffusing species which causes net
diffusion of particles toward the other, larger peak (top middle of Fig. 7).
In a third step, chemical reactions revert the particles back into the slow
state $u$, thus letting the larger peak grow (“particle incorporation”, the
opposite of the release process at the first peak; right part of Fig. 7). In
the next paragraphs, we explain why particles are released from the smaller
and incorporated into the larger peak. This leads to the criterion for
coarsening given in Eq. (4) in Ref. [46].

Consider the weakly perturbed initial configuration where, starting from a
symmetric configuration of equal masses, a small amount of mass $\delta M/2$
is redistributed from the left to the right half peak. This creates a mass
imbalance $\delta M$ between the two half peaks (see Fig. 7).111111Note that
the peak mass $M$ as defined in Eq. (16) decreases/increases by $\pm\delta M$
at the left/right peak because we analyze only half peaks. This mass imbalance
implies that reactive turnover balance, Eq. (18), can no longer be fulfilled
for the same cytosolic densities at the two peaks. At the smaller peak, the
stationary cytosolic concentration ${v\approx{\eta_{\mathrm{stat}}}(M-\delta
M)}$ increases compared to the uniform concentration
${v\approx{\eta_{\mathrm{stat}}}(M)}$ of the symmetric steady state
[${\partial_{M}{\eta_{\mathrm{stat}}}<0}$, see Eqs. (20), (21)]. Thus, the
local relaxation of the smaller peak toward the shifted stationary state
corresponding to its reduced peak mass drives the detachment of particles from
the membrane to increase the local cytosolic density toward
${v\approx{\eta_{\mathrm{stat}}}(M-\delta M)>{\eta_{\mathrm{stat}}}(M)}$. The
same but reverse process occurs at the larger peak: Particles attach onto the
membrane and deplete the cytosolic density toward
${v\approx{\eta_{\mathrm{stat}}}(M+\delta M)<{\eta_{\mathrm{stat}}}(M)}$. The
resulting cytosolic gradient couples both processes by diffusive transport of
the released particles from the smaller toward the larger peak, inducing the
positive feedback that drives the mass-competition instability. Taken
together, the characteristic timescale $1/\sigma$ of the mass-competition
instability is given by the sum of the timescales of reactive conversion
$1/\sigma_{\mathrm{R}}$ at the peaks (particle release/incorporation) and
diffusive redistribution $1/\sigma_{\mathrm{D}}$ between the peaks. Next, we
determine the rates $\sigma_{\mathrm{R}}$ and $\sigma_{\mathrm{D}}$.

In the limit ${D_{v}\to\infty}$, diffusive transport in the cytosol becomes
instantaneous and the only dynamic processes are the release and incorporation
processes (${1/\sigma_{\mathrm{D}}\to 0}$). The dynamics of release and
incorporation follows from integration of Eq. (1a) (${\varepsilon=0}$) over,
say, the left half peak $[-\Lambda/2,0]$. One can use that the gradients of
the $u$ profile in the plateau are small and that the peak mainly forms on the
membrane such that $-\partial_{t}\delta
M/2\approx\partial_{t}\int_{-\Lambda/2}^{0}\mathrm{d}x\,u$. The factor 1/2
accounts for the fact that we consider only a half-peak. It follows that [cf.
reactive turnover balance, Eq. (18), and Fig. 7]

$\partial_{t}\frac{-\delta
M}{2}\approx\int_{-\frac{\Lambda}{2}}^{0}\mathrm{d}x\,\tilde{f}(\rho(x),{\eta_{\mathrm{stat}}}(M))\,.$
(25)

The value of the mass-redistribution potential stays constant at its
unperturbed value ${\eta_{\mathrm{stat}}}(M)$ of the symmetric steady state
because the cytosolic density is uniform in the limit of instantaneous
redistribution and any inflow of particles at the smaller peak must be
balanced by an outflow of particles at the larger peak (antisymmetry of the
eigenmode associated with mass competition, see Fig. 6). The reactive turnover
integral (right-hand side) can be evaluated approximately using that the
plateau regions of the pattern rapidly relax toward reactive equilibrium
${\tilde{f}\approx 0}$ (see Secs. V, VI.1). Thus, the integral is dominated by
contributions from the peak, and one can approximate the integral by the
product of half the peak width $\ell_{\mathrm{int}}$ and the (linearized)
average net reactive flux ${\bar{f}_{\mathrm{peak}}(\delta
M)=\langle\tilde{f}_{\eta}\cdot[{\eta_{\mathrm{stat}}}(M)-{\eta_{\mathrm{stat}}}(M{-}\delta
M)]\rangle_{\mathrm{int}}}$, introducing the shorthand notation
$\tilde{f}_{\eta}=\partial_{\eta}\tilde{f}$. Taken together, this yields

$\displaystyle\partial_{t}\frac{-\delta M}{2}$
$\displaystyle\approx\ell_{\mathrm{int}}\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}\Big{[}{\eta_{\mathrm{stat}}}(M)-{\eta_{\mathrm{stat}}}(M{-}\delta
M)\Big{]}$
$\displaystyle=\ell_{\mathrm{int}}\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}\,\partial_{M}{\eta_{\mathrm{stat}}}\,\delta
M\,,$ (26)

where the second line applies for the linear regime where $\delta M$ is small.
This gives the following estimate for the rate $\sigma_{\mathrm{R}}$ of the
release and incorporation processes

$\sigma_{\mathrm{R}}\approx-2\ell_{\mathrm{int}}\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}\partial_{M}{\eta_{\mathrm{stat}}}\,.$
(27)

The mathematical analysis given in Appendix E yields explicit expressions for
the ‘half-peak’ width or interface width $\ell_{\mathrm{int}}$ as well as the
average conversion rate $\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}$ [see
Eqs. (94), (95)]. The resulting expression for $\sigma_{\mathrm{R}}$ is exact
in the sense of a singular-perturbation analysis.

If the diffusive transport between the peaks is not instantaneous, its rate
$\sigma_{\mathrm{D}}$ is determined by the strength of the gradient
$\delta\eta(x)$ resulting from the shifts in the cytosolic density at the
peaks. In the limit of fast reactive conversion ${1/\sigma_{\mathrm{R}}\to
0}$, the only dynamic process is the diffusive mass transport between the
peaks, and the cytosolic density at the peaks directly relaxes toward its
equilibrium values ${{\eta_{\mathrm{stat}}}(M\pm\delta M)}$ (QSS
approximation). Thus, the mass transport between the peaks is approximately
given by [integrate the continuity equation, Eq. (3) (${\varepsilon=0}$), over
the domain half $[-\Lambda/2,0]$]

$\displaystyle\partial_{t}\frac{-\delta M}{2}$ $\displaystyle\approx
D_{v}\frac{{\eta_{\mathrm{stat}}}(M{+}\delta
M)-{\eta_{\mathrm{stat}}}(M{-}\delta M)}{\Lambda}$
$\displaystyle\approx\frac{2D_{v}}{\Lambda}(\partial_{M}{\eta_{\mathrm{stat}}})\delta
M\,,$ (28)

where the second approximation again holds in the linear regime for small mass
differences $\delta M$. The rate $\sigma_{\mathrm{D}}$ of diffusive mass
transport between the peaks can then be read off as

$\sigma_{\mathrm{D}}\approx-\frac{4D_{v}}{\Lambda}\,\partial_{M}{\eta_{\mathrm{stat}}}\,.$
(29)

Taken together, the above analysis explains that mass competition proceeds in
three subsequent steps. Its characteristic timescale is the sum of the
characteristic timescales of the “reactive-conversion steps”,
$1/\sigma_{\mathrm{R}}$, (particle release and incorporation) and diffusive
mass transport, $1/\sigma_{\mathrm{D}}$, between the competing domains.

The rates of the coarsening modes. — Indeed, the singular-perturbation
analysis given in Appendix E shows that the growth rate $\sigma$ of the mass-
competition instability in peak and mesa patterns has the structure

$\frac{1}{\sigma}=\frac{1}{\sigma_{\mathrm{D}}}+\frac{1}{\sigma_{\mathrm{R}}}\,.$
(30)

Please note that the mathematical derivation of this relationship and the
following expressions for each individual rate does not use the assumption
${D_{u}\ll D_{v}}$ we have used above to simplify the interpretation (see
Appendix E).

For peak competition (cf. Fig. 2), the singular perturbation analysis yields
the rates expected from the above heuristic arguments (see Appendix E)

$\displaystyle\sigma_{\mathrm{D}}$
$\displaystyle=-\frac{4D_{v}}{\Lambda}\,\partial_{M}{\eta_{\mathrm{stat}}}\,,$
(31a) $\displaystyle\sigma_{\mathrm{R}}$
$\displaystyle=-\frac{2\ell_{\mathrm{int}}\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}}{1+d}\,\partial_{M}{\eta_{\mathrm{stat}}}\,.$
(31b)

For the two scenarios of mesa competition and coalescence, we obtain (see
Appendix E)

$\displaystyle\sigma_{\mathrm{D}}^{\pm}$
$\displaystyle=-\frac{4D_{v}}{\xi_{\mp}\Lambda}\,\partial_{M}^{\pm}{\eta_{\mathrm{stat}}}\,,$
(32a) $\displaystyle\sigma_{\mathrm{R}}^{\pm}$
$\displaystyle=-\frac{2\ell_{\mathrm{int}}\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}}{1+d}\,\,\partial_{M}^{\pm}{\eta_{\mathrm{stat}}}\,,$
(32b)

where ${\xi_{\pm}=2L_{\pm}/\Lambda}$ are the relative plateau widths. The
superscript $+$ indicates the mode driven by competition of the high-density
plateaus [mesa competition; Fig. 2(a)], while $-$ indicates the mode driven by
competition of the low-density plateaus [mesa coalescence; Fig. 2(b)]. We
introduced the notation
${\partial_{M}^{\pm}{\eta_{\mathrm{stat}}}=\pm\frac{1}{2(\rho_{+}-\rho_{-})}\partial_{L_{\pm}}{\eta_{\mathrm{stat}}}}$
which describes the change $\delta{\eta_{\mathrm{stat}}}$ of the stationary
mass-redistribution potential only due to the change in length ${\delta
L_{\pm}=\pm\delta M/(\rho_{+}-\rho_{-})}$ of the high-density (low-density)
plateau, respectively. The rates for mesa competition and coalescence depend
on the distinct, modified derivatives
$\partial_{M}^{\pm}{\eta_{\mathrm{stat}}}$ since for the competition
(coalescence) mode only the high-density (low-density) plateaus change in
length, respectively (see Fig. 2). The other plateau in each scenario merely
shifts as a whole and does not change in length during the mass-competition
process. Because peak coalescence also depends on the length change of the
lower plateaus, it is analogous to the case of mesa coalescence and is
described by the same rate $\sigma^{-}_{\mathrm{R,D}}$.121212For peak
coalescence, the rate expressions $\sigma_{\mathrm{D,R}}^{-}$ only hold
qualitatively because the gradient $\delta\eta(x)$ between the competing (low-
density) domains lies within the peak and not within a plateau. Therefore, it
cannot simply be approximated by a linear gradient within the sharp-interface
limit (cf. Appendix E and Fig. 16).

The condition for uninterrupted coarsening. — In the above heuristic
discussion of the mass-competition process, we have argued that the mass-
competition instability is driven by the _decrease_ of the stationary mass-
redistribution potential ${\eta_{\mathrm{stat}}}(M)$ as a function of the
domain mass $M$. Indeed, both rates $\sigma_{\mathrm{D}}$ and
$\sigma_{\mathrm{R}}$ are proportional to the derivative
$\partial_{M}{\eta_{\mathrm{stat}}}$. Thus, the singular perturbation analysis
recovers the condition that uninterrupted coarsening occurs if
${\partial_{M}{\eta_{\mathrm{stat}}}<0}$ holds for all stationary peaks that
are stable as isolated elementary patterns. Similarly, mesa patterns undergo
uninterrupted coarsening if ${\partial_{M}^{+}{\eta_{\mathrm{stat}}}<0}$ or
${\partial_{M}^{-}{\eta_{\mathrm{stat}}}<0}$ if fulfilled for all stable
stationary mesas. Then, these stable stationary domains are all destabilized
by mass exchange with neighboring domains: The growth rate $\sigma$ is always
positive and any mass difference $\delta M$ between two domains grows
exponentially by $\partial_{t}\delta M=\sigma\delta M$, i.e., one has $\delta
M\sim\operatorname{e}^{\sigma t}$. This stability condition was previously
derived mathematically [30, 45]. In Ref. [46] it is shown that the criterion
${\partial_{M}{\eta_{\mathrm{stat}}}<0}$ (and
${\partial_{M}^{\pm}{\eta_{\mathrm{stat}}}<0}$) is generically fulfilled in
strictly mass-conserving two-component reaction–diffusion systems, implying
that these systems always exhibit uninterrupted coarsening [see also Eqs.
(20), (21) Sec. V]. In mass-conserving systems with more than two components,
it is possible for $\partial_{M}{\eta_{\mathrm{stat}}}$ to become positive,
resulting in interrupted coarsening. For example, interrupted coarsening is
observed for specific three-component mass-conserving systems where the total
density of the three components is conserved [6, 105, 85, 80]. For two-
component systems, the coarsening processes can be brought to a halt by
introducing weak production and degradation terms that break mass
conservation. This latter scenario may serve as a prototype for the analysis
of interrupted coarsening in more-component models and will be presented in
Sec. VIII below.

Diffusion- and reaction-limited coarsening. — The growth rate of the mass-
competition instability exhibits two distinct limits, depending on whether the
diffusive mass transport or the reactive conversion step is rate-limiting. In
the first case we find ${\sigma\to\sigma_{\mathrm{D}}}$ while in the second
case we have ${\sigma\to\sigma_{\mathrm{R}}}$. Therefore, we term
$\sigma_{\mathrm{D}}$ the _diffusion-limited_ and $\sigma_{\mathrm{R}}$ the
_reaction-limited_ growth rate (cf. Ref. [46]). We have actually applied these
two limiting cases in the above heuristic discussion to determine the two
growth rates $\sigma_{\mathrm{D,R}}$. From the explicit expression for the
growth rates $\sigma$, $\sigma_{\mathrm{D}}$, and $\sigma_{\mathrm{R}}$, Eqs.
(30)–(32), one infers that the crossover between these regimes is located at

$\frac{D_{v}}{\Lambda}\sim\ell_{\mathrm{int}}\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}\,,$
(33)

i.e., where the timescale of diffusive redistribution between the peaks/mesas,
$1/\sigma_{\mathrm{D}}$, is comparable to the timescale of reactive conversion
at the peak or mesa interface, $1/\sigma_{\mathrm{R}}$.

####  VI.2.1 The diffusion-limited regime

In the diffusion-limited regime $\sigma\to\sigma_{\mathrm{D}}$ one recovers
the result Eq. (3) from Ref. [46].131313The expressions differ by a factor of
$2$ because here we consider the competition of two half peaks or half mesas.
There, the mass-redistribution potential $\eta$ at each peak (or mesa) was
approximated to be in a quasi-steady state (QSS), i.e.,
${\eta|_{\mathrm{peaks}}\approx{\eta_{\mathrm{stat}}}(M\mp\delta M)}$. This
assumption holds in the diffusion-limited regime considered here because the
reactive conversion at the peaks and the ensuing relaxation of $\eta$ toward
its QSS are fast. In this regime, the growth rate of the mass-competition
instability is determined solely by the amplitude of the gradients in $\eta$
that are induced by a mass difference between the peaks/mesas. During
coarsening, this regime will always be reached at late times because the peak
or interface separation increases further and further, and mass redistribution
becomes slower over larger distances. Hence, the asymptotic long-time behavior
of the coarsening process only depends on $\sigma_{\mathrm{D}}$.

The diffusion-limited regime manifests itself in the generalized CH equation,
in which the mass-redistribution potential $\eta$ does not have a time
evolution of its own, but instantly adapts to the density profile $\rho$ (see
Eq. (9b)). Consequently, the only dynamic process is mass redistribution, and
$\sigma_{\mathrm{D}}$ gives the total growth rate of the mass competition
instability in the generalized CH equation. The derivation of the growth rate
$\sigma=\sigma_{\mathrm{D}}$ of the mass-competition instability in the
generalized CH equation proceeds analogously to Appendix E.

####  VI.2.2 The reaction-limited regime

If, in contrast, diffusive mass redistribution between domains is fast
compared to the local reactive conversion, we have
${\sigma\to\sigma_{\mathrm{R}}}$. Due to fast mass redistribution through the
cytosol, the species $v$ acts as a global pool for mass exchange and the
coarsening rate is set by the release of mass into the pool by the shrinking
peak/mesa and the intake of mass from the pool by the growing peak/mesa.

The reaction-limited regime is reached in the shadow limit ${D_{v}\to\infty}$
such that $\sigma_{\mathrm{R}}$ is the growth rate expression found in cAC
systems. In contrast to CH equations, which describe the dynamics based solely
on mass redistribution, cAC systems describe only the reactive dynamics of
peak growth and shrinking (“interface-controlled kinetics”).

##  VII Scaling analysis of the coarsening dynamics

Having calculated the growth rate $\sigma$ of the mass-competition
instability, we now apply a scaling argument [51, 106, 46] to argue how the
coarsening law, i.e., the time evolution of the average wavelength
$\langle\Lambda\rangle(t)$ of the pattern on a large domain, can be obtained
from $\sigma$. For specificity, we again focus on peak patterns [cf. Fig.
1(a)].

An extended pattern consisting of many peaks can, in principle, show both peak
competition and peak coalescence. Also, mass can be exchanged not only between
nearest neighbor peaks but also between peaks that are further apart. Out of
these coarsening scenarios the fastest competition (or coalescence) process
will lead to the fastest collapse of domains and thus drive the coarsening
process. For peak patterns, we can focus on the mass-competition instability
between two neighboring peaks.141414Mass competition is faster over shorter
distances. Moreover, in the shadow limit (cAC system) mass exchange is
instantaneously mediated via a “global pool”. Thus, the competition rate is
(approximately) independent of the spatial arrangement of the peaks. Refer to
Appendix E for a detailed discussion. The coalescence mode can be neglected
because its rate ${\sigma\sim\partial_{M}^{-}{\eta_{\mathrm{stat}}}}$ is
exponentially suppressed with the plateau length. The change
$\partial_{M}{\eta_{\mathrm{stat}}}$ of the stationary mass-redistribution
potential induced by the height change of the peaks—which induces mass
transport in the peak-competition scenario—falls off much more slowly [cf. Eq.
(21), in agreement with the observed phenomenology [cf. Fig. 1(a)].

To find the time evolution of the average peak separation
$\langle\Lambda\rangle(t)$, we need to determine how fast the total number of
peaks $N(t)$ changes. For this, we first determine the characteristic collapse
timescale $t_{\mathrm{col}}$ for the smaller of two neighboring peaks with
masses $M_{1}$ and $M_{2}$. Let us denote the average mass by
${\bar{M}=(M_{1}+M_{2})/2}$ and the mass difference by ${2\,\delta
M=M_{1}-M_{2}}$ (assuming ${\delta M>0}$). The symmetric state ${M_{1}=M_{2}}$
is a stationary state, and we approximate the time evolution of the mass
difference $\delta M$ by the linear linearized growth law for small mass
differences. Thus, we have

$\partial_{t}\delta M\approx\sigma(\bar{M},\Lambda)\delta M\,,$ (34)

where the growth rate $\sigma$ of the mass-competition instability depends on
the average mass $\bar{M}$ and the spatial separation $\Lambda$ of the peaks
[see Eqs. (31)]. One can then estimate the collapse timescale
$t_{\mathrm{col}}$ as the time until all the mass is transferred from the
smaller to the larger peak, i.e., ${0\approx\bar{M}-\delta
M\operatorname{e}^{\sigma(\bar{M},\Lambda)\,t_{\mathrm{col}}}}$. This yields

$t_{\mathrm{col}}\approx\frac{\log\left(\bar{M}/\delta
M\right)}{\sigma(\bar{M},\Lambda)}\,.$ (35)

Thus, the collapse timescale depends on on the relative size of the initial
mass difference $\delta M/\bar{M}$ and the rate $\sigma$ of mass-competition.
On a large domain showing many simultaneous collapse events, the magnitude of
the mass differences $\delta M$ of the different, competing pairs of peaks can
be estimated by the standard deviation of the peak-size distribution. Within
the scaling hypothesis common to coarsening [49, 50, 52], the shape of the
peak-size distribution is invariant under the coarsening dynamics and only
scales with the average peak mass. Under this scaling hypothesis, one can thus
assume that $\log\left(\bar{M}/\delta M\right)$ is constant on average during
the coarsening process. Moreover, the scaling hypothesis suggests that the
initial mass $\bar{M}$ and the peak separation $\Lambda$, which enter the rate
of mass competition $\sigma$ as arguments, scale with the average peak mass
$\langle M\rangle$ and average peak separation $\langle\Lambda\rangle$,
respectively. Consequently, the average rate of a peak collapse is
proportional to the rate of mass competition ${\sim\sigma(\mu\langle
M\rangle(t),\nu\langle\Lambda\rangle(t))}$ during the coarsening process,
where $\mu$ and $\nu$ are scaling amplitudes.

Using that the average collapse rate per peak is determined by $\sigma$, the
total number of peaks $N(t)$ in a large system of many peaks evolves as

$\partial_{t}N\sim-\sigma\big{(}\mu\langle
M\rangle(t),\nu\langle\Lambda\rangle(t)\big{)}\,N\,.$ (36)

We are still missing the connection between the total peak number $N(t)$, the
average peak mass $\langle M\rangle(t)$, and the average peak separation
$\langle\Lambda\rangle(t)$ to write this evolution equation in a closed form.
This relationship is given by the mass-conservation constraint: The total mass
in the system has to be distributed between the peaks and the background
plateau $\rho_{-}$ from which the peaks rise. This gives the scaling ${\langle
M\rangle\approx\left(\bar{\rho}-\rho_{-}\right)\langle\Lambda\rangle\sim
N^{-1}}$. For simplicity, we set ${\rho_{-}\approx 0}$.151515At late times
during coarsening, we can neglect any changes of the mass in the low-density
plateau from which the peaks rise (see Sec. VI.1).  With this mass-
conservation constraint, we arrive at the closed evolution equation for the
typical pattern length scale $\langle\Lambda\rangle(t)$,

$\partial_{t}\langle\Lambda\rangle\sim\sigma\big{(}\mu\bar{\rho}\langle\Lambda\rangle,\nu\langle\Lambda\rangle\big{)}\,\langle\Lambda\rangle\,.$
(37)

This relation shows that the functional form of the growth rate $\sigma$ of
the mass-competition instability completely determines the coarsening law.

Finally, one can reduce the relationship between $\sigma$ and the length-scale
evolution $\langle\Lambda\rangle(t)$ given in Eq. (37) to an algebraic
relation if the collapse rate $\sigma$ decreases sufficiently strongly as the
peaks grow during the coarsening process [cf. Eqs. (20), (21)]. If the
decrease in the collapse rate is sufficiently strong, the growth of the
average length scale is limited by the duration of the last collapse events.
Then, the time $t$ which is necessary to increase the average peak separation
to $\langle\Lambda\rangle(t)$, scales with the duration of the peak collapses
at size ${\langle
M\rangle(t)\approx(\bar{\rho}-\rho_{-})\langle\Lambda\rangle(t)}$, and one
finds the coarsening law by inversion of the scaling relation

$t\sim\sigma\big{(}\mu\bar{\rho}\langle\Lambda\rangle(t),\nu\langle\Lambda\rangle(t)\big{)}^{-1}\,.$
(38)

This asymptotic scaling can be verified explicitly by integration of the
dynamic equation, Eq. (37), for power-law (peak patterns) or exponential
suppression (mesa patterns) of $\sigma$ with $\langle\Lambda\rangle$.

In summary, we employed a scaling argument [51, 106, 46] to show that the
elementary motive of coarsening—mass competition between neighboring domains
with rate $\sigma$—determines the macroscopic evolution of the average pattern
length scale $\langle\Lambda\rangle(t)$. In particular, for a one-dimensional
system, the scaling analysis results in power-law coarsening for peaks and
logarithmically slow coarsening in mesa-forming systems, independent of
whether these are CH, cAC or 2cMcRD systems. The scaling relation is tested
numerically in Ref. [46]. This comparison also shows that the crossover from
peak to mesa patterns in systems with highly asymmetric $\mathsf{N}$-shaped
nullclines is faithfully predicted [46].

##  VIII Mass competition in the presence of source terms

Building on the understanding of coarsening in 2cMcRD systems, in Ref. [46] a
simple criterion was found to determine the length scale
$\Lambda_{\mathrm{stop}}(\varepsilon)$ of interrupted coarsening in 2cRD
systems with source terms of strength ${0<\varepsilon\ll 1}$. Because the
mass-competition instability is weak, and coarsening slow for large peak or
mesa masses, one expects that only weak source terms are necessary to suppress
the instability and interrupt the coarsening process (cf. Sec. II.1). To give
a basis of interrupted coarsening in 2cRD systems beyond the QSS approximation
used in Ref. [46], we discuss in this section the full growth rate of the
mass-competition instability under the influence of weak source terms. The
detailed mathematical analysis using singular perturbation theory can be found
in Appendix G. We restrict ourselves here to stating the expressions thus
obtained and focus on heuristic considerations that explain their mathematical
structure. We then discuss when the simple criterion for coarsening arrest
based on the QSS approximation for the mass-redistribution potential $\eta$ in
the peak or interface regions applies.

###  VIII.1 The growth rate and stability threshold

To derive the growth rate, we employ the same assumptions as in the mass-
conserving case, i.e., we use the sharp-interface approximation and assume
that mass competition is slow compared to the local relaxation of elementary
stationary patterns (see Sec. VI.1). We analyze the contribution from the weak
source terms perturbatively to first order in the source strength
$\varepsilon$. To analyze the regime in which the mass-competition instability
competes with the stabilizing effect of the source terms, one chooses
$\varepsilon$ such that the stabilization rate
$\varepsilon\sigma_{\mathrm{S}}$ (to be specified below) is of the same order
as the mass-competition rate $\sigma$. For peak patterns, one has to
additionally assume that the profile of the stationary mass-redistribution
potential is approximately constant at the peak (for details see Appendices F,
G). While in the mass-conserving case ${\eta_{\mathrm{stat}}}$ is strictly
constant, this is no longer the case under the influence of weak source terms
[see Eq. (15a) and Appendix F]. Importantly, this assumption does not imply
that the mass-redistribution potential $\eta_{\mathrm{stat}}^{\varepsilon}$ is
uniform _on the scale of the wavelength_ $\Lambda$, i.e., that diffusion on
the scale of the wavelength is fast compared to the reactive timescales. That
is, we do not assume a well-mixed cytosolic reservoir (shadow limit).

The growth rate of the mass-competition instability under the influence of
weak source terms for peak and mesa patterns then has the form (see Appendix
G)

$\sigma^{\varepsilon}\approx\frac{\sigma_{\mathrm{R}}}{\sigma_{\mathrm{D}}+\sigma_{\mathrm{R}}}\,\big{[}\sigma^{\varepsilon}_{\mathrm{D}}-\varepsilon\sigma_{\mathrm{S}}\big{]}\,.$
(39)

For mesa patterns, $\sigma_{\mathrm{D}}$ and $\sigma_{\mathrm{R}}$ have to be
replaced by $\sigma^{\pm}_{\mathrm{D}}$ and $\sigma^{\pm}_{\mathrm{R}}$ that
describe mesa competition or mesa coalescence, respectively [cf. Eqs. (32)].
The rate $\sigma^{\varepsilon}$ consists of two contributions: The first is a
generalization of the growth rate $\sigma$ found in mass-conserving systems
[note that Eq. (30) implies
${\sigma=\sigma_{\mathrm{D}}\sigma_{\mathrm{R}}/(\sigma_{\mathrm{D}}+\sigma_{\mathrm{R}})}$],
while the second, negative term is genuinely new and reflects a stabilization
process due to the added source terms. The influence of the stabilizing
process is quantified by the “stabilization rate” $\sigma_{\mathrm{S}}$. For
mesa competition (superscript $+$) and mesa coalescence (superscript $-$), it
is given by

$\sigma_{\mathrm{S}}^{\pm}=\frac{\left|s_{\mathrm{tot}}^{\pm}\right|}{\Delta\rho}\,,$
(40)

where
${s_{\mathrm{tot}}^{\pm}=s_{\mathrm{tot}}(\rho_{\pm},\eta_{\mathrm{stat}}^{\infty})}$
and ${\Delta\rho=\rho_{+}-\rho_{-}}$. For peak patterns, we find the
stabilization rate as

$\sigma_{\mathrm{S}}=|\langle\partial_{\rho}s_{\mathrm{tot}}\rangle_{\mathrm{int}}|\,.$
(41)

The rate $\sigma^{\varepsilon}_{\mathrm{D}}$ constitutes a generalization of
the diffusion-limited mass-competition rate $\sigma_{\mathrm{D}}$ [cf. Eqs.
(31a), (32a)],

$\sigma^{\varepsilon}_{\mathrm{D}}=-\frac{4D_{v}}{\Lambda}\,\partial_{M}\big{(}{\eta_{\mathrm{stat}}}+\varepsilon\,\delta\eta_{\mathrm{stat}}^{\varepsilon}\big{)}\,,$
(42)

and implicitly depends on the source strength $\varepsilon$. The additional
term $\varepsilon\,\delta\eta_{\mathrm{stat}}^{\varepsilon}$ accounts for the
shift of the _stationary_ mass-redistribution potential at the stationary peak
or interface due to the weak source terms. It is given by (see Appendix F)

$\delta\eta_{\mathrm{stat}}^{\varepsilon}=-\frac{\langle
s_{1}+d\,s_{2}\rangle_{\mathrm{int}}}{\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}}\,.$
(43)

Notably, for mesa patterns, the shift in the stationary mass-redistribution
potential is independent of the mass $M$,
(${\partial_{M}^{\pm}\delta\eta_{\mathrm{stat}}^{\varepsilon}=0}$) because the
shift, Eq. (43), only depends on the interface profile, which does not change
with the mesa $M$ within the sharp-interface approximation (the interface is
only shifted, see Sec. V). In other words,
$\delta\eta_{\mathrm{stat}}^{\varepsilon}$ does not contribute to gradients in
$\eta$ between mesas with different masses. This implies that the contribution
to the growth rate $\sigma^{\varepsilon}$ that stems from mass competition
[first term in Eq. (39)] is unchanged compared to the mass-conserving case for
mesa patterns
(${\sigma_{\mathrm{D}}^{\varepsilon}=\sigma^{\pm}_{\mathrm{D}}}$).

Diffusion- and reaction-limited regimes. — As in the mass-conserving case, we
distinguish a diffusion- and a reaction-limited regime for the growth rate
$\sigma^{\varepsilon}$. In the diffusion-limited regime
[${D_{v}\ll\Lambda\ell_{\mathrm{int}}\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}}$,
see Eq. (33)], the growth rate $\sigma^{\varepsilon}$ simplifies to (we reason
in Appendix G why $\sigma^{\varepsilon}_{\mathrm{D}}$ may be replaced by
$\sigma_{\mathrm{D}}$)

$\sigma^{\varepsilon}\approx\sigma_{\mathrm{D}}-\varepsilon\sigma_{\mathrm{S}}\,,$
(44)

where one again substitutes ${\sigma_{\mathrm{D}}=\sigma_{\mathrm{D}}^{\pm}}$
in the case of mesa competition and mesa coalescence, respectively. In this
limiting regime, the rates of mass transport and source production or
degradation simply add up because single domains and the value of the mass-
redistribution potential at the domains can be approximated by a QSS. The mass
of the quasi-stationary domains then slowly changes by mass transport between
these and the production or degradation of particles within the domains.

In contrast, in the reaction-limited regime
(${D_{v}\gg\Lambda\ell_{\mathrm{int}}\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}}$),
the shift $\varepsilon\,\delta\eta_{\mathrm{stat}}^{\varepsilon}$ becomes
relevant and one has to distinguish between peak and mesa patterns (a scaling
argument is given at the end of Appendix G). The growth rate
$\sigma^{\varepsilon}$ for mesa patterns takes the form

$\sigma^{\varepsilon}\approx\sigma_{\mathrm{R}}^{\pm}-\varepsilon\frac{\sigma_{\mathrm{R}}^{\pm}}{\sigma_{\mathrm{D}}^{\pm}}\sigma_{\mathrm{S}}\,.$
(45)

This corresponds to the result in the diffusion-limited regime but scaled by
the ratio of the reaction- and diffusion-limited rates of mass competition.
The scaling of the magnitude of the rate accounts for the reactive limitation
of the rate for the mass change of single domains while the shift
$\varepsilon\,\delta\eta_{\mathrm{stat}}^{\varepsilon}$ does not contribute
for mesa patterns
(${\partial_{M}\delta\eta_{\mathrm{stat}}^{\varepsilon}\approx 0}$). At first
order in $\varepsilon$, these growth rates for mesa patterns in the reaction-
limited regime agree with the growth rates derived by McKay and Kolokolnikov
[59] when evaluating their expressions for nearest-neighbor
competition.161616For comparison, we identify their reaction terms $f$ and $g$
with our terms ${f+\varepsilon s_{1}}$ and ${-f+\varepsilon s_{2}}$, the
diffusion constants ${\varepsilon^{2}=D_{u}}$, ${D=D_{v}}$, and set
${\tau=1}$. Then one obtains ${g_{\pm}=\varepsilon s_{\mathrm{tot}}^{\pm}}$
and
${\alpha_{\pm}=-\Delta\rho\,\partial_{M}^{\pm}{\eta_{\mathrm{stat}}}[1+\mathcal{O}\\!\left(\varepsilon\right)]}$.

For peak patterns in the reaction-limited regime, we find instead

$\sigma^{\varepsilon}\approx\sigma_{\mathrm{R}}-2\ell_{\mathrm{int}}\langle\tilde{f}_{\eta}\rangle_{\mathrm{int}}\,\varepsilon\,\partial_{M}\delta\eta_{\mathrm{stat}}^{\varepsilon}-\varepsilon\frac{\sigma_{\mathrm{R}}}{\sigma_{\mathrm{D}}}\sigma_{\mathrm{S}}\,,$
(46)

where the middle term is due to the shift of the stationary mass-
redistribution potential at the peak. The last term is negligible in
comparison to the second if the source terms fulfill ${s_{1}\gtrsim s_{2}}$
[see Eq. (43)]. We will explain the structure of the growth rates in the
following section and also discuss what causes the shift
$\varepsilon\,\delta\eta_{\mathrm{stat}}^{\varepsilon}$. Before, the stability
threshold is analyzed.

Stability threshold. — For a sufficiently large source strength or large
wavelength, $\sigma^{\varepsilon}$ may become negative as the stabilization
rate $\varepsilon\sigma_{\mathrm{S}}$ increases or the rate of destabilizing
mass transport $\sigma_{\mathrm{D}}^{\varepsilon}$ decreases, respectively. A
negative growth rate $\sigma^{\mathrm{\varepsilon}}$ indicates that the mass-
competition instability is suppressed, and a mass difference $\delta M$
between the domains decreases exponentially in time, ${\delta
M\sim\operatorname{e}^{\sigma^{\varepsilon}t}}$. Because the coarsening
process increases the pattern wavelength $\Lambda$ due to the collapse or
coalescence of peaks or mesas, it is expected to interrupt at the stability
threshold ${\sigma^{\varepsilon}(\Lambda,\varepsilon)=0}$. This behavior of
patterns with several peaks or mesas is analyzed in Sec. VIII.3. Using Eqs.
(39), (40), the criterion ${\sigma^{\varepsilon}(\Lambda,\varepsilon)=0}$
yields for mesa patterns the critical source strength

$\varepsilon_{\mathrm{stop}}(\Lambda)=\frac{\sigma^{\pm}_{\mathrm{D}}(\Lambda)}{\left|s_{\rho}^{\pm}\right|/\Delta\rho}\,.$
(47)

For peak patterns, the stability threshold follows from Eqs. (39), (42), (41)
as

$\varepsilon_{\mathrm{stop}}(\Lambda)=\frac{\sigma_{\mathrm{D}}(\Lambda)}{-\langle\partial_{\rho}s_{\rho}\rangle_{\mathrm{int}}+\frac{4D_{v}}{\Lambda}\partial_{M}\delta\eta_{\mathrm{stat}}^{\varepsilon}}\,,$
(48)

where the dependence on $\Lambda$ appears implicitly in the various terms. The
critical wavelength $\Lambda_{\mathrm{stop}}(\varepsilon)$ is found by
inverting the above relations. In general, it is not possible to perform this
inversion analytically.

The growth rate Eq. (44) in the diffusion-limited regime as well as the
threshold of interrupted coarsening for mesa patterns, Eq. (47), recover the
expressions found based on the QSS approximation in Ref. [46]. However, the
mathematical structure of the general growth rate Eq. (39) is quite intriguing
and demands further analysis.

###  VIII.2 The mathematical structure of the growth rates

In the following, we will analyze how the different terms of the growth rate
$\sigma^{\varepsilon}$ arise [Eq. (39)]. We will find that the rate
$\sigma_{\mathrm{S}}$ accounts for net production and degradation at the peaks
or mesas. We then show that the source terms determine, together with mass
transport in the cytosol, whether a peak or mesa accumulates or loses mass.
The rate $\sigma_{\mathrm{R}}$ of reactive conversion (particle release and
incorporation) will be shown to only affect how fast the accumulated mass can
be incorporated into a peak or mesa (or how fast mass is released from the
peak/mesa). Thus, the rate $\sigma_{\mathrm{R}}$ only affects the magnitude of
the growth rate $\sigma^{\varepsilon}$, and enters $\sigma^{\varepsilon}$ only
as a prefactor. To begin, we first consider the effect of the source terms
alone and explain the stabilization rate $-\varepsilon\sigma_{\mathrm{S}}$.

Stabilization rate due to weak source terms. — The relaxation rate
$-\varepsilon\sigma_{\mathrm{S}}$ describes the direct stabilizing effect of
the source terms. In the case of two competing half-mesas (cf. Fig. 2, mesa
competition/coalescence and also peak coalescence), a shift in the interface
positions (the peak) by $\delta\ell$ due to mass transfer from one mesa to the
other only changes the length of the outer plateaus, i.e., the high-density
plateaus for mesa competition and the low-density plateaus for mesa/peak
coalescence. For specificity, we now analyze mesa competition. In the high-
density plateaus, the source terms lead to net degradation, i.e.,
${s_{\mathrm{tot}}(\rho_{+},{\eta_{\mathrm{stat}}})=s_{\mathrm{tot}}^{+}<0}$.171717As
shown in Appendix F, single (half-)mesas are stable only if
${s_{\mathrm{tot}}(\rho_{+},{\eta_{\mathrm{stat}}})<0<s_{\mathrm{tot}}(\rho_{-},{\eta_{\mathrm{stat}}})}$.
Similarly, in peak patterns, stability mandates that degradation increases for
larger peaks, i.e.,
${\langle\partial_{\rho}s_{\mathrm{tot}}\rangle_{\mathrm{int}}<0}$. Thus, an
increase of mass ${\delta M=\Delta\rho\,\delta\ell}$ in one (half-)mesa there
leads to additional degradation ${\varepsilon
s_{\mathrm{tot}}^{+}\delta\ell=\partial_{t}^{\mathrm{source}}\delta M<0}$, and
to additional production ${-\partial_{t}^{\mathrm{source}}\delta M>0}$ in the
other (half-)mesa. Consequently, production and degradation together result in
a relaxation of the initial perturbation and drive the interfaces back to
their symmetric rest position with the rate
${-\varepsilon\sigma_{\mathrm{S}}=-\varepsilon|s_{\mathrm{tot}}^{+}|/\Delta\rho}$
(respectively,
${-\varepsilon\sigma_{\mathrm{S}}=-\varepsilon|s_{\mathrm{tot}}^{-}|/\Delta\rho}$
for mesa/peak coalescence). Analogously, in the case of peak competition,
additional mass $\delta M$ increases the peak size and leads to increased
degradation

$\displaystyle\partial_{t}^{\mathrm{source}}\delta M$
$\displaystyle\approx\varepsilon\int_{0}^{\frac{\Lambda}{2}}\mathrm{d}x\,s_{\mathrm{tot}}(\rho_{\mathrm{stat}}+2\delta
M\partial_{M}\rho_{\mathrm{stat}},$ $\displaystyle\hskip
73.97733pt{\eta_{\mathrm{stat}}}+2\delta
M\partial_{M}{\eta_{\mathrm{stat}}})\,.$ (49)

Using that $\partial_{M}{\eta_{\mathrm{stat}}}$ becomes negligible for
sufficiently large peaks [cf. Eqs. (20), (21) and Appendix F] one arrives at
the simplified mass evolution

$\displaystyle\partial_{t}^{\mathrm{source}}\delta M$ $\displaystyle\approx
2\delta
M\varepsilon\int_{0}^{\frac{\Lambda}{2}}\mathrm{d}x\,\left(\partial_{\rho}s_{\mathrm{tot}}\right)\,2\partial_{M}\rho_{\mathrm{stat}}$
$\displaystyle=\varepsilon\langle\partial_{\rho}s_{\mathrm{tot}}\rangle_{\mathrm{int}}\delta
M\,,$ (50)

where the last equality follows from the definition of the interface average
$\langle\cdot\rangle_{\mathrm{int}}$ [see Eq. (94)]. Because
${\langle\partial_{\rho}s_{\mathrm{tot}}\rangle_{\mathrm{int}}<0}$ is mandated
by the stability of a single peak (see Appendix D), particle production and
degradation at the peaks stabilize the symmetric peak configuration. The