Domain-of-Dependence Stabilization in DG Methods

• Domain-of-Dependence Stabilization is a cut-cell technique for discontinuous Galerkin discretizations that redistributes mass and flux to maintain the correct numerical domain of dependence.
• It overcomes small-cell CFL restrictions by extending inflow data and bypassing tiny cut-cell interiors, thus preserving stability and consistency.
• The method has been rigorously validated for scalar advection, acoustic and wave models, and SBP formulations, demonstrating accuracy and energy stability.

Domain-of-Dependence Stabilization is a cut-cell stabilization technique for discontinuous Galerkin discretizations of hyperbolic problems. Its defining purpose is to remove the small-cell time-step restriction that appears when explicit schemes are applied on cut-cell meshes, while retaining a numerical domain of dependence that matches the physical domain of dependence associated with a time step chosen from the background mesh rather than from arbitrarily small cut cells (Birke et al., 2023). In the cut-cell literature, the method is formulated by adding localized stabilization terms to the semi-discrete DG operator, initially for scalar advection, then for higher-order advection, linear symmetric hyperbolic systems, the two-dimensional acoustic wave equation, and, more recently, energy-preserving wave and summation-by-parts formulations for linear kinetic models (Streitbürger et al., 2023).

1. Small cut cells and the domain-of-dependence mismatch

The method arises from the small-cell CFL problem on cut-cell meshes. A domain $\Omega$ is embedded in a structured background mesh $\widehat{\mathcal M}_h$; intersecting the physical boundary or an internal cut with background cells produces cut cells whose area or volume fraction can be arbitrarily small. For explicit hyperbolic discretizations, a standard CFL condition depends on the smallest element size, so a single tiny cut cell can force an impractically small time step even when the rest of the mesh is of size $h$ (Birke et al., 2023).

The central observation is that the obstruction is not only geometric but also causal. For a hyperbolic equation, the physical domain of dependence of a point over one time step is determined by characteristic propagation. On a uniform mesh, an explicit DG stencil and a CFL based on $h$ produce a numerical domain of dependence that is compatible with that physical picture. On a cut-cell mesh, if $\Delta t$ is still chosen from the background mesh, a wave can physically traverse a tiny cut cell and influence a downwind neighbor in one step, but an unstabilized scheme only couples nearest neighbors through the tiny cell itself. The cut cell then becomes a numerical bottleneck: too much mass or wave content is forced through a degree of freedom with volume $|E_{\text{cut}}|\ll h^2$, and instability follows (Birke et al., 2023).

Domain-of-Dependence stabilization addresses this mismatch by enlarging the discrete domain of dependence around small cut cells. In the scalar advection formulation, the method adds extra flux terms that directly transfer information between inflow neighbors and outflow neighbors of a small cut cell, thereby bypassing the cut-cell interior as the sole conduit. In the one-dimensional fully discrete analysis, the same idea is described as redistributing mass and hence mass matrix contributions over a local neighborhood of the small cut cell, so that the norm of the discrete operator scales with the background cell size $\Delta x$, not with the cut-cell factor $\alpha$ (Petri et al., 7 Aug 2025).

2. DG formulation and the local DoD mechanism

The foundational setting is a semi-discrete DG method. For linear symmetric hyperbolic systems in two space dimensions,

$u_t + A u_x + B u_y = 0,$

with piecewise constant DG unknowns, the unstabilized semi-discretization has the form

$$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$

DoD stabilization augments this by a bilinear form $\widehat{\mathcal M}_h$0 supported only on a set $\widehat{\mathcal M}_h$1 of small cut cells: $\widehat{\mathcal M}_h$2 In the original two-dimensional scalar setting, $\widehat{\mathcal M}_h$3 coupled a triangular cut cell with exactly one inflow face and one outflow face by using an extension of the inflow-neighbor state beyond its own cell (Birke et al., 2023).

For higher-order advection, the method acquires two distinct components,

$\widehat{\mathcal M}_h$4

where $\widehat{\mathcal M}_h$5 is an interface term and $\widehat{\mathcal M}_h$6 is a volume term. The interface term transfers information from the inflow neighbor to the outflow edge of the cut cell through an extension operator $\widehat{\mathcal M}_h$7, while the volume term penalizes the discrepancy between the cut-cell polynomial and the inflow-neighbor extension inside the cut cell itself (Streitbürger et al., 2023). In one dimension, the same decomposition appears as a flux stabilization term

$\widehat{\mathcal M}_h$8

and a volume stabilization term

$\widehat{\mathcal M}_h$9

with the explicit interpretation that the first routes a fraction of the inflow directly from the upstream cell to the downwind cell and the second regularizes gradients in the small-cell region (Petri et al., 7 Aug 2025).

A central algebraic device is the extension operator. For a cell $h$0, $h$1 or $h$2 extends the polynomial defined on $h$3 to a polynomial on a larger region, often the whole domain. In the scalar advection constructions, this lets the scheme evaluate the inflow-neighbor polynomial on the outflow face of the cut cell. In later wave formulations, reflected extensions $h$4 are introduced to incorporate reflecting boundary conditions by mirroring the velocity component normal to a wall (Birke et al., 27 Jan 2026).

The amount of stabilization is controlled by a cell-dependent parameter. In higher-order advection this is written through a capacity

$h$5

so $h$6 is near $h$7 for a very small cut cell and near $h$8 for a cell that can sustain the chosen explicit time step without special treatment (Streitbürger et al., 2023). In the energy-preserving wave formulation, the analogous quantity is the capacity

$h$9

which encodes the same background-mesh-versus-cut-cell balance (Birke et al., 27 Jan 2026).

3. From scalar advection to systems, acoustics, and energy-preserving wave schemes

A major generalization concerns cut cells with multiple inflow and outflow faces. Earlier DoD work in two dimensions effectively required a triangular cut cell with exactly one inflow face, one outflow face, and one no-flow face. The extension in (Birke et al., 2023) allows triangular cut cells with arbitrary flow direction, so that a small cut cell may have two inflow and one outflow face, or two outflow and one inflow face. For a face $h$0 with outward normal $h$1, the flux matrix is

$h$2

with positive and negative parts

$h$3

The generalized stabilization on a cut cell $h$4 is then written as

$h$5

where the weights $h$6 distribute inflow over all outflow faces (Birke et al., 2023).

These weights are constrained by a partition-of-unity condition and a flux-balance condition,

$h$7

and

$h$8

together with symmetry and positive-semidefiniteness of $h$9. For simultaneously diagonalizable systems on triangular cut cells, a specific choice is

$\Delta t$0

interpreted as a normalized inflow fraction (Birke et al., 2023).

The acoustic wave equation introduced a different difficulty. In the two-dimensional first-order acoustic system, the matrices

$\Delta t$1

do not commute, so the matrix products that were symmetric in the simultaneously diagonalizable setting are no longer symmetric. The acoustic DoD formulation therefore replaces $\Delta t$2 by its symmetrized part and adds a corrective dissipative term based on the negative part of that symmetrized matrix, scaled by a parameter $\Delta t$3. This modification is used to recover semi-discrete $\Delta t$4-stability for the two-dimensional acoustic wave equation on cut-cell meshes (Birke et al., 2023).

A later wave formulation recasts the construction in a more abstract and explicitly energy-preserving form. For the linear wave equation on cut-cell meshes, the stabilization is built from propagation forms $\Delta t$5, $\Delta t$6, and $\Delta t$7, together with extension and reflected-extension operators. The forms satisfy balance identities that redistribute central numerical fluxes through a small cell while preserving either energy conservation or energy dissipation, depending on whether the underlying DG flux uses a central or dissipative split (Birke et al., 27 Jan 2026). This suggests a unifying viewpoint: DoD stabilization can be interpreted as a local flux-redistribution calculus designed to preserve the structural invariant of the base discretization.

4. Stability theory

The earliest rigorous results established semi-discrete $\Delta t$8-stability for piecewise constants. For the generalized multi-face formulation on simultaneously diagonalizable systems, the semi-discrete scheme

$\Delta t$9

satisfies

$|E_{\text{cut}}|\ll h^2$0

under homogeneous boundary conditions and the stated assumptions on $|E_{\text{cut}}|\ll h^2$1 and $|E_{\text{cut}}|\ll h^2$2. The proof uses positivity of the upwind DG form and a cut-cellwise decomposition of $|E_{\text{cut}}|\ll h^2$3 into negative semidefinite terms that are compensated by the positive face contributions of the DG operator (Birke et al., 2023).

For higher-order advection in two dimensions, the semi-discrete $|E_{\text{cut}}|\ll h^2$4-stability theorem shows that, for a ramp geometry with constant velocity $|E_{\text{cut}}|\ll h^2$5 parallel to the ramp and compactly supported solution away from inflow and outflow boundaries,

$|E_{\text{cut}}|\ll h^2$6

for all $|E_{\text{cut}}|\ll h^2$7. The proof exhibits an energy identity in which the standard upwind DG form produces positive jump terms, while the DoD terms generate an additional positive “extended jump” between the inflow and outflow neighbors across the cut cell (Streitbürger et al., 2023).

A crucial later development is fully discrete stability. For one-dimensional linear advection with DoD stabilization, the semi-discrete system $|E_{\text{cut}}|\ll h^2$8 satisfies an operator norm estimate

$|E_{\text{cut}}|\ll h^2$9

where $\Delta x$0 depends only on polynomial degree $\Delta x$1 and node type, and there is no dependence on the small cut-cell factor $\Delta x$2. Combined with semiboundedness of $\Delta x$3 and strong-stability-preserving Runge–Kutta theory, this yields a CFL-like time-step restriction that does not depend on $\Delta x$4, so the fully discrete scheme satisfies

$\Delta x$5

under a background-mesh CFL (Petri et al., 7 Aug 2025).

The summation-by-parts extension to the telegraph equation places DoD stabilization inside a broader energy framework. There, the central DoD operator is shown to be a periodic SBP operator,

$\Delta x$6

and a symmetrized upwind DoD pair satisfies

$\Delta x$7

These properties imply semidiscrete energy stability for the telegraph system, and, with appropriate IMEX Runge–Kutta splitting, the fully discrete scheme is asymptotic preserving with respect to the heat-equation limit as $\Delta x$8 (Petri et al., 9 Jan 2026).

For the energy-preserving wave formulation, the stabilized semi-discrete scheme

$\Delta x$9

satisfies

$\alpha$0

Hence the method is energy conservative when the base discretization is central and energy dissipative when it uses dissipative numerical fluxes (Birke et al., 27 Jan 2026).

5. Accuracy, consistency, and numerical evidence

For piecewise constants on general cut-cell meshes, the underlying DG scheme is formally first order, and the generalized DoD stabilization preserves consistency because the added terms vanish for constant solutions and behave as higher-order residuals in smooth regions. Numerical experiments for linear simultaneously diagonalizable systems confirm approximately first-order convergence in both $\alpha$1 and $\alpha$2, and all solution values on small cut cells remain within the range of the initial data, with no over- or undershoots reported near cut cells (Birke et al., 2023).

For higher-order advection in two dimensions, numerical convergence tests indicate orders of $\alpha$3 in the $\alpha$4 norm and between $\alpha$5 and $\alpha$6 in the $\alpha$7 norm. The reported experiments use smooth advection parallel to a ramp, explicit SSP Runge–Kutta time stepping, and a time step

$\alpha$8

chosen from the background mesh size $\alpha$9, not from the smallest cut cell (Streitbürger et al., 2023).

For the two-dimensional acoustic wave equation, the numerical study uses $u_t + A u_x + B u_y = 0,$0, final time $u_t + A u_x + B u_y = 0,$1, background grids ranging from $u_t + A u_x + B u_y = 0,$2 to $u_t + A u_x + B u_y = 0,$3, and cut-cell volume fractions in $u_t + A u_x + B u_y = 0,$4. The results show first-order convergence in $u_t + A u_x + B u_y = 0,$5 for all components. In $u_t + A u_x + B u_y = 0,$6, the pressure exhibits nearly first-order convergence, while the velocity components are more sensitive to the corrective parameter $u_t + A u_x + B u_y = 0,$7; larger $u_t + A u_x + B u_y = 0,$8 reduces absolute errors and improves convergence rates but can introduce mild wiggles (Birke et al., 2023).

The fully discrete advection analysis also includes numerical CFL studies. In one dimension, the sharp CFL number is essentially independent of $u_t + A u_x + B u_y = 0,$9 for first order, and for higher orders the small-$$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$0 regime remains close to the background method. In two-dimensional channel advection, optimized choices of the tuning parameter $$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$1 produce stable CFL numbers around $$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$2–$$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$3 for $$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$4, $$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$5–$$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$6 for $$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$7, and $$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$8–$$(\partial_t u_h(t), v_h)_{L^2(\Omega)} + a_h^{\mathrm{upw}}(u_h(t),v_h) + l_h(v_h) = 0.$$9 for $\widehat{\mathcal M}_h$00, over wall angles $\widehat{\mathcal M}_h$01, again without further restriction from the small cut cells (Petri et al., 7 Aug 2025).

The telegraph-equation SBP formulation supplies additional numerical evidence of structural robustness. Convergence tests with five small cut cells and $\widehat{\mathcal M}_h$02 show observed order approximately $\widehat{\mathcal M}_h$03 for alternating upwind fluxes. In the heat-equation limit, the condition number

$\widehat{\mathcal M}_h$04

is reported as $\widehat{\mathcal M}_h$05 without DoD on cut-cell meshes and $\widehat{\mathcal M}_h$06 with DoD-SBP, while background-mesh values are $\widehat{\mathcal M}_h$07. Implicit midpoint simulations are stable with DoD and unstable without it (Petri et al., 9 Jan 2026).

Consistency for arbitrary polynomial degree was established later. The key result is that, for sufficiently smooth exact solutions $\widehat{\mathcal M}_h$08, the DoD stabilization satisfies

$\widehat{\mathcal M}_h$09

for every stabilized small cell $\widehat{\mathcal M}_h$10, both for linear advection and for the linear wave system. The argument relies on extending the cellwise extension operators from the discrete DG space to the sum of discrete and continuous spaces and proving that, on the intersection of these spaces, the extension operators reduce to the identity. This is the analytical step needed for a refined high-order error analysis beyond the previously treated $\widehat{\mathcal M}_h$11 case (Birke et al., 11 Mar 2026).

6. Scope, limitations, and current directions

The present theory is broad but not yet complete. The generalized multi-face $\widehat{\mathcal M}_h$12-stability result for linear systems assumes constant-coefficient symmetric hyperbolic systems with simultaneously diagonalizable flux matrices and is formulated for piecewise constants on triangular cut cells in two dimensions (Birke et al., 2023). The acoustic-wave extension overcomes simultaneous diagonalizability, but its 2023 presentation remains limited to lowest-order DG and straight-line cuts producing triangular cut cells; small cut cells are also assumed not to be neighbors of one another (Birke et al., 2023).

Higher-order theory remains the main analytical frontier. Semi-discrete $\widehat{\mathcal M}_h$13-stability has been obtained for arbitrary polynomial degree in two-dimensional advection (Streitbürger et al., 2023), fully discrete stability with a background-mesh CFL has been obtained in one-dimensional linear advection (Petri et al., 7 Aug 2025), and consistency for arbitrary polynomial degree has now been proved for both advection and the linear wave equation (Birke et al., 11 Mar 2026). A plausible implication is that a full high-order a priori error theory is now structurally within reach, but that step itself is not yet part of the cited results.

There are also known high-order difficulties. In the fully discrete analysis, the extension operator becomes increasingly ill-conditioned as $\widehat{\mathcal M}_h$14 grows, and the constants in the operator norm estimate can grow rapidly with $\widehat{\mathcal M}_h$15 and with the size of the extension region. The proposed mitigation is a cut-cell-dependent tuning parameter,

$\widehat{\mathcal M}_h$16

with $\widehat{\mathcal M}_h$17 chosen by a min–max operator-norm optimization; this is verified numerically in one and two dimensions (Petri et al., 7 Aug 2025).

More recent work broadens the conceptual reach of the method. The SBP-based telegraph construction shows that DoD stabilization can be combined with asymptotic-preserving IMEX time integration and with diffusion-limit structure (Petri et al., 9 Jan 2026). The energy-preserving wave construction shows that central and dissipative DoD variants can be designed to preserve the corresponding invariant of the underlying DG spatial discretization exactly, even on cut-cell meshes with reflecting walls (Birke et al., 27 Jan 2026). This suggests that DoD stabilization is no longer only a remedy for scalar advection but a general local operator design principle for hyperbolic cut-cell discretizations.

In contemporary usage, Domain-of-Dependence Stabilization therefore denotes a family of DG-compatible cut-cell techniques whose common mechanism is to alter local flux and volume couplings around small cut cells so that the discrete domain of dependence reflects the physical one at a background-mesh time step. The defining outputs of that mechanism are conservation or controlled dissipation, semi-discrete or fully discrete stability independent of arbitrarily small cut cells, and consistency with the base DG method across increasing levels of algebraic and geometric complexity (Birke et al., 2023).