Wavelength selection by interrupted coarsening in reaction-diffusion systems

Wavelength selection in reaction-diffusion systems can be understood as a coarsening process that is interrupted by counteracting processes at a certain wavelength. We show that coarsening in mass-conserving reaction-diffusion systems is driven by self-amplifying mass transport between neighboring high-density domains. We derive a general coarsening criterion and show that coarsening in two-component systems is generically uninterrupted. The theory is then generalized to study interrupted coarsening and anti-coarsening due to weakly-broken mass conservation, providing a general path to analyze wavelength selection in pattern formation far from equilibrium.

Wavelength selection in reaction-diffusion systems can be understood as a coarsening process that is interrupted by counteracting processes at a certain wavelength. We show that coarsening in mass-conserving reaction-diffusion systems is driven by self-amplifying mass transport between neighboring high-density domains. We derive a general coarsening criterion and show that coarsening in two-component systems is generically uninterrupted. The theory is then generalized to study interrupted coarsening and anti-coarsening due to weakly-broken mass conservation, providing a general path to analyze wavelength selection in pattern formation far from equilibrium.
Coarsening-the growth of the typical length scale of patterns-is observed in a wide range of non-equilibrium systems, including cell polarization [1][2][3], chemotactic cells [4,5], granular media [6], and active matter systems [7][8][9][10][11]. While coarsening is well understood as minimization of the free energy for systems relaxing to thermal equilibrium (e.g. for phase separation in binary mixtures), this reasoning is generally not applicable for non-equilibrium systems whose steady states break detailed balance. For one-component systems, a theory for coarsening based on a multiple-scale analysis has been developed [12,13], but generalizations to multi-component systems have remained elusive.
In this Letter, we propose that wavelength selection in reaction-diffusion systems can be understood as a coarsening process that is interrupted by counteracting processes at a certain wavelength. We start by developing a theory for the coarsening dynamics in massconserving reaction-diffusion (McRD) systems with two components (2cMcRD), building on the recently introduced local equilibria theory [14,15]. These systems serve as paradigmatic models for intracellular pattern formation [2,3,[16][17][18][19][20], and are used as phenomenological models for a wide range of systems including precipitation patterns [21], granular media [6], and braided polymers [22]. It has long been speculated that 2cM-cRD systems generically exhibit uninterrupted coarsening [1,3,18,23]. For specific mathematical forms of the reaction kinetics, an effective free energy (Lyapunov function) can be constructed such that uninterrupted coarsening follows from a minimization argument [23,24]. However, a systematic understanding of the coarsening dynamics for generic reaction kinetics is missing, largely owing to a lack of insight into the underlying physical processes. Whether coarsening goes to completion in all 2cMcRD systems remains unclear.
Here, we show that coarsening is driven by positive feedback in the competition for mass, derive a simple and quantitative description of coarsening dynamics, and explain why coarsening is generically uninterrupted in 2cM-cRD systems. As they are grounded in a recently intro-duced phase-space analysis for 2cMcRD systems [14], our results are independent of the specific mathematical form of the reaction kinetics. Generalizing this phase-space analysis, we elucidate and quantify the physical mechanisms underlying interrupted coarsening and domain splitting in the presence of weak source terms (weakly broken mass conservation). Since our approach builds on studying the spatial redistribution of a (nearly) conserved quantity, we expect that it can be generalized beyond two-component reaction-diffusion systems; for instance, to systems with more components, to extensions of the Cahn-Hilliard equation [25,26] and to hydrodynamic models for active matter systems [8,11,[27][28][29].
In the (strictly) mass-conserving case, the general form of a 2cMcRD system with components u and v is on a domain Ω, where the averageρ = |Ω| −1 Ω dxρ(x, t) of the total density ρ = u+v is conserved. For specificity, we choose D u < D v . The time evolution of the conserved density ρ is given by [1,2,14,22] ∂ t ρ(x, t) = D v ∇ 2 η(x, t) (2) with the mass-redistribution potential defined by η := v + (D u /D v )u; the corresponding dynamical equation for η(x, t) is given in the SM. For stationary patterns [ũ(x),ṽ(x)], the mass-redistribution potential must be spatially uniform, η(x) = η stat . Based on this one can analyze 2cMcRD systems in the (u, v) phase plane [14]: There, stationary patterns are constrained to a linear subspace, v + (D u /D v )u = η stat , called flux-balance subspace (FBS); see Fig. 1b. The intercept η stat is determined by the balance of the spatially integrated reactive flows (total turnover balance), corresponding (approximately) to a balance of areas (shaded in red in Fig. 1b (a) Illustration of a stationary peak with peak mass M . Increasing the mass to M + δM leads to an increases the peak amplitude toû + δû. (b) Representation of the stationary peak in phase space (thick blue line), which is constrained to the FBS (dashed blue line). The FBSoffset ηstat(M ) is determined by a balance of total reactive turnovers (areas shaded in red). For a peak with increased mass M + δM , and thus increased peak amplitude δû, the FBS shifts downwards δηstat until total turnover balance is restored (balance of green-shaded areas). (c) After a perturbation of two identical stationary peaks, the gradient in the mass-redistribution potential η (orange line) drives masstransport between the peaks (orange arrow) such that the larger (smaller) peak grows (shrinks) further (blue arrows). peaks. A peak forms when the maximum density does not saturate in a high-density plateau (Fig. 1a, compare Fig. 2a) [14] [30]. The relationship between the total mass M of a stationary peak/mesa and the FBS-offset η stat (M ), as determined by total-turnover balance, will turn out to be central for the coarsening dynamics. We begin the analysis with peak patterns and then generalize the results to mesas.
A mass-competition instability drives coarsening. Coarsening requires the transport of mass between peaks. Because mass transport is diffusive, it is fastest on the shortest length scales; hence, the dominant process is competition for mass between neighboring peaks (Fig. 1a). Thus, as an elementary case, we study two peaks in a 'box' with no-flux boundary conditions. Consider a situation ("coarsening limit") where the peaks are well separated, such that diffusive transport is limiting. We can then approximate the peaks to be in (regional) quasi-steady state (QSS), such that η = η stat (M ) at a given peak with total mass M . This is akin to the local thermal equilibrium approximation in Lifshitz-Slyozov-Wagner theory for Ostwald ripening [31,32]. A technical discussion on the validity of this QSS approximation (based on multi-scale analysis and renormalization group theory [33]) will be published elsewhere.
Starting from two identical, stationary peaks, each with total mass M 0 , the dynamics of the mass difference between them (M R,L = M 0 ± δM )-obtained by integration of Eq. (2) over a single peak-is determined by the η-gradients in the plateau between them (indicated by the orange arrow in Fig. 1c). Using QSS at each peak separately, the mass-redistribution potential at the peaks is given by η R,L = η stat ±(∂ M η stat | M0 ) δM . The resulting gradient in η is linear because diffusive relaxation within the plateau is fast compared to the change in η R,L caused by the flux of mass from one peak to another (see SM for details). For a given peak separation Λ, this approximation determines the dynamics of mass redistribution The subscript D the diffusion-limited regime. If the growth rate σ D is positive, an instability driven by positive feedback in competition for mass results in coarsening. Hence, the condition for uninterrupted coarsening reads i.e. that η stat (M ) is a strictly monotonically decreasing function for all stable stationary single-peak solutions. This recovers a previous, mathematically derived coarsening condition [1,2]. Importantly, the analysis presented here gives insight into the underlying physical mechanism and shows that not only the criterion for coarsening, but the entire temporal evolution of coarsening is determined by ∂ M η stat via Eq. (3) [34]. We learn that the functional dependence of the mass-redistribution potential on the peak mass, η stat (M ), plays a role analogous to the functional dependence of the chemical potential on the droplet size that drives Ostwald ripening. Moreover, the mass-competition instability driving coarsening is a manifestation of the mass-redistribution instability (introduced in Ref. [14]), on the level of elementary patterns (peaks/mesas). Generic coarsening laws for mass-conserving systems. A simple geometric argument shows that the coarsening criterion, Eq. (4), generically holds for both peak and mesa patterns, regardless of the specific mathematical form of the reaction term f (u, v). Consider a single stationary peak with mass M (see Fig. 1a) and its representation in phase space, the blue line in Fig. 1b. Add an amount δM of mass and hold η stat fixed for the moment (for the sake of argument). The additional mass will increase the peak amplitudeû (Fig. 1b), causing the reactive turnover to the right of u 0 to increase (Fig. 1b). The resulting imbalance of total turnover entails a net reactive flow that shifts the flux-balance subspace downwards, i.e. lowers η stat , to restore total turnover balance. Therefore, η stat (M ) is generically a monotonically decreasing function. (In the SM, we present a more rigorous argument to show that ∂ M η stat < 0 holds for stable peaks.) (a) Illustration of the peak to mesa transition as the total mass M is increased. (b) The function ηstat(M ) obtained by numerical continuation of the stationary solutions for the reaction kinetics fex. Crossover from power law for peak patterns (amplitude not saturated) to exponential approach to η ∞ stat for mesa patterns. (c) Coarsening dynamics from finite element simulations for fex (black circles; mean peak distance averaged over four independent runs started from random initial conditions; parameters: Du = 1, Dv = 10 4 ,ρ = 1.5 and system size |Ω| = 2 × 10 5 ). The red line shows the analytic prediction based on σ D from ηstat(M ), shown in (b), via Eq. (3). After an initial transient, power-law coarsening Λ ∼ t 3/8 for peaks is observed, which flattens into logarithmic coarsening for mesas.
As an example, consider a reaction kinetics with where the first and second terms may, for instance, describe protein recruitment and first-order enzymatic detachment, respectively. A simple scaling argument [35] yields a power-law relation η stat (M ) ∼ M −α , where the exponent depends on the specific reaction kinetics (for the example here α = 2 /3); see Fig. 2b. In a large system containing multiple peaks, the average peak separation Λ is linked to the characteristic peak mass by M = (ρ − ρ − ) Λ , where ρ − is the total density in the low density plateau between the peaks, and · denotes an average over the entire system. As peaks collapse, with a typical time given by the inverse growth rate of the masscompetition instability t ∼ σ −1 D , the average peak separation Λ will increase. Combining σ D ∼ − ∂ M η stat / Λ with ∂ M η stat ∼ M −α−1 ∼ (ρ Λ ) −α−1 yields powerlaw coarsening with Λ (t) ∼ t 1/(2+α) ; see Fig. 2c and Fig. S2. Moreover, using appropriate scaling amplitudes, the coarsening trajectories for different average masses ρ can be collapsed onto a single master curve obtained from ∂ M η stat (see SM, Fig. S1).
As peaks collapse, those remaining grow in mass and height. When the density at the peak maximum saturates in a high-density plateau (corresponding to a FBS-NC intersection point in phase space), a mesa pattern starts to form (Fig. 2a). For such mesas, somewhat more subtle arguments show that η stat (M ) remains a monotonically decreasing function (see SM). In essence, changing M shifts the interface positions and thus changes the width of a mesa's plateau. As the exponential tails in the density profile approaching the limiting plateau u ± (η ∞ stat ), the concentration maximum increases only exponentially slowly as the mesa width w grows; here, we define η ∞ stat as the limit of η stat for the stationary pattern on an infinite domain (see SM). As a result, η stat (M ) approaches η ∞ stat exponentially slowly (see inset in Fig. 2b).
Using the same scaling arguments as for peaks, one obtains a logarithmic coarsening law for all mesa patterns, as in the one-dimensional Cahn-Hilliard model [36]. For the concrete example f ex , we find excellent agreement between finite-element simulations and Λ (t) obtained from η stat (M ) by these scaling arguments (see Fig. 2c).
The limit of large D v . For D v → ∞, mass redistribution by v-diffusion becomes instantaneous, such that the reactive conversion between u and v, which drives the growth/shrinking of mesas or peaks, becomes limiting. In this reaction-limited case, we find σ R ≈ (∂ M η stat ) int f v int , where · int denotes the average over the interface region (see SM for details and numerical verification). Comparing with Eq. (3) shows that the coarsening criterion Eq. (4) holds in both regimes, and the crossover from diffusion-to reaction-limited coarsening occurs at Weakly broken mass conservation. We now incorporate slow production and degradation processes in form of source terms s 1,2 into the dynamics, with a small, dimensionless source strength ε. The time evolution of the total density then contains the total source s := s 1 + s 2 Hence, the average mass ρ is no longer a control parameter but a time-dependent variable that is determined indirectly by a balance of production and degradation (in short: source balance). In phase space, there are now two reactive nullclines, one each for u and v, which both converge to f = 0 for ε → 0. Their intersection point(s) determine(s) the homogeneous steady state (HSS) ρ hss that balances the total source term. Consider a mesa pattern. Along the plateaus, the spatial gradients induced by slow productiondegradation (ε small) are shallow, such that the dynamics is (approximately) slaved to the nullcline f = 0 (see Fig. 3a) justifying a local equilibrium approximation On the short scale of the interface width, the weak source term is negligible, such that the densities at each interface are constrained to a flux-balance subspace. Then the 'half lengths' L ± of the upper and lower plateaus (see Fig. 3a) are determined from source balance 0 ≈ s * (ρ + )L + + s * (ρ − )L − together with L + + L − = Λ/2. Here, we approximate the plateau densities by ρ ± + O(ε) and neglect contributions from the interface ( int Λ).
We are now in a position to generalize the phase-space analysis introduced in Ref. [14] and analyze interrupted coarsening [25,27,37] and peak/mesa splitting [38-40].
(i ) Peak/mesa splitting. Consider the fully coarsened state for ε = 0 and add a small source term such that s * (ρ + ) < 0 and s * (ρ − ) > 0 (i.e. ρ − < ρ hss < ρ + , see Fig. 3a) [41]. The upper plateau is depressed by net degradation and is refilled by inflow from the interfaces that connect to the lower plateau where net production prevails. The longer the plateaus (and the larger ε), the more they curve towards ρ hss . Since ρ − < ρ hss < ρ + , ρ(x) will eventually enter the interval of lateral instability [ρ − lat , ρ + lat ] (where ∂ ρ η * < 0), triggering a regional lateral instability resulting in splitting of the mesa (see Fig. 3a and Movie 2) [42]. A simple approximation for the threshold wavelength Λ split (ε) where this happens is derived in the SM. Comparison with numerical simulations shows excellent agreement (see Fig. 3c).
In two and more dimensions, additional mechanisms, such as shape instabilities, can lead to domain splitting [43][44][45]. Studying such instabilities building on the phase space analysis presented here is an interesting avenue for future research.
(ii ) Interrupted coarsening. Intuitively, production and degradation can counteract the mass-competition instability. To determine the corresponding length scale Λ stop where coarsening arrests, we consider the stability of two neighboring, symmetric mesas. A perturbation that moves a small amount of mass from one mesa to the other (Fig. 3b) has two effects: First, it shifts the mass-redistribution potential at the interfaces, leading to mass transport that further amplifies the perturbation with rate σ D (Λ)δM as in the strictly massconserving situation; cf. Eq. (3). Second, the changed lengths δL = δM/(ρ + − ρ − ) of the two mesas result in net production (degradation) in the shorter (longer) mesa with rate ε|s * (ρ outer )|δL (indicated by the purple arrows in Fig. 3b). Here ρ outer denotes the total density of the outer plateau (the inner plateau shifts as a whole and does not change in length, see Fig. 3). Together, the balance of both processes determines Λ stop (see SM for details) As a concrete example, we apply Eq. (7) to the "Brusselator" model [46] (f = u 2 v − u, s = p − u), and find excellent agreement with numerics ( Fig. 3c). Notably, the simple estimate given by Eq. (7) generalizes a previous, mathematically obtained results [37,47]. Our analysis shows that the mechanisms underlying mesa splitting and interrupted coarsening are distinct. Notably, the length scale where coarsening stops is much smaller than the length scale where mesas/peaks split (see Fig. 3c). This implies that there are stable periodic patterns for a large, continuous range of wavelengths (multistability), as was shown previously for the 'Brus-selator' [37,38,46]. Similarly, multistability of wavelengths was recently found in a hydrodynamic model for flocking [28]. Interestingly, a unique length scale is selected once noise is accounted for [48]. Noise-driven wavelength selection was also observed in an "active Model B" [49]. It would be interesting to study whether this phenomenon is also found in reaction-diffusion systems.
Outlook. We have combined the physical mechanism of mass transport with phase-space geometry to provide a new avenue for studying length-scale selection in pattern-forming systems. We expect that our approach can be generalized to systems with more than two components, higher spatial dimensions and also beyond reaction-diffusion systems. In particular, conserved densities (particle numbers) are a generic feature of many active matter systems in which coarsening and lengthscale selection ("micro-phase separation") are of growing interest [5, 7-11, 27-29, 45, 49-51].
This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)project-id 201269156 -Collaborative Research Center (SFB) 1032 -Project B2.
FB and HW contributed equally to this work.
stat , for sufficiently large M and Dv Du. Following a similar scaling analysis for other reaction terms, e.g. with different nonlinearities in the recruitment term, yields other exponents α. [41] In steady state, net degradation in high-density regions (ρ > ρ hss ) and net production in low-density regions (ρ > ρ hss ) must balance.

Basic setup, examples, and numerics
In this section, we provide the basic setup and establish the notation for the following analysis. In particular, we establish the notion of "elementary patterns" and their defining equations. Moreover, we give an overview of the models used to exemplify our results in the subsequent sections and give a brief explanation of the numerical methods employed.
Sections 2-6, regard mass-conserving reaction-diffusion systems with two components (2cMcRD): on a domain Ω with periodic or no-flux boundary conditions, conserving the average In Sec. 7, we generalize our analysis to systems with weakly broken mass conservation by including weak source terms on the RHS of Eq. (1); see Eq. (6) in the main text. We restrict our analysis mostly to one spatial dimension and address generalization to more spatial dimensions briefly in Sec. 5.4.
Let us denote the stationary solutions of Eq.
Thanks to diffusive flux balance, extrema inũ andṽ always coincide. Therefore, any stationary pattern can be dissected, by inserting no-flux boundaries at its extrema, into monotonic pieces that we will refer to as elementary or single-interface patterns. In particular, periodic stationary patterns with a wavelength Λ can be constructed by translation and reflection of the elementary pattern on the domain [0, Λ/2] with no-flux boundary conditions. Multiplying Eq. (2) by ∂ xũ (x) and integrating over the entire domain of an elementary pattern (from extremum to extremum) yields the total turnover balance condition which determines η stat .

Tranformation to (ρ, η) variables
Using the variables ρ = u + v and η = du + v, the 2cMcRD dynamics Eq. (1) can be cast as with the reaction termf For stationary patterns [ρ(x),η(x)], Eq. (5a) together with no-flux or periodic boundary conditions mandates thatη(x) = η stat be spatially uniform. This is the diffusive fluxbalance condition. From Eq. (5b) the stationary density profile of an elementary pattern,ρ(x), is then determined by subject to the total density constraint

Example models
Throughout the SI, we use four different models to illustrate our results, compare analytics to numerics, and relate our work to previous literature. Each model and its specific purpose are briefly described below. Because the models serve only illustrative purposes, we do not specify units.
Simple recruitment and enzyme-driven detachment. To demonstrate the peak-mesa crossover during coarsening (see Fig. 2b,c in the main text), we use a model motivated by intracellular protein kinetics where u, v ≥ 0 represent protein concentration on the membrane and in the cytosol, respectively. The first term represents self-recruitment to the membrane and the second term accounts for enzyme-driven membrane detachment (Michaelis-Menten kinetics).
Ref. [1] with all rate constants set to unity.) Otsuji's "Model II". A reaction term that exhibits only peak patterns and no crossover to mesas for large Λ was introduced in Ref. [2], where it was studied as a minimal toy model in the context of cell polarization and coarsening: with u, v ∈ R. The nullcline of this model becomes asymptotically parallel to the FBS. This becomes immediately evident when f is transformed to (ρ, η) variables: which yields the nullcline segments ρ = 0 and η * (ρ) = 1/ρ. Hence, there is no third FBS-NC intersection point (for ρ ≥ 0) such that peaks grow unlimited (never cross over to mesas). * This property is useful to demonstrate the power law coarsening of peaks, since the power law is not cut off at large masses by the crossover to mesa patterns (see Sec. 3).
Brusselator. To exemplify our results for nearly mass-conserving dynamics, we use a non-dimensionalized form of the Brusselator, a paradigmatic, phenomenological model for pattern formation [3]: with u, v ≥ 0. We decompose the reaction kinetics into a mass-conserving 'core' term f = u 2 v − u and a source term εs = ε (p − u). We will refer to the mass-conserving version (ε = 0) as 'Brusselator core'. Since the stationary elementary pattern of the 'Brusselator core' on the infinite line has a simple analytic form, closed-form expressions can be obtained for its coarsening rate (see Sec. 5).
For the Brusselator, mesa-splitting and the stability of periodic patterns in the limit of small ε were previously studied mathematically in Refs. [4][5][6]. We will compare our results to these previous works in Sec. 7.3 below.
Cubic model. Arguably, the most elementary nonlinearity is the cubiĉ where ρ, η ∈ R can take both positive and negative values, i.e. ρ should not be thought off as a density but rather as a generic order parameter. This "reaction term" is constructed such that the stationary patterns, determined by Eq. (7), do not depend on the ratio of the diffusion constants, specifically not on D v . This property is useful to illustrate the crossover from diffusion-to reaction-limited coarsening as D v is increased (see Sec. 6). Note that casting Eq. (13) in terms of u and v yields a "reaction term" that depends on the ratio of the diffusion constants d.
Further, a nearly mass-conserving version of the simple cubic model with a source term εs = ε(p − ρ) is used to illustrate mesa splitting (see Movie 2).

Numerical methods
Finite element simulations. Performed in COMSOL Multiphysics (Version 5.4) [7]. Setup files for COMSOL simulations and Python scripts used to evaluate the simulation data are available at https://github.com/f-brauns/2cMcRD-coarsening.
Numerical continuation and numerical stability analysis. To numerically find stationary patterns and their linear stability we used a finite differences discretization (second-order central differences scheme). The equations for the stationary patterns (Eq. (2) under the constraint (3) in the case of strictly mass-conserving systems, Eqs. (12) in the case of the nearly mass-conserving Brusselator), become sets of algebraic equations that are solved with a Newton algorithm. In addition, we implemented a pseudo-arclength continuation method (see e.g. Ch. 4 in Ref. [8]) to follow solution branches as a function of a free parameter. This is particularly useful to obtain the relation η stat (M ), from which the coarsening law can be calculated (cf. Eq. (3) in the main text). Linear stability of the stationary patterns, e.g. the growth rate of the competition instability, was obtained by numerically solving the eigenvalue problem obtained by linearizing the discretized dynamics around the stationary solutions.

Gradients in η are linear between peaks/interfaces
In the derivation of Eq. (3) for the dynamics of δM in the main text, we argue that there is a linear gradient in η in the plateau region between the peaks (see Fig. S1). This approximation is justified as follows. First, gradients in the plateau are shallow (by definition) such that the local reactions there are much faster than diffusion. Therefore the slow dynamics in the plateau is purely diffusive (we can set f ≈ 0 in Eq. (5b) on the slow timescale). Moreover, the amount of mass redistributed as the plateau relaxes to its linear QSS profile is much smaller than the mass transported from one peak to another during coarsening. This ensures that relaxation of the η-profile to a linear gradient in the plateau is fast compared to the coarsening process it drives.
A more detailed analysis and a precise criterion for the validity of this approximation will be presented in a forthcoming technical paper.

Scaling relations for the coarsening law
The collapse of a peak due to competition for mass with larger, neighboring peaks is determined by the growth rate of the mass-competition instability. In the main text, we derived this growth rate in the diffusion-limited regime: The inverse growth rate σ D (M, Λ) −1 determines the timescale for a peak of mass M to collapse in the diffusion limited regime, with its distance to neighboring peaks given by Λ.
On a large spatial domain containing many peaks, one expects a distribution of peak masses where the average mass of the single peaks is approximatelyρ Λ because for large peaks, ρ − Λ is small compared to the peak mass. The mass-competition instability causes the smallest peaks to collapse the fastest, while the largest peaks grow. It is, therefore, not necessarily the mean peak massρ Λ that determines the typical collapse time. However, during coarsening the distribution of peak masses typically evolves to a shape whose only characteristic scale is the mean peak massρ Λ (t) (cf. the droplet size distribution in Ostwald ripening). Assuming such a scaling behavior we expect the mean collapse time to scale with σ D (µρ Λ , Λ ) −1 , where we introduce the scaling amplitude µ to account for the distribution of peak masses.
Next, we use that −∂ M η stat is a rapidly decreasing function of M (power law for peaks, see main text; exponential for mesas, see Sec. 5.1), such that the collapse time increases as peaks become larger. Therefore, the time t that passes until an average peak massρ Λ (t) is reached during coarsening is dominated by the duration of the most recent collapse σ −1 D . Thus, we set where the scaling amplitude ν translates the dominating collapse timescale (the duration of the recent collapse events) into the evolution time t of the coarsening process. Writing out σ D explicitly by Eq. (14), we then have the master curve that implicitly defines the coarsening law Λ (t). In the following, we use this relation to analyze two examples and compare with numerical simulations.
Example 1: Recruitment-detachment model. In the main text we discussed the recruitment-detachment model f ex (see Fig. 2b,c in the main text). We simulated the full reaction-diffusion system on a large domain (see series of snapshots in Fig. S2) and extract the time evolution of the average peak separation Λ (t) (black trajectories in Fig. S3a). Simulations were performed for average massesρ = 1.5, 3, and 6. To compare to the scaling arguments above we calculated the relation ∂ M η stat by numerical continuation of an elementary stationary pattern, * and fitted the scaling amplitudes ν and µ in Eq. (16) for each value ofρ separately (see red dashed line in Fig. S3a). After an initial transient, the simulations follow the predicted coarsening dynamics Eq. (15) (see red, dashed lines in Fig. S3). In particular, we find that the scaling amplitude µ ≈ 1.3 is identical for all three average masses. This affirms our assumption that the average peak mass and the average peak distance alone are the characteristic properties of the coarsening process. Finally, using the fitted scaling amplitudes, all three curves can be collapsed onto a single master curve (see Fig. S3b).
Example 2: Otsuji's 'Model II'. As a second example, we performed finite element simulations of 'Model II' from Otsuji et al. [2]; cf. Eq. (11). In this model, peaks can grow to arbitrary mass because the nullcline becomes asymptotically parallel to the FBS, such that there is no third FBS-NC intersection point at large density. Correspondingly, coarsening follows a power-law up to arbitrarily large length scales (and respectively peak masses). The coarsening exponent can be derived from η stat (M ) ∼ M −α as discussed in the main text. Notably, for this model the single-peak solution on the infinite domain is known analytically [2] ρ ∞ (x) =ρ(η stat ) sech Using this, one finds η stat ∝ 1/M which yields α = 1. Using the relation derived in the main text, this results in the coarsening law Λ (t) ∼ t 1/(2+α) = t 1/3 . Numerical simulations show excellent agreement with this analytically predicted coarsening law (see Fig. S4). * Analytic expressions for ∂ M ηstat are only available for idealized peaks and mesas, but not in the crossover regime; see Fig. 2b in the main text.

Coarsening of peaks is always uninterrupted
Consider a stable, symmetric peak with mass M 0 on a domain of size Λ, for which ∂ M η stat | M0 < 0 holds. Then, it follows that patterns constructed from such elementary peaks will initially coarsen due to the mass-competition instability. The question is whether coarsening always remains uninterrupted, i.e. whether ∂ M η stat remains negative as small peaks disappear and the mass of the remaining peaks increases. In other words, does ∂ M η stat < 0 continue to hold along the entire branch of stable stationary peaks, that is, for the family of solutions to Eq. (2), parameterized by η stat ? In the following we show that this is indeed the case, assuming far-separated peaks (Λ int ). Moreover, we demand that for the initial peak ∂ρη stat < 0 holds in addition to ∂ M η stat < 0. These two derivatives are closely related viaρ = ρ − + M/Λ. As in the scaling analysis in section 3, note that ρ − Λ is negligible compared to the amplitude of large peaks such that this additional, technical assumption is typically not restrictive.
For far-separated peaks, the properties of individual peaks are the same as for a single peak on the infinite line (up to corrections exponentially small in their distance Λ). Asymptotically far away from the peak position (x → ±∞), the pattern profile fulfills ∂ xũ (x) → 0. Therefore, the stationary profile equation (2) implies lim x→±∞ũ (x) = u − (η stat ). With these boundary conditions, the boundary value problem Eq. (2) has a unique solutionũ(x) for each given η stat . * As a consequence, there is a unique peak mass M = dx (ρ − ρ − ) for each given η stat . Hence, for sufficiently well-behaved reaction terms f , the function M (η stat ) is continuous, such that the slope ∂ ηstat M can only change sign in extrema. However, extrema of M (η stat ) are folds of its inverse η stat (M )-the curve we are ultimately interested in.
We argue that already a single (elementary) peak becomes unstable before reaching the fold along this curve. If the peak mass M were a control parameter, a fold in η stat (M ) would directly imply that an eigenvalue of the linearized dynamics close to the stationary single-peak pattern crosses through zero (and therefore a change from stability to instability of the single peak). However, in a closed system (with reflective or periodic boundary conditions) with size Λ, the conserved average massρ = ρ − +M/Λ is a control parameter. As we will see, a fold bifurcation in the curve η stat (ρ) occurs before the fold in η stat (M ). Taking the derivative of M = Λ(ρ − ρ − ) with respect to η stat yields ∂ ηstatρ = ∂ ηstat ρ − at the fold in M where ∂ ηstat M = 0. Lateral stability of the plateau mandates ∂ ηstat ρ − > 0 [1]. Therefore, one has ∂ ηstatρ > 0 at the fold in M . Since we demanded ∂ ηstatρ < 0 for the initial peak, we must have passed an extremal point inρ(η stat ) (i.e. fold in η stat (ρ)) before reaching the fold in M (η stat ). Sinceρ is a control parameter, and therefore the fold in η stat (ρ) a bifurcation where an eigenvalue crosses zero, this implies that the elementary peak becomes unstable before reaching an extremum in M (η stat ). As a result, ∂ M η stat < 0 will hold along the entire branch of stable peaks such that coarsening will never interrupt.

Mesa patterns
In this section, the goal is to derive the rate of the mass-competition instability for mesa patterns. To this end, the equivalent of ∂ M η stat for peaks has to be constructed for mesas. Subsequently, we discuss the resulting growth rates of the mass-competition modes that drive coarsening and show that coarsening is generic for mesas.

The elementary mesa pattern
Consider the elementary pattern consisting of one interface connecting a low-and a high-density plateau on the infinite line (green profile in Fig. S5a, green line connecting u − and u + in Fig. S5b), determined by equation Eq. (2). The plateau densities are given by the FBS-NC intersection points u ± := u ± (η ∞ stat ) = lim x→±∞ũ ∞ (x), such that total turnover balance, Eq. (4), reads which determines the FBS-offset η ∞ stat . Linearization of Eq. (2) around the plateau densities for the stationary pattern, withf u | u± := ∂ uf (u ± , η ∞ stat ). This shows that the profile approaches the plateau densities u ± through exponential tails for x → ±∞ as: where 2 ± = −D u /f u | u± is the length scale of exponential relaxation to the plateaus (note thatf u | u± < 0 because the plateaus are laterally stable [1]).
For example, for the Brusselator core (Eq. (12) with ε = 0) the stationary singleinterface solution with one interface at x = 0 on the infinite line is given bỹ where η ∞ stat = 3 d/2 and = ± = √ D u (recall that we set the timescale to unity; introducing a timescale k −1 explicitly via the reaction term, f → kf , yields ± = D u /k). From this, one obtains the amplitudes a ± = 2/d. Consider now a stationary, single-interface pattern on a finite-sized domain with no-flux boundary conditions at distances L ± from the inflection point respectively (see Fig. S5). This elementary pattern is the basic building block from which periodic patterns with wavelength Λ = 2(L − + L + ) are constructed. Neglecting contributions from the interface (sharp interface approximation), the average mass is given by Solving for L ± yields Our goal is now to find the mass-redistribution potential η stat (L − , L + ) on the finite domain. To this end, we again use total turnover balance Eq. (4), which now reads because on the length scales L ± the pattern has not fully approached the plateau densities u ± . To find the boundary concentrationsũ(L ± ), we use the linearized profile equation (19) again-this time with no-flux boundary conditions at L ± -which yields for the tails of the pattern in the upper ("+") and lower ("−") plateau, respectively. The amplitudes δu ± describe the deviation from the plateaus at the boundaries,ũ(L ± ) = u ± ∓ δu ± , and can be found by asymptotic matching to the tails on the infinite domain Eq. (20): Substituting expression (26) into Eq. (24) and approximating to lowest order in δu ± and δη stat = η stat (L − , L + ) − η ∞ stat we find The three terms of this equation have a simple graphical interpretation. In phase space, total turnover balance corresponds to a balance of the areas enclosed by the FBS and NC (see Fig. S5, red-shaded area). The first term of Eq. (27) then corresponds to the right purple triangle and the second one to the left purple triangle whose areas grow and shrink with δu ± respectively. The third term corresponds to the green stripe whose area is proportional to δη stat . Solving Eq. (27) for δη stat and taking the derivative with respect to L ± yields This describes the change in the stationary mass-redistribution potential if the length of the high-or low-density plateau is changed.

Two coarsening modes for mesa patterns in 1D
The process underlying the coarsening of mesas is the shifting of their interfaces. There are two distinct elementary modes of coarsening that differ in how neighboring In the specific setting of a pattern with two interfaces in a box with no-flux boundary conditions there are two distinct stationary patterns, see dark blue profiles in Fig. S6 where the vertical gray lines indicate the reflective boundaries. For each of these patterns, only one of the two coarsening modes is compatible with mass conservation, respectively; namely the one that changes the lengths of the outer plateaus and shifts the inner plateau. This allows one to study the two modes separately.
(i ) When the high-densities plateaus are located at the reflective (no-flux) boundaries, separated by a trough in the center (Fig. S6a), the coarsening shifts the trough and changes the length of the two high-density plateaus each containing mass M ≈ (ρ + − ρ − )L + . Therefore, the change in η stat driving the coarsening in this situation follows from the change of the high-density plateau alone, as the width of the trough does not change. Using the chain rule to translate the mass change into a length change of the high-density plateau, we define the derivative that determines the coarsening rate for the mode changing the width of mesas.
(ii ) In the opposite case of low-density outer plateaus (Fig. S6b), the high-density (inner) plateau shifts as a whole, and only the low-density (outer) plateaus change in length. The corresponding change of the mass-redistribution potential at the interfaces is given by which determines the coarsening rate of the mode changing only the width of troughs.
We are now in a position to write down the growth rates σ ± D of these two modes of coarsening in the diffusion-limited regime. By the same reasoning that lead us to Eq. (3) in the main text, i.e. that the gradient in the mass-redistribution potential η is linear between the interfaces, we obtain (Recall that the plateaus' half-lengths, L ± (Λ), are determined by the average massρ via Eq. (23).) For a pattern with two interfaces, Eq. (30) highlights that the "interaction" or "overlap" of the inner exponential tails forming the trough between the two high-density plateaus does not contribute to coarsening in the leading order. Instead, the cause for coarsening is the gradient in η caused by the length change of the outer exponential tails at the no-flux boundaries of the box.
For a pattern with more than two interfaces, as illustrated by the light blue lines beyond the vertical dashed lines in Fig. S6, both coarsening modes contribute simultaneously. However, the growth rates of the two modes are typically vastly different because the upper and lower plateaus will, in general, not have (on average) similar effective lengths L ± / ± and the corresponding rates are exponentially small in these effective plateau lengths L ± / ± . Only ifρ is tuned to the "symmetric" situation L − / − ≈ L + / + , i.e. ifρ ≈ ( − ρ − + + ρ + )/( − + + ), both modes have similar growth rates and contribute equally to coarsening.
Example: Brusselator core. For the Brusselator core, one hasf u u± = −1 and With ρ + − ρ − = (1 − d) 2/d, one then finds Substituted into Eq. (30) yields for the growth rates of the two coarsening modes. To verify this analytic result, we performed linear stability analysis numerically for stationary patterns with two interfaces on a domain with no-flux boundary conditions (see Fig. S7). Varying the average densityρ changes the relative length of the upper and lower plateau. This increases the rate for the one mode while the other rate decreases.

Mesa coarsening is uninterrupted
To show that the coarsening is uninterrupted for mesas in general we need to show that ∂ ± M η stat < 0. From Eq. (28) and using the chain rule ∂ ± M η stat = ±(ρ + − ρ − ) ∂ L± η stat , it follows that ∂ ± M η stat < 0 holds iff u | u± /F η < 0. Moreover, the stability of plateaus mandatesf u u± < 0 [1]. Thus, we are left to show that F η > 0. We will argue that this condition is necessary for elementary (single-interface) patterns to be stable.
Consider a perturbation to the mass-redistribution potential δη. In the interface region, δη will quickly homogenize by fast diffusion (cf. Eq. (5b)). The dynamics of δη int close to the interface is then driven by the imbalance of reactive flows: * Hence, if F η < 0, the perturbation δη int will grow, destabilizing the interface. Conversely, stability of the interface requires F η > 0. A a mathematical proof for this stability condition in the limit D u D v can be found in Ref. [6] (Lemma 3.3). Taken together, the above arguments show that ∂ ± M η stat < 0, i.e. coarsening is generic for all mesa patterns.

Mesas in 2D
The Laplace operator in polar coordinates is ∇ 2 = ∂ 2 r + 1 r ∂ r +∂ 2 θ . Hence, for a spherical droplet of density u + in a surrounding density u − < u + total turnover balance reads For R int (sharp interface approximation), we can approximateũ(r) ≈ũ 1D (r − R), whereũ 1D (x) is the stationary solution on the infinite line with one interface at x = 0. With η stat = η ∞ stat + δη stat , we then obtain The resulting proportionality δη stat ∝ R −1 , is identical to functional dependence of the chemical potential on the droplet radius that drives Ostwald ripening. Therefore, the results from Lifshitz-Slyozov-Wagner (LSW) theory [11,12] also hold for 2cMcRD systems in two and more spatial dimensions. In particular, mesa-droplets will follow t 1/3 coarsening. Peaks, in contrast,-which exhibit power-law coarsening with a system-dependent exponent already in 1D-will have system dependent coarsening exponents in higher spatial dimensions as well. (Details will be presented in a forthcoming publication.)

Reaction-limited regime
So far, we have considered the case where the diffusive transport of mass between peaks/interfaces is the limiting process, such that η is determined by total turnover balance in quasi-steady state on the scale of the interface width int . This is no longer the case when the diffusive flux of mass between peaks/interfaces ∼ D v /Λ is faster than the reactive flows locally driving growth and shrinking of peaks (resp. interface movement of mesas). In this reaction-limited regime, v(x, t) (and hence η, since η = v for D u /D v → 0) can be approximated as globally uniform such that Eq. (1) maps to a generalized form of the conserved Allen-Cahn equation [13] ∂ t u( This type of dynamics has previously been used as as phenomenological model for phase separation in granular media [14] and for intracellular pattern formation wherē v(t) is the protein concentration in an effectively well-mixed cytosol [15][16][17]. It is worth emphasizing that whether such an approximation is justified cannot be decided a priori as it depends on the relevant length-and time-scales in the phenomenon under investigation. This is often not trivial and, in general, requires to study the full dynamics, Eq. (1), first. For coarsening, the crossover from the diffusion-to the reaction-limited regime takes place at D v /Λ ∼ f η int int as we show in the following. Other examples where cytosolic concentration gradients play an important role, such that a well-mixed assumption for the cytosol is not warranted, are geometry sensing [18,19] and oscillatory patterns in cells and in vitro [20][21][22]. How can we understand coarsening in the reaction-limited regime? Above we showed that for stationary peaks/mesas of mass M the mass-redistribution potential is determined by total turnover balance which determines the relation η stat (M ). Moreover, we argued that η stat (M ) is a strictly decreasing function. Therefore, a spatially uniform η =η cannot simultaneously balance total turnover for mesas of different width (resp. peaks of different height). Hence, there will be net (residual) reactive flows leading to movement of the mesa-interfaces (resp. growth/shrinking of peaks) in the u-profile. As a consequence, the timescale of reactive conversion between u and v determines the growth rate σ R of the competition instability while mass redistribution throughv(t) is instantaneous.
As in the diffusion-limited case, we study the mass-competition instability in the elementary scenario of two (half-)mesas on the domain [−Λ/2, Λ/2] with no-flux boundary conditions. Let us start with the observation that since such a stationary state is, by construction, symmetric around the origin, all eigenmodes of the linearized dynamics are even and odd functions. The mode corresponding to the mass-competition instability, must be an odd function as it transports mass from one interface to the other. It follows that for this mode, the spatially uniform perturbation to the massredistribution potential must vanish, δη(t) = 0. The dynamics of the mass-competition instability is therefore fully contained in the shifting positions of the interfaces, δx, as the mesas gain/lose mass.
As we argued above, the dynamics in the reaction-limited regime is determined by turnover imbalance. Therefore it is useful to introduce the turnover integral With this, total turnover balance (cf. Eq. (24)) which determinesη = η stat in the stationary state, reads F (L − , L + , η stat ) = 0. The shift of the interface positions δx, changes the length of the outer plateaus as the inner plateau shifts as a whole (see Fig. S6a). For specificity, let the outer plateau be a high-density plateau. Then the turnover imbalance at the right/left interface is For small δx, and hence small turnover imbalance, the ensuing dynamics can be approximated as a slow shift of the stationary interface profile, analogously to slowly propagating fronts in bistable media [10]. Using ∂ t u(x, t) ≈ −(∂ t δx) ∂ xũ (x), the propagation velocity of the right interface is given by where I = dx (∂ xũ ) 2 . With the ansatz δx ∝ e σ + R t this yields to a growth rate σ + R ≈ I −1 ∂ L+ F for the mass-competition instability in the reaction-limited regime. (As before, the superindex "+" indicates that the upper plateau(s) change in width while the lower plateau(s) shift.) We now use Eq. (27) to express ∂ L± F in terms of ∂ ± M η stat . Along the branch of stationary patterns, we have F L − , L + , η stat (L − , L + ) = 0 (cf. Eq. (27)). Taking the derivative with respect to L + along the branch yields and thus where we applied the definition of ∂ ± M , Eq. (29). Taken together, we obtain Let us briefly discuss the physics implied by this expression and compare with its analog in the diffusion-limited case. Recall that the factor ∂ ± M η stat also appears in σ ± D in the diffusion-limited regime (see Eq. (30)). There, we showed that gradients in η between peaks/mesas with different mass M drive coarsening, while each peak/mesa is in quasi-steady state. In contrast, in the reaction-limited regime, the factor F η ∂ ± M η stat represents the imbalance in turnover caused by perturbations to the peak/mesa mass M that drive η stat (M ) away from the globally uniformη.
To obtain the expression given in the main text, we make two further approximations in Eq. (43). First, the slope ∂ xũ in the interface can be estimated as (u Second, we approximate f η in the interval [u − , u + ] by its average over the interface f η ≈ f η int which yields F η ≈ (u + − u − ) f η int . Together, these approximations yield the expression given in the main text: Comparing to the growth rate σ ± D in the diffusion-limited regime, we find the crossover at D v /Λ ∼ f η int int . Further using the approximation int ≈ D u / f u int for the interface width [1], shows that the crossover depends on both diffusion constants, the wavelength Λ, and the reactive timescales in a non-trivial way.
Notably in his analysis of Ostwald ripening [12], Wagner also considered a "reactionlimited" case, in addition to the diffusion-limited process typically associated with Ostwald ripening. Wagner considers droplets that grow by attachment of molecules from the surrounding bulk solution. He postulates slow attachment at the droplet surface with a "Stoffübergangsgeschwindigkeit" (mass conversion rate) k. In our 2cM-cRD system, the droplet is formed in u while the faster diffusing v represents the bulk solution. Hence, the net rate of conversion between v and u in the interface region (ρ + − ρ − )F η /I ≈ int f v int corresponds to the rate k in Wagner's work. By this correspondence, Wagner's results directly carry over to 2cMcRD systems in the reaction-limited regime and can be used to study coarsening of these systems in higher spatial dimensions.
Example: cubic model. The stationary pattern on the infinite line is given bỹ In the diffusion-limited regime, Eq. (30) yields In the reaction-limited regime, Eq. (43) yields The crossover σ D = σ R takes place at D v /D u = 3L ∓ / ± . In Fig. S8, we compare these analytic results with numerical linear stability analysis based on a finite difference discretization of the spatial domain.

Weakly broken mass conservation
In the presence of a weak source term, the time evolution of the total density is given by (cf. Eq. (6) in the main text) To find the stationary, single-interface pattern on domain [0, Λ/2] (see Fig. S9), we make a perturbation ansatz in ε: whereρ 0 (x; L − , L + ) denotes the stationary pattern of the mass-conserving system, with plateaus of lengths L ± , where L + + L − = Λ/2.
Recall that in the strictly mass-conserving system, L ± were determined by the conserved average massρ. In the presence of a source, the average mass is no longer a control parameter. Instead, the source has to be balanced in steady state, that is, the integral of Eq. (47) over the entire domain has to vanish. Using the local equilibrium approximation s ≈ s * (ρ) introduced in the main text, and neglecting contributions from the interface (sharp interface approximation), source balance mandates that Deviations from the plateaus ρ ± due to the source terms are of order O(ε) and therefore do not contribute to L ± to lowest order: In summary, to lowest order in ε, the effect of the source is to single out an interface position (that is, the relative lengths of the high-and low-density plateaus), that balances the source. Equation (50) has a solution with L ± > 0 only if s * (ρ ± ) has opposite signs for ρ ± . Moreover, stability of an elementary (single-interface) pattern requires that s * (ρ + ) < 0 (and thus s * (ρ − ) > 0). To see why, suppose s * (ρ + ) > 0 instead. Then a perturbation increasing the length of the upper plateau (and respectively shortening the lower plateau) will result in net production, leading to an increase of the average density and thereby a further increase of the upper plateau length. This feedback loop destabilizes elementary patterns. Since all stationary patterns can be decomposed into such elementary patterns, no stable patterns exist if s * (ρ + ) > 0. Solving Eq. (50) together with L − + L + = Λ/2, one finds where we defined the relative plateau lengths ξ ± , which fulfill ξ − + ξ + = 1.
As a concrete example, for the Brusselator, Eq. (12), one has Figure S9. Illustration of a stationary, single-mesa pattern in the presence of a weak source term. (Left) Stationary profilesρε(x) andηε(x), with the respective high-and low-density plateaus denoted by superindices ±. The stationary pattern for ε = 0 is shown as a dashed line. The reflected profile is shown, such that the 'center' of the high density plateau is located at x = 0.

Mesa splitting
To numerically investigate mesa splitting, we performed finite element simulations where the source rate ε is increased adiabatically (see also Refs. [5,23]). For specificity, we discuss the splitting of a high-density plateau, below. Splitting of a low-density plateau, i.e. 'insertion' of a new mesa in a trough, can be analyzed analogously.
To gain some intuition into the mesa-splitting process, we first performed finite element simulations, starting from a fully phase separated state, slowly increasing ε (see Movie 2 and Fig. S9). The interface remains in a flux balance subspace, while the plateaus are slaved to the nullcline η * (ρ) along the two laterally stable branches (∂ ρ η * > 0). As the source strength increases, the minimum (resp. maximum) of the upper (resp. lower) plateau approach ρ hss . Eventually, one of the plateaus enters the laterally unstable regime where the nullcline slope is negative (∂ ρ η * < 0). This triggers a regional lateral instability in the region in space centered around the respective extremum of the curved plateau. This regional instability initiates the formation of a dip (resp. peak) that continues to grow until it reaches the stable nullcline branch where it saturates and forms a new plateau. Finally the two "daughter mesas" slowly move until a new stationary state is reached.
In addition to the finite element simulation, we performed numerical continuation of the stationary mesa patterns for the parameter ε. Continuation reveals that the branch of stationary patterns undergoes a fold bifurcation as the mesa splits (see Fig. S10). In the vicinity of the fold point, the formation of the 'dimple' in the center of the high-density plateau follows a 'universal' evolution that was previously analyzed in Ref. [5]. Here, we only estimate the critical source strength ε split , using the intuition that splitting is driven by a regional lateral instability. For specificity, we derive the splitting condition for the high-density plateau. Analogous arguments can be used for the low-density plateau.
As explained above (and in the main text), splitting is triggered if the density profile of the upper plateau,ρ + ε (x), drops below ρ + lat . Since the plateau is slaved to the nullcline η * (ρ), this is equivalent toη + ε (x) dropping below η * (ρ + lat ). (Actually, in the vicinity of ρ + lat , the plateau density starts to deviate from the nullcline. This becomes relevant when one studies the splitting process in more detail (see e.g. Ref. [5]) but can be neglected for the lowest order approximation we attempt here.) To find the stationary η-profile in the high-density plateau,η + ε (x; ε, Λ), we need to solve Eq. (47) in steady state. Using local equilibrium approximation we find, to lowest order in ε, To first order in ε, this is solved by a quadratic profilẽ where η * (ρ + ) = η ∞ stat is the flux-balance subspace position of the interface (see Fig. S9).
The critical source strength for splitting, ε split (Λ), is determined by the conditioñ η + ε (0; ε + split , Λ) = η * (ρ + lat ), which yields . (55a) Following analogous steps for the low density plateau, one finds (55b) Figure 3d in the main text shows the curve ε + split (Λ) (green line) together with the location of fold points (black squares) obtained by numerical continuation (cf. Fig. S10). We find good agreement for large Λ (resp. small ε), while a deviation becomes obvious for small Λ. The main reasons for this are that corrections to the sharp interface approximation (of order ∼ int /Λ) and higher order terms in Eq. (53) become important as the critical value ε ± split (Λ) increases. Better approximations could be obtained by taking these higher order terms into account.
As we briefly noted in the main text, splitting of droplets in two or more spatial dimensions can be driven by additional shape instabilities [23][24][25] that are not present in 1D. The mechanism driving these shape instabilities is potentially different from the regional lateral instability identified here. Generalizing the phase-space analysis presented here to account for shape instabilities is a promising avenue for future research. Droplet splitting is particularly interesting in the context of membraneless organelles that form by protein phase separation [26,27].

Interrupted coarsening
Let us now turn to interrupted coarsening. To that end, we again study the stability of a the stationary pattern on the domain [−Λ/2, Λ/2], obtained by reflecting an elementary (single-interface) mesa at x = 0. For specificity, we present the details for the case that the outer plateaus are at high density ρ + .
As we argued in the analysis of mesa coarsening above (Sec. 5.2), a perturbation transferring a small amount of mass δM from one side to the other (M R,L = M 0 ±δM ), results in a shift of the interfaces by δx ≈ δM/(ρ + − ρ − ). To find the dynamics of δM , we integrate Eq. (47) over the 'half box' [−δx, Λ/2], where the inner boundary is shifted by −δx to account for the fact that the inner plateau shifts as a whole and only the outer plateaus change in length. In the source term, u =ũ+δu and v =ṽ+δv denote the concentration profiles perturbed by the mass transfer δM from the left to the right mesa. As we argued above, this perturbation results in a shift of both interfaces by −δx, such that δu ≈ −(∂ xũ ) δx, where ∂ xũ is the "translation" (or Goldstone) mode; analogously δv = −(∂ xṽ ) δx. To approximate the source term, we use the local equilibrium approximation s ≈ s * (ρ) and neglect contributions from the interface (sharp interface approximation, cf. Eq. (49)): To approximate the flux term, we use the growth rate of mass-competition instability from the strictly mass conserving situation: ∂ This means that, to lowest order in ε, the growth rate of the mass-competition mode is simply shifted by ε s * (ρ + )/(ρ + − ρ − ). Recall that s * (ρ + ) is negative, so the source has a stabilizing effect. The growth rate vanishes for wavelengths larger than Λ + stop given by By an analogous calculation for the case that the outer plateau is at ρ − , one obtains Recall that s * (ρ + ) < 0 and s * (ρ − ) > 0, so both equations above can be combined to Eq. (7) in the main text. In a large system, containing many mesas, coarsening stops once both modes are stable, i.e. if the wavelength Λ exceeds max(Λ − stop , Λ + stop ). Finally, using Eqs. (28) and (30) we can write equations (59a,b) as where ξ ± , the relative width of the plateaus, are determined by the source term via Eq. (51) and all other quantities are determined by the single-interface pattern of the strictly mass-conserving system (ε = 0) on the infinite line (cf. Sec. 5.1). Implementing this for the Brusselator, yields where = √ D u and the relative plateau lengths ξ ± are given by Eq. (52). When ξ + > ξ − , then Λ − stop > Λ + stop and vice versa. In Fig. 3d in the main text, we plot the analytically obtained Λ stop (ε), and find excellent agreement with the critical wavelength (where the dominant eigenvalue vanishes) obtained from numerical linear stability analysis. Kolokolnikov et al. (2006). In Ref. [4], the Brusselator in the limit of slow production/degradation was analyzed using a mathematical approach that is conceptually very different from our analysis. This warrants a direct comparison of the results.

Comparison to the results by Kolokolnikov et al.
Casting the stability criterion Eq. (5.3) from Ref. [4] in our notation yields In contrast to our criterion, this does not explicitly distinguish the two different competition modes (see Fig. S6). Notably, for plateaus of equal relative width, ξ + = ξ − = 1/2, Eq. (62) and our criteria Eq. (61) agree exactly. For plateaus of unequal width, ξ + = ξ − , the two competition modes become stable above different wavelengths Λ ± stop , and the pattern is stable (i.e. coarsening stops) if Λ > max(Λ − stop , Λ + stop ). This inequality is approximated by the criterion Eq. (62) from Ref. [4] (see Fig. S11). This hints at a deeper relationship between the different approaches used to obtain the stability conditions. The mathematical reasoning used in Ref. [4] has recently been generalized to a mass-conserving system with three components [23]. This might be a good starting point for future generalizations of the framework presented here.
McKay & . In Ref. [6], the stability of mesa-type solutions is analyzed for general two-component reaction-diffusion systems in the limit D u 1, D v 1 and with 0 < τ O(D −1/2 u ). The systems with (weakly) broken mass conservation that we analyzed in the present manuscript (see Eq. (6) in the main text) can be cast in a closely related form ∂ t u = D u ∂ 2 x u +f (u, η) + εs 1 (u, η) (64a) wheres(u, η) = s(u, η−du). Importantly, we expect that the additional term (1−d)∂ t u is negligible during the slow coarsening dynamics in the diffusion limited regime. Thus, via the relations w ↔ η,f ↔ f , εs ↔ g, and setting τ = 1 we can map the results from Ref.
[6] to our system. Notably, the stability criterion as stated in "Principal Result 3.1" in Ref. [6] is identical to our result Eq. (60) under this mapping. This is striking, because the approximations in Ref. [4] require D v 1, while we derived Eq. (60) in the diffusion limited regime, which requires D v not too large (see Sec. 6). In other words, our criterion generalizes the previous result to the diffusion limited regime. However, the growth rates of the instability away from the critical point do not agree (Eq. (56) in Ref. [6] vs Eq. (58) here). This discrepancy, and the surprising agreement of the stability criteria merit a detailed investigation in future work. Kolokolnikov et al. [4] +-mode (analytic) -mode (analytic) +-mode (numeric) -mode (numeric) Figure S11. We set the production rate p = ξ+ 2/d (cf. Eq. (52)), which allows us to control the relative plateau lengths directly by tuning ξ+ as a control parameter. The red, dashed line shows the solution to Eq. (62) from Ref. [4]. The blue and green, solid lines show the solutions to Eq. (61) for the two competition modes. The mode is stable above and unstable below the respective curve. Symbols indicate the curves along which the numerically determined stationary pattern changes stability (numerical LSA using finite element discretization). Parameters: Du = 1, Dv = 10, ε = 10 −4 .