Chemically Active Mullins–Sekerka Models

• Chemically active Mullins–Sekerka models are sharp-interface, reaction–diffusion formulations that incorporate bulk and interfacial chemical reactions to drive droplet growth, division, and pattern formation.
• They couple governing equations, including Gibbs–Thomson and normal flux conditions, to describe the evolution and stability of interfaces under nonlinear reaction kinetics and capillarity.
• Numerical schemes using parametric finite element methods validate the models by capturing complex behaviors such as droplet splitting, multilayer shell formation, and interface instabilities.

Chemically active Mullins–Sekerka models are reaction-diffusion-driven sharp-interface formulations that extend the classical Mullins–Sekerka problem by incorporating both bulk and interfacial chemical reactions. These extensions are motivated by scenarios such as the growth, division, and compartmentalization of active droplets, with direct conceptual relevance for protocell models and the organization of chemical compartments in living systems. The interplay between diffusive transport, reaction kinetics, capillarity, and interface motion enables a rich variety of dynamical phenomena, including steady growth, interface instabilities, splitting, droplet division, and multilayer (shell) patterns (Garcke et al., 16 Jan 2026).

1. Sharp-Interface Formulation and Model Overview

The chemically active Mullins–Sekerka model describes a time-dependent partition of a domain $\Omega \subset \mathbb{R}^d$ $(d=2,3)$ into a droplet region $\Omega^+(t)$ and its complement $\Omega^-(t)$, separated by a moving hypersurface $\Sigma(t)$. The governing equations are:

• Bulk reaction–diffusion (quasi-static):

$$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$

with $m_\pm > 0$ (diffusivities), $\rho_\pm > 0$ (sink rates), $S_+ < 0$, $S_- > 0$ (source/sink).

• Gibbs–Thomson condition at $(d=2,3)$0:

$(d=2,3)$1

where $(d=2,3)$2 is the mean curvature, $(d=2,3)$3 (capillary length).

• Normal flux jump and interfacial reaction:

$(d=2,3)$4

$(d=2,3)$5 is the normal velocity of the interface, $(d=2,3)$6 is an interfacial source or sink.

• Boundary and dissipation: $(d=2,3)$7 on $(d=2,3)$8. Under suitable conditions, total interface length (2D) or area (3D) dissipates monotonically [(Garcke et al., 16 Jan 2026), (2.1a–e)].

2. Non-Dimensionalization and Reduced System

Characteristic scales are introduced: $(d=2,3)$9 (length), $\Omega^+(t)$0 (time), $\Omega^+(t)$1 (potential). Dimensionless parameters:

• $\Omega^+(t)$2
• $\Omega^+(t)$3
• $\Omega^+(t)$4
• $\Omega^+(t)$5
• $\Omega^+(t)$6

The reduced, dimensionless sharp-interface system is:

$\Omega^+(t)$7

Key dynamical parameters (e.g., $\Omega^+(t)$8) are generically $\Omega^+(t)$9 after rescaling [(Garcke et al., 16 Jan 2026), (2.3a–e)].

3. Analytical Structure: Radial and Planar Steady States, Stability

Radially Symmetric Solutions

For a spherical droplet of radius $\Omega^-(t)$0 in $\Omega^-(t)$1, the problem reduces to matching radial solutions (using modified Bessel functions) in $\Omega^-(t)$2 and $\Omega^-(t)$3, yielding:

$\Omega^-(t)$4

where $\Omega^-(t)$5 is an explicit but nontrivial function of $\Omega^-(t)$6, model parameters, and combinations of Bessel functions [(Garcke et al., 16 Jan 2026), (3.11)]. Real solutions $\Omega^-(t)$7 of $\Omega^-(t)$8 correspond to steady radii. Typically, for suitable parameters, there exist two roots, the smaller being radially unstable, the larger stable.

Linear Stability Analysis

Spherical or planar interfaces are subjected to infinitesimal harmonic perturbations (e.g., in spherical harmonics $\Omega^-(t)$9 or Fourier $\Sigma(t)$0). For the spherical case,

