Diffuse Domain Method (DDM)

• Diffuse Domain Method (DDM) is an unfitted domain-embedding technique that replaces sharp interfaces with a smooth, diffuse layer for handling complex geometries.
• It reformulates partial differential equations onto a simple computational domain using weighted integrals and asymptotic analysis, yielding controlled error rates linked to the layer thickness.
• DDM accommodates various boundary conditions and problem types—including Neumann, Robin, and Dirichlet—and extends to transmission, moving interfaces, and multiphysics, supported by rigorous error estimates.

Searching arXiv for recent and foundational papers on diffuse domain methods to ground the article in the current literature. The diffuse domain method (DDM) is a class of unfitted domain-embedding techniques for partial differential equations on geometrically complicated domains. Its central operation is to replace a sharp physical domain $D$ or interface $\partial D$ by a thin diffuse layer of thickness $\varepsilon$, represented by a phase-field or smoothed characteristic function, and to reformulate the original sharp-interface problem on a larger, simple computational domain such as a rectangle or box. In this reformulation, bulk integrals are weighted by an approximation of $\chi_D$, while boundary or interface terms are converted into volumetric terms concentrated in the diffuse layer; this makes the method compatible with standard Cartesian finite differences, rectangular finite elements, and related solvers on fixed meshes (Burger et al., 2014, Hao et al., 3 Apr 2025, Xu, 13 May 2026).

1. Geometric embedding and diffuse reformulation

In the DDM literature, the physical domain is typically represented by a signed or oriented distance function. For a bounded domain $D\subset\mathbb R^d$, one common definition is

$$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$

A diffuse indicator is then introduced through a sigmoidal profile. One formulation uses

$$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$

with examples such as $S(t)=\tanh(t)$ or $S(s)=\tanh(3s)$. As $\epsilon\to 0$, $\partial D$0 converges formally to the indicator $\partial D$1, while the transition from $\partial D$2 to $\partial D$3 is confined to a layer of thickness $\partial D$4 (Hao et al., 3 Apr 2025, Xu, 20 Oct 2025). A closely related formulation uses

$\partial D$5

where $\partial D$6 is the signed distance to the boundary, negative inside the physical domain (Lervåg et al., 2014).

The diffuse reformulation rests on two standard replacements. Volume integrals over $\partial D$7 are approximated by weighted bulk integrals,

$\partial D$8

and boundary or interface integrals are approximated by layer integrals involving the diffuse surface measure,

$\partial D$9

Equivalently, in phase-field notation, $\varepsilon$0 or $\varepsilon$1 can act as regularizations of the surface delta distribution $\varepsilon$2 (Hao et al., 3 Apr 2025, Lervåg et al., 2014).

For the second-order parabolic problem

$\varepsilon$3

the sharp weak form

$\varepsilon$4

is replaced by the diffuse variational problem

$\varepsilon$5

Because $\varepsilon$6 vanishes outside $\varepsilon$7, the same problem can be extended to a fixed rectangular domain $\varepsilon$8, which is the main practical reason DDM avoids body-fitted meshes (Hao et al., 3 Apr 2025).

A parallel strong-form presentation appears in asymptotic analyses of elliptic problems. For Poisson-type equations, the extended-domain form is written as

$\varepsilon$9

or, for reaction-diffusion,

$\chi_D$0

where $\chi_D$1 is a diffuse approximation of the original boundary condition (Lervåg et al., 2014). This same structural pattern reappears in later DDM variants for transport, phase-field, fluid, and interface systems.

2. Functional-analytic framework and weighted Sobolev structure

A defining mathematical feature of modern DDM analysis is the use of weighted Sobolev spaces induced by the diffuse indicator. For $\chi_D$2,

$\chi_D$3

and

$\chi_D$4

In the Hilbert case, one writes $\chi_D$5, with weighted inner product $\chi_D$6 (Burger et al., 2014, Hao et al., 3 Apr 2025).

This weighted setting is not merely notational. Near the diffuse boundary, the weight degenerates in a controlled way, and the natural coercivity, trace, and embedding estimates are uniform in $\chi_D$7. Foundational results include the weighted trace estimate

$\chi_D$8

the weighted Sobolev embedding

$\chi_D$9

and Poincaré–Friedrichs-type inequalities of the form

$D\subset\mathbb R^d$0

