Theory of Elastic Microphase Separation

Elastic microphase separation refers to equilibrium patterns that form by phase separation in elastic gels. Recent experiments revealed a continuous phase transition from the homogeneous phase to a regularly patterned phase, whose period decreased for stiffer systems. We here propose a model that captures these observations. The model combines a continuous field of the elastic component to describe phase separation with nonlocal elasticity theory to capture the gel's microstructure. Analytical approximations unveil that the pattern period is determined by the geometric mean between the elasto-capillary length and a microscopic length scale of the gel. Our theory highlights the importance of nonlocal elasticity in soft matter systems, reveals the mechanism of elastic microphase separation, and will improve the engineering of such systems.

Phase separation in elastic media is a ubiquitous phenomenon, which is relevant in synthetic systems to control micro-patterning [1][2][3] and in biological cells, where droplets are embedded in the elastic cytoskeleton or chromatin [4][5][6].While biological systems are typically dynamic and involve active processes, the simpler synthetic systems can exhibit regular stable structures.These patterns harbor potential for metamaterials and structural color, particularly since they are easier to produce and manipulate than alternatives like self-assembly by block co-polymers [7] or chemical cross-linking [8].In these applications, it is crucial to control the length scale, the quality, and the stability of the pattern.Such control might be possible in a recent experiment, which found stable equilibrium patterns [1].However, the underlying mechanism for this elastic microphase separation is unclear, complicating further optimization.
The elastic microphase separation experiment [1] proceeds in two steps (Fig. 1A): First, a PDMS gel is soaked in oil at high temperatures for tens of hours until the system is equilibrated.When the temperature is lowered in the second step, the sample develops bicontinuous structures, reminiscent of spinodal decomposition.However, in contrast to spinodal decomposition, the length scale of the structure does not coarsen but stays arrested at roughly one to ten micrometers, depending on the gel's stiffness.Interestingly, this transition is reversible and the pattern disappears upon reheating, suggesting a continuous phase transition.Moreover, the resulting pattern is independent of the cooling rate, in contrast to other experiments on similar materials [3,9].Consequently, the experiments should be explainable by an equilibrium theory that captures elastic deformations in PDMS due to oil droplets formed by phase separation.
In this paper, we propose a theoretical model that explains the experimental observations [1].The model combines the continuous density field of the elastic component, which naturally describe phase separation, with a nonlocal elasticity theory to capture the microstructure of the gel [10][11][12][13][14].This approach allows us to capture the continuous phase transition to a patterned phase and predict its equilibrium period.

II. RESULTS
To explain the experimental results [1] using an equilibrium theory, we define a free energy comprising contributions from elastic deformation as well as entropic and enthalpic contributions that can induce phase separation.While elastic deformations are naturally described by the strain tensor ϵ(X) defined in a reference frame X, phase separation is typically described by the volume fraction field ϕ(x) of the elastic component in the lab frame x.For simplicity, we will focus on one-dimensional systems in this paper, where volume conservation connects the scalar strain ϵ to the fraction ϕ in the reference frame, where ϕ 0 denotes the fraction in the relaxed homogeneous initial state.The fraction ϕ(x) in the lab frame then follows from the coordinate transform dx/dX = ϵ(X)+1.This connection between strain ϵ and volume fraction ϕ permits a theory in terms of only one scalar field.

A. Local elasticity models cannot explain equilibrium pattern
We start by investigating a broad class of elastic models, where the elastic energy density is a function of strain ϵ.Since ϵ can also be expressed in terms of the volume fraction ϕ (Eq.1), the free energy of the system reads where k B is Boltzmann's constant, T is the constant absolute temperature of the system, and ν is a relevant molecular volume, e.g., of the solvent molecules.Here, f 0 captures the elastic energy density as well as molecular interactions and translational entropy associated with ordinary phase separation, while the second term proportional to the interfacial parameter κ penalizes volume fraction gradients and gives rise to surface tension.Eq. 2 is identical to basic models of phase separation without elastic contributions [15,16].Such models exhibit phase separation and subsequent coarsening to minimize interfacial costs, known as Ostwald ripening [17].While adding local elasticity alters f 0 (ϕ), functions that minimize F local can only have a single interface [18] and equilibrium patterns with finite length scales are thus impossible.We show in the Appendix A that this result generalizes to higher dimensions.Taken together, field theories based on local elasticity, including sophisticated non-linear finite strain models, cannot explain the equilibrium length scales observed in experiments.

