Surface-tension-driven coarsening in mass-conserved reaction-diffusion systems

Mass conservation in chemical species appears in a broad class of reaction-diffusion systems (RDs) and is known to bring about coarsening of the pattern in chemical concentration. Recent theoretical studies on RDs with mass conservation (MCRDs) reported that the interfacial curvature between two states contributes to the coarsening process, reminiscent of phase separation phenomena. However, since MCRDs do not presuppose a variational principle, it is largely unknown whether description of surface tension is operative or not. In this study, we numerically and theoretically explore the coarsening process of patterns in MCRDs in two and three dimensions. We identify the parameter regions where the homogeneous steady state becomes stable, unstable, and metastable. In the unstable region, pattern formation is triggered by usual Turing instability, whereas in the metastable region, nucleation-growth-type pattern formation is observed. In the later stage, spherical droplet patterns are observed in both regions, where they obey a relation similar to the Young-Laplace law and coarsen following the evaporation-condensation mechanism. These results demonstrate that in the presence of a conserved variable, a physical quantity similar to surface tension is relevant to MCRDs, which provides new insight into molecular self-assembly driven by chemical reactions.


I. INTRODUCTION
Reaction-diffusion systems (RDs) are one of the most generic mathematical frameworks that give rise to spatiotemporal patterns and have been applied to various phenomena, ranging from physics, chemistry, biology, geology, to ecology [1][2][3][4]. The Turing mechanism explains how a homogenous initial state is destabilized to develop a spatial structure; infinitesimally small fluctuations with a particular wave number q m selectively grows exponentially. Notably, despite this analysis being valid only for every early linear time regime, the length scale obtained from this analysis m (≡ 2π/q m ) usually provides a good estimation of the characteristic length scale of final steady-state patterns, such as stripes or hexagonal dots in two-dimensional (2D) systems.
Mass-conserved reaction-diffusion systems (MCRDs) do not follow the above trend. MCRDs were originally developed to model molecular localization of membranebounded proteins inside cells, specifically Rho-family GTPases, which regulate cell polarity [5][6][7][8][9][10][11][12][13][14] (see also a recent review article Ref. [15]). In such systems, it was found that the characteristic length scale of the patterns easily grows beyond m after a linear time regime and often reaches a length scale comparable to the system size L. In the context of biomolecular self-assembly, this feature indicates that proteins accumulate at one site, not spreading over multiple sites, and is crucial to determining the unique directionality of cells. It was also discussed that the resulting scalability of the pattern against * c-tateno@g.ecc.u-tokyo.ac.jp the system size has biological significance in morphogenesis [16][17][18][19][20][21].
Pattern formation of MCRDs in one dimension has been intensively investigated in previous research. Two types of MCRDs were mainly studied: in the first case, referred to as the "Turing type" [11], multiple peak-like domains, illustrated in the right panel of Fig. 1a, appear from a homogeneous initial state via the Turing mechanism [5,6]; the second case is the "wave-pinning type" [7], where the system is bistable with high and low chemical density states, and by exposing a finite (not infinitesimally small) perturbation to the homogeneous state, mesa-like concentration profiles emerge as illustrated in the left panel of Fig. 1a. Both systems exhibit coarsening, i.e., smaller peaks and mesas shrink and disappear while larger ones grow. For the Turing type, both the height and width of the peak increase, whereas for the wavepinning type, only the width of the mesa grows while keeping the height (see Fig. 1a). For a long time scale, the system reaches a single isolated domain in the Turing type. On the other side, for the wave-pinning type, the decrease in the number of domains significantly slows down at some point, resulting in the apparent coexistence of multiple domains.
Despite the difference between the Turing and wavepinning types mentioned above, recent studies revealed that these differences would be rather superficial [11,12,14]. Chiou et al. [11] found that the peak-like profile in Turing type MCRDs is owing to the unsaturation of the peak height. That is, the peak-like pattern is transient and it will eventually become mesa-like for a sufficiently large system size and long simulation time. Their study also discussed that the difference in the peak heights arXiv:2010.03900v1 [cond-mat.stat-mech] 8 Oct 2020 causes a difference in "recruitment power", which enhances the flux of the chemical component from lower peaks to higher ones (see ∇p introduced in Sec. II). This drives faster coarsening of peaks seen in the Turing type than the wave-pinning type. Furthermore, Turing instability can occur in the wave-pinning type MCRDs by choosing appropriate parameter values [10][11][12]. Taken together, it is understood that there is quantitative, but no essential difference between long term dynamics of two cases. We note that the coexistence of multiple peaks can also be induced by a slight violation of mass conservation [12,14] or the introduction of negative feedback [12].
However, since earlier studies have conducted intensive analyses only in 1D systems, there remains a fundamental untested component in understanding the dynamics of MCRDs. For higher dimensions, the interfacial morphology of patterns may be crucial to the pattern formation dynamics of MCRDs. Indeed, Refs. [11,12,14] pointed out that the curvature of the interface affects the coarsening of the patterns in the late stage. This observation raises a possible coarsening mechanism governed by "surface tension" in MCRDs, similar to the phase separation phenomena [22]. However, in MCRDs the chemical reaction violates the detailed balance condition in general, and the existence of the variational functional, which is prerequisite for the theory of phase separation, is not obvious. Thus, it is unclear whether the physical description of the surface tension could be rationalized.
Motivated by the question raised above, this study explores the influence of curved interface on pattern formation dynamics in MCRDs. We perform extensive numerical simulations using multiple GPUs for both 2D and 3D systems, and the obtained dynamics are analyzed using methods developed for phase separation phenomena, together with technique of dynamical systems theory. By investigating the similarities and differences between MCRDs and phase separation phenomena, we explore whether mathematical structures such as energy variational principles are also inherent in MCRDs.

