Multi-Phase Field Model Overview

• Multi-phase field model is a diffuse-interface method representing multiple coexisting phases using smooth order parameters to capture complex interfaces.
• It couples interfacial energetics, bulk thermodynamics, and multiphysics such as mechanics, heat, and chemical reactions to simulate diverse phenomena.
• Modern implementations use variational formulations and advanced numerical strategies to handle topology changes and robust phase interactions.

Multi-phase field models are diffuse-interface formulations for systems with multiple coexisting phases, domains, fluids, grains, or cells. Rather than tracking sharp interfaces explicitly, they represent state by smooth order parameters, most classically as simplex-constrained phase fractions such as $\phi=(\phi_1,\dots,\phi_M)$ with $\sum_{\alpha=1}^M \phi_\alpha=1$ and $\phi_\alpha\ge0$, but also as cell-resolved fields $\phi_i(\mathbf{x},t)$ in tissues, as dual fields for coupled transition-and-fracture processes, or as the dichotomy-based vector $\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$ that avoids the simplex constraint altogether (Gräser et al., 2016, Monfared et al., 7 Mar 2025, Zhang et al., 9 Nov 2025). Across these variants, the central theme is the same: interfacial energetics, bulk thermodynamics, kinetics, and couplings to mechanics, heat transport, hydrodynamics, chemistry, or activity are encoded in a continuum free-energy or entropy structure that remains applicable to topology change, multiple junctions, and complex geometries (Toth et al., 2015, Dunbar et al., 2018).

1. Order-parameter representations and state spaces

A defining feature of the multi-phase field model is its choice of state variables. In classical multiphase formulations, the unknowns are local fractions of the competing phases. Representative notations include $\phi_\alpha$, $u_i$, and $p_i$, typically subject to a local sum constraint, for example

$\sum_{i=1}^N u_i(\mathbf r,t)=1$

or

$\sum_{i=1}^{N} p_i = 1,$

with $\sum_{\alpha=1}^M \phi_\alpha=1$0 on the Gibbs simplex [(Gräser et al., 2016); (Toth et al., 2015); (Pogorelov et al., 2013)]. In incompressible multiphase flow, the phase variables are often written as volume-fraction contrasts $\sum_{\alpha=1}^M \phi_\alpha=1$1, with $\sum_{\alpha=1}^M \phi_\alpha=1$2 and $\sum_{\alpha=1}^M \phi_\alpha=1$3 (Huang et al., 2020). In several solids models, the phase fields represent variants, grains, or precipitate states rather than volume fractions in a mixture, but the simplex structure remains central (Safi et al., 3 Mar 2025, Chatterjee et al., 2023).

That classical construction is not universal. A distinct line of work replaces the simplex-constrained $\sum_{\alpha=1}^M \phi_\alpha=1$4-tuple by $\sum_{\alpha=1}^M \phi_\alpha=1$5 independent scalar fields on the cube $\sum_{\alpha=1}^M \phi_\alpha=1$6, using a dichotomic or nested binary encoding of bulk phases. In that DBPF framework, the simplex constraint disappears, while mechanic consistency, energetic consistency, algebraic consistency, and dynamic consistency are imposed through interpolation rules for the surface-tension functions $\sum_{\alpha=1}^M \phi_\alpha=1$7 (Zhang et al., 9 Nov 2025). This addresses a longstanding modeling difficulty identified in earlier simplex-based theories: the construction of nonlinear potentials and the enforcement of consistency for many phases (Zhang et al., 9 Nov 2025, Toth et al., 2015).

A second major departure appears in biological tissue models, where each cell is represented by its own smooth phase field. In that setting, $\sum_{\alpha=1}^M \phi_\alpha=1$8 changes from “inside” to “outside” for cell $\sum_{\alpha=1}^M \phi_\alpha=1$9, so the model is a continuum–discrete hybrid: interfaces remain diffuse, but cells remain individually resolved (Monfared et al., 7 Mar 2025). Related monolayer models use one scalar field $\phi_\alpha\ge0$0 per cell, with $\phi_\alpha\ge0$1 inside and $\phi_\alpha\ge0$2 or $\phi_\alpha\ge0$3 outside depending on the convention, and exploit this representation to capture cell shape change, contact rearrangements, T1 transitions, and rosette formation without explicit vertex tracking (Wenzel et al., 2021, Chiang et al., 2024, Graham et al., 26 Aug 2025).