B. Microscopic picture suggests nonlocal elasticity theory
Why do standard elastic theories fail to explain the observed patterns?One answer is that only the interfacial parameters κ carries dimensions of length in Eq. 2, so on dimensional grounds we cannot expect another length scale beyond the interfacial width to emerge.While the interfacial width is typically governed by molecular sizes (∼ 1 nm), realistic elastic meshes exhibit additional length scales like the mesh size (∼ 10 nm [14,19,20]) and correlation lengths of spatial inhomogeneities (∼ 100 nm [21][22][23]).Since the last two quantities are comparable to the pattern length scale (several 100 nm to several µm [1]), we hypothesize that a characteristic length of the mesh is key for explaining the observed patterns.
If microscopic lengths of the elastic mesh are relevant, local elastic theories are insufficient [10][11][12][13][14].This is because moving a particular crosslink transmits forces to connected crosslinks in the vicinity (see Fig. 1B), implying stresses are no longer local, and the associated elastic energy cannot be expressed as a function of the strain.Instead, the stress on a particular crosslink is now given by a sum over all connected crosslinks, whose contributions decay with distance X in the reference frame [10,13].In a continuous field theory, this nonlocal averaging is expressed as a convolution operation [10,13].Using a simple linear elastic model for the local stress Eϵ, with elastic modulus E, we obtain the nonlocal stress where we choose a Gaussian convolution kernel [13,24], with a characteristic length ξ, which quantifies the microscopic length of the gel [10,13,24].This nonlocal model can also be derived more rigorously, either generically (see Appendix B) or from a more explicit microscopic model [10,25].Note that the convolution is performed in the reference frame since the topology of the network, governing which crosslinks interact with each other, is determined in this unperturbed state.
The elastic energy density of the system is then given by the product of strain and nonlocal stress, so the free energy of the entire system reads where F local now only captures the contributions associated with phase separation.To capture the essence of phase separation, we consider a simple Flory-Huggins model for the local free energy density [26][27][28], where 1 − ϕ is the solvent fraction.Here, the first two terms capture entropic contributions, while the last term describes the interaction between elastic and solvent components, quantified by the Flory-Huggins parameter χ.
Taken together, Eqs.1-6 define the free energy F nonlocal as a functional of the fraction ϕ of the elastic component.

C. Nonlocal elasticity enables equilibrium patterns
We start by analyzing equilibrium states of the model by determining profiles ϕ(x) that minimize F nonlocal using a numerical scheme described in the Appendix C. Beside typical macroscopic phase separation, we also find periodic patterns for some parameter sets; see Fig. 1C and Fig. S1.In soft systems (small stiffness E), dilute regions, corresponding to solvent droplets, alternate with dense regions, where the elastic mesh is hardly strained (ϵ ≪ 1).In contrast, a harmonic profile emerges for stiff systems (large E).Taken together, the nonlocal elastic theory supports periodic patterns that qualitatively resemble the experimentally observed ones [1].
To understand when periodic patterns form, we next investigate the simple case where components can freely exchange with a surrounding reservoir kept at fixed exchange chemical potential µ.This situation allows solvent molecules to rush in and out of the system, adjusting Grand-canonical phase diagrams reveal patterned phase.(A) Phase diagram as a function of the chemical potential µ and the interaction strength χ for E = 0.01 kBT /ν.Homogeneous phases (region H) coexist on the brown line between the critical point of phase separation (black disk) and the triple point (gray disk), while the patterned phase (region P) coexists with the homogeneous phase on the blue-brown-dashed line.(B) Phase diagram as a function of µ and χ for E = 0.2 kBT /ν.The binodal line separating the homogeneous and patterned phase exhibits either a first-order transition (blue-brown-dashed line) or a continuous transition (red dotted line with associated critical points marked by red disks; details in the Appendix D). (A-B) Model parameters are ϕ0 = 1 and κ = 0.05 ξ 2 .
the average fraction φ of the elastic component.Fig. 2 shows two phase diagrams of this grand-canonical ensemble at different stiffnesses E. In the soft system (left panel), the phase diagram mostly resembles that of ordinary phase separation: For weak interactions (χ < 2), we find only a homogeneous phase and µ simply controls φ.In contrast, above the critical point at χ ≈ 2 (black disk), we observe a first-order phase transition (brown line) between a dilute phase (µ ≲ 0) and a dense phase (µ ≳ 0).However, at even stronger interactions (χ ≳ 3.3), an additional patterned phase (denoted by P) emerges, where the periodic patterns exhibit the lowest free energy.The line of the first-order phase transitions between the patterned phase and the dilute or dense homogeneous phase (blue-brown-dashed curves) meet the line of the phase transition between the two homogeneous states at the triple point (gray disk), where these three states coexist.
The grand-canonical phase diagram of soft systems (left panel of Fig. 2) qualitatively resembles simple pressure-temperature phase diagrams, e.g., of water.Assuming that the chemical potential µ plays the role of pressure and that the interaction χ is negatively correlated with temperature, the dilute and dense homogeneous phases respectively correspond to the gas and liquid phases.They become indistinguishable at the critical point at low interaction strength (corresponding to high temperatures).In contrast, the patterned phase, with its periodic microstructure, resembles the solid phase.
The general form of the grand-canonical phase diagram persists for stiff system (right panel of Fig. 2), although the parameter region of the patterned phase is much larger.However, the first-order transition between the dilute and dense homogeneous phases disappears together with the normal critical point of phase separation.Instead, we now find a continuous phase transition (dotted red line) between the homogeneous and the patterned phases, which we will discuss in more detail below.Taken together, these phase diagrams suggest that stable patterned phases emerge for sufficiently large stiffness E and interaction χ for intermediated φ.
The grand-canonical ensemble that we discussed so far is suitable when the time scale of an experiment is long compared to the time scale of particle exchange with the reservoir.In the experiments [1], the initial swelling takes place over tens of hours with a measurable increase in size and mass, indicating that solvent soaks the sample until it is equilibrated with the surrounding bath.In contrast, the temperature quench, during which the patterned phase is observed, takes place on a time scale of minutes without the solvent bath.This suggests that this process is better described by a closed system.