II. SYSTEM DESCRIPTION
We consider a system composed of two chemical species, U and V , and denote their concentration fields as u(t, r) and v(t, r), respectively. These species are interconvertible by chemical reactions U → V and V → U , the rates of which are k u u and k v v, respectively. Note that these two reactions do not have to be in reverse relationship and the system may not satisfy detailed balance. Both the transition rates k u and k v can be dependent on the density states u and v, as a result of positive feedback regulation in the chemical reaction network. Moreover, both U and V migrate over space by thermal diffusion. We model the time-evolution equations for u and v as the following reaction-diffusion equations: Here, D u and D v are the diffusion coefficients of U and V , respectively, and we assume D u < D v . f is the reaction term and written as f = k u u − k v v. By adding the above two equations, we obtain where s = u + v and p = D u u + D v v. Eq. 3 indicates that s is a conserved variable without material sources. The relevance of adopting MCRDs to describe cellular processes can be rationalized as follows: (i) Many subcellular processes, such as cell polarity formation, are sufficiently rapid (∼1 min) compared with the production and degradation rates of the molecules (>10 min). Such processes are usually conducted only through the activation and deactivation of molecules, such as phosphorylation and dephosphorylation. The ratio of molecules between "on" and "off" states changes; however, the total mass of the molecule is conserved throughout the processes. (ii) Molecular interaction contains positive feedback loops that regulate the switching rate between molecular states, e.g., the acceleration of phosphorylation is found in subcellular processes. (iii) Molecular states are often associated with an affinity to the cellular membrane; thus, molecules in different states are likely to attach to or detach from the membrane. As a result, molecular diffusion coefficients are significantly different between the molecular states, being small on the cell membrane but large in the cytosol. As these processes are basic, many previous models, explicitly or implicitly, contain chemical components that are conserved through dynamics. In Ref. [5], it was shown that in a model with these conditions, a uniform distribution of molecules can be destabilized in a Turing-like manner, leading to an accumulation of molecules at one site.
To examine the generality of the results presented below, we employ two models, following earlier works on MCRDs, both of which violate the detailed balance condition. Model I employs k u = ca/(b + u 2 ) and k v = c [5,6]. For Model II, k u = α and k v = k 0 +βu 2 /(K 2 +u 2 )) are employed [7]. The former and latter models correspond to Turing type and wave-pinning type, respectively, by the criteria mentioned in the Introduction. Because these models share many aspects in their dynamical trends, we primarily describe the results of Model I in the main text, and make comments when there are remarkable differences in the results between the two models. The results from Model II are provided in Appendix D.
In the results presented below, the statistical data are based on simulations performed on 4096 2 and 1024 3 grids for 2D and 3D systems, and have normalized units, described in Appendix B. After normalization the governing equations Eq. 1 and 2 have two free parameters: one FIG. 1. Steady state analysis of MCRD systems. a. Two typical density profiles in MCRDs. u1 and u2 denote u at the bistable points. Depending on whether u in the domain is saturated to the stable points or not, the density profile of u show "mesa-like" (wave-pinning type, left panel) or "peak-like" (Turing type, right panel) behavior. b. Dispersion relation in the homogeneous state. The largest growth rate λ is shown for various parameter values indicated by the labels corresponding to those in panel e. For large scale (q → 0), λ follows λ = −(Dvfu − Dufv)/(fu − fv)q 2 + O(q 4 ), and the sign of the prefactor of q 2 determines the stability of the homogeneous state. c. Mapping of the steady stable points on the u-v space for two different parameter values of T . The gray curve is the nullcline for the reaction term f . The blue and green lines represent p∞ = Duu + Dvv for T = 0.02 and 0.002, respectively. The blue and green filled circles represent the stable fixed points. Note that f is independent of the choice of the model parameters in our normalization (see Appendix B). d. Function profiles of Ψ for T = 0.02 (blue) and 0.002 (green) at p = p∞. e. Phase diagram. The blue curve represents the border of the unstable region where the Turing mechanism is operative. The black and gray curves represent the set of stable points, u1 and u2, respectively. We denote the region outside of the two curves as "stable" and regions that are neither stable nor unstable as "metastable".
is associated with the reaction term and is defined as R = a/b and α/β for Model I and II, respectively; the other is D = D v /D u , associated with the diffusion constants. As we will show below, the stability of the homogeneous state is determined by the product of these parameters T −1 ≡ RD, and the initial chemical concentrations that determine the conserved quantity in the system.

