Immersed Interface Method (IIM)

• Immersed Interface Method (IIM) is a family of unfitted discretizations that sharply enforces discontinuities using correction terms derived from interface jump conditions.
• It applies to elliptic equations, incompressible flows, and fluid–structure interactions by accurately handling solution and flux jumps on simple, typically Cartesian, meshes.
• Recent advancements integrate high-order compact formulations and scalable solver technologies, achieving precise and robust performance on complex, multi-dimensional problems.

The immersed interface method (IIM) is a family of unfitted, sharp-interface discretizations for partial differential equations with embedded boundaries or material interfaces. Its defining idea is to retain a simple background mesh—typically Cartesian or otherwise interface-independent—while incorporating discontinuities through correction terms derived from interface jump conditions, rather than through body-fitted remeshing or the regularized spreading used in classical immersed boundary formulations. In the elliptic setting, the interface may carry discontinuities in the solution and flux; in incompressible flow and fluid–structure interaction, it induces jumps in pressure and velocity derivatives generated by singular interfacial forces. Across these settings, the common objective is accurate resolution of interface physics on simple meshes, with local modifications confined to stencils or approximation spaces that intersect the interface (Kolahdouz et al., 2018, Singhal et al., 2021).

1. Defining formulation and jump structure

In the elliptic interface literature represented here, a standard model problem is

$$\nabla \cdot (\beta \nabla u) = f \quad \text{in } \Omega,$$

with piecewise-constant or discontinuous coefficient $\beta$, and interface conditions

$$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$

This formulation captures both solution jumps and flux jumps across an embedded interface without requiring the computational mesh to conform to $\Gamma_M$ (Gabbard et al., 28 Mar 2025). A more general two-dimensional elliptic setting treated in higher-order compact IIM work includes singular sources and discontinuous coefficients,

$$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$

with interface conditions such as

$[u]=\hat S, \qquad [\beta u_n]=\sigma,$

where the embedded interface splits the domain into $\Omega^-$ and $\Omega^+$ (Singhal et al., 2021).

For incompressible flow, the interface carries stress discontinuities rather than a direct elliptic jump alone. In the discrete-surface and fluid–structure interaction formulations, the governing jump conditions take the form

$\llbracket \mathbf{u} \rrbracket = 0,$

$$\llbracket p \rrbracket = -\jmath^{-1}\mathbf{F}\cdot \mathbf{n},$$

$\beta$0

so that normal force components determine the pressure jump and tangential components determine jumps in velocity gradients or shear stress. This separation is central to IIM formulations for viscous incompressible flow, discrete surfaces, and sharp-interface fluid–structure interaction (Facci et al., 2024).

A recurring point in this literature is that IIM is not defined by a particular PDE class, but by a treatment principle: the mesh stays simple, the interface is embedded, and the discontinuity is imposed sharply through jump-aware corrections. This suggests that “immersed interface” denotes a discretization philosophy that can be instantiated in elliptic, Navier–Stokes, lattice Boltzmann, and higher-order multiphysics solvers, provided the relevant jump structure can be expressed in usable discrete form (Kolahdouz et al., 2018).

2. Discrete mechanism: irregular points, polynomial extension, and stencil correction

Classical IIM behavior is local. Away from the interface, one uses standard finite differences or standard finite elements; only at irregular points—grid points or elements whose stencil is cut by the interface—does the method deviate from the underlying background discretization. In finite-difference formulations, this typically means correcting derivatives or Laplacians by Taylor expansions involving jump terms, so that the interface is felt only by cut stencils. In discrete-surface formulations, the same principle appears through projected jump data and correction terms inserted into momentum or Poisson operators near cut grid lines (Kolahdouz et al., 2018).

A particularly explicit reinterpretation appears in the vorticity–velocity Navier–Stokes formulation, where the immersed interface method is recast as a polynomial extrapolation scheme. There, one introduces control points, affected points, and extrapolated values inside the obstacle. For a Dirichlet wall, the jump correction becomes equivalent to evaluating an interpolating polynomial built from nearby fluid-side values and boundary data, then applying standard conservative finite differences. This perspective makes IIM closely related to ghost-cell and sharp-interface ideas, while preserving the formal language of jump corrections (Gabbard et al., 2021).

High-order immersed Poisson solvers use the same structure in a different notation. Near a boundary control point $\beta$1, a local interpolant $\beta$2 is constructed from physical-side data, and the discrete operator is applied to the extended field

$\beta$3

For Neumann data, the boundary value is reconstructed from a discrete normal derivative formula; for material interfaces, one solves a local $\beta$4 system combining the solution jump and flux jump to obtain $\beta$5. This is a canonical IIM pattern: missing or crossed stencil values are not guessed diffusely, but reconstructed from interface conditions and local polynomial structure (Gabbard et al., 28 Mar 2025).

The same locality governs higher-order compact formulations. In the higher-order compact explicit-jump immersed interface method, regular points use a nine-point compact stencil, while irregular points are treated by expanding crossed stencil values about nearby regular points and adding explicit jump corrections at interface crossing locations. The central design requirement is not merely correctness at the crossing, but preservation of the same compact nine-point structure throughout the domain, with only the right-hand side modified by explicit corrections (Singhal et al., 2021).

3. Accuracy, compactness, and high-order developments

A long-standing issue in the IIM family is the trade-off between sharp enforcement, stencil compactness, and global order. The higher-order compact explicit-jump immersed interface method explicitly positions itself against two earlier reference points: the standard IIM of LeVeque–Li, which is second-order and uses local Taylor expansions around the interface, and the explicit-jump IIM of Wiegmann and Bube, which is explicit but does not retain the same compact nine-point high-order structure throughout the domain. In contrast, HEJIIM blends a higher-order compact nine-point stencil with explicit jump treatment and chooses $\beta$6 in the Taylor/jump expansions so that irregular points also achieve $\beta$7 accuracy. For the circular-interface Poisson problem, the reported infinity-norm error decreases from $\beta$8 at $\beta$9 to $$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$0 at $$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$1, with observed rates of convergence around $$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$2–$$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$3 (Singhal et al., 2021).

High-order sharp immersed methodology has since been extended from elliptic solves to full incompressible Navier–Stokes time integration. A fourth-order immersed solver combines a Runge–Kutta-based projection method with a new fifth-order immersed discretization of the advection term and high-order interface-corrected finite differences for the remaining operators. The reported result is fourth-order convergence of both velocity and pressure in the infinity norm, both inside the domain and on immersed boundaries, for stationary and moving-boundary manufactured problems. The same work also emphasizes that maintaining high-order pressure accuracy requires a pseudo-pressure Neumann boundary condition that includes boundary acceleration terms, because the more common projection boundary condition becomes only high-order in time for steady boundary velocities (Ji et al., 20 Aug 2025).

For elliptic problems on complex domains, high-order immersed discretization has also been coupled to scalable linear algebra. A matrix-free Krylov solver preconditioned by a low-order Shortley–Weller multigrid method supports fourth- and sixth-order immersed finite-difference discretizations of Laplace and Poisson equations with Dirichlet, Neumann, and interface jump conditions. The analysis identifies a boundary truncation error of order $$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$4, so that preserving $$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$5-th order global accuracy requires $$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$6 for the boundary/interface interpolants. This line of work shows that high-order IIM is no longer limited to small verification problems; it has been demonstrated on complex 3D domains and multiresolution adaptive grids (Gabbard et al., 28 Mar 2025).

These developments clarify a common misconception: sharp-interface treatment does not inherently restrict IIM to second-order accuracy. The literature here shows second-order, fourth-order, fifth-order advection discretization within fourth-order flow solvers, and sixth-order elliptic convergence, with the caveat that the local reconstruction, jump representation, and linear solver must all be designed coherently for the target order (Singhal et al., 2021).

4. Relation to immersed finite elements and other unfitted interface methods

IIM and immersed finite element methods share the same nonconforming geometric premise but differ in where the interface information is encoded. Classical IIM modifies discrete equations or finite-difference stencils near the interface. Immersed finite element methods modify the local approximation space on interface elements while keeping the global mesh interface-independent. The 3D trilinear immersed finite element construction is explicit on this point: on non-interface Cartesian elements one uses $$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$7, while on interface elements one uses piecewise trilinear polynomials tied together by approximate jump conditions across an approximating plane. The resulting local space still has dimension $$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$8, preserves standard nodal/Lagrange degrees of freedom, and is unisolvent regardless of interface location and coefficient contrast (Guo et al., 2019).

The partially penalized 3D trilinear immersed finite element method extends this idea to rigorous a priori analysis on cuboidal unfitted meshes. It establishes trace and inverse inequalities for trilinear IFE functions on interface elements with arbitrary cutting configurations, proves first-order energy-norm and second-order $$[u] = j_0(s), \qquad [\beta \partial_n u] = j_1(s) \quad \text{on } \Gamma_M.$$9-norm error estimates, and shows that only a small fraction of elements are interface elements, with the fraction scaling like $\Gamma_M$0. Here the interface physics enters through locally modified basis functions and penalties on interface faces, not through finite-difference stencil corrections, but the governing ambition—accurate interface resolution on unfitted Cartesian meshes—remains the same (Guo et al., 2020).

In one dimension, immersed finite element analysis goes further and shows superconvergence properties that mirror standard fitted finite elements. For elliptic interface problems with a single discontinuity point, superconvergence occurs at roots of generalized orthogonal polynomials that satisfy both orthogonality and interface jump conditions. On the interface element, the discrete solution is superconvergent at roots of generalized Lobatto and generalized Legendre polynomials; for the pure diffusion case, the IFE solution is exact at all mesh nodes. This result is significant because it demonstrates that unfitted interface handling need not destroy fine-scale structure in the error distribution (Cao et al., 2015).

The same “immersed finite element realization” of interface methodology now reaches beyond second-order elliptic equations. An immersed $\Gamma_M$1 interior penalty method for biharmonic interface problems constructs high-order local IFE spaces in a least-squares sense, proves unisolvency and partition of unity, and incorporates the spaces into a modified $\Gamma_M$2 interior penalty scheme with additional penalties on interface segments. Numerical experiments confirm optimal convergence in $\Gamma_M$3, $\Gamma_M$4, and $\Gamma_M$5 norms for several interface geometries (Chen et al., 16 Sep 2025).

This comparison matters because IIM is sometimes used informally to describe any unfitted interface method. The more precise distinction in the present literature is that IIM usually denotes stencil correction in finite differences or related sharp-interface Cartesian schemes, whereas IFE modifies the finite element basis itself. The two approaches are conceptually adjacent but not identical (Guo et al., 2019).

5. Incompressible flow, discrete surfaces, and fluid–structure interaction

In incompressible flow, IIM has developed into a sharp-interface alternative to classical immersed boundary coupling. A major step was the introduction of an immersed interface method for discrete $\Gamma_M$6 surfaces, designed to work with standard nodal Lagrangian finite element surface meshes rather than smooth analytic parameterizations. In that framework, only the lowest-order jump conditions for pressure and velocity gradient are imposed, yet the method reports global second-order convergence together with nearly second-order local convergence in Eulerian velocity, between first- and second-order global convergence in Eulerian pressure together with first-order local convergence, second-order local convergence in interfacial displacement and velocity, and first-order local convergence in fluid traction. The same method was demonstrated on an anatomical model of the inferior vena cava, showing that sharp-interface treatment can be combined with complex, imaging-derived or FE-generated geometries without analytic geometry reconstruction (Kolahdouz et al., 2018).

A complementary direction formulates IIM directly in vorticity–velocity variables for two-dimensional Navier–Stokes flow with multiple bodies, nonconvex obstacles, and outflow boundaries. There the vorticity transport equation is discretized conservatively and the velocity reconstruction is reduced to a scalar Poisson problem with immersed corrections. Conservative differencing throughout leads to exact enforcement of a discrete Kelvin’s theorem, which is the key to multiply connected domains and circulation evolution around multiple bodies. The full solver achieves second-order spatial accuracy and third-order temporal accuracy, and its conservative structure supports force, pressure, and shear recovery on immersed surfaces without body-fitted meshes (Gabbard et al., 2021).

The IIM philosophy has even been transferred to lattice Boltzmann methods. In the immersed interface–lattice Boltzmann method, the central derivation is a jump condition for the particle distribution functions,

$\Gamma_M$7

which sharply resolves the normal-force-induced discontinuity, while tangential force is still treated by immersed boundary spreading. This hybrid approach produces higher-order accuracy than IB-LBM in selected benchmarks, much smaller velocity error, and markedly improved volume conservation. In a dynamic thin elastic interface test, the reported total area loss at the finest grid is around $\Gamma_M$8 for II-LBM, compared with $\Gamma_M$9 for standard IB-LBM and $$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$0 for IB-LBM with Guo forcing (Qin et al., 2020).

Taken together, these developments show that the role of IIM in flow computation is broader than boundary imposition alone. It serves as a coupling mechanism between Eulerian solvers and Lagrangian geometries, a means of obtaining sharp pressure and traction fields, and a way to preserve physically important circulation and stress-discontinuity structure that diffuse-interface methods tend to smear (Kolahdouz et al., 2018).

6. Sharp features, pressure robustness, mesh-ratio instability, and near contact

Recent work has focused less on whether IIM can resolve a smooth interface and more on where it fails when the geometry or coupling becomes difficult. One such difficulty is the treatment of corners and edges on immersed surfaces. In early discrete-surface implementations, jump conditions were projected into a continuous Galerkin finite element space on the interface. That works well for smooth geometries, but when the surface normal is genuinely discontinuous, the pressure jump and shear-jump data inherit those discontinuities, and projection into a continuous basis causes $$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$1 errors near sharp features. A discontinuous Galerkin representation of the jump fields avoids this forced smoothing. For smooth cylinders and spheres, CG and DG give comparable accuracy; for sharp geometries such as square cylinders, triangular appendages, and cubes, DG gives substantially better accuracy and a maximum stable time step that is essentially insensitive to geometric sharpness, whereas CG becomes increasingly restrictive as the angle becomes more acute (Facci et al., 2024).

A related issue is pressure robustness on discrete $$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$2 triangulations. Earlier discrete-surface IIM formulations used piecewise-constant element normals to decompose interfacial force into normal and tangential parts. This was already substantially more accurate than regularized immersed boundary methods in shear-dominated settings, but it struggled under pressure loading because the pressure jump depends directly on the surface normal. A recent remedy reconstructs a continuous normal field before projecting the jump conditions, using either an $$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$3 projection of the discontinuous normal field into a continuous finite element space or inverse centroid-distance weighted vertex normals with linear interpolation. Numerical experiments show that these reconstructed normal fields reduce leakage by up to six orders of magnitude across a range of pressures (Facci et al., 7 Mar 2026).

Robustness issues also arise in the Eulerian–Lagrangian coupling itself. In a discrete-surface IIM/ILE framework for fluid–structure interaction, stability previously required the local Lagrangian element size $$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$4 to exceed the local Eulerian grid spacing $$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$5, expressed as $$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$6. When this mesh factor ratio is violated, the force spreading operator develops a null space and parts of the interface motion become uncontrolled. A Tikhonov-inspired stabilization of the velocity restriction operator,

$$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$7

adds surface-gradient penalization and broadens the practical regime to $$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$8 in the reported tests, without meaningful change in the flow dynamics when $$(\beta u_x)_x+(\beta u_y)_y + \kappa(x,y) u=f(x,y)+\sigma \delta\{(x-x^{*})(y-y^{*})\},$$9 is chosen appropriately (Sun et al., 4 Apr 2025).

The near-contact regime exposes another limitation of single-interface assumptions. In thin gaps, a single interpolation stencil may intersect two immersed boundaries. An enhanced IIM for near-contact incompressible flows introduces a bilinear velocity interpolation operator that superposes jump corrections from multiple nearby interfaces in the same stencil. The method uses local ray casting to detect intersections with the second interface and evaluates additional correction terms there. For shearing plates and concentric or eccentric cylinders, the two-correction method remains accurate when the gap is as small as $[u]=\hat S, \qquad [\beta u_n]=\sigma,$0, whereas one-correction methods lose accuracy as $[u]=\hat S, \qquad [\beta u_n]=\sigma,$1 (Facci et al., 29 Oct 2025).

These examples reveal that the practical frontier of IIM is no longer basic interface capture but robust sharp-interface enforcement under discontinuous normals, acute corners, extreme mesh anisotropy, and multiple-interface occupancy of a single stencil.

7. Solver technology, scalability, and current scope

Because sharp-interface corrections usually destroy the symmetry and translation invariance of the background operator, scalable linear algebra has become a defining part of modern IIM. High-order immersed Poisson solvers now use matrix-free GMRES or FGMRES together with a low-order geometric Shortley–Weller multigrid preconditioner, deliberately separating high-order accuracy in the outer operator from geometric robustness in the preconditioner. The resulting systems are generally nonsymmetric, especially in the presence of immersed boundaries, so standard conjugate gradients are not appropriate. This architecture has enabled high-order immersed methods on large 3D domains and multiresolution adaptive grids (Gabbard et al., 28 Mar 2025).

The reported computational scale is notable. On adaptive 3D problems, matrix-free high-order immersed solvers have been demonstrated with about $[u]=\hat S, \qquad [\beta u_n]=\sigma,$2 billion unknowns on $[u]=\hat S, \qquad [\beta u_n]=\sigma,$3 cores, with weak scaling of about $[u]=\hat S, \qquad [\beta u_n]=\sigma,$4. The same framework supports complex geometries such as perturbed spheres and rings of linked tori, and is explicitly presented as a route toward high-order immersed treatments of linear and non-linear elasticity, incompressible Navier–Stokes equations, and fluid–structure interactions (Gabbard et al., 28 Mar 2025).

At the application level, the scope of IIM now includes elliptic interface problems with singular sources and discontinuous coefficients, bluff-body flows, multiple immersed bodies, prescribed moving boundaries, conjugate heat transfer with multiple immersed solids, discrete-surface fluid–structure interaction, lattice Boltzmann coupling, and anatomical internal flows. The literature summarized here consistently emphasizes two advantages: avoidance of body-fitted meshing and sharp resolution of discontinuities that are physically present in pressure, flux, shear, or distribution functions (Singhal et al., 2021, Ji et al., 20 Aug 2025).

A plausible implication is that IIM should now be viewed less as a single method than as a toolbox of jump-aware unfitted discretizations. Within that toolbox, the major design choices are where the interface information enters—stencil correction, extrapolation, basis modification, projection of jumps, or reconstruction of normals—and how that local correction is stabilized, accelerated, and coupled to the background solver.