D. Patterned and homogeneous phases coexist in closed systems
In the closed system, corresponding to a canonical ensemble, the average fraction φ of elastic components, and thus also the average fraction of solvent, is fixed.In this situation, we find that multiple different phases can coexist in the same system; see Fig. 3.This is again rem-iniscent of phase separation, where the common-tangent construction reveals the fractions in coexisting homogeneous states.Indeed, we find exactly this behavior in soft systems (left panel of Fig. 3A), where a dilute and dense phase coexist for fractions between the two vertical dotted lines, while the free energy of the patterned phase (blue line) is always larger and thus unfavorable.The picture changes for larger stiffness (right panel of Fig. 3A), where the patterned phase has lower energy and we can construct two separate common tangents, which respectively connect the dilute and dense homogeneous phase with the patterned phase.Analogously to phase separation, we thus expect situations in which a patterned phase coexists with a homogeneous phase (when φ is in the region marked with H+P or P+H).Fig. 3B corroborates this picture and shows various coexisting phases as a function of the stiffness E and the interaction strength χ.Taken together, the main additional feature of the canonical phase diagrams is the coexistence of multiple phases, which was only possible exactly at the phase transition in the grand-canonical phase diagram.

E. Higher stiffness and interaction strength stabilize patterned phase
The canonical phase diagrams shown in Fig. 3B are complex, but they generally preserve three crucial aspects of the grand-canonical phase diagram shown in Fig. 2: Higher stiffness (i) slightly favors the homogeneous phases, (ii) greatly expands the parameter region of the patterned phase, and (iii) induces a continuous phase transition.The first point is illustrated by the binodal line of the homogeneous phase (thick brown lines and red dotted lines), which moves up with increasing stiffness E, implying that larger interaction strengths χ are necessary to stabilize inhomogeneous systems.Inside the binodal line the system exhibits various behaviors, which can be categorized by χ.At a critical value χ * , the patterned phase (blue disk) coexists with the dilute and dense homogeneous phase (brown disks), and the associated tie line corresponds to the triple point in Fig. 2. For weaker interactions (χ < χ * ), we mostly observe coexistence of a dilute and dense homogeneous phase (region H+H), which corresponds to normal phase separation.For stronger interactions (χ > χ * ), the system exhibits the patterned phase, either exclusively (colored region) or in coexistence with a homogeneous phase (regions H+P and P+H).Larger stiffness E lowers the critical value χ * , thus expanding the parameter region where the patterned phase exists.Eventually, for sufficiently large E, χ * approaches the critical point of the binodal (gray point), a tiny region with patterned phase appears, and part of the binodal line becomes a continuous phase transition (red dotted line), reproducing the behavior predicted by the grand-canonical phase diagram of stiff systems (right panel in Fig. 2).The influence of stiffness E and interaction strength χ becomes even more apparent in the three-dimensional phase diagram shown in Fig. 3C: With increasing E, the χ associated with the critical point of phase separation (black line) increases slightly, whereas the states of three-phase coexistence (blue line and brown lines) shift to lower χ.All lines meet at the tricritical point (black sphere) for E ≈ 0.037 k B T /ν, φ ≈ 0.54, and χ ≈ 2.14.Increasing E further, a part of the binodal line exhibits a continuous phase transition, which expands with larger E. The phase diagram thus summarizes three main aspects of our model: First, the binodal line of phase separation, which is only weakly affected by E, determines whether the system can exhibit nonhomogeneous states.Second, if the system can be inhomogeneous, the stiffness E determines at what value of χ patterned phases emerge.Third, for sufficiently large E, these patterned phases form immediately due to the continuous phase transition.