A. Phase diagram of MCRDs
Before investigating pattern formation dynamics in MCRDs, we begin with identifying their parameter dependencies by analyzing the initial homogeneous and eventual bistable solutions.
We first consider the stability of the homogeneous steady solution (u 0 , v 0 ). u 0 and v 0 are determined by the conditions f (u 0 , v 0 ) = 0 and s 0 = u 0 + v 0 , where s 0 ≡ sdr/V is the mean molecular concentration conserved through the dynamics (see Eq. 3), and V is the volume of the system. We suppose that the fixed point (u 0 , v 0 ) is linearly stable when the diffusion terms are absent in Eqs. 1 and 2 (otherwise the system runs to infinity). This is conditioned by f u − f v > 0, where f u and f v represent ∂f /∂u and ∂f /∂v, respectively. In the presence of the diffusion terms, the homogeneous state becomes linearly unstable against an infinitesimally small perturbation with wave number q between 0 < q 2 < is satisfied. This condition corresponds to the Turing instability. Since the denominator is positive as mentioned above, Turing instability condition is simply written as For this condition to hold for Model I, T < T c ≡ 1/8 is necessary (see Appendix B). The dispersion relation obtained in the analysis is shown in Fig. 1b, where λ(q) represents the growth rate of the perturbation with wave number q. Note that, owing to mass conservation, λ(q) always starts from zero in MCRDs unlike typical RDs. Next, we explore the dynamic behavior in the long time regime. The most significant feature in MCRDs is the existence of a conserved variable s, and the transport of s is driven by the spatial gradient of a scalar variable p as indicated by Eq. 3. At the steady state, p satisfies ∇ 2 p = 0. This implies that p becomes a constant value p ∞ in the absence of material source. Owing to the constant nature of p, (u, v) is constrained on a line p ∞ = D u u+D v v which has three intersection points with the nullcline f = 0 (see Fig. 1c). Two of them, indicated by (u 1 , v 1 ) and (u 2 , v 2 ) (u 1 < u 2 ) in the figure, are locally stable points, to which (u, v) tends to relax. Therefore, it is expected that the system behaves similar to a bistable system, where the bistability is induced by the diffusion rate asymmetry between components U and V .
We then consider a one-dimensional solution u(x, t) where two stable states u 1 and u 2 are connected via a single interface. The solution expected for the steady state is of the form u(x − ct); it either translates with a constant speed c while maintaining the concentration profile, or stays at the same positions (i.e., c = 0) [23]. According to earlier studies [5][6][7][8][9][10][11][12][13][14], the latter static solution is typical in MCRDs (see also Appendix C). We obtain this solution when is satisfied, where Ψ(u, p) is defined by Eq. 5 determines the value of p ∞ (see Appendix C for derivation of the above relation). u 1 and u 2 corresponds to stable solutions for dΨ(u, p ∞ )/du = 0. In higher dimensions, p ∞ obtained above is for a solution u(r) where two stable states are segregated by a flat interface. We will see that p takes larger value than p ∞ for curved interface, which is responsible for material transport among domains and dominates long-time coarsening dynamics. Examples of Ψ(u) at p = p ∞ are shown in Fig. 1d. The double-minimum functional shape of Ψ(u) seen in the figure, together with the functional form shown in Eq. 6, reminds us of energy functional in phase separation dynamics [24]. Based on this analogy, we can construct the phase diagram of the MCRD, as shown in Fig. 1e. In the figure, the horizontal and vertical axes represent the composition of u at the initial state and the dimensionless parameter T = 1/RD, respectively. The region enclosed by the blue curve is obtained by the condition corresponding to the spinodal line in the phase separation. Interestingly, this condition coincides with that for Turing instability shown in Eq. 4. Hereafter we refer to this region as the unstable region. The brown and black curves represent two values of u that minimize Ψ(u), obtained by solving dΨ/du = 0. These are interpreted as a binodal line. We call the region between the spinodal and the binodal lines as metastable region, where homogeneous states are stable in a deterministic sense; however, inclusion of the density fluctuation to MCRDs provokes pattern formation.
In the following sections, we will investigate the pattern formation dynamics of MCRDs in more detail, separately in the unstable and metastable regions. Through the analysis, we will test the analogy between MCRDs and phase separation phenomena. Before proceeding to the next section, we make two remarks: first, the topology of phase diagrams depends on the specific form of the reaction term [13]. In Appendix D, Fig. 5b, we show a phase diagram for Model II, which does not have a parabolic shape unlike in Fig. 1e. Second, a relation equivalent to Eq. 5 was derived in a recent paper by Brauns et al. [13] in which the reaction-term dependence of the phase diagrams is discussed in detail. Both their and our approaches do not require any approximations, and thus the resulting phase diagram (Fig. 1e) improve on those presented in previous studies (see, e.g., [10,12]).

