---
title: Layered Dyadic Green's Function Analysis
url: https://www.emergentmind.com/topics/layered-dyadic-green-s-function
type: topic
---

# Layered Dyadic Green's Function Analysis

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 [2601.00709]. 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 [1902.08250].

## 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}$ [1902.08250]. 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 [1902.08250]. 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'$ [2601.00709].

In homogeneous space, the electric-field dyadic Green’s function is
\[
\bar G^E(\mathbf r,\mathbf r')=\left(\bar I+\frac{1}{k^2}\nabla\nabla\right)\frac{e^{ik|\mathbf r-\mathbf r'|}}{4\pi|\mathbf r-\mathbf r'|},
\]
or equivalently
\[
\bar G^E(\mathbf r,\mathbf r')=\frac{1}{k^2}\nabla\times\nabla\times\left(\frac{e^{ik|\mathbf r-\mathbf r'|}}{4\pi|\mathbf r-\mathbf r'|}\bar I\right),
\]
with $k=\omega\sqrt{\epsilon\mu}$ [1310.4241]. 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 [1610.04479].

At each interface one imposes continuity of tangential electric and magnetic fields and continuity of the normal components of $D=\epsilon E$ and $B=\mu H$. In dyadic form one convenient statement is
\[
\mathbf n\times G_E=0,\qquad \mathbf n\cdot(\epsilon G_E)=0,\qquad \mathbf n\times G_H=0,\qquad \mathbf n\cdot(\mu G_H)=0,
\]
with $\mathbf n=\hat z$ for planar layering [2601.00709]. 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 [2008.01047].

## 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 $\lambda$ and vertical propagation constants
\[
\gamma(\lambda;k)=\sqrt{k^2-\lambda^2},\qquad \kappa(\lambda;k)=\sqrt{\lambda^2-k^2},
\]
with branch choice such that $\operatorname{Re}\sqrt{\lambda^2-k^2}\ge 0$ for $|\lambda|>k$ and $\sqrt{\lambda^2-k^2}=i\sqrt{k^2-\lambda^2}$ with nonnegative imaginary part for $|\lambda|<k$ [1902.08250]. The general scalar layered kernel treated in the translated fast multipole framework has the form
\[
G(x,y;x_0,y_0)=\int_{-\infty}^{\infty} e^{-\kappa(\lambda;k)(y+d)} e^{i\lambda x} e^{\pm \kappa(\lambda;k_0)y_0} e^{-i\lambda x_0}\frac{\sigma(\lambda)}{4\pi \kappa(\lambda;k)}\,d\lambda,
\]
where $\sigma(\lambda)$ is an interface density determined by matching conditions and converges to a constant as $|\lambda|\to\infty$ [1902.08250].

In three-dimensional electromagnetics, the dyadic Green’s function admits the spectral representation
\[
G(\mathbf r,\mathbf r')=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{dk_x\,dk_y}{(2\pi)^2}
e^{ik_x(x-x')+ik_y(y-y')}\big[G^{TE}(k_\rho;z,z')+G^{TM}(k_\rho;z,z')\big],
\]
where $k_\rho=\sqrt{k_x^2+k_y^2}$ and the TE/TM pieces are built from polarization projectors and multilayer propagators [2601.00709]. A related formulation writes
\[
\bar{\bar G}(\rho,z;\rho',z')=\int \frac{d^2k_\rho}{(2\pi)^2}e^{ik_\rho\cdot(\rho-\rho')}
\big[\bar{\bar U}_{dir}+\bar{\bar U}_{sc}\big],
\]
with
\[
\bar{\bar U}_{sc}=\frac{1}{2\gamma_\ell}\big[R_s(k_\rho)\bar{\bar P}_s+R_p(k_\rho)\bar{\bar P}_p\big]e^{i\gamma_\ell(z+z')}+\text{transmitted terms},
\]
where $\gamma_\ell=\sqrt{k_\ell^2-k_\rho^2}$, $\operatorname{Im}\gamma_\ell\ge 0$, $\bar{\bar P}_s=\hat s\hat s$, and $\bar{\bar P}_p=\hat p\hat p$ [1902.08250].

This spectral decomposition separates propagating and evanescent sectors. In the scalar formulation, $|\lambda|<k$ yields oscillatory propagating plane waves, while $|\lambda|>k$ yields exponentially decaying evanescent contributions associated with surface and lateral waves [1902.08250]. 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
\[
\hat u=\frac{k_x}{k_\rho}\hat x+\frac{k_y}{k_\rho}\hat y,\qquad
\hat v=\hat z\times \hat u,
\]
the transformed fields are expressed through scalar spectral amplitudes $V_e,V_h,I_e,I_h$ [2601.00709]. In each homogeneous layer the Fourier-transformed Maxwell system decouples into scalar Helmholtz equations for the TM quantity $I_e$ and the TE quantity $V_h$:
\[
\frac{\partial^2 I_e}{\partial z^2}+k_z^2 I_e=k_\rho \hat J_z-\partial_z \hat J_u,\qquad
\frac{\partial^2 V_h}{\partial z^2}+k_z^2 V_h=-\frac{k^2}{\omega\epsilon}\hat J_v,
\]
with $k_z=\sqrt{k^2-k_\rho^2}$ [2601.00709]. The associated scalar Green’s functions $\hat G_1,\hat G_2,\hat G_3$ satisfy layered Helmholtz problems with interface conditions
\[
[\hat G_1]=0,\qquad \left[\frac{1}{\mu}\partial_z \hat G_1\right]=0,
\]
\[
[\hat G_2]=0,\qquad \left[\frac{1}{\epsilon}\partial_z \hat G_2\right]=0,
\]
and $\hat G_3=-\partial_{z'}\hat G_2$ [2601.00709].

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 $G_A$ satisfying
\[
\nabla^2 G_A(\mathbf r,\mathbf r')+k^2G_A(\mathbf r,\mathbf r')=\frac{1}{\omega}\delta(\mathbf r-\mathbf r')I,
\]
and reconstructs
\[
G_E=-\omega\left(I+\frac{\nabla\nabla}{k^2}\right)G_A,\qquad
G_H=\frac{1}{\mu}\nabla\times G_A
\]
under the Lorenz gauge [2601.00709]. The matrix basis $J_1,\dots,J_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 $b_1,b_2,b_3$, with
\[
b_1=\frac{1}{\omega}\hat G_1,\qquad
b_2=\frac{1}{\omega\mu_{\ell'}}\hat G_2,\qquad
b_3=\frac{1}{\omega\mu_{\ell'}}\hat G_3
\]
[2601.00709]. 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
\[
\mathfrak{R}^0=\operatorname{span}_{\mathbb F_0}\{J_1,\ldots,J_5\},
\]
with coefficients independent of the azimuthal angle $\alpha$ [2008.01047]. 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 [2601.00709].

## 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 $a$ and $b$,
\[
r_s=\frac{\mu_b k_{z,a}-\mu_a k_{z,b}}{\mu_b k_{z,a}+\mu_a k_{z,b}},\qquad
t_s=\frac{2\mu_b k_{z,a}}{\mu_b k_{z,a}+\mu_a k_{z,b}},
\]
\[
r_p=\frac{\epsilon_b k_{z,a}-\epsilon_a k_{z,b}}{\epsilon_b k_{z,a}+\epsilon_a k_{z,b}},\qquad
t_p=\frac{2\epsilon_b k_{z,a}}{\epsilon_b k_{z,a}+\epsilon_a k_{z,b}}
\]
[2601.00709]. In multilayers, generalized reflection and transmission coefficients are obtained recursively. For a three-layer slab, one representative generalized reflection coefficient is
\[
\bar R_{12}=\frac{R_{12}+R_{23}e^{2ik_{z2}d}}{1+R_{12}R_{23}e^{2ik_{z2}d}},
\]
with analogous generalized transmission coefficients [1610.04479].

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 $\sigma(\omega)$ 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 [2104.03581]. 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 [2104.03581].

For layered topological insulators, the constitutive relations include the axion term
\[
D=\epsilon_0\epsilon_j(\omega)E+\bar\alpha_j(\omega)B,\qquad
H=\frac{1}{\mu_0\mu_j(\omega)}B-\bar\alpha_j(\omega)E,
\]
with
\[
\bar\alpha_j(\omega)=\frac{\alpha}{\pi}\frac{\theta_j(\omega)}{\mu_0 c},
\]
so the bulk propagation operator is unchanged for constant $\theta$, but interface jumps in $\theta$ modify the boundary conditions and cause TE/TM mixing on reflection and transmission [1509.03012]. In that case, the single-interface reflection operator becomes a $2\times 2$ 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 [1509.03012].

## 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 $g^{R}_{1,5},\dots,g^{R}_{1,8}$, and all transmitted components in the adjacent layer use five integrals $g^{T}_{2,5},\dots,g^{T}_{2,9}$ [1610.04479]. 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 $10^{-10}$ in representative cases [1610.04479].

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 [1902.08250]. 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 [1902.08250]. 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 [1902.08250].

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 [2507.18491]. 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 [2507.18491]. 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 $O(p^4)$ to $O(p^2)$ [2507.18491]. Numerical experiments demonstrate $\mathcal O(N\log N)$ complexity and rapid convergence for low-frequency electromagnetic sources in 3-D layered media [2507.18491].

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 [2111.09089]. There, the spectral variable is tilted and then Cayley-mapped so that the transformed Green’s function approaches finite constants at both $\eta\to 0$ and $|\eta|\to\infty$, enabling rational approximation with equal numerator and denominator degrees [2111.09089]. 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
\[
\frac{|J_{p+1}(k_0 r)\Phi_{p+1}(x,y)|}{|J_p(k_0 r)\Phi_p(x,y)|}\approx \frac{r}{\rho},
\qquad
\frac{|J_{p+1}(k\tilde r)\Psi_{p+1}(x_0,y_0)|}{|J_p(k\tilde r)\Psi_p(x_0,y_0)|}\approx \frac{\tilde r}{\tilde \rho},
\]
show exponential decay of multipole and local terms, just as in free space, once modified distances are used [1902.08250]. In three-dimensional dyadic formulations, the reaction field can be written as a sum
\[
G\approx \sum_\kappa I_\kappa(\rho;z,z')\,M_\kappa,
\]
isolating angular matrices from radial Sommerfeld integrals and thereby facilitating steepest-descent or stationary-phase analysis [2601.00709]. Space-wave contributions decay like $1/\sqrt{\rho}$ in two-dimensional cylindrical far-field asymptotics, while surface or leaky-wave contributions arise from pole residues when present [2601.00709].

Low-frequency behavior is a longstanding issue for electric-field dyadic Green’s functions because the homogeneous-space formula contains factors of $1/k^2$. The generalized Debye source approach was proposed precisely to avoid this low-frequency breakdown in layered media [1310.4241]. In that formulation, Maxwell fields are represented by scalar densities $r,q$ and tangential vector fields $\mathbf J,\mathbf K$, and in the static limit the electric and magnetic channels decouple rather than suffering catastrophic cancellation [1310.4241]. The paper states that the resulting layered spectral system has $\omega$-weighted off-diagonal couplings and block-diagonalizes as $\omega\to 0$ [1310.4241]. 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 [2601.00709]. In 2D elastodynamics, the Green’s tensor
\[
G(x,z;x_s,z_s;\omega)=
\begin{bmatrix}
G_{xx} & G_{xz}\\
G_{zx} & G_{zz}
\end{bmatrix}
\]
for P–SV motion and the scalar $G_{yy}$ 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 [2603.19080]. 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 [2601.00709].

Source: https://www.emergentmind.com/topics/layered-dyadic-green-s-function