F. Continuous phase transition explains experimental measurements
The continuous phase transition that we identified at sufficiently large stiffness E implies that the system can change continuously from a homogeneous phase to a patterned phase when the interaction strength χ is increased (corresponding to cooling).Indeed, the amplitude of the predicted pattern vanishes near the transition (right panel of Fig. 3B), while the length scale stays finite (left panel of Fig. 3B).This behavior is not expected for phase separating systems, where first-order transitions are typically, which are associated with a jump in observables (see gray dots in Fig. 4A for an example).
Pattern length scale exhibits scaling laws.Length scale L as a function of stiffness E (panel A) and interfacial parameter κ (panel B) for various parameters.Putative scaling laws are indicated and the prediction by Eq. 9 is shown for ϕ0 = 1, φ = 0.5, χ = 4, and γ ≈ kBT κ 1/2 /ν (green line).
The continuous phase transition was already hypothesized for the experiments [1], based on a lack of hysteresis and a continuous change of the contrast measured by light intensity.To connect to experiments, we mimic the contrast using the square of the amplitude of the optimal volume fraction profile.Fig. 4A and the right panel of Fig. 3B show that the contrast changes continuously from zero when the interaction strength χ is increased for sufficiently stiff systems.Moreover, Fig. 4B shows that the associated pattern length scale changes only slightly, consistent with the experiments.Note that deviations in the form of the curves could stem from thermal fluctuations, finite resolution in the experiment, and also deviations in model details.

G. Stiffness and interfacial cost control pattern length scale
We next use the numerical minimization of the free energy F nonlocal to analyze how the length scale L of the patterned phase depends on parameters.Fig. 5 shows that L decreases with larger stiffness E and increases with the interfacial cost parameterized by κ.The data in Fig. 5A suggests the scaling L/ξ ∝ E −1/2 over a significant parameter range, which matches the experimental observations [1].Moreover, Fig. 5B suggests L/ξ ∝ ξ −1/2 κ 1/4 , which has not been measured experimentally.Taken together, the two scaling laws suggest that the equilibrium length scale emerges from a competition between elastic and interfacial energy.
The two scaling laws emerge qualitatively from a simple estimate of the elastic and interfacial energies: Since shorter patterns have more interfaces, the interfacial energy per unit length is proportional to γL −1 , with surface Approximate model explains scaling laws.(A) Example for a volume fraction profile (pink lines) and the corresponding piecewise approximation (dotted gray lines) in the reference (top) and lab frame (bottom).(B) Derivatives of the average energy density (in units of kBT ξ/ν) as a function of the pattern period L. Shown are data from full numerics (dots), numerics for the piecewise profile (solid lines), and asymptotic functions (dashed lines) for the elastic (gray) and negative interfacial energy (violet).The stable length L corresponds to the crossing point of the elastic (black) and the interfacial terms (violet).Model parameters are E = 0.02 kBT /ν, κ = 0.05 ξ 2 , ϕ0 = 1, φ = 0.5, and χ = 4. tension γ ∝ κ 1/2 [15].In contrast, the elastic energy of a single period originates from stretching a part of material from initial length ξ to final length L, resulting in an elastic energy density proportional to ELξ −1 .Minimizing the sum of these two energy densities with respect to L results in L/ξ ∝ ξ −1/2 E −1/2 κ 1/4 , which explains the observed scalings qualitatively.

H. Approximate model predicts length scale
To understand the origin of the length scale L in more detail, we consider the limit of strong phase separation, where the interfacial width is small compared to L; see Fig. 1C.We thus approximate the volume fraction profile ϕ(x) of the elastic component by a periodic step function with fixed fractions ϕ − and ϕ + ; see dotted lines in Fig. 6A.Material conservation implies that the relative size of these regions is dictated by the average fraction φ in the swollen state, so we can only vary the period L of the profile.The stable period L then corresponds to the L that minimizes F nonlocal given by Eq. 5, implying F ′ nonlocal (L) = 0. Since changing L does not affect the local free energy f 0 , we only investigate the average free energy of the interface, fint ( L) ≈ 2γ L−1 , and the average elastic free energy, is the period in the reference frame.Fig. 6B shows the derivatives of these contributions with respect to L, indicating that they sum to zero for L = L.We show in the Appendix E that indicating three regimes bounded by Fig. 6B shows that this approximation of ∂ L fel captures the main features of the full numerical data.Fig. 6B suggests that stable patterns are mainly possible in the gray region (L min < L < L max ), which we interpret further below.In this region, we use Eq. 7 to solve ∂ L fel + ∂ L fint = 0 for L, resulting in L ≈ (8π) consistent with numerical results; see transparent green lines in Fig. 5.This expression shows that the stable period L is governed by the geometric mean of the elastocapillary length γ/E and the microscopic length ξ.Moreover, L increases with a larger average fraction φ of the elastic component, i.e., less swelling.In contrast, the fraction ϕ + has only a weak influence since it is close to 1 in the case of strong phase separation, implying that the interaction strength χ affects L only weakly.