B. Pattern formation in the unstable region
We first study the unstable region. To focus on similarities and differences between pattern formation in MCRDs and phase separation phenomena, in the following section, the time evolution of the conserved variable s is investigated. Fig 2a shows 2D pattern formation process of s for the parameter set A, indicated in Fig. 1e. An initial pattern with a characteristic length spontaneously appears everywhere in space (t = 50). Then, elongated and branched domains are relaxed to spherical shape (t = 200). Using linear stability analysis, the deviation of s from the initial homogeneous density s 0 , δs = s − s 0 , is expected to obey δs(q, t) ∼ δs(q, 0) exp(−λ(q)t), where q is the wave number. λ is the maximum growth rate of the perturbation with the wave number q (see also Fig. 1b), which can be estimated from the simulation re- is the structure factor at time t and the bracket represents the average over the angle. If R q (t 1 , t 2 ) is independent of time, R q is to be equivalent to λ(q). Fig 2b shows the growth rate measured by the simulation results, where R q at different times collapse onto λ(q) estimated by linear analysis (see red curve) for a short time. This indicates that the observed pattern formation seen is driven by the Turing instability.
In the later stage (t ≥ 400), spherical dropletlike domains coarsen over time. More specifically, the smaller droplets tend to shrink and eventually disappear, whereas the larger droplets tend to grow. This trend is reminiscent of the evaporation-condensation (Lifshitz-Slyozov-Wagner, or LSW) mechanism in phase separation phenomena as pointed out in Refs. [11,12,14]. To qualitatively characterize this growth behavior, we track how the number of droplets changes over time. Fig 2c shows the number density of the droplets n(t), where n shows a power-law decay for the long time regime. The estimated power-law exponent is approximately 0.60 (see red line in Fig 2c). Because the volume (area) occupied by a single droplet is roughly n −1 , the characteristic length scale at time t is estimated as = n −1/d (d being spatial dimension, i.e., d = 2), we accordingly obtain ∝ t −0.30 . This power-law exponent, often called the growth exponent, takes a different value depending on the coarsening mechanism. The value of 0.30 observed here is close to the growth exponent known for the LSW mechanism (1/3). We perform the same analysis for a 3D simulation at the same density (i.e., for parameter set A in Fig. 1e) and the obtained growth exponent is 0.29. Furthermore, the obtained prefactors of the coarsening law are 9.6 and 9.8 for 2D and 3D, respectively. These results are consistent with the fact that coarsening law by the LSW mechanism is independent of dimensionality for d ≥ 2. In phase separation phenomena, the power-law decay in the characteristic length is a consequence of a selfsimilar growth in the patterns. To examine whether this is the case for our problem, we perform a dynamic scaling analysis for S q . In Fig 2d, we show that g ≡ (t) −d S q (t) as a function of q. When a self-similar growth is the case, all the structure factors at different time should be mapped onto a unique mater curve. In this figure, we can observe that g for small q collapses onto a single powerlaw function after the formation of spherical droplets (t ≥ 400), indicating that the distance between centerof-mass positions of neighboring droplets are scaled by a unique length . On the other side, for large q, S q slowly shifts toward the higher-q side for intermediate time regime (e.g., t ∼ 3200), meaning that the density profile of droplets is not completely scaled by in this time regime. This is because s inside some droplets is not saturated to the stable fixed point s 2 . We mention that such breakdown of self-similarity originating from unsaturation is not specific to MCRDs, but is also observed in the phase separation phenomena when the shape of double-well potential is highly asymmetric (for example, see Refs. [25][26][27]).
For low-q, S q is proportional to q 4 , which is often observed in phase transition phenomena with a conserved order parameter (see, e.g., [28][29][30]). According to the literature, this trend appears in the case where surface tension plays a central role in the coarsening process, rather than thermal fluctuation [31] and phase separating structure is statistically isotropic [32].

