Invariant Imbedding Method

• Invariant imbedding is a method that recasts global boundary-value problems as a sequence of initial-value evolution equations with known start conditions.
• It applies broadly—from stratified electromagnetics and multiple scattering to T-matrix light scattering and continuous-depth neural network formulations.
• The technique improves numerical efficiency and accuracy by incrementally building media or network depth, thus simplifying complex matching problems.

Invariant imbedding method is a reformulation in which a boundary-value or terminal-value problem is recast as evolution with respect to a growth variable—typically slab thickness $L$, optical depth $\tau$, network depth $p$, or shell radius $r$—so that reflection/transmission operators, adjoints, or scattering matrices satisfy initial-value equations as the system is enlarged. In stratified electromagnetics it advances reflection and transmission matrices as a medium is built up from zero thickness (Kim et al., 2016); in multiple scattering it evolves the Fourier reflection function in optical depth (Kawabata, 2016); in continuous-depth learning it treats depth as an explicit variable and reduces the usual forward–backward optimal-control problem to forward-facing initial value problems (Corbett et al., 2022); and in light scattering by faceted particles it grows the T-matrix through concentric spherical shells governed by a matrix Riccati equation (Wu et al., 19 Jun 2026).

1. Core formulation

The defining operation is to replace a fixed global problem by a nested family of problems. For stratified media, one considers a slab of thickness $L$ and asks how the scattering data change when an infinitesimal layer is appended. For continuous-depth networks, one varies the left endpoint $p$ of the interval $[p,q]$ while holding the input $x$ fixed and comparing the nested systems $[p+\Delta,q]\subset[p,q]$. For the invariant-imbedding T-matrix method, one grows a particle outward from an inner core through concentric spherical shells and evolves the T-matrix in the shell radius (Kim et al., 2016, Corbett et al., 2022, Wu et al., 19 Jun 2026).

In this sense, the method is not tied to a single physical observable. The evolved quantity may be a reflection matrix $r(L)$, a transmission matrix $\tau$0, a reflection function $\tau$1, an imbedded adjoint $\tau$2, or a T-matrix $\tau$3. What remains invariant depends on the application: in continuous-depth learning the input $\tau$4 is kept invariant while depth varies; in layered-wave problems the lower structure provides the initial data from which the full system is built; in IITM the inner core supplies the starting matrix $\tau$5 (Corbett et al., 2022, Kim et al., 2016, Wu et al., 19 Jun 2026).

A plausible implication is that invariant imbedding is best regarded as a structural change of independent variable rather than a single numerical scheme. The common feature is the conversion of a global matching problem into evolution equations whose initial data are known at zero thickness, zero added depth, or a simple core configuration.

2. Canonical evolution equations

In stratified bi-isotropic media, the reduced field vector is

$\tau$6

and the coupled wave equation can be written as

$\tau$7

with

$\tau$8

Invariant imbedding yields differential equations in the slab thickness $\tau$9 for the reflection and transmission matrices. In the normalized form used for arbitrarily inhomogeneous stratified bi-isotropic media,

$p$0

$p$1

with Fresnel-type initial values at $p$2 (Kim et al., 2016).

For scalar multiple scattering in a plane-parallel atmosphere, the Fourier coefficient of the reflection function satisfies

$p$3

where $p$4 contains direct single scattering and terms involving the current reflection function, so the equation is nonlinear and implicit in $p$5 (Kawabata, 2016).

For faceted-particle light scattering, IITM evolves the T-matrix in radius through

$p$6

with shell-coupling matrix

$p$7

The common pattern is that growth in $p$8, $p$9, or $r$0 induces a closed evolution law for a reduced scattering object rather than for the full field everywhere (Wu et al., 19 Jun 2026).

3. Stratified electromagnetics and bi-isotropic media

Bi-isotropic media are the most general form of linear isotropic media in the papers considered, with constitutive relations

$r$1

Chiral media correspond to $r$2, and Tellegen media to $r$3. In such media, $r$4 and $r$5 polarizations are coupled, so invariant imbedding is formulated as a two-channel matrix scattering problem rather than as two decoupled scalar problems (Kim et al., 2016).

The generalized formulation for arbitrarily inhomogeneous stratified bi-isotropic media computes not only $r$6 and $r$7 but also the internal field matrix $r$8. This permits direct evaluation of reflectance, transmittance, absorptance, and field distributions in absorbing, evanescent, and strongly inhomogeneous regimes. The method was verified on a uniform chiral slab, on surface-wave excitation in a Tellegen/metal bilayer, and on mode conversion in inhomogeneous Tellegen media. A key conclusion is that, unlike ordinary isotropic media, surface-wave excitation and mode conversion occur for both $r$9 and $L$0 waves in bi-isotropic media (Kim et al., 2016).

For the interface between a metal and a general bi-isotropic medium in the Kretschmann configuration, invariant imbedding was used together with an analytical surface-wave dispersion relation. The latter yields a complex propagation constant $L$1, from which the effective index

$L$2

and propagation length

$L$3

follow. The numerical observables computed by invariant imbedding include absorptance, cross-polarized reflectance, and internal field profiles. The agreement between the invariant-imbedding calculation and the dispersion relation is reported as perfect (Kim et al., 2016).

The resulting polarization physics differs qualitatively from the ordinary metal/dielectric case. In chiral media, the effective index is an increasing function of the chirality index, whereas in Tellegen media it is a decreasing function of the Tellegen parameter. The propagation length in both cases increases substantially as either parameter increases. In the chiral case, sufficiently large chirality produces quasi surface waves; in the Tellegen case, the propagation length diverges when the effective index goes to zero (Kim et al., 2016).

4. Multiple scattering, stiffness, and hybrid fast invariant imbedding

In plane-parallel atmospheric multiple scattering, invariant imbedding is applied to the reflection function

$L$4

for an atmosphere consisting of $L$5 homogeneous slabs above a Lambert surface. The challenge is that the imbedding equation is stiff for large

$L$6

which occurs for highly slanted directions. Standard explicit ODE solvers then become inefficient (Kawabata, 2016).

The fast invariant imbedding method of Sato et al. (1977), as described in the hybrid algorithm paper, addresses this by approximating the source term locally by low-order Lagrange polynomials, integrating the exponential integrating-factor form analytically over each step, and solving the resulting implicit nonlinear equation by successive iteration. The first step uses a first-order Lagrange polynomial through two points; subsequent steps use a quadratic Lagrange polynomial through three points. If convergence fails within the maximum number of iterations $L$7, the step is reduced by a factor $L$8; step growth is controlled by $L$9. An optional termination criterion is based on the derivative magnitude $p$0 (Kawabata, 2016).

The hybrid method combines fast invariant imbedding with doubling–adding. The reflection function of the lowermost homogeneous slab is computed by the doubling–adding method, including the Lambert ground through an adding correction, and this reflection function then supplies the initial value for fast invariant imbedding as the remaining slabs are added upward. The stated computational consequence is that the execution speed of the hybrid method is no slower than one half of that of the doubling–adding method even in the most unsuitable cases for fast invariant imbedding, while for some cases it is approximately four times faster (Kawabata, 2016).

This division of labor is methodologically specific. Doubling–adding is used where thick lower layers are most costly for pure invariant imbedding, whereas invariant imbedding is used where repeated layer addition becomes favorable. The paper further reports that the efficiency of the hybrid method increases rapidly with the number of atmospheric slabs and the optical thickness of each slab (Kawabata, 2016).

5. Invariant-imbedding T-matrix method for faceted particles

The invariant-imbedding T-matrix method is a standard tool for light scattering by large, sharply faceted, non-axisymmetric particles such as atmospheric ice crystals and mineral dust. Its attraction is that the T-matrix is obtained not from a single boundary integral, as in the extended boundary condition method, but by growing the particle outward from an inner core through concentric spherical shells. This formulation is robust at large size parameter and for sharp-edged shapes where surface-integral formulations can become ill-conditioned (Wu et al., 19 Jun 2026).

The limiting numerical issue is “staircasing.” On each shell, the dielectric contrast is piecewise constant inside and outside the particle, so the quadratures over $p$1, $p$2, and $p$3 encounter non-smooth integrands wherever the shell crosses sharp boundaries. The paper identifies a single geometric origin for this non-smoothness: tangencies between the integration sphere and particle faces or edges. These tangencies generate three classes of singular behavior: jumps, kinks, and half-integer branches. In the radial direction the local branch law is

$p$4

with

$p$5

The key claim is that the leading non-analytic exponent is fixed by local contact geometry (Wu et al., 19 Jun 2026).

The boundary-conformal correction is direction-specific. In azimuth, boxcar-like cross-sections are integrated using closed-form azimuthal coefficients, so for a hexagonal prism the azimuthal integration becomes exact. In zenith, the interval is split at analytically known breakpoints, including

$p$6

and the endpoint substitution

$p$7

is applied when a square-root branch is present. The identity

$p$8

shows that the inscribed-cylinder crossing is analytic in $p$9, not in $[p,q]$0. The paper therefore states that the zenithal crossing is a square-root branch rather than a kink, so interval-splitting alone yields only $[p,q]$1 Gauss convergence, whereas the substitution restores spectral convergence (Wu et al., 19 Jun 2026).

In the radial step, the Riccati equation is lifted to a linear system,

$[p,q]$2

and advanced by fourth-order Runge–Kutta together with the Möbius update

$[p,q]$3

Without singularity removal, a branch $[p,q]$4 limits the effective order to $[p,q]$5. Applying $[p,q]$6 absorbs the half-integer branch into an integer power of $[p,q]$7, restoring the designed fourth order of the RK4 lift (Wu et al., 19 Jun 2026).

For a hexagonal prism, the reported asymptotic behavior is exact azimuthal integration, spectral zenithal convergence after the square-root substitution, and fourth-order radial convergence with the substituted Möbius lift; the unsubstituted radial lift is limited to about $[p,q]$8. For the solid hexagonal bullet, the same mechanism persists even though tilted faces add more tangency loci: zenithal panel splitting alone again stalls at $[p,q]$9, the square-root substitution restores spectral convergence, and the radial Möbius lift with substitution again reaches fourth order. The convergence orders are verified size-independent up to $x$0; what grows with size is the resolution needed to reach the asymptotic regime, not the order (Wu et al., 19 Jun 2026).

6. Continuous-depth learning and optimal control

In continuous-depth neural networks, invariant imbedding is used to make network depth explicit. Starting from the Neural ODE

$x$1

the solution is rewritten as $x$2, so that each left endpoint $x$3 defines a network on $x$4. The input $x$5 remains invariant while depth varies. The central forward relation proved in the paper is

$x$6

This treats depth as a dynamical variable and links differentiation with respect to depth to differentiation with respect to input (Corbett et al., 2022).

For the general nonlinear vector-valued Bolza problem with running loss $x$7 and terminal loss $x$8,

$x$9

the invariant-imbedding formulation yields an imbedded adjoint satisfying

$[p+\Delta,q]\subset[p,q]$0

with terminal condition

$[p+\Delta,q]\subset[p,q]$1

The paper also derives a first-order optimality condition,

$[p+\Delta,q]\subset[p,q]$2

The stated consequence is that the usual forward ODE plus backward adjoint boundary-value problem is converted into a system of forward initial value problems in the depth variable $[p+\Delta,q]\subset[p,q]$3 (Corbett et al., 2022).

The discrete realization is the Invariant Imbedding Network (InImNet). The algorithm initializes at the trivial deepest network,

$[p+\Delta,q]\subset[p,q]$4

and then moves to shallower depths by Euler-type updates. The paper presents the resulting architectures as discrete implementations of the invariant-imbedding formulation and as resources comparable to imbedded residual neural networks. Experiments on supervised learning and time series prediction are reported as competitive; on the bouncing-balls problem the paper gives about $[p+\Delta,q]\subset[p,q]$5 s/epoch for InImNet at $[p+\Delta,q]\subset[p,q]$6 versus about $[p+\Delta,q]\subset[p,q]$7 s/epoch for the shooting-based baseline. The major practical limitation is the appearance of nested Jacobians such as $[p+\Delta,q]\subset[p,q]$8, which can be memory-intensive and numerically rough in high-dimensional or high-depth regimes (Corbett et al., 2022).

Across these applications, invariant imbedding serves as a domain-independent mechanism for exposing the dependence of a solution on added extent. In layered-wave problems that extent is thickness or optical depth; in IITM it is shell radius; in continuous-depth learning it is network depth. The technical content varies sharply by domain, but the methodological constant is the same: build the system incrementally, derive evolution equations for the reduced observable of interest, and exploit the resulting initial-value structure.