I. Patterned phase is governed by reference state
Finally, we use the approximate model to understand when the patterned phase emerges.Here, it proves useful to interpret Eq. 8 in the reference frame, where the convolution of the nonlocal elastic energy takes place.Defining the length L 0 = φ ϕ0 L in the reference frame and the associated fraction α 0 = ϕ− ϕ0 (ϕ + − φ)/(ϕ + − ϕ − ) occupied by the solvent droplet (Fig. 6A), we find where the numerical pre-factors are very close to one.The first condition (L 0 ≳ ξ) suggests that two solvent droplets need to be separated by more than ξ in the reference frame since L 0 roughly estimates their separation; see Fig. 6A.If droplets were closer, they would feel each other's deformations, which is apparently unfavorable.In the extreme case (L 0 ≪ ξ), the average elastic energy is almost constant, essentially because short-ranged variations are averaged by the comparatively large nonlocal kernel.In contrast, the second condition implies that the droplet size in the reference frame (α 0 L 0 ) must be smaller than the microscopic length scale ξ.Assuming ξ corresponds to the mesh correlation length, this suggests that the droplet can at most deform the correlated part of the mesh, which might correspond to large soft regions in natural meshes.If droplets were larger (α 0 L 0 ≫ ξ), nonlocal features would only be relevant at interfaces, so the system would behave as if it had only local elasticity and coarsen indefinitely.This analysis highlights that the existence of the periodic pattern depends on the reference frame, while its length scale L also depends on the different stretch of the dilute and dense region; see Fig. 6A.This observation suggest an intuitive explanation for the influence of the interaction χ: Assuming that ϕ − and ϕ + correspond to equilibrium volume fractions and φ = 1  2 for simplicity, we find α 0 ∝ ϕ − , which decreases with larger χ.Consequently, the lower bound L min is unaffected, while L max increases, consistent with our observation that the patterned phase forms easier at higher χ and the scaling law given by Eq. 9 holds for broader parameter range with higher interaction strength (Fig. 5).