C. Pattern formation in the metastable region
Next, we focus on pattern formation dynamics in the metastable region. We confirm from simulations that the homogeneous state is destabilized for parameter sets in the unstable region (A and B in Fig. 1e), but is stable in the metastable region (C and D), which is consistent with the Turing mechanism. However, by choosing appropriate non-homogeneous initial conditions (while the mean density is kept constant in the metastable region), we find cases where systems relax to a bistable spatial structure. To see a typical kinetics in the metastable region, in this section we investigate how patterns develop in the presence of density fluctuation. We incorporate intrinsic fluctuations originating from both diffusion and chemical reactions, referring to Ref. [33], see Appendix A for details. Time series of the maximum value of s, smax, for four independent simulations using parameter set D. The sudden increase in smax at random time indicates that the formation of the droplet is a stochastic activation process. The system size is chosen to be 162 2 , which is small enough to observe a single droplet formation in homogeneously mixed space, and large enough compared to the size of the nucleus. c. Nucleation rate under various initial conditions. d. Temporal changes in the number of droplets per unit volume n(t) for parameter set (u0 = 0.70, T = 0.002). The red line represents a power-law function with exponent 0.60.
down. This is because molecules in the region surrounding the droplets are almost exhausted, where s is in the vicinity of the minimum stable point. For a droplet to grow further, it is necessary to transport material from other droplets. Indeed, at the later stage (t = 940), the growth of the droplets is accompanied with shrinkage and eventual disappearance of small droplets. These pattern formation dynamics are similar to those in the nucleation-growth type phase separation of pure fluids, multi-component liquid mixtures, and alloys [34], but differs from those in so-called stochastic Turing patterns which are induced by stochastic resonance [33,[35][36][37][38][39].
As an indicator of the formation of the nucleus, we monitor the time evolution of a maximum density value of s max . The results obtained by four independent simulations are shown in Fig. 3b. We can see that s max initially fluctuates around a particular value and suddenly starts to increase. The onset time t inc clearly indicates the birth of nucleus. From the onset time, we compute the nucleation rate I via I = N/ t inc V , where N is the number of nuclei at time t inc , V is the volume (area), and the bracket represents the ensemble average. 16 independent simulations were performed to determine I under various initial conditions. The results are shown in Fig. 3c, where the vertical axis is displayed using a logarithmic scale. The nucleation rate I decreases significantly for a slight decrease in the initial density u 0 . Such a steep change in the nucleation rate for a small change of environment is widely observed in nucleation-type phase separation phenomena [34].
Although the driving force of pattern formation in the early stage is entirely different between the unstable and the metastable regions, we observe that the coarsening behavior in the late stage are very similar for these regions; small droplets evaporate while large ones grow. This similarity is quantitatively confirmed by plotting n(t), the temporal change in the number density of droplets (Fig. 3d). We find the same trend ∝ t −0.30 as in Fig. 2c, further supporting that the LSW mechanism can be applied to MCRDs.

D. Coarsening process in the late stage
Here, we investigate the evaporation-condensation (LSW) process observed in both the unstable and metastable regions in more detail. Since the time evolution of the conserved variable s proceeds following the gradient of p = D u u + D v v (see Eq. 3), the variable p is responsible for the transport of s among the droplets. Figs. 4a and b show an enlarged image of Fig. 2a at t = 3200 and the corresponding spatial map of p, respectively. We can see that smaller droplets have a higher p value. In this sense, we may call p as "chemical potential" of s, on the analogy of phase separation phenomena [22].
To verify the validity of this physical interpretation in a  Fig. 4a). When a section is located outside the droplets, we use a red cross symbol ('background' in the figure). We determine whether a section is inside or outside of the droplets by whether s at that section is above or below a threshold value s th = 2.2, which is slightly higher than the smaller stable point s 1 . In this figure, we can see that points in the same droplet are distributed on a line with constant p, and smaller droplets have higher p. This result indicates that "chemical potential" p is uniquely determined by the morphology of the droplet, which is suggestive of the existence of the surface tension. We also confirm similar results for the 3D simulation and the 2D simulation of Model II (see Fig. 4f and Appendix D, Figs. 5 and 6, respectively).
To analyze these results in further detail, in Fig 4d we show the simultaneous distribution of droplet size, S drop = (s − s 1 )dr (the integral is taken within each droplet) and the value of p in the droplets. The data are sampled at three different points in time (t = 800, 3200, and 12800). The results clearly show that there is a one-to-one correspondence between S drop and p. Furthermore, for the large S drop region we find the relation Fig. 4d). This asymptotic behavior can be understood as follows: when S drop is large, s in the droplet saturates to the larger stable point s 2 and is connected to the outside of the droplet via the interface, the width of which is almost constant and independent of the droplet size. (see the bottom panel of Fig 4e). In this situation, S drop can be approximated as S drop ∼ (s 2 − s 1 )V drop , where V drop is the volume (area) of the droplet. Since the droplets have a spherical shape, by denoting the radius as R, we can obtain the relationship (p − p ∞ ) ∝ 1/R. This is the same form as the Young-Laplace equation, where the proportionality constant of the above relation is expected to be associated with the surface tension in MCRDs. Note that p coincides with p ∞ for a large curvature limit (R → ∞), a situation in which two stable states are segregated by a flat interface, as discussed In Sec. III A. The same relationship is confirmed for the 3D simulation (see Fig 4g) and a 2D simulation with a different reaction term (Model II, see Appendix D, Fig 7c). When S drop is small, s in the droplets is far from s 2 and the above approximation does not hold (see the top panel of Fig 4e), leading to a deviation of p−p ∞ from ∝ S The reason why the density profile of the droplets changes depending on the droplet size is that initial density s 0 is very far from the stable fixed point s 2 ; thus, it takes long time for s to saturate to s 2 while retaining the mass conservation law for s. We note again that such a trend is not a feature specific to MCRDs, but is also seen in phase separation phenomena with a highly asymmetric shape in the free energy form [25][26][27].

E. Reduced mass transport equation
Finally, we discuss the physical origins of "surface tension" in MCRDs. This can be translated into a question of whether or not the variable p that we call "chemical potential" can be written in some variational form. The governing equation of p is Eq. 3 and 7 provide a closed set of partial differential equations with respect to s and p, which is equivalent to the original equation for MCRDs (Eq. 1 and 2). Since variable s is a conserved quantity, the time evolution of s is driven only by the diffusive transportation of molecules (see Eq. 3) and is expected to proceed slowly compared with that of p (Eq. 7). Indeed, in the simulations we observed that the time derivative of p quickly decays in the early linear time regime and ∂p/∂t is close to 0, where the terms on the right-hand side of Eq. 7 are balanced. Thus, p is subjective to s, i.e., we can regard p as a functional of s. The validity of this description can be confirmed by a one-to-one correspondence between p and s, as shown in Figs. 4d and g. In addition, since p at steady states becomes spatially constant, the spatial profile of p is gradual in the coarsening regime (see Fig. 4b).
Thus, ∇ 2 p can be regarded as a small perturbative term and we can expand p as where we neglect the spatial derivatives higher than the forth order (see Appendix E for detailed derivation and precise meaning of the above expression). P (s) provides the leading approximation, which is given by a steadystate solution of p. A(s) is also a function of s that is determined by bistable steady stable solution. Substituting Eq. 8 into Eq. 3, we obtain a mass transport equation in the following form: where κ(s) and ζ(s) are defined as κ(s) ≡ −A(s)dP/ds and ζ(s) ≡ A(s)d 2 P (s)/ds 2 . This equation for mass transport satisfies generic conditions such as mass conservation and spatial inversion symmetry of the system. Notably, Eq. 9 has the same form as the coarse-grained equation obtained for the system showing motility induced phase separation (MIPS), for which Solon et al. showed that the system can be mapped into a variational form [40]. This indicates that Eq. 9 also follows the same variational principle. At the lowest order, s obeys the non-linear diffusion equation ∂s/∂t ∇ [(dP/ds)∇s], indicating that the sign of "diffusion coefficient" dP (s)/ds determines the stability of the homogeneous solution.
In general, one can prove the relationship , and dP/ds < 0 is nothing but Turing instability condition (see also the growth rate λ(q) in the caption of Fig. 1b). For Model I, we can analytically evaluate P (s) and obtain the same stability condition T < T c = 1/8 as before (Appendix F). These consistencies support the validity of the above reduction.