$\Sigma(t)$1

where $\Sigma(t)$2 is an explicit function of $\Sigma(t)$3, parameters, and combinations of Bessel functions, and $\Sigma(t)$4 is the perturbation mode number [(Garcke et al., 16 Jan 2026), (3.24), (4.15)]. Positive $\Sigma(t)$5 signals instability to mode $\Sigma(t)$6. A similar analysis for planar fronts yields a dispersion relation

$\Sigma(t)$7

recovering a chemically generalized Mullins–Sekerka spectrum [(Garcke et al., 16 Jan 2026), (4.24)].

4. Multilayer (“Shell”) Solutions and Their Stability

Multilayer configurations correspond to multiple interfaces—e.g., concentric shells (radial) or parallel planar fronts. The system reduces to a finite-dimensional system for the interface locations $\Sigma(t)$8:

$\Sigma(t)$9

For the planar double-front case, the evolution has a strict gradient-flow structure with a convex energy:

$$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$0

The unique symmetric minimizer yields stable paired fronts. Generalization to radial shells leads to analogous ODEs with stability determined by the eigenvalues of the linearized matrix [(Garcke et al., 16 Jan 2026), (4.6), (4.25a–b), (4.30)].

5. Numerical Scheme: Parametric Finite Element Method

The numerical method utilizes a parametric finite element approach:

• Interface representation: $$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$1 is discretized as a polyhedral mesh.
• Bulk discretization: Standard $$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$2 FEM on a (generally unfitted) triangulation.
• Coupled advance: At each timestep, the following linear system is solved, including curvature-driven interface updates and reaction/bulk source terms:

$$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$3

• Curvature and normal update: Coupled linear conditions for $$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$4, interface curvature $$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$5, and position $$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$6.
• Adaptive refinement is employed near $$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$7; topological changes (splitting, merging) are robustly detected and remeshed following [Benninghoff–Garcke ’17].
• Solvers: Direct (UMFPACK) in 2D or Schur-preconditioned CG in 3D. The scheme is unconditionally stable and energy-dissipative [(Garcke et al., 16 Jan 2026), (5.1a–c)].

6. Representative Computational Results

Computational studies validate the sharp-interface theory and reveal a spectrum of dynamical behaviors:

• Single planar fronts track the ODEs for $$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$8 exactly.
• Sinusoidal perturbations display mode selection: final patterns correspond to the fastest growing $$-m_\pm \Delta \mu_\pm = S_\pm - \rho_\pm \mu_\pm \quad \text{in } \Omega^\pm(t)$$9.
• Planar multilayer (strip) evolutions match the $m_\pm > 0$0 ODE structure and select the leading unstable Fourier mode.
• Droplet evolution in $m_\pm > 0$1 quantitatively traces the radial ODE.
• Radial shells converge to double-layer steady states if parameters are in the stable regime.
• Droplets subject to instability split into $m_\pm > 0$2-lobes for $m_\pm > 0$3, matching the dispersion relation.
• In 3D, an off-center “cigar” droplet can exhibit up to 23 pinch-offs, resulting in 24 daughter droplets.
• Vesicle models with curvature, surface viscosity, and surfactant coupling demonstrate shell and pearling instabilities concordant with linear theory [(Garcke et al., 16 Jan 2026), Figs. 5.1–5.14].

7. Significance and Applications

Chemically active Mullins–Sekerka models capture phenomena inaccessible to classical interface dynamics, including chemically driven droplet division, compartmentalization, and multilayer shell formation. These behaviors are considered relevant for realistic protocell models and may underpin fundamental mechanisms in the emergence and organization of life, where compartmentalization and growth-division cycles play a key organizing role (Garcke et al., 16 Jan 2026). The explicit coupling of sharp-interface motion, nonlinear chemical kinetics, and topological transition algorithms enables rigorous quantitative analysis and prediction of these complex pattern-forming dynamics.