III. DISCUSSION
We propose a theory that explains the experimentally observed elastic microphase separation [1] based on nonlocal elasticity, which captures aspects of the microscopic gel structure.Within this theory, regular periodic patterns appear for sufficiently strong phase separation (large enough χ) and stiffness E, while surface tension γ opposes the trend.Essentially, solvent droplets inflate a region of the elastic mesh of the size of the microscopic length ξ.The pattern period L then results from a balance of elastic and interfacial energies, so that L scales as the geometric mean between ξ and the elasto-capillary length γ/E.In contrast, the interaction strength χ, leading to phase separation in the first place, affects L only weakly, but it determines whether the patterned phase is stable, similar to ordinary phase separation.However, the normal first-order transition between the homogeneous and heterogeneous phase (at the binodal line) can now also exhibit a continuous phase transition.Consequently, the patterned phase can appear with arbitrarily small amplitude in a reversible process.
Our model captures the main features of the experiment [1], including the continuous phase transition leading to reversible dynamics.Moreover, it explains that the pattern length scale L is independent of the cooling rate, only weakly affected by the final temperature, and decreases with stiffness E. Importantly, our model predicts that a structural length ξ of the mesh is essential for the emergence of the observed L. Our numerics indicate that L can be an order of magnitude larger than ξ, suggesting that ξ could relate to observed correlation lengths of the order of a few hundred nanometers [21].Since ξ is small compared to the distance between droplets (see Eq. 10), the nonlocal effects of elasticity do not affect droplet positioning.Furthermore, we found that a coexisting homogeneous phase does not affect the free energy of the patterned phase strongly (see Appendix C), suggesting that the two phases can be interspersed, which would contribute to irregularity of the droplet placement in real systems.In contrast, the observed variation in droplet size [1] likely originates from local heterogeneity in material properties, like ξ, E, and γ.Taken together, our theory makes clear predictions that could be tested experimentally.
To capture the mesh's microstructure, we employ nonlocal elasticity [10][11][12][13][14] based on a convolution of the stress field, which is similar to theories used in fracture mechanics [29].Our work complements related theories, which either modeled pores explicitly [20,[30][31][32][33] or resorted to particle-based methods [34,35].Nonlocality is generally responsible for the emergence of microstructure in multiple physical systems, such as the Ohta-Kawasaki model [36], phase separation with electrostatic interac-tion [37], and also nonlocal elasticity [38,39], e.g., to study polymeric materials [40,41].In contrast to the first two theories, we use a convolution in the reference frame, capturing the microscopic topology of the elastic mesh.More generally, the convolution kernel given by Eq. 4 can be interpreted as a Green's function of a diffusion process in the reference frame, suggesting that the nonlocal elasticity is similar to the damage field introduced in fracture mechanics [42].
We analyzed our model in the simple case of one dimension to highlight fundamental properties, but to capture experimental details, including various morphologies, we need to generalize the model to higher dimensions, which will require a tensorial convolution kernel [11].Moreover, we might require more realistic models of phase separation (including different molecular sizes and higher-order interactions terms) and elasticity (involving finite extensibility, viscoelasticity [43], as well as plastic deformation, like fracture [44,45] and cavitation, which can lead to regular droplet patterns [46]).Finally, experimental systems exhibit heterogeneities in key model parameters including ξ, E, and γ, which will contribute to uncertainty and might even induce large scale rearrangements [9,47].Such extended theories will allow us to compare the full pair correlation and scattering functions to experiments, shedding light on how we can manipulate this pattern forming system to control microstructures.where δ i,N is the Kronecker delta.Note that the last two equations are simply the incompressibility and the volume conservation of the elastic component, respectively.The constraint term F con [ϕ i , u, ζ, η] does not contribute to the free energy when the constraints are satisfied.
To minimize the free energy, we follow two steps: First we require that the phase separated structure is periodic, and describe the unit cell by a group of parameters θ [48].For a given θ, Eqs.A3 can be solved to obtain the free energy minimum with fixed unit cell where the symbols with asterisk denote the solution of Eqs.A3.The minimum of the free energy can then be obtained by optimizing F * with respect to the period θ.Note that F * is not a functional, but a function of θ instead.
Follow the steps in [48], we find that the derivative of F * (θ) can be obtained from the partial derivative of the free energy functional with respect to the period θ while keeping the shape of all the spatial functions unchanged, Assuming the total volume of the system V is constant, the total free energy can be replaced by the average free Note that Eq.A5 does not impose any assumptions on the exact form of the terms.For any local volume-fractiondependent term F frac [ϕ i ], e.g., the Flory-Huggins free energy, the partial derivative with respect to θ vanishes, since the variables ϕ i are dimensionless and have no explicit dependence on the pattern length scale.For any local elasticity, including nonlinear ones with large deformations, the elastic energy F el takes a local form of the deformation gradient tensor ∇ X u, which is also dimensionless and has no explicit dependence on the length scale.Therefore, the partial derivative of the local elastic energy vanishes as well, leading to Free energy exhibits a minimum when varying pattern length scale.Average free energy density fnonlocal = F nonlocal /V in units of kBT /ν as a function of the length scale L of the periodic pattern obtained numerically.Model parameters are E = 0.02kBT /ν, κ = 0.05ξ 2 , ϕ0 = 1, φ = 0.5, and χ = 4.
Here, term ∂ fint /∂θ * does not vanish since the interfacial term usually depends on the gradient of the volume fraction fields ∇ϕ i , which has the dimension of length −1 and contains explicit dependence on the length scale.
The period θ typically contains not only the size but also the shape of the unit cell, including the angles between the base vectors of the unit cell.Keeping the shape of the unit cell unchanged and denoting its size by L, we have which is usually negative since the interfacial energy prefers larger structure size.Specifically, for the interfacial energy fint = k B T ν where κ i quantifies the interfacial cost for component i, we have which is negative.Consequently, Eq.A6 states that with local elasticity, even if a periodic patterned structure is formed, it still cannot be the equilibrium state.Instead, Ostwald ripening is inevitable, since the free energy can always be lowered by coarsening, and the equilibrium length scale diverges.With nonlocal elasticity, the elastic energy will no longer be a local function of the deformation gradient tensor, implying two non-vanishing terms in Eq.A5, which compete with each other, leading to equilibrium pattern length scale (Fig. S1).
The conclusion drawn from the free energy derivative above, although not a strict proof, can be understand in a more intuitive way: Imagine a periodic patterned structure with average free energy f , which is first scaled to be slightly larger, and then relaxed.With local elasticity, the free energy decreases during scaling since only the interfacial term is affected.Since the free energy must not be higher after relaxation due to the variational principle, the new free energy is generally lower than f .This process can be repeated until the length diverges, implying local elasticity cannot explain the finite equilibrium length scale in a continuous field theory.
the elastic energy F el is a functional of the displacement field u(X), where X is the coordinates in the reference frame.Expanding F el to the second order, we have where Φ i are the derivatives of order i of F el .Since the constant part is irrelevant and the expansion is made near the relaxed state, the Φ 0 and Φ 1 terms can be dropped.Consequently, the lowest order approximation of F el reads Since the elastic energy must be invariant under constant shifts of the displacement field u, Extracting the diagonal element from Φ(X, X ′ ), and integrating both side with respect to X ′ and using Eq.B3, we find Note that the function Ψ is arbitrary and not subject to a constraint similar to Eq. B3, in contrast to Φ. Inserting Eq.B4 into Eq.B2 and using Eq.B3 again, we find The equation above has a very clear physical picture: The system contains multiple springs, with both ends of each spring tied to X and X ′ in the reference frame, and stretched by u(X) − u(X ′ ) after the deformation of the elastic component.The elastic energy of the system is just the total potential energy of all the springs, while the function Ψ is related to the stiffness of the springs.Defining the strain ϵ = du/dX, this can be written as where Π γ (X) is the box function centered at 0 with width γ.Assuming no explicit dependence on position, we have Note that this does not require a homogeneous deformation of the elastic component, but only that the property of the continuum is homogeneous.Defining the nonlocal stress σ nonlocal (X) = dX ′ ϵ(X)c(X − X ′ ), Eq.B7 becomes which is the elastic energy term in Eq. 5 of the main text.The nonlocal elasticity given by Eq.B10 can also be obtained from a more explicit model, such as a full network model of polymer network [25].Assuming that the network is formed by multiple Gaussian chains, and each polymer where J = du/dX + 1.
In addition to Eq. C3, other self-consistent equations for numerical calculation can be obtained by minimizing the sum of the four energy terms in Eq.C1 and Eqs.C4-C6, which reads Note that we transform all the coordinates to the reference frame, where the convolution exhibits a simple form and can be calculated efficiently with fast Fourier transforms (FFT).In practice, we use the periodic alternative u * (X) = u(X) − (ϕ 0 / φ − 1)X as the free variable since the displacement field u(X) is not periodic.Using u, w el and w s as the main variables, the following iteration scheme solves Eqs.C7, To improve numerical stability, a simple mixing method is used in most cases, where the differences between the new fields and the old ones are partially accepted, Here λ u and λ w are two empirical constants, which usually take value smaller than 0.1.To accelerate the convergence, Anderson mixing is also used every few steps [49].The variable-cell method is used to simultaneously optimize the period L of the structure during iteration.L is evolved in the direction of lowering the free energy [48,50], FIG. S3.Continuous phase transition is identified in two different ways.Phase diagram as a function of the average fraction φ of the elastic component and interaction strength χ for two stiffnesses E. The brown and black dashed lines mark the spinodal curves considering the microphase separation and the macrophase separation, obtained by evaluating Eq.D8 and Eq.D4, respectively.The binodal curves are obtained from the full numerics.The red dotted lines mark the window of continuous phase transition, obtained by considering the inequality in Eq.D12.
σ nonlocal with J, and transform the integral into the reference frame.The average free energy density then reads Since we have J = ϕ 0 / φ for the homogeneous phase average volume fraction φ, we find To test the stability, we first take the second-order derivative of the fhomo.with respect to φ, The spinodal of macrophase separation between two homogeneous phase corresponds to f ′′ homo.( φ) = 0, yielding see the black dashed lines in Fig. S3.
To consider microphase separation, we perturb J from a constant value, J = (ϕ 0 / φ) 1 + a cos(qX) , where a is the amplitude of the perturbation and q is the associated wave number.Evaluating Eq.D1, taking the second-order derivative of f with respect to a, and taking the limit a → 0, we find Stability of the homogeneous state requires ∂ 2 f /∂a 2 a=0 ≥ 0 for all q, implying χ holds for all q.χ u assumes its minimum at q = q * with q if (νE/k B T )ξ 2 ϕ 3 0 > 16κ φ5 , leading to the spinodal of microphase separation, see the brown dashed lines in Fig. S3.Comparing this line to the binodal curve obtained from the full numerics, we find that they overlap in a parameter windows, which hints at the existence of a continuous phase transition (Fig. S3).Note that if (νE/k B T )ξ 2 ϕ 3 0 > 16κ φ5 does not hold then χ u takes its minimum at q → +0, where χ u turns to χ macro and the macrophase spinodal given by Eq.D4 is recovered.