IV. CONCLUSION
In summary, we find that the pattern dynamics in MCRDs, starting from a uniform state towards an eventual single isolated domain, are classified into two types, similar to phase separation phenomena (see also Ref. [13]). One is triggered by the absolute instability of the initial uniform state and the other goes through a nucleation and growth process. This classification is systematically addressed via a phase diagram, which can be constructed by the derivative of a free-energy like function introduced by Eq. 6. Furthermore, we performed large-scale numerical simulations in both two and three dimensions. We confirmed that in the late stage, droplet patterns appearing both in stable and metastable regions obey a relation similar to the Yong-Laplace equation, and their coarsening process is limited by the evaporation-condensation mechanism with growth exponent 1/3. These results suggest that in the presence of conserved variables, a surface-tension-like quantity generally emerges in RDs. By a reduced description of the systems focusing on the dynamically slow mode associated with mass conservation, we show that the origins of "surface tension" can be understood from chemical potential like functional p[s] following a variational form.
MCRDs were originally introduced and have been discussed as models for molecular localization such as membrane-bounded GTPases which are responsible for the formation of cell polarity [5][6][7][8][9][10][11][12][13][14] and neuron dendrites [41], and phosphoinositide lipids [42]. Recently, MCRDs have been also applied to various systems such as oscillatory motion in min-protein [43] and chemical turbulence [44], by which the importance of MRCDs for pattern formation in biological systems is further recognized. Although the details of chemical reactions are significantly different among the systems, MCRDs show universal dynamics and can be relevant for many existing systems. One interesting possible application of MCRDs is for liquid-liquid phase separation observed inside cells [45][46][47][48]. Liquid droplets found in cells contain many chemical components that regulate each other; thus, some droplets may be formed by chemical reactions consuming cellular energy rather than non-molecular specific physical interactions among proteins. MCRDs can provide a natural framework to model and understand the general mechanism of such phenomena.
Although MCRDs exhibiting pattern formations do not satisfy detailed balance, they show very similar behavior with phase separation driven by the thermodynamic variational principle. In this sense, phase separation observed in MCRDs can be regarded as a nonequilibrium extension of the typical phase separation phenomena. It is then interesting to compare MCRDs with MIPS, another non-equilibrium extension of phase separation, where self-propelled particles with repulsive interaction spontaneously form spatially localized assemblies [49,50]. In the systems showing MIPS, transport of molecules is supposed to be active, while in MCRDs transport itself is by simple diffusion and activity is introduced into chemical reactions. Despite this difference, both out-of-equilibrium systems show similar dynamics of phase separation, indicating that one could find a theoretical principle to understand both phenomena in an integrated manner. For MIPS, Solon et al. derived a time evolution equation for particle density, and based on the equation they proposed a generalized variational principle, as a generalization of the variational principle of equilibrium system [40]. For MCDRs, the reduction equation Eq. 9 that we derived in the present study is of the same form as that derived for the particle density in MIPS. Therefore, the generalized variational principle could be applicable to MCRDs. Although an energetic interpretation of the extended variational functionals is still lacking at present, the study of this common mathematical structure is one of directions to be explored in the future.

ACKNOWLEDGMENTS
Numerical simulations were performed on the SGI ICE XA/UV hybrid system at the ISSP at the University of Tokyo. This study was supported by JSPS KAK-ENHI (JP20K14424 and JP18057992) and JST CREST JPMJCR1923, Japan.

