Diffuse-Interface Models (DIMs)

• Diffuse-interface models are frameworks that describe multi-phase systems by representing interfaces as smooth regions of finite thickness using phase-field variables.
• They employ free-energy functionals with higher-order derivative terms to regulate interfacial energy and ensure numerical convergence across different codimensions.
• Applications include multi-phase flows, fracture mechanics, and material microstructure evolution, offering versatile tools for simulating complex physical processes.

A diffuse-interface model (DIM) is a mathematical framework for describing the evolution of multi-phase systems in which interfaces between phases are represented as regions of finite thickness, rather than as sharp boundaries. DIMs are formulated by introducing one or more continuous phase-field variables whose spatial gradients or higher-order derivatives encode interfacial energy and structure. Critical applications include multi-phase flows, fracture, elastic breakdown, vesicle dynamics, reactive flows, and material microstructure evolution. DIMs are characterized by free-energy functionals penalizing spatial inhomogeneity, regularized through higher-order or nonlinear gradient terms as required by the nature of the dispersed phase.

1. Free-Energy Functionals and Governing Equations

DIMs are grounded in variational principles, typically through a free-energy functional of the form: $$F[\phi] = \int_\Omega W(\phi,\nabla\phi, \ldots)\,d\Omega$$ where $\phi$ is the phase field (or a vector of fields for multi-component or multi-phase systems), and $W$ encapsulates bulk and interfacial energetics. The classical Ginzburg–Landau (Cahn–Hilliard) functional includes terms such as $|\nabla\phi|^2$ and a double-well potential $\psi(\phi)$, producing diffuse interfaces whose thickness and surface energy are controlled by the coefficients and scaling parameters.

Extensions to higher codimension inclusions (e.g., lines in 3D, points in 2D) necessitate inclusion of higher derivative terms such as $(\Delta\phi)^2$ or nonlinear gradients $|\nabla\phi|^p$ $(p>2)$ to ensure well-posedness and physical locality of the interface (Zipunova et al., 2020). The Euler–Lagrange equations derived from such functionals define phase-field kinetics and equilibrium, often resulting in coupled systems involving mechanics (Navier–Stokes, elasticity), transport, and (electro-/magneto-)statics.

2. Regularity, Codimension Effects, and Sobolev Embedding

A central theoretical insight is that the regularity imposed by DIMs depends critically on the codimension of the dispersed inclusion. For codimension-1 objects (surfaces in 3D), the classical $|\nabla\phi|^2$-based energy ensures trace properties strong enough to impose boundary conditions on interfaces. For codimension-2 inclusions (curves in 3D, lines in 2D), the function $\phi$ must either belong to $W^{2,2}(\Omega)$ (via a biharmonic penalty, i.e., $\int(\Delta\phi)^2d\Omega$) or to $W^{1,p}(\Omega)$ with $p>2$ (e.g., $\int|\nabla\phi|^p d\Omega$ with $p\geq4$). This arises from Sobolev embedding theorems and capacity theory: a 1D curve in 3D has zero $W^{1,2}$-capacity, but positive $W^{2,2}$- and $W^{1,4}$-capacity (Zipunova et al., 2020).

Table: Required regularity of phase-field energy terms by codimension (in 3D):

Inclusion type	Energy term required	Sobolev space
Surface (codim 1)	$|\nabla\phi|^2$	$W^{1,2}(\Omega)$
Curve (codim 2)	$(\Delta\phi)^2$ or $|\nabla\phi|^p$ ($p>2$)	$W^{2,2}$ or $W^{1,p}$ ($p>2$)

If higher regularity is not enforced, minimizers fail to resolve the interface: under mesh refinement, the phase-field narrows to a spike, and boundary conditions on the inclusion cannot be imposed in a physically or mathematically meaningful way (Zipunova et al., 2020).

3. Model Construction and Corrected Energy Densities