These results underpin well-posedness of diffuse variational problems for elliptic and parabolic PDEs and make it possible to compare diffuse and sharp solutions in norms natural to the weighted geometry (Burger et al., 2014, Xu, 20 Oct 2025).

A second core analytic ingredient is a family of diffuse-versus-sharp integral consistency estimates. Representative examples are

$D\subset\mathbb R^d$1

and, for suitable $D\subset\mathbb R^d$2 vanishing on $D\subset\mathbb R^d$3,

$D\subset\mathbb R^d$4

For elliptic problems, analogous estimates for diffuse volume and boundary integrals yield $D\subset\mathbb R^d$5 or $D\subset\mathbb R^d$6 consistency depending on available derivatives (Hao et al., 3 Apr 2025, Burger et al., 2014).

This framework explains why many rigorous DDM results are stated first in weighted norms. Unweighted conclusions on the original domain are often recovered by restriction, but the weighted setting is the native analytic environment of the diffuse approximation.

3. Boundary conditions, boundary representations, and accuracy mechanisms

A recurring theme in DDM research is that the treatment of boundary conditions is structurally decisive. The geometric embedding by itself does not determine the asymptotic order; the diffuse representation of the boundary operator does.

For Neumann and Robin conditions in elliptic DDMs, two classical approximations are

$D\subset\mathbb R^d$7

for Neumann, and

$D\subset\mathbb R^d$8

for Robin problems. Matched asymptotic analysis shows that for Neumann problems both BC1 and BC2 are second-order accurate in $D\subset\mathbb R^d$9, while for Robin problems BC1 is second-order but BC2 is only first-order because the first correction $$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$0 is generally nonzero: $$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$1 That classification is one of the central asymptotic results in the DDM literature (Lervåg et al., 2014).

The same paper proposed corrected Robin variants,

$$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$2

designed to cancel the $$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$3 defect and restore second-order accuracy. More local corrections, such as BC1M1 and BC1M2, were also derived, though some of these may compromise ellipticity outside the physical domain (Lervåg et al., 2014).

Dirichlet conditions are more delicate. Standard DDM formulations include

$$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$4

$$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$5

$$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$6

Matched asymptotics for Poisson problems with Dirichlet data show that DDM2 is generally $$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$7, whereas DDM1 and DDM3 are generally $$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$8 when the boundary normal derivative is nonzero. Modified forms using the first-order normal correction

$$d_D(\bm x)=(\bm x,D)-(\bm x,\mathbb R^n\setminus D), \qquad D=\{d_D<0\}.$$9

yield

$$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$0

asymptotically when the thinner internal layer is resolved, while mDDM2 remains first-order (Yu et al., 2018).

A later line of work reformulated Dirichlet DDMs through mixed and Nitsche-type constructions. New methods such as Mix0DDM and Mix1DDM derive the diffuse formulation from the mixed relation $$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$1, thereby converting the essential Dirichlet condition into a natural one in the mixed system. Nitsche-based variants NDDM and NSDDM were also constructed, with coercivity proofs for the new and key existing approximations. Numerical experiments reported that NSDDM almost always gave the best $$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$2 errors, whereas Mix0DDM gave the best $$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$3 errors (Benfield et al., 29 Sep 2025).

Boundary/data type	Representative diffuse form	Principal conclusion
Neumann	$$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$4, $$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$5	Both asymptotically second-order in the elliptic matched-asymptotic analysis
Robin	$$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$6, $$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$7	BC1 second-order; BC2 first-order unless corrected
Dirichlet	Penalty terms using $$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$8 or $$\varphi^\epsilon = S\!\left(-\frac{d_D}{\epsilon}\right), \qquad \omega^\epsilon=\frac{1+\varphi^\epsilon}{2},$$9	Standard methods are low-order; modified and Nitsche/mixed forms improve behavior

A common misconception is that DDM has a single accuracy order determined only by $S(t)=\tanh(t)$0. The literature does not support that view. Boundary representation, penalty scaling, two-sided versus one-sided embedding, and the norm in which error is measured all materially affect the observed and provable order (Lervåg et al., 2014, Yu et al., 2018, Benfield et al., 29 Sep 2025).

4. Discretization, DDFE formulations, and rigorous error estimates

Practical DDM solvers usually combine the diffuse reformulation with a standard discretization on a simple background mesh. A representative recent construction for parabolic equations is the diffuse domain finite element method (DDFE): the transformed weighted problem is posed on a fixed rectangle