Appendix A: MCRD with fluctuations
The mass-conserved reaction-diffusion equation including fluctuations is given as follows: Here D u and D v are the diffusion coefficients of chemical species, U and V , and we assume D u < D v . f is the reaction term. j u and j v are fluctuations accompanying with diffusion of U and V , respectively, satisfying the below statistics, η r is fluctuation accompanying with reaction and satisfies In our study, we use a reaction term of the form In this case, Appendix B: Non-dimensionalization By denoting the time and space units as τ and , respectively, and taking the units for u, v and f as u 0 , v 0 and f 0 , respectively, we can rewrite Eq. A1 and A2 as, where we define R = v 0 /u 0 , D = D v /D u , σ = 1/ v 0 d , = √ D u τ , and τ = v 0 /f 0 .Ã represents scaled variable of A. The statistics of the noise terms are then given as j u,α (t,r)j u,β (t ,r ) = 2ũδ(t −t )δ(r − r )δ α,β , j v,α (t,r)j v,β (t ,r ) = 2ṽδ(t −t )δ(r −r )δ α,β and η r (t,r)η r (t ,r ) = (f + +f − )δ(t −t )δ(r −r )δ α,β . From the above equations, we can immediately obtain the time evolution equation for the conserved variables =ũ + Rṽ as,∂ wherep =ũ+RDṽ. Here we use u 0 as a common unit for s and p. The last term on the right-hand side represents the diffusion noise that satisfies j p,α (t,r)j p,β (t ,r ) = 2pδ(t −t )δ(r −r )δ α,β . For Model I (k u = ca/(b + u 2 ), k v = c), we choose u 0 = √ b, v 0 = a/ √ b and f 0 = ac/ √ b, and the resulting reaction term isf =ũ/(1+ũ 2 )−ṽ. For Model II (k u = α, k v = k 0 + βu 2 /(K 2 + u 2 )), we set u 0 = K, v 0 = αK/β, f 0 = αK and k 0 = 0, and the resulting reaction term is f =ũ 2ṽ /(1 +ũ 2 ) −ũ.
In the absence of noise terms (σ = 0), bistable points in the steady-state solution, (ũ 1 ,ṽ 1 ) and (ũ 2 ,ṽ 2 ), are obtained using two equalitiesf = 0 andũ + RDṽ = p ∞ (see Fig. 1c), wherep ∞ is determined by Eq. 5 (see also Appendix C). Becausef is parameter-free in our normalization, the existence of such a solution depends solely on RD ≡ 1/T as shown in Fig. 1e. Note that if we allow k 0 to be varied freely for Model II, another parameter needs to be considered in addition to T .
By the non-dimensionalization described above, the condition for the Turing instability (Eq. 4) is rewritten as fũ − Tfṽ < 0. For Model I, T < 1/8 (≡ T c ) is necessary for this inequality to be satisfied and T c is the value of T at the vertex of a parabolic function shown in Fig. 1e. Such an upper bound does not exist for Model II (see Fig. 5b).
All the data presented in the figures are normalized by the units introduced here.
Appendix C: Derivation of Eq. 5 We consider a solution for MCRDs in one dimension, where (u, v) approaches (u 1 , v 1 ) and (u 2 , v 2 ) for x → −∞ and +∞, respectively. In a steady condition, Eq. 3 becomes ∇ 2 p = 0 and the general solution is given by p = Ax + B, where A and B are constants. If A is nonzero, p (and thus, at least one of u or v) diverges for x → ±∞, which is unphysical. Thus, p has a constant value, which we denote as p ∞ . As discussed in the main text, we consider a solution that propagates with a constant velocity c [23]. Substituting u(x, t) = u(x − ct) into Eq. 1, we obtain the following relationship where z = x − ct is a moving coordinate at speed c. By multiplying the above equation by du/dz and taking integral over the entire space, we obtain The first term on the right-hand side vanishes because du/dz = 0 for z = ±∞. Therefore, c is related to u(z) as follows: Here, Ψ is an integral of f with respect to u introduced in Eq 6. At the stationary state where c = 0, two stable states are related by Ψ(u 2 , p ∞ ) = Ψ(u 1 , p ∞ ). This equality determines a position-independent value of p = p ∞ . Note also that u 1 and u 2 are two minima of Ψ(u, p ∞ ), In the main text, p ∞ is evaluated by assuming the periodic boundary conditions for which the system size L is sufficiently large. In the main text, we investigated the pattern formation process for Model I. To examine model dependence of the results obtained there, in this section we analyze Model II and compare it with Model I. By performing the steady state analysis explained in the main text, we obtain results corresponding to Fig. 1c and e, presented in Fig. 5a  the reaction term f (u, v) = 0 has two branches that extend to infinitely large v, and always crosses with the line p ∞ = D u u + D v v at two points for any large value T (see Fig. 5a). Note that we have set k 0 = 0 here, and drawing the full phase diagram for Model II requires an additional axis representing k 0 . Below we present the results of coarsening dynamics in 2D at point (u, T ) = (0.125, 1) (see the red cross symbol in Fig. 5b). Fig. 6a shows the temporal change of pattern (the density field of s), where droplet-like patterns form in the early time (t = 200) and coarsen over time. Fig. 7a shows the temporal change of the number of the droplets n(t), which asymptotically approaches a powerlaw function with exponent −0.60 (red line). Thus, the growth rate is 0.60/d = 0.30 (d = 2), which is coincident with the results for Model I (Fig. 2c and 3d). Fig. 7b shows the scaled structure factor, where we can see that the structure factors at different time are scaled for both small and large q after the system reaches the coarsening regime where n(t) shows a power-law decay. In Fig. 2d, we have seen that the scalability for large q in intermediate time regime does not hold for Model I. This difference is ascribed to whether the density inside the droplets is saturated to the stable fixed point s 2 (see Fig 1a). For the current parameter choice for Model II, the values of the bistable steady solution s 1 and s 2 are not so far from that of the initial homogeneous state. The density inside the droplets easily saturates to the stable points, and as a result, the system shows a self-similar coarsening. Fig. 6b shows the map of p at the same spatial region and time as Fig. 6a, where we can see a trend similar as Fig. 4b. This similarity, together with the growth exponent 0.30, indicates that the evaporation-condensation mechanism also works in Model II. Moreover, in Fig. 6c, one-to-one correspondence between p and the droplet size by a power-law with exponent −1/d is also confirmed as in Model I (see Fig. 4d and g). These results indicate that a surface-tension-like quantity generally exists in MCRDs without depending on details of the reaction terms, and plays a crucial role in the coarsening process of the patterns.

Appendix E: Derivation of reduction equation
Here we provide derivation of a reduction form Eq. 8 by taking account of the slow dynamical mode associated with mass conservation and spatial inversion symmetry.