Higher-order analysis identifies the continuous phase transition
Eq. D8 defines the spinodal curve of the microphase separation in the ϕ-χ plane, which implies that the homogeneous phase is (meta-)stable below the curve, and unstable above the curve.No information about stability is provided right on the spinodal curve, since Eq.D6 is not a sufficient condition for stability when χ = χ u .In order to test the stability on the spinodal, we approximate J as J = (ϕ 0 / φ) 1 + a q * cos(q * X) + q̸ =q * a q cos(qX) and perturb the average volume fraction φ as φ + δ φ at the same time.The average free energy density can then be written as On the spinodal line of the homogeneous phase, we have while all second-order cross derivatives vanish.Here •| 0 indicates the value at the homogeneous state on the spinodal.
Expanding the free energy up to fourth-order, we find where omitted terms are either zero or absorbed in the remainder.The stability of the homogeneous phase on the spinodal then requires 1 24 Solving this inequality for φ numerically, we obtain the window of φ with a continuous phase transition.Combined with Eq.D8, the phase boundary of the continuous phase transition can be obtained, which is marked as red dotted lines or red surface in Fig. S3, Fig. 2 and Fig. 3.In fact, in all these phase diagrams, the continuous transitions are verified with both the overlapping of spinodal and binodal, as well as the higher-order stability analysis.