Many contemporary formulations are also explicitly multi-field rather than merely multi-phase. Frozen-soil modeling couples a freezing field $\phi_\alpha\ge0$4 and a damage field $\phi_\alpha\ge0$5 so that homogeneous pore freezing and segregated ice-lens growth remain distinct processes (Suh et al., 2021). Fluid-driven fracture in porous media similarly uses a saturation field $\phi_\alpha\ge0$6 for the invading-fluid/defending-fluid interface and a damage field $\phi_\alpha\ge0$7 for fracture (Guével et al., 2023). Hyperelastic multiphase solids with Eulerian interfaces couple a deformation $\phi_\alpha\ge0$8 to a phase/composition field $\phi_\alpha\ge0$9 defined on the deformed configuration $\phi_i(\mathbf{x},t)$0 rather than the reference domain (Grandi et al., 2020). This suggests that “multi-phase field model” denotes a family of diffuse-interface constructions rather than a single canonical PDE system.

2. Free-energy and entropy structures

The governing structure of a multi-phase field model is usually variational. In the non-isothermal multi-phase Penrose–Fife system, the model is built from a generalized entropy functional

$\phi_i(\mathbf{x},t)$1

where $\phi_i(\mathbf{x},t)$2 is a multi-well obstacle potential and the thermodynamic state couples phase fractions to inverse temperature $\phi_i(\mathbf{x},t)$3 (Gräser et al., 2016). In consistent multiphase theories for interface-driven dynamics, the analogous object is a free energy of the form

$\phi_i(\mathbf{x},t)$4

with $\phi_i(\mathbf{x},t)$5 designed so that binary interfaces remain equilibrium solutions and absent phases are not generated deterministically (Toth et al., 2015).

The interfacial part can be constructed in markedly different ways. Pairwise formulations write

$\phi_i(\mathbf{x},t)$6

with interface energy densities based on gradient terms and double-obstacle potentials (Schiedung et al., 2017). The N-dimensional extension of the Folch–Plapp program instead constructs barrier functions $\phi_i(\mathbf{x},t)$7 satisfying flatness and stability conditions on every dual interface $\phi_i(\mathbf{x},t)$8, specifically so that “ghost” phases do not appear on binary interfaces (Pogorelov et al., 2013). In DBPF models, the free energy is built by recursive interpolation of binary Ginzburg–Landau energies,

$\phi_i(\mathbf{x},t)$9

with the coupling encoded entirely through the interpolation functions $\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$0 (Zhang et al., 9 Nov 2025).

Mechanical and multiphysics variants enrich the free energy rather than replacing it. Surface-tension-induced elasticity uses

$\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$1

with an interfacial density

$\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$2

so that surface energy depends on tangential deformation and acts as a localized interfacial stress (Schiedung et al., 2017). Hyperelastic multiphase solids with Eulerian interfaces use

$\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$3

thereby measuring interfacial area in the deformed configuration rather than the reference configuration (Grandi et al., 2020). Tissue models add double-well interfacial terms, volume penalties, cell–cell repulsion, and cell–cell or cell–substrate adhesion to a cell-wise free energy $\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$4 (Monfared et al., 7 Mar 2025).

Thermodynamic consistency is a central criterion across the literature. The Penrose–Fife system proves entropy monotonicity in the no-source/no-flux setting, $\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$5 (Gräser et al., 2016). XMPF is constructed to ensure $\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$6 and non-negative entropy production (Toth et al., 2015). The DBPF model satisfies

$\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$7

(Zhang et al., 9 Nov 2025). Consistent incompressible multiphase flow adds a full kinetic-plus-free-energy law, while multi-phase surfactant flow yields a closed-system free energy that is nonincreasing under viscous and diffusive dissipation (Huang et al., 2020, Dunbar et al., 2018).

3. Evolution equations, constraints, and consistency requirements

The kinetic equations vary with the physical transport mechanism. Nonconserved phase evolution often appears in Allen–Cahn or Ginzburg–Landau form. Examples include the elastic capillarity model,

$\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$8

and many solids models in which order parameters evolve by variational relaxation (Schiedung et al., 2017, Safi et al., 3 Mar 2025). Conserved dynamics, by contrast, appear in capillarity-driven interface diffusion,

$\mathbf{\Phi}^N=(\phi_1,\dots,\phi_{N-1})$9

and in Cahn–Hilliard-type DBPF systems,

$\phi_\alpha$0

(Schiedung et al., 2017, Zhang et al., 9 Nov 2025). Biological tissue models typically combine advection and relaxation,

$\phi_\alpha$1

with overdamped translational force balance for cell motion (Monfared et al., 7 Mar 2025).

Consistency conditions became a major theoretical issue in generalized multiphase formulations. XMPF formalized a set of requirements including the local sum constraint, formal indistinguishability, equilibrium consistency, thermodynamic consistency, reducibility from $\phi_\alpha$2 to $\phi_\alpha$3 phases, no spurious phase generation, and independent pairwise interfacial and kinetic data (Toth et al., 2015). The N-dimensional Folch–Plapp extension imposed flatness and stability on dual interfaces and at pure-phase vertices: $\phi_\alpha$4 specifically to eliminate “ghost” phases (Pogorelov et al., 2013). The DBPF formulation restated this agenda as mechanic consistency, energetic consistency, algebraic consistency, and dynamic consistency, now without a simplex constraint (Zhang et al., 9 Nov 2025).

Sharp-interface asymptotics provide the principal validation of these constructions. The conserved surface/phase-boundary diffusion model recovers the classical law

$\phi_\alpha$5

in the sharp-interface limit and reproduces the von Neumann relation at triple junctions (Schiedung et al., 2017). The DBPF ternary asymptotics recover curvature laws on each interface and the Neumann triangle condition

$\phi_\alpha$6

(Zhang et al., 9 Nov 2025). The general $\phi_\alpha$7-phase model in $\phi_\alpha$8-dimensional phase-field space verifies Young’s law for three and four phases with high accuracy (Pogorelov et al., 2013). Surfactant-coupled multiphase flow derives a sharp moving-boundary problem with continuity of chemical potential at triple junctions and Young’s law modified by surfactant-dependent surface tensions (Dunbar et al., 2018).

A further generalization concerns kinetics itself. In ferroelectrics, the classical Allen–Cahn assumption implies linear interface kinetics at small driving force. The multi-phase-field model with general kinetics replaces that restriction by a phase-fraction evolution law containing prescribed monotone odd functions $\phi_\alpha$9 and $u_i$0, so that 180° and 90° domain walls propagate with arbitrarily chosen nonlinear kinetic relations (Guin et al., 2022). This is not a change in interfacial regularization so much as a change in what is encoded as the sharp-interface limit.

4. Couplings to heat, mechanics, hydrodynamics, chemistry, and activity

Thermal coupling is explicit in non-isothermal phase-change models. In the multi-phase Penrose–Fife system, the unknowns are the multi-phase field and the positive inverse temperature field $u_i$1, and the weak formulation couples phase evolution and heat transport through the terms $u_i$2 and $u_i$3 (Gräser et al., 2016). The coupling is physically transparent: phase fractions respond to thermal state and latent heat, while temperature responds to phase change through latent heat release or absorption (Gräser et al., 2016). Frozen-soil modeling uses an analogous but distinct two-field strategy, with a freezing order parameter $u_i$4 for water–ice transition and a damage field $u_i$5 for fracture, so that homogeneous freezing and ice-lens growth remain separate mechanisms in a coupled thermo-hydro-mechanical setting (Suh et al., 2021).

Mechanical coupling appears in several forms. Surface-tension-induced elasticity derives a total stress

$u_i$6

so that surface or interface energy acts as a tangential interfacial stress whose divergence yields the curvature force (Schiedung et al., 2017). Coherently stressed three-phase solids couple phase evolution to elasticity through a partial rank-one homogenization scheme that enforces static and kinematic compatibilities in diffuse interfacial regions, in contrast to simpler Voigt-Taylor interpolation (Chatterjee et al., 2023). T$u_i$7 precipitate evolution in Al–Cu–Li alloys combines CALPHAD-based chemistry, anisotropic interfacial energy, variant-specific Allen–Cahn kinetics, Cahn–Hilliard solute diffusion, and an FFT-based spectral elasticity solve, concluding that anisotropy in interfacial energy and linear reaction rate dominates shape evolution while elastic effects are minor under the conditions studied (Safi et al., 3 Mar 2025).

Hydrodynamic coupling is especially intricate in incompressible multiphase flow. A key result is that the momentum equation should not use $u_i$8 as the inertial flux when diffusive phase transport is present; instead it should use the consistent mass flux

$u_i$9

so that

$p_i$0

(Huang et al., 2020). This enforces consistency of mass conservation and of mass and momentum transport, yielding Galilean invariance and kinetic-energy consistency (Huang et al., 2020). Surfactant-coupled multiphase flow extends this picture by adding a diffuse surfactant balance, surfactant-dependent surface tension, Marangoni forcing, and triple-junction transport under local chemical equilibrium (Dunbar et al., 2018). Porous-media fracture couples Darcy flow, a Cahn–Hilliard saturation field, and phase-field damage so that capillary/Korteweg stresses, permeability evolution, and fracture nucleation interact within a single thermodynamic framework (Guével et al., 2023).

Biological applications add active matter ingredients rather than classical constitutive couplings. Tissue models separate translational overdamped motion from interface relaxation and include cell cortex tension, volume constraint, cell–cell repulsion, cell–cell adhesion, cell–substrate interaction, self-propulsion, polarity alignment, nematic contractility, and shape-based active stress (Monfared et al., 7 Mar 2025). Comparative studies of collective migration show that random-orientation, elongation-based, polar, and nematic activity prescriptions can all produce solid-to-liquid transitions, defect dynamics, and active-flow patterns, but with qualitatively different rosette statistics, defect motion, and confinement response (Wenzel et al., 2021). A related monolayer theory derives an effective 2D model from a 3D thin-layer description and shows how intercellular friction emerges from a viscous stress as a pairwise relative-velocity coupling (Chiang et al., 2024). The same framework has been extended to stochastic junctional adhesion, where pairwise adhesion coefficients follow an Ornstein–Uhlenbeck process and drive T1 rearrangements, fluidization, and diffusive cell motion (Graham et al., 26 Aug 2025).

5. Numerical approximation and computational strategies

Because the governing equations span variational inequalities, conserved and nonconserved gradient flows, incompressible hydrodynamics, and nonlinear elasticity, numerical strategy is a major part of the field. In the Penrose–Fife system, Rothe’s method is used: the model is discretized implicitly in time and with linear finite elements in space, producing a stepwise saddle-point problem in the multi-phase field and inverse temperature (Gräser et al., 2016). The resulting algebraic variational inequality is reformulated as a dual minimization problem in the temperature variable and solved by a non-smooth Schur–Newton method. Each nonlinear iteration evaluates the constrained phase minimization by truncated non-smooth Newton multigrid and solves the generalized Schur-complement step by multigrid with a Vanka-type smoother or by preconditioned GMRES (Gräser et al., 2016). The reported behavior is effectively mesh independent, and adaptive refinement concentrates nodes around phase boundaries and thermal gradients (Gräser et al., 2016).

Finite-difference and spectral approaches are equally prominent. Surface/phase-boundary diffusion is implemented in OpenPhase with a 27-point stencil for the Laplacian, a 3-point first-order central finite difference scheme for gradients and divergences, and explicit Euler time stepping; because the double-obstacle potential is non-smooth at $p_i$1, a bent-cable model is used to interpolate derivatives smoothly near the minima (Schiedung et al., 2017). The T$p_i$2 precipitate model uses semi-implicit Fourier-spectral updates for both logarithmic composition variables and Allen–Cahn order parameters, together with an FFT-based spectral elasticity solver (Safi et al., 3 Mar 2025). The multi-component phase field crystal model uses semi-implicit Fourier-space time stepping for density and concentration fields (Ofori-Opoku et al., 2012).

Recent multiphase formulations often aim for linearity, decoupling, and provable stability. The DBPF paper develops a mobility operator splitting framework and, for the ternary case, a second-order, linear, decoupled, energy-stable scheme based on three Strang-composed substeps and a modified Crank–Nicolson-type linearization (Zhang et al., 9 Nov 2025). The consistent and conservative multiphase incompressible flow scheme instead focuses on preserving structure after discretization: a gradient-based phase selection procedure maintains the summation constraint and suppresses fictitious phases in convection; balanced-force and conservative surface-force discretizations preserve, respectively, low spurious currents and exact discrete momentum conservation; and discrete consistency theorems ensure mass conservation, reduction consistency, and the use of exactly the same discrete transport in the phase and momentum equations (Huang et al., 2020).

Large coupled systems are frequently implemented in general-purpose PDE frameworks. The frozen-soil model uses Taylor–Hood elements for $p_i$3, linear elements for $p_i$4, implicit backward Euler time integration, and an operator-split staggered strategy in FEniCS with PETSc (Suh et al., 2021). The Darcy–Cahn–Hilliard fracture model is implemented in MOOSE with RACCOON, uses backward Euler, Newton–Raphson, and a primal-dual active set for damage irreversibility (Guével et al., 2023). The mechanically compatible three-phase solid model is solved in MOOSE with CALPHAD-precomputed thermodynamic and kinetic data (Chatterjee et al., 2023). Composite-laminate fracture is implemented in Abaqus UEL with a block-diagonal Newton solve for displacement and the two damage fields (Kumar et al., 10 Mar 2026). The 3D-to-2D monolayer reduction uses finite differences, fourth-order central differences for gradients, a nine-point Laplacian, a third-order upwind advection scheme, LAPACK for the velocity matrix, and OpenMP parallelization (Chiang et al., 2024).

6. Application domains, limitations, and current directions

The application range of multi-phase field models is unusually broad. In materials processing, they have been used for grain growth and liquid phase crystallization of silicon thin films (Gräser et al., 2016), thermal grooving and annealing of multi-nano-clusters on deformable free surfaces (Schiedung et al., 2017), surface-tension-induced stress in spherical heterogeneities, thin plates, ellipsoids, and sintered structures (Schiedung et al., 2017), multicomponent phase separation and grain coarsening in polycrystals (Toth et al., 2015), structural transformations and precipitation in multicomponent alloys at the phase-field-crystal level (Ofori-Opoku et al., 2012), coherently stressed three-phase intermetallic microstructures in Ni–Al and Al–Cr–Ni (Chatterjee et al., 2023), stoichiometric precipitate growth in Al–Cu–Li (Safi et al., 3 Mar 2025), ferroelectric domain-wall motion with nonlinear kinetics (Guin et al., 2022), and intralaminar failure of fiber-reinforced composite laminates using separate fiber and inter-fiber phase fields (Kumar et al., 10 Mar 2026). In fluids and porous media, the same general methodology supports incompressible multiphase flow at large density ratios (Huang et al., 2020), surfactant transport through triple junctions (Dunbar et al., 2018), fluid-driven fracture (Guével et al., 2023), and ice-lens growth and thaw in frozen soils (Suh et al., 2021).

Biological uses are equally extensive. Reviews of dense, soft multicellular systems describe applications to embryogenesis, morphogenesis, wound repair, tumor invasion, cell extrusion, collective migration, heterogeneous cell populations, and confined systems (Monfared et al., 7 Mar 2025). Monolayer studies reproduce solid-to-liquid transitions, cell-shape variability, nematic properties, vorticity correlations, and confinement-dependent flows across several activity prescriptions (Wenzel et al., 2021). Thin-layer reductions provide a derivation of the commonly used 2D monolayer model from a 3D continuum description on a substrate and show that passive deformation stresses and cell–cell friction both tend to solidify the tissue (Chiang et al., 2024). Junctional-fluctuation models demonstrate that stochastic pairwise adhesions can fluidize an epithelial sheet and produce a non-monotonic dependence of the effective diffusion coefficient on the persistence time of adhesion fluctuations (Graham et al., 26 Aug 2025).

Several recurring limitations are also explicit in the literature. Some capillarity–elasticity models assume constant surface/interface energy and linear elasticity, and diffuse-interface deviations from sharp-interface benchmarks become noticeable when the characteristic size approaches the interface width (Schiedung et al., 2017). XMPF emphasizes that a fully general proof for all higher-order equilibria remains incomplete, even though the formulation satisfies the main consistency criteria and performs well numerically (Toth et al., 2015). The 2D laminate-fracture model does not include delamination, plasticity, or through-thickness effects (Kumar et al., 10 Mar 2026). Biological reviews identify open problems in biochemical signaling, nuclear mechanics, poroelasticity, water transport, interstitial flow, consistent treatment of active, thermal, and biochemical noise, and scaling to large 3D tissues and tumors (Monfared et al., 7 Mar 2025).

Current research directions therefore combine formal consistency, richer constitutive coupling, and computational scalability. One trajectory relaxes the classical simplex constraint through dichotomic encodings (Zhang et al., 9 Nov 2025). Another strengthens mechanical compatibility in interfacial regions through homogenization schemes tailored to coherent multiphase solids (Chatterjee et al., 2023). A third pushes toward full 3D biological and fluid-mechanical couplings, mechanochemical feedback, and data-integrated calibration (Monfared et al., 7 Mar 2025). A plausible implication is that the modern multi-phase field model is less a single method than a common thermodynamic language for diffuse-interface problems in which interfacial geometry, multiple junctions, and strongly coupled bulk physics must be resolved simultaneously.