Basic property of MCRDs
We consider the equations for MCRDs closed with p = D u u + D v v, s = u + v (Eqs. 7 and 3 in the main text). Volume of the system V = L d (d being the spatial dimension of the system) is sufficiently large and a periodic boundary condition is assumed.
is the conserved quantity of the system. As discussed in the main text p(t, r) becomes constant at the eventual stationary state p = P (s 0 ), which depends on the conserved quantity s 0 . For given s 0 , s at the stationary state is the solution of We denote this stationary solution as s(r) = S(r, s 0 ). The stationary solution (p, s) = (P (s 0 ), S(r, s 0 )) is parameterized by s 0 and constitutes a center manifold in the space of (p, s). The derivative of the stationary solution with respect to s 0 , is the tangential vector of the manifold, and is a zeroeigenmode of a liner operator L, Lη s = 0, where L is derived from Eqs. 7 and 3, The conjugate vector of η s is η s = (0, 1) for which the following relationship is shown.
Here, L † is the conjugate operator of L, and the scalar product between a(r) and b(r) is defined by All the relationships presented above are general consequences of mass conservation. Note that there exist other zero-eigenmodes of L, such as η x = (0, ∂ x S(r, s 0 )) that originates from the spatial translation invariance of the system along x-axis. Similarly, η y = (0, ∂ y S(r, s 0 )), etc., are zero-eigenmodes in multi-dimension. They are irrelevant to the following analysis.

Dynamics on the center manifold
We derive a single variable equation from Eqs. 7 and 3 to describe the coarsening process of droplets based on the center manifold projection reduction method [52,53].
As p and s in late stage are supposed to be quasi-static, we put an ansatz that they are expressed in the following form: p(t, r) = P (c) + ρ p (r, c) , (E8) s(t, r) = S(r, c) + ρ s (r, c). (E9) In each equation, the first term describes the dynamics along the center manifold, and the second is the correction that is orthogonal to it. Therefore, ρ ≡ (ρ p , ρ s ) is constrained by η s , ρ = 0 and η s , ∂ t ρ = 0. c is "coarse-grained mass density", defined as Here,Ṽ is a volume element centered at X, taken to be larger than the interfacial width but smaller than a typical center-of-mass distance between droplets. X is a coordinate to represent corresponding large-scale spatial variation of c. Then we employ the multi-scale variable method [52,53] in which r and X are treated as independent variables, and ∇ is replaced by ∇ + ∇ X ; ∇ X is defined by derivatives with respect to X and is a small parameter representing gradual spatial variation of c. By substituting Eqs. E8 and E9 into Eqs. 7 and 3, one obtains the following equations where N is the non-linear part of f defined by By taking scalar product between η s and Eq. E11, and using the relations η s , ∂ t ρ = 0, η s , Lρ = L † η s , ρ = 0, and η s , ∇a = 0 for an arbitrary function a(r), we obtain Eq. E14 provides the mass transport equation in this system. The first term is the lowest order estimation of ∂c/∂t. For evaluating higher term, our aim below is to analyze the behavior of ρ p . We expand ρ by and evaluate ρ p for higher orders. ρ = ρ (1) + 2 ρ (2) + · · · . (E16) Here, each ρ (i) is supposed to be orthogonal to η s , i.e., η s , ρ (i) = η s , ∂ t ρ (i) = 0. Note that ∂ t ρ = (∂ρ/∂c)(∂c/∂t) = O( 3 ) and does not contribute to B up to the second-order analysis performed below.
At the first order, Eq. E11 leads into The solution of the linear equation is given by where η s is the homogeneous solution given by Eq. E3, while η (1) is the particular solution for given g (1) in Eq. E17. a (1) is constant with respect to r, and is determined to satisfy the orthogonality condition η s , ρ (1) = 0. The second line of Eq. E17, ∇ 2 ρ (1) p = 0 indicates that ρ (1) p is independent of r. By taking into account the first line of Eq. E17, we find that the solution is ascribed to the homogeneous one in Eq. E18, i.e., a (1) ∝ ρ (1) p . Notice that η s and η (1) have opposite parity against spatial reflection r → −r, i.e., η s (S(r), r) → η s (S(−r), r) and η (1) (S(r), r) → −η (1) (S(−r), r), because g (1) consists of only the first order derivative about r. The inner product between the functions with opposite parity vanishes; therefore, a (1) = − η s , η (1) = 0. Thus, B(c) in Eq. E14 is zero at this order. This conclusion is simply understood by the fact that the system is invariant against spatial reflection, so that the equation should not include odd terms about spatial derivative ∇ X , excluding odd orders of . We can also determine ρ (1) s (r, c) from Eq. E17 in the same way, although its explicit form is not necessary below.
The solution of ρ (2) p has the form ρ The higher term of B(c) is O( 4 ), compatible with system invariance against spatial reflection. By substituting B(c) into Eq. E14, we reach the fol-lowing equation for c: which describes the mass transport process between droplets. The second term corresponds to the perturbative term ρ p at the lowest (second) order. Numerical data show that p exhibits gradual spatial variation, compatible to center-of-mass distance among the droplets (see Fig. 4b and 6b). Thus, at the scale of our interest, we may replace ∇ X and c by ∇ and s, respectively, and write a reduction form for p as follows. A is the function of s determined by steady stable solution S. However, it is usually difficult to obtain their exact expression for given MCRDs. On the other side, P (s) can be rather easily obtained since P coincides with p for a homogeneous steady solution (u 0 , v 0 ). In other words, P satisfies P For Model I, u 0 is a solution for u 3 0 − su 2 0 + (a + b)u 0 − sb = 0. When the left-hand side of the cubic equation is monotonic for u 0 , a unique (real) solution u 0 exists for any s, and thus we can analytically determine the specific form of P . In Fig. 8, we show the functional profiles P for various parameter sets. The functional profile of P changes from a monotonically increasing function to nonmonotonic one as dimensionless parameter T crosses the critical value T c = 1/8. Here we fix R = 1/4 and vary D. In the figures, we can see that P monotonically increases for T > Tc = 1/8 whereas for T < Tc, the region s exists where P = P (s) has three solutions.