Appendix E: Approximate model and asymptotic solutions
To understand the scaling law of the length scale L, we assume sharp interfaces and approximate the volume fraction profile ϕ(x) as a box function, ϕ(x) = ϕ + + (ϕ − − ϕ + )Π αL (x − L/2) , (E1) within one period x ∈ [0, L).Here α = (ϕ + − φ)/(ϕ + − ϕ − ) is the fraction of the solvent-rich region relative to the period L, φ is the average volume fraction of the elastic component in the deformed state, and ϕ − and ϕ + are the minimum and the maximum value of the volume fraction profile, respectively.Converting the profile to the reference frame and making use of the relationship between strain and the volume fraction given by Eq. 1, we find where L0 = ( φ/ϕ 0 ) L and α 0 = (ϕ − / φ)α are the period and relative droplet size in the reference frame, respectively.
To evaluate the elastic energy, we first consider the case where the period is much larger than the microscopic length scale ( L ≫ ξ).In this case, we can safely ignore the interference between the neighboring periods since the Gaussian convolution kernel decays exponentially with distance.The elastic energy density thus reads fel = E 2 where fel,0 is a term with no dependence on L, with ϵ 0 = ϕ 0 /ϕ + − 1. Expanding Eq.E3 around L → +0, we have while expanding around L → +∞ leads to Next, we consider the case L ≪ ξ, where the elastic energy can be derived from the physical picture.Since the period is much smaller than the microscopic length scale ξ, the convolution of the strain field simply gives the average strain which is not affected by the period L. Evaluating Eq. 5 and converting it to energy density, we obtain Combining the results in Eq.E5, Eq.E6 and Eq.E7, we find fel ≈ where the boundary values L min and L max will be estimated later.Differentiating fel with respect to L, we have Since the derivatives of the average free energy density govern the equilibrium length scale (see Eq. A10), we determine L max by balancing the last two terms of the derivatives of fel , resulting in In contrast, the derivatives are all constants in the first two regimes of Eq.E9, so we cannot estimate the boundary in the same way.We thus balance fel directly in the first two regimes to get Converting the two bounds L min and L max to the reference frame yields Eq. 10 in the main text.
For completeness, we here also present the generic expression of fel of the approximated model with the ϑ function This integral can be evaluated numerically; see the black line in Fig. 6C for its derivative.Note that all orders of derivative with respect to L at L → +0 are zero, so the free energy is almost independent of L, consistent with our picture to derive Eq.E7.

FIG. 1 .
FIG. 1. Nonlocal elasticity yields regular equilibrium patterns.(A) Schematic picture of the experiment [1]: A relaxed elastic gel is swollen in a solvent at high temperature; After cooling, a regular pattern emerges.(B) Schematic of a polymer network displaying the displacement (opaque chains) from the reference state (transparent) after the central crosslink has been moved by u.We model the forces transmitted along the network using a nonlocal convolution kernel (blue density) of size ξ.(C) Equilibrium states for various stiffnesses E and interaction parameters χ for ϕ0 = 1, φ = 0.5, and κ = 0.05 ξ 2 .

FIG. 3 .
FIG. 3. Closed systems exhibit phase coexistence.(A) Schematic free energy of homogeneous and patterned phases with common-tangent construction (thin gray lines) for two stiffnesses E. Fig. S2 shows corresponding numerical results.(B) Phase diagram as a function of the average fraction φ of the elastic component and interaction strength χ for various E.Only the homogeneous phase (region H) is stable outside the binodal (brown line; black disk marks critical point) with a continuous phase transition at the red dotted part.Only the patterned phase (region P) is stable inside the blue lines with color codes indicating length scale and amplitude in the left and right column, respectively.Two indicated phases (H+P, P+H, H+H) coexist in other regions.The triple point corresponds to the tie line (thin gray line), where fractions φ of coexisting homogeneous and patterned phases are marked by brown and blue disks, respectively.(C) Phase diagram as a function of φ, χ, and E. The binodal of the homogeneous phase (brown surface) and the patterned phase (blue surface) overlap in the continuous phase transition (red surface).The critical points in panel B now correspond to critical lines, which all merge in the tricritical point (large black disk).A rotating version of the diagram is available as a movie.(A-C) Model parameters are ϕ0 = 1 and κ = 0.05 ξ 2 .
Continuous phase transition recovers experimental measurements.Squared amplitude (panel A) and length scale (panel B) of periodic patterns as a function of in- teraction strength χ for various parameters indicated in panel B, ϕ0 = 1, and κ = 0.05 ξ 2 .The amplitude indicates a continuous (colored data) and first-order (gray data) transition.