Immersed Boundary Smooth Extension (IBSE)

• IBSE is a high-order numerical method that smoothly extends the solution and its derivatives to enforce global Ck regularity on Cartesian grids.
• It achieves fourth-order accuracy for elliptic and third-order for Neumann problems, with proven convergence in both stationary and moving boundary conditions.
• IBSE utilizes spectral differentiation, specialized spread/interpolation operators, and a Schur-complement strategy for efficient, high-accuracy simulations.

The Immersed Boundary Smooth Extension (IBSE) method is a high-order numerical technique for solving partial differential equations (PDEs) in arbitrary smooth domains using a Cartesian computational grid. The core principle is the construction of a solution that possesses global $C^k$ regularity by smoothly extending the unknown field from the physical domain to the full computational domain, a strategy that overcomes the low-order limitations of classical immersed boundary approaches, particularly near boundaries. IBSE methods have demonstrated fourth-order accuracy for elliptic equations and third-order for Neumann problems, and have been successfully applied to time-dependent, nonlinear, and fluid dynamics problems, including incompressible Navier–Stokes and Stefan-type interface evolution.

1. Mathematical Foundations and Core Formulation

Given a physical domain $\Omega$ with smooth boundary $\Gamma$ embedded in a simple computational domain $C$ (e.g., a torus $\mathbb{T}^d$), the extension region is $E = C \setminus \Omega$. The IBSE-$k$ method constructs an extension $\eta$ of the unknown solution $u$ such that $\eta$ agrees with $u$ and its first $k$ normal derivatives on $\Gamma$, enforcing $u \in C^k(C)$.

For Dirichlet Poisson problems: $$\Delta u = f \quad \text{in}~\Omega, \quad u = g \quad \text{on}~\Gamma$$ the IBSE formulation is

$$\Delta u_e - \chi_E \Delta \eta = \chi_\Omega f \quad \text{in}~C$$

with

$\mathcal{H}^k \eta + T_k F = 0 \quad \text{in}~E$

and interface matching constraints ensuring normal-derivative continuity up to order $k$: $$R_k^*(\eta - u_e) = 0 ~ \text{on}~\Gamma, \quad S_{(0)}^* u_e = g ~\text{on}~\Gamma$$ where $\mathcal{H}^k = \Delta^{k+1} + (-1)^{k+1} \Theta$ is an elliptic smoothing operator. The spread $S_{(j)}$ and interpolation $S_{(j)}^*$ operators are defined via regularized delta-functions and their normal derivatives, forming composite operators $T_k$ and $T_k^*$.

The solution $u_e |_\Omega$ recovers the original PDE with enforced boundary regularity, enabling the use of spectral methods and inducing high-order pointwise convergence.

2. Extension Operators, Discretization, and Boundary Treatments

Spread and interpolation operators are central to IBSE:

• $$S_{(j)} F(x) = (-1)^j \int_{\Gamma} F(\alpha) \partial^j_n \delta(x - X(\alpha))\, dX(\alpha)$$
• $$S_{(j)}^* \xi(\alpha) = (-1)^j \int_C \xi(x) \partial^j_n \delta(x - X(\alpha))\, dx$$

These operators provide block coupling between Eulerian grid functions and Lagrangian boundary data, supporting higher-order derivative matching.

Spatial discretization is typically performed on $N \times N$ uniform grids, with periodic Fourier spectral differentiation ensuring diagonalizable operators. The regularized delta-function is constructed for $C^3$ regularity and 16$\Delta x$ support, resulting in fourth-order accuracy for $j = 0...3$. The conditioning of the high-order operator $\mathcal{H}^k$ is controlled by the parameter $\Theta \sim (1/(N \Delta x))^{2(k+1)}$.

Temporal integration for parabolic or nonlinear problems employs implicit-explicit (IMEX) schemes or high-order BDF methods (e.g., BDF4), leveraging the elliptic solvers for efficient time stepping.

3. Application to Fluid Dynamics and Moving Boundaries

The IBSE method generalizes to incompressible Stokes and Navier–Stokes systems by coupling velocity and pressure extensions with divergence constraints: $$(\alpha I - \Delta) U + \nabla P - \chi_E \left[ (\alpha I - \Delta) \xi_u + \nabla \xi_p \right] = \chi_\Omega f$$

$$\nabla \cdot U - \chi_E (\nabla \cdot \xi_u) = \chi_\Omega f_p$$

Boundary matching is enforced for up to $k$ normal derivatives of velocity and $k-1$ for pressure, ensuring global regularity:

• $R_k^* \xi_u = R_k^* U$ on $\Gamma$
• $T_{k-1}^* \xi_p = T_{k-1}^* P$ on $\Gamma$
• $S^* U = g$ on $\Gamma$