A general corrected DIM free-energy density for codimension-2 inclusions takes the form: $$\pi(\Phi, \phi, \nabla\phi, \Delta\phi) = -\frac12\epsilon[\phi]|\nabla\Phi|^2 + \frac{\Gamma}{l^2}(1-f(\phi)) + \frac{\Gamma}{4}|\nabla\phi|^2 + \alpha \frac{\Gamma l^2}{8}(\Delta\phi)^2 + \beta \frac{\Gamma l^{p-2}}{p}|\nabla\phi|^p,$$ where $\alpha,\,\beta\geq0$, $p\geq4$, $f(\phi)$ is an interpolation function (typically a cubic/quartic with $f(0)=0,\ f(1)=1$), $\Gamma$ is a line (or higher-codimension) energy density, and $l$ is the interface thickness (Zipunova et al., 2020). Nonzero values of $\alpha$ or $\beta$ (with $\alpha+\beta>0$) are required to prevent collapse of the diffuse region for codimension-2 inclusions.

The corresponding Euler–Lagrange equations include higher-order spatial derivatives, yielding, in equilibrium, a fourth-order semilinear elliptic PDE for the phase field. Kinetics are typically modeled by Allen–Cahn or Cahn–Hilliard evolution, where the regularization terms ensure well-posedness and mesh-convergence of diffuse interfaces even for high codimension objects.

4. Numerical Implementation and Convergence

Numerical schemes for DIMs may employ finite differences, finite elements, or spectral methods, depending on the context. In the axisymmetric case, the governing equations reduce to nonlinear ODEs for the phase field, which can be efficiently solved via Newton's method (Zipunova et al., 2020). Crucially, computational experiments demonstrate that uncorrected models ($\alpha=\beta=0$) are not mesh-convergent—solutions degenerate to singular spikes as the grid is refined. When higher-regularity terms are present (e.g., $\alpha=0.1$, $\beta=0.01$, $p=4$), the phase field remains monotone, the diffuse boundary is well resolved, and the mesh-converged solution retains the expected finite thickness $O(l)$.

A parametric study in the regularization parameters reveals a smooth transition between models regularized purely by a biharmonic term (linear, $\alpha>0,\beta=0$) and models regularized by a nonlinear $p$-Laplacian term ($\alpha=0,\beta>0$), offering flexibility for linear vs. nonlinear numerical solution strategies.

5. Physical and Mathematical Implications

The codimension dependence inherent to DIMs has strong consequences for the design of phase-field models in physics and engineering. For classical examples—interfaces (codimension-1) such as cracks or phase boundaries in fracture mechanics, or phase separation in alloys—the standard gradient-squared regularization is sufficient. For problems involving line singularities (e.g., dielectric breakdown channels, 1D cracks in 2D), models must include either a $(\Delta\phi)^2$ term or a $|\nabla\phi|^p$ term ($p>2$) to retain the ability to impose physically meaningful boundary conditions and to avoid collapse of the diffuse region (Zipunova et al., 2020).

This constraint is also visible in the sharp-interface (Γ-convergence) limit: only enriched energies (with higher-order or higher-power derivatives) produce sharp-limit functionals whose interfacial (surface) terms are proportional to the Hausdorff measure $H^{n-2}(\Gamma)$. From the standpoint of capacity theory, objects of codimension-$k$ in $n$ dimensions require control over higher-order Sobolev norms for the phase-field variable, with the regularizing energy scaling determined by the codimension.

6. Guidelines for Model Selection and Practical Trade-offs

For applications requiring diffuse interface representations of objects with codimension $>1$, the choice between biharmonic (linear, $(\Delta\phi)^2$) and p-Laplacian (nonlinear, $|\nabla\phi|^p,\ p>2$) regularization is dictated by solver capabilities and desired phase-field smoothness (Zipunova et al., 2020). The well-posedness and convergence of the model, as well as the physical interpretability of imposed boundary data, depend entirely on these enriched energy terms.

In practice:

• Codimension-1: Use classic Ginzburg–Landau energy with $|\nabla\phi|^2$.
• Codimension-2: Introduce $(\Delta\phi)^2$ and/or $|\nabla\phi|^p$, $p\geq4$.
• Choose regularization strategy in accordance with numerical solvers (linear vs. nonlinear).
• Ensure that the mesh resolves the interface thickness $l$; this is a physical and numerical parameter.

This framework generalizes DIMs beyond standard applications, enabling well-posed and convergent modeling of high-codimension inclusions across physics, materials science, and engineering contexts (Zipunova et al., 2020).