$S(t)=\tanh(t)$1

discretized by a tensor-product finite element space of continuous piecewise multilinear basis functions,

$S(t)=\tanh(t)$2

on a quasi-uniform rectangular grid, together with BDF2 time stepping

$S(t)=\tanh(t)$3

The fully discrete recurrence is written as

$S(t)=\tanh(t)$4

The weight $S(t)=\tanh(t)$5 appears directly in the mass and stiffness forms, so the finite element method is applied to the diffuse weighted PDE rather than to the original sharp-domain problem (Hao et al., 3 Apr 2025).

For second-order linear parabolic equations with Neumann boundary conditions on irregular domains, rigorous error analysis separates three contributions: $S(t)=\tanh(t)$6 The diffuse-interface modeling error satisfies

$S(t)=\tanh(t)$7

the semidiscrete FE error satisfies

$S(t)=\tanh(t)$8

and, under the CFL-type condition

$S(t)=\tanh(t)$9

the temporal error satisfies

$S(s)=\tanh(3s)$0

Consequently, the fully discrete DDFE estimate is

$S(s)=\tanh(3s)$1

This is one of the clearest decompositions of modeling, spatial, and temporal error currently available for parabolic DDMs (Hao et al., 3 Apr 2025).

For semilinear parabolic equations with Neumann data, a later weighted-space analysis established the same optimal diffuse-interface rates

$S(s)=\tanh(3s)$2

under regularity assumptions on $S(s)=\tanh(3s)$3, $S(s)=\tanh(3s)$4, and the nonlinearity $S(s)=\tanh(3s)$5, together with a local weighted Lipschitz estimate for $S(s)=\tanh(3s)$6 (Xu, 20 Oct 2025).

For elliptic problems, the rigorous theory is older and more norm-dependent. In the weighted $S(s)=\tanh(3s)$7-framework for second-order elliptic boundary value problems, Robin and Neumann DDMs satisfy low-regularity estimates of order

$S(s)=\tanh(3s)$8

while in the smooth-data regime the weighted $S(s)=\tanh(3s)$9 error improves to $\epsilon\to 0$0. Dirichlet data are treated by a Robin penalty with parameter $\epsilon\to 0$1, leading to weaker estimates such as

$\epsilon\to 0$2

and, in the smooth case with $\epsilon\to 0$3,

$\epsilon\to 0$4

These results explain why elliptic Dirichlet DDMs are often analytically and numerically less straightforward than Neumann or Robin formulations (Burger et al., 2014).

5. Transmission problems, moving interfaces, and multiphysics generalizations

The DDM has expanded well beyond one-sided scalar elliptic and parabolic boundary value problems. A major direction concerns transmission problems across internal interfaces. In the two-sided formulation

$\epsilon\to 0$5

with

$\epsilon\to 0$6

matched asymptotic analysis for a one-dimensional Robin-type transmission problem shows that the diffuse approximation is generically exactly first-order in $\epsilon\to 0$7 when $\epsilon\to 0$8 is fixed. The first correction $\epsilon\to 0$9 is explicitly nonzero, for example through the jump relation

$\partial D$00

This result corrected a possible expectation, inherited from one-sided asymptotics, that second-order behavior would persist under practical two-sided stabilization (Luong et al., 2024).

A complementary variational approach studies the diffuse energy rather than only the PDE. For elliptic transmission with internal interface $\partial D$01, one diffuse energy is

$\partial D$02

For the multidimensional Neumann transmission case $\partial D$03, the energies $\partial D$04 were shown to $\partial D$05-converge in strong $\partial D$06 to the sharp-interface energy $\partial D$07, and minimizers converge strongly in $\partial D$08 and, up to subsequences, in $\partial D$09. Robin transmission with $\partial D$10 remains fully open in general dimension, although the one-dimensional case is covered (Luong et al., 23 Apr 2025).

Moving-domain and interface-evolution problems have also been recast in diffuse form. In the vertical-growth model for stacked 2D materials, a scalar phase field $\partial D$11 represents substrate, first layer, and second layer. Diffuse characteristic functions $\partial D$12, $\partial D$13 localize adatom diffusion equations to the appropriate terraces, while kinetic attachment laws are encoded by $\partial D$14-weighted source terms. The phase field itself evolves by a Cahn–Hilliard-type law so that $\partial D$15 approximates curvature. For the one-layer prototype problem, matched asymptotics showed overall second-order convergence in $\partial D$16, and numerical tests supported second-order behavior for concentration and interface velocity (Guo et al., 2019).

In multiphase and fluid problems, DDM has been used to avoid explicit enforcement of complicated wall conditions. A diffuse-domain, consistent and conservative multiphase model rewrites the reduction-consistent Cahn–Hilliard and Navier–Stokes equations on a regular box using a smooth characteristic function $\partial D$17, with wall wettability encoded by $\partial D$18-localized source terms and no-slip enforced by the penalty

$\partial D$19

The resulting DD-CC equations were coupled to a Hermite-moment lattice Boltzmann discretization (Liu et al., 2023). In a different direction, a chemotaxis-fluid diffuse-domain model for sessile drops embedded the droplet into a rectangle, used $\partial D$20-weighted conservative equations for cell density, oxygen concentration, and incompressible flow, and formally recovered mixed sharp-interface conditions, including

$\partial D$21

as $\partial D$22 (Wang et al., 2022).

A recent unifying development formulates DDMs by Onsager’s variational principle. In that framework, sharp-surface free-energy and dissipation functionals are embedded through a diffuse surface density $\partial D$23, with typical choices

$\partial D$24

and the governing equations are derived from the Rayleighian rather than from term-by-term extension of the sharp PDE. This variational organization is intended to systematize DDM construction for moving surfaces, surface transport, and fluid-structure interaction (Xu, 13 May 2026).

6. Software implementations, practical use, and open directions

The methodological maturity of DDM has led to software frameworks intended to expose multiple variants through a common interface. A recent example is ddfem, a Python package in which both the original PDE and the transformed diffuse-domain PDE are represented in UFL. The package provides geometry definition through signed distance functions and constructive solid geometry, a collection of DDM transformers such as DDM1, Mix0DDM, NDDM, and NSDDM, and a new mechanism for combining distinct boundary conditions on disjoint boundary segments through normalized weights

$\partial D$25

This construction was introduced to prevent artificial over-enforcement where multiple component phase fields overlap near corners or CSG artifacts (Benfield et al., 22 Jul 2025).

In practice, several implementation patterns recur across the literature. The interface thickness $\partial D$26 is usually chosen proportional to the mesh size, the computational box is kept sufficiently large that its artificial outer boundary does not contaminate the solution near the physical domain, and data defined on $\partial D$27 or $\partial D$28 are extended off the boundary along approximate normals. Standard rectangular finite elements, central finite differences, BDF2 or Crank–Nicolson time stepping, multigrid solvers, adaptive Cartesian refinement near $\partial D$29, and projection methods for incompressible flow all appear repeatedly (Hao et al., 3 Apr 2025, Yu et al., 2018, Wang et al., 2022).

Several limitations remain stable across problem classes. Many rigorous results assume $\partial D$30 is at least $\partial D$31, and transmission $\partial D$32-convergence results may require $\partial D$33 interfaces (Luong et al., 23 Apr 2025). Strong regularity is often imposed on the exact solution, extensions, and forcing; for parabolic problems one commonly assumes bounds such as

$\partial D$34

(Hao et al., 3 Apr 2025). Dirichlet formulations remain more sensitive than Neumann or Robin formulations, and some higher-order corrections improve asymptotics at the cost of ellipticity or conditioning (Lervåg et al., 2014, Yu et al., 2018). In several papers, the deepest analytic development still concerns diffuse-interface consistency and weighted-space stability rather than the full combined $\partial D$35-error theory (Hao et al., 3 Apr 2025, Xu, 20 Oct 2025).

Open questions are correspondingly structural rather than cosmetic. They include multidimensional Robin transmission in the variational $\partial D$36-convergence setting, sharper norm equivalences between weighted diffuse norms and sharp-domain norms, rigorous asymptotic analysis for newer mixed and Nitsche-type Dirichlet formulations, and systematic variational formulations for moving-surface models with coupled surface and bulk physics (Luong et al., 23 Apr 2025, Benfield et al., 29 Sep 2025, Xu, 13 May 2026). Taken together, these developments suggest that DDM is best understood not as a single scheme, but as a family of diffuse-interface embeddings whose analytic and numerical properties depend crucially on how geometry, boundary operators, constitutive structure, and discretization are combined.