For moving-boundary problems, such as Stefan-type phase transitions coupled to convection (Huang et al., 2020), the IBSE framework is integrated with a $\theta$–$L$ interface-evolution scheme. The curve $\Gamma$ is represented by tangent angle $\theta(\alpha,t)$ and length $L(t)$, evolving via surface dynamics and regularization: $$\theta_t = \frac{1}{L} \left[ \beta/Pe\, \partial^2 c / \partial n \partial \alpha + V_s\,\theta_\alpha \right] + \frac{\varepsilon}{L^2} \theta_{\alpha\alpha}$$ Bulk advection-diffusion and streamfunction-vorticity equations are solved in IBSE form each timestep, with explicit Schur-complement and GMRES-based solvers maintaining efficiency as boundaries evolve.

4. Algorithmic Workflow and Computational Complexity

IBSE employs a Schur-complement strategy to reduce the coupled PDE-ODE system to a small dense block for interfacial unknowns—Lagrange multipliers enforcing boundary and regularity conditions. The high-level workflow at each timestep includes:

1. Evaluation of needed boundary derivatives, interface evolution via advanced BDF/IMEX schemes.
2. Construction/refresh of spread/interpolation operators.
3. Schur-complement matrix re-factorization if boundary geometry has significantly changed.
4. Solution of the IBSE system for bulk fields using SC preconditioned GMRES.
5. Recovery of gradient and normal derivatives for subsequent time steps.
6. Assembly of advective terms for nonlinear coupling.

The complexity per solve is $O(N \log N)$ for 2D or $O(N^{4/3})$ for 3D FFT-based spectral methods; setup cost for the Schur-complement matrix is $O((n_\Gamma)^{3})$. For stationary domains, LU factorization of the SC is precomputed; for moving boundaries, iterative or hierarchical updating strategies are necessary.

Minimal geometric data—Cartesian boundary points and normals—are required, and the method adapts directly to both stationary and moving domain problems with smooth boundaries.

5. Convergence Analysis and Empirical Results

Rigorous convergence studies confirm precise high-order accuracy:

• For pure elliptic (Poisson, Dirichlet), IBSE-$k$ achieves $O(\Delta x^{k+1})$ in $L^\infty$ norm for the solution field.
• Neumann problems lose one order, yielding $O(\Delta x^{k})$.
• In 2D, IBSE-3 demonstrates fourth-order spatial convergence (Dirichlet) and third-order (Neumann) (Stein et al., 2015).
• For incompressible flows, IBSE-2 delivers third-order for velocity and second-order for pressure/stress tensor components (Stein et al., 2016).
• In Stefan-type dissolution coupled to fluid flow, third-order temporal and spatial convergence for interface variables (pure Stefan), and global second-order accuracy for the coupled Navier–Stokes–Stefan system (Huang et al., 2020).

Numerical benchmarks include analytic Frank-disk melting, Mullins–Sekerka instability, flow past obstacles, and pattern-formation in dissolution, demonstrating qualitative and quantitative matches to expected morphological phenomena.

6. Extensions, Limitations, and Practical Guidelines

IBSE methods generalize to broader classes of PDEs, including nonlinear reaction-diffusion, coupled eigenvalue problems, and time-dependent convection-dominated systems. The extension order $k$ governs achievable accuracy; for double precision, $k = 2, 3$ is practical before operator conditioning degrades. The delta-function regularization must be sufficiently smooth and wide (e.g., $C^3$ with 16$\Delta x$ support). Extension to three dimensions hinges on accurate quadrature for boundary integration (e.g., spherical harmonics or high-order panels).

A plausible implication is that for moving-boundary or highly nonlinear domains, IBSE's cost and setup can be minimized by deploying fast, updateable solvers for the Schur-complement block. The use of arbitrarily high $k$ is theoretically possible, contingent on construction of compact-support delta functions or Fourier-type kernels with requisite regularity.

7. Relation to Other High-Order Extension Methods

The Smooth Forcing Extension (SFE) method (Qadeer et al., 2020) implements similar principles via explicit Fourier continuation for forcing extension and NUFFT-based interpolation at non-uniform boundary nodes. SFE achieves arbitrary spatial orders by matching normal derivatives of the inhomogeneous term, with extension accuracy tied directly to the regularity of the basis and order of interpolation. Both IBSE and SFE avoid grid generation and boundary-fitted mesh complexities, relying instead on global extensions and small dense systems for success. Numerical tests in SFE indicate sub-geometric convergence in highly regular 1D/2D settings and robust stability over extended time integration.

The IBSE and SFE frameworks represent significant developments in high-order PDE solvers on complex domains, advancing the efficiency and accuracy of immersed boundary paradigms and enabling new classes of problems—fluid–structure interaction, interface evolution, and spectrally accurate eigenvalue computations—previously impeded by boundary regularity constraints.