Layered Dyadic Green's Function Analysis

• Layered dyadic Green’s functions are the dyadic fundamental solutions to Maxwell’s equations in stratified media characterized by piecewise constant material properties.
• They employ spectral, Sommerfeld-integral representations to separate free-space contributions from layered corrections and support both TE/TM and matrix-basis formulations.
• Interface continuity conditions, Fresnel coefficients, and specialized computational strategies enable accurate field predictions and efficient fast multipole translations in complex media.

Searching arXiv for the cited layered-media dyadic Green's function papers and closely related work. Layered dyadic Green’s function denotes the dyadic fundamental solution of Maxwell’s equations in a stratified medium whose constitutive parameters are piecewise constant by layer and whose interfaces are planar. In such media, the dyadic Green’s tensor encodes the electric- and magnetic-field response to a point source while enforcing continuity of tangential fields and the appropriate normal flux conditions at every interface. In homogeneous space the dyadic Green’s tensor reduces to differential operators acting on the scalar Helmholtz kernel, whereas in layered media it acquires a spectral, Sommerfeld-integral structure in which TE/TM polarization channels, reflection/transmission coefficients, and vertical propagation factors replace the single free-space kernel (Yuan et al., 2 Jan 2026). The same layered structure underlies both classical TE/TM formulations and matrix-basis or operator formulations, and it also supports fast algorithms, asymptotic analysis, and extensions to anisotropic, graphene-coated, and topological layered systems (Cho et al., 2019).

1. Definition and governing equations

In a planar, horizontally invariant layered medium, layers are indexed by $\ell$ and separated by interfaces at fixed depths such as $z=d_\ell$ or $y=y_\ell$, depending on the coordinate system adopted. Each layer has piecewise constant material parameters $(\epsilon_\ell,\mu_\ell)$ and wavenumber $k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$ (Cho et al., 2019). For Maxwell’s equations in three dimensions, the dyadic Green’s tensor $\bar{\bar G}$ satisfies

$$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$

with tangential field continuity across interfaces (Cho et al., 2019). A closely related formulation introduces the Maxwell DGF pair $(G_E,G_H)$ through

$$\nabla\times G_E(\mathbf r,\mathbf r')=-\omega\mu_\ell G_H(\mathbf r,\mathbf r'),\qquad \nabla\times G_H(\mathbf r,\mathbf r')=\omega\epsilon_\ell G_E(\mathbf r,\mathbf r')-\delta(\mathbf r-\mathbf r')\frac{1}{\omega\mu_{\ell'}}\hat t,$$

for a Hertzian current dipole source located in layer $\ell'$ (Yuan et al., 2 Jan 2026).

In homogeneous space, the electric-field dyadic Green’s function is

$z=d_\ell$0

or equivalently

$z=d_\ell$1

with $z=d_\ell$2 (O'Neil, 2013). In layered media, the total field is typically decomposed into a free-space part and a reaction, scattered, reflected, or transmitted part. This separation is operationally important because the singular free-space term can be treated analytically while the layered correction is represented by Sommerfeld-type spectral integrals (Cho et al., 2016).

At each interface one imposes continuity of tangential electric and magnetic fields and continuity of the normal components of $z=d_\ell$3 and $z=d_\ell$4. In dyadic form one convenient statement is

$z=d_\ell$5

with $z=d_\ell$6 for planar layering (Yuan et al., 2 Jan 2026). In isotropic layered media these conditions decouple under TE/TM decomposition; in other formulations they are recast as scalar jump conditions for auxiliary layered Helmholtz Green’s functions (Zhang et al., 2020).

2. Spectral and Sommerfeld representations

The defining analytic structure of layered dyadic Green’s functions is their representation in the transverse Fourier domain. In two dimensions, the scalar layered kernel is written with spectral variable $z=d_\ell$7 and vertical propagation constants

$z=d_\ell$8

with branch choice such that $z=d_\ell$9 for $y=y_\ell$0 and $y=y_\ell$1 with nonnegative imaginary part for $y=y_\ell$2 (Cho et al., 2019). The general scalar layered kernel treated in the translated fast multipole framework has the form

$y=y_\ell$3

where $y=y_\ell$4 is an interface density determined by matching conditions and converges to a constant as $y=y_\ell$5 (Cho et al., 2019).

In three-dimensional electromagnetics, the dyadic Green’s function admits the spectral representation

$y=y_\ell$6

where $y=y_\ell$7 and the TE/TM pieces are built from polarization projectors and multilayer propagators (Yuan et al., 2 Jan 2026). A related formulation writes

$y=y_\ell$8

with

$y=y_\ell$9

where $(\epsilon_\ell,\mu_\ell)$0, $(\epsilon_\ell,\mu_\ell)$1, $(\epsilon_\ell,\mu_\ell)$2, and $(\epsilon_\ell,\mu_\ell)$3 (Cho et al., 2019).

This spectral decomposition separates propagating and evanescent sectors. In the scalar formulation, $(\epsilon_\ell,\mu_\ell)$4 yields oscillatory propagating plane waves, while $(\epsilon_\ell,\mu_\ell)$5 yields exponentially decaying evanescent contributions associated with surface and lateral waves (Cho et al., 2019). A plausible implication is that the same separation governs numerical stiffness in dyadic calculations, since the dyadic formulations inherit the scalar Sommerfeld structure channel by channel.

3. TE/TM decomposition and equivalent matrix-basis formulations

The classical formulation is based on TE/TM decomposition. Introducing the horizontal unit vectors

$(\epsilon_\ell,\mu_\ell)$6

the transformed fields are expressed through scalar spectral amplitudes $(\epsilon_\ell,\mu_\ell)$7 (Yuan et al., 2 Jan 2026). In each homogeneous layer the Fourier-transformed Maxwell system decouples into scalar Helmholtz equations for the TM quantity $(\epsilon_\ell,\mu_\ell)$8 and the TE quantity $(\epsilon_\ell,\mu_\ell)$9: $k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$0 with $k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$1 (Yuan et al., 2 Jan 2026). The associated scalar Green’s functions $k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$2 satisfy layered Helmholtz problems with interface conditions

$k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$3

$k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$4

and $k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$5 (Yuan et al., 2 Jan 2026).

The vector-potential matrix-basis formulation arrives at the same layered dyadic Green’s function through a different algebraic route. There one solves for a dyadic vector potential $k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$6 satisfying

$k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$7

and reconstructs

$k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$8

under the Lorenz gauge (Yuan et al., 2 Jan 2026). The matrix basis $k_\ell=\omega\sqrt{\epsilon_\ell\mu_\ell}$9 is designed so that differential operators act through a simple multiplication table and the spectral dyadics can be expressed in terms of three scalar coefficient functions $\bar{\bar G}$0, with

$\bar{\bar G}$1

(Yuan et al., 2 Jan 2026). The paper shows term-by-term equivalence between the TE/TM and matrix-basis formulations.

A related 3×3 matrix-basis framework for Maxwell and elastic layered Green’s functions states that all spectral dyadic blocks live in the “real-like” subspace

$\bar{\bar G}$2

with coefficients independent of the azimuthal angle $\bar{\bar G}$3 (Zhang et al., 2020). This suggests that the matrix basis is not merely algebraic bookkeeping; it isolates rotational symmetry and separates non-symmetric angular factors from radially symmetric spectral densities, which is central both for asymptotics and for fast multipole translation operators (Yuan et al., 2 Jan 2026).

4. Interface physics, Fresnel coefficients, and layered variants

For isotropic interfaces, spectral reflection and transmission coefficients retain the Fresnel form. At a single interface between layers $\bar{\bar G}$4 and $\bar{\bar G}$5,

$\bar{\bar G}$6

$\bar{\bar G}$7

(Yuan et al., 2 Jan 2026). In multilayers, generalized reflection and transmission coefficients are obtained recursively. For a three-layer slab, one representative generalized reflection coefficient is

$\bar{\bar G}$8

with analogous generalized transmission coefficients (Cho et al., 2016).

The layered dyadic Green’s function has been specialized to a number of physically important interface types. For graphene–dielectric stacks, graphene is modeled by an isotropic scalar surface conductivity $\bar{\bar G}$9 and enters only through impedance boundary conditions at interfaces, which modify the TE/TM reflection and transmission coefficients and lead to surface-plasmon-polariton and waveguide poles (Raad et al., 2021). In that setting, the scattered Green’s function in each layer is expanded in cylindrical vector wave functions with unknown coefficients, and recurrence relations are derived for arbitrary source and field layers via Kronecker-delta bookkeeping (Raad et al., 2021).

For layered topological insulators, the constitutive relations include the axion term

$$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$0

with

$$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$1

so the bulk propagation operator is unchanged for constant $$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$2, but interface jumps in $$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$3 modify the boundary conditions and cause TE/TM mixing on reflection and transmission (Crosse et al., 2015). In that case, the single-interface reflection operator becomes a $$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$4 matrix in the TE/TM basis rather than a diagonal pair of Fresnel coefficients. This directly contradicts the common simplification that TE and TM channels always decouple in planar media; they decouple in isotropic layered dielectrics, but not when the interface law itself mixes the channels (Crosse et al., 2015).

5. Computational formulations and fast algorithms

A major computational challenge is the efficient evaluation of Sommerfeld integrals and the construction of scalable solvers. One strategy is to derive concise spatial-domain dyadic formulas from spectral TE/TM expressions. In the half-space case, the electric-field dyadic Green’s function can be written so that all reflected components in the source layer are expressed using only four Sommerfeld integrals $$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$5, and all transmitted components in the adjacent layer use five integrals $$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$6 (Cho et al., 2016). This reduction is achieved by Bessel identities and by regrouping spectral terms so that the singular free-space dyadic is separated from the regular layered correction. Interface continuity errors reported in two- and three-layer validations are as small as $$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$7 in representative cases (Cho et al., 2016).

A different computational direction is the translated fast multipole method. For scalar layered Green’s functions, the key idea is not to compress the original convolution matrix directly, but to apply a translation so that the source dependence is factored into the same multipole coefficients as in the free-space FMM (Cho et al., 2019). In that framework, the source-to-multipole and multipole-to-multipole operators are unchanged from free space, the local-to-local operator is also unchanged, and only the multipole-to-local translation requires layered Sommerfeld integrals (Cho et al., 2019). Asymptotic analysis shows that the layered multipole and local expansions decay with the same geometric ratios as in free space when measured using modified distances that include layer-induced vertical shifts (Cho et al., 2019).

The 2025 extension to Maxwell’s equations in 3-D layered media uses the magnetic vector potential under the Lorenz gauge and represents the dyadic Green’s function through three scalar layered Helmholtz Green’s functions (Yuan et al., 24 Jul 2025). It introduces equivalent polarization images for sources and effective locations for targets so that multiple and local expansions are governed by the actual transmission distance of different reaction-field components (Yuan et al., 24 Jul 2025). To accelerate multipole-to-local translations, the method employs a Chebyshev polynomial expansion of associated Legendre functions, reducing the number of required Sommerfeld-type integrals in M2L tabulation from $$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$8 to $$\nabla\times\nabla\times \bar{\bar G}(\mathbf r,\mathbf r')-k_0^2\epsilon(\mathbf r)\mu(\mathbf r)\bar{\bar G}(\mathbf r,\mathbf r')=\bar{\bar I}\,\delta(\mathbf r-\mathbf r'),$$9 (Yuan et al., 24 Jul 2025). Numerical experiments demonstrate $(G_E,G_H)$0 complexity and rapid convergence for low-frequency electromagnetic sources in 3-D layered media (Yuan et al., 24 Jul 2025).

For direct Green-function evaluation in difficult regimes such as thick lossy media, negative-parameter media, and dominant branch-cut contributions, Padé–Fourier approximation in a conformally mapped spectral plane has been proposed (Hadad, 2021). There, the spectral variable is tilted and then Cayley-mapped so that the transformed Green’s function approaches finite constants at both $(G_E,G_H)$1 and $(G_E,G_H)$2, enabling rational approximation with equal numerator and denominator degrees (Hadad, 2021). This suggests a complementary acceleration strategy to contour deformation and discrete complex image methods when layered spectra exhibit strong pole and branch-cut structure.

6. Asymptotics, low-frequency behavior, and broader extensions

Large-distance and large-order asymptotics play a central role in interpreting layered dyadic Green’s functions. In the scalar translated-FMM analysis, the large-order ratios

$(G_E,G_H)$3

show exponential decay of multipole and local terms, just as in free space, once modified distances are used (Cho et al., 2019). In three-dimensional dyadic formulations, the reaction field can be written as a sum

$(G_E,G_H)$4

isolating angular matrices from radial Sommerfeld integrals and thereby facilitating steepest-descent or stationary-phase analysis (Yuan et al., 2 Jan 2026). Space-wave contributions decay like $(G_E,G_H)$5 in two-dimensional cylindrical far-field asymptotics, while surface or leaky-wave contributions arise from pole residues when present (Yuan et al., 2 Jan 2026).

Low-frequency behavior is a longstanding issue for electric-field dyadic Green’s functions because the homogeneous-space formula contains factors of $(G_E,G_H)$6. The generalized Debye source approach was proposed precisely to avoid this low-frequency breakdown in layered media (O'Neil, 2013). In that formulation, Maxwell fields are represented by scalar densities $(G_E,G_H)$7 and tangential vector fields $(G_E,G_H)$8, and in the static limit the electric and magnetic channels decouple rather than suffering catastrophic cancellation (O'Neil, 2013). The paper states that the resulting layered spectral system has $(G_E,G_H)$9-weighted off-diagonal couplings and block-diagonalizes as $$\nabla\times G_E(\mathbf r,\mathbf r')=-\omega\mu_\ell G_H(\mathbf r,\mathbf r'),\qquad \nabla\times G_H(\mathbf r,\mathbf r')=\omega\epsilon_\ell G_E(\mathbf r,\mathbf r')-\delta(\mathbf r-\mathbf r')\frac{1}{\omega\mu_{\ell'}}\hat t,$$0 (O'Neil, 2013). This addresses a common misconception that low-frequency instability is an unavoidable property of layered electromagnetics; the instability is tied to a particular dyadic representation, not to the layered problem itself.

Layered dyadic Green’s functions also extend beyond electromagnetics. The matrix-basis approach has been carried over to elastic wave equations in layered media, where Maxwell TE/TM decomposition is mirrored by S/P decomposition of the elastic dyadic (Yuan et al., 2 Jan 2026). In 2D elastodynamics, the Green’s tensor

$$\nabla\times G_E(\mathbf r,\mathbf r')=-\omega\mu_\ell G_H(\mathbf r,\mathbf r'),\qquad \nabla\times G_H(\mathbf r,\mathbf r')=\omega\epsilon_\ell G_E(\mathbf r,\mathbf r')-\delta(\mathbf r-\mathbf r')\frac{1}{\omega\mu_{\ell'}}\hat t,$$1

for P–SV motion and the scalar $$\nabla\times G_E(\mathbf r,\mathbf r')=-\omega\mu_\ell G_H(\mathbf r,\mathbf r'),\qquad \nabla\times G_H(\mathbf r,\mathbf r')=\omega\epsilon_\ell G_E(\mathbf r,\mathbf r')-\delta(\mathbf r-\mathbf r')\frac{1}{\omega\mu_{\ell'}}\hat t,$$2 for SH motion have been compressed by Greedy Tucker Approximation with PGD-type alternating least squares, yielding substantial reductions in memory requirements relative to full-order models (Farooq et al., 19 Mar 2026). This broader use reinforces that “layered dyadic Green’s function” is not a narrow Maxwellian construct but a general spectral-response framework for stratified wave systems.

In contemporary usage, the subject encompasses exact spectral representations, TE/TM or matrix-basis decompositions, asymptotic and low-frequency reformulations, and fast hierarchical algorithms. Across these variants, the invariant core is the same: a dyadic response operator whose layered structure is encoded through scalar or polarization-resolved spectral factors, interface matching, and physically chosen branches of the vertical wavenumber (Yuan et al., 2 Jan 2026).