GF_NLayer

N-layer Sommerfeld dyadic Green’s function for planar dielectric stacks.

The implementation follows the Tomas/Welsch convention used in local/derivation/multi_layer_GF.tex. Layers are ordered from bottom to top. The first and last layers are normally semi-infinite; finite internal layers carry a physical thickness. Source and observer are assumed to lie in the same layer, with their z-coordinates measured upward from that layer’s lower interface.

class mqed.Dyadic_GF.GF_NLayer.LayerSpec(epsilon: complex, thickness_m: float | None, name: str = 'layer')

One layer in a bottom-to-top planar stack.

Parameters:

• epsilon – Complex relative permittivity at the active frequency.

• thickness_m – Layer thickness in meters. Use None for semi-infinite boundary layers.

• name – Human-readable layer label.

__delattr__(name)

Implement delattr(self, name).

__eq__(other)

Return self==value.

__init__(epsilon: complex, thickness_m: float | None, name: str = 'layer') → None

__repr__()

Return repr(self).

__setattr__(name, value)

Implement setattr(self, name, value).

class mqed.Dyadic_GF.GF_NLayer.NLayerGreenFunction(layers: Sequence[LayerSpec], source_layer: int, omega: float, qmax: float | None = None, epsabs: float = 1e-09, epsrel: float = 1e-09, limit: int = 400, split_propagating: bool = False, integration_method: str = 'direct', dcim_q_start: float = 0.0, dcim_q_stop: float | None = None, dcim_sample_count: int = 128, dcim_image_count: int = 16, hybrid_direct_q_stop: float | None = None, hybrid_tail_q_stop: float | None = None, hybrid_sample_count: int | None = None, hybrid_image_count: int | None = None, hybrid_validation_rtol: float = 0.05, hybrid_validate_tail: bool = True, hybrid_fallback_to_direct: bool = True, fixed_grid_points: int = 2048)

Dyadic Green’s function for source/observer in one layer of an N-layer stack.

The scattered tensor is assembled from seven Sommerfeld integrals. All multilayer physics enters through the generalized upward/downward Fresnel coefficients r_plus and r_minus computed by recursive Airy formulas.

For a source and observer in layer j, the total Green tensor is

\[\mathbf G(\mathbf r,\mathbf r_0,\omega) = \mathbf G_0^{(j)}(\mathbf r,\mathbf r_0,\omega) + \mathbf G_{\mathrm{sc}}^{(j)}(\rho,\phi,z,z_0,\omega).\]

The scattered term has the angular-spectrum form

\[\mathbf G_{\mathrm{sc}}^{(j)} = \frac{i}{4\pi} \left(\mathbf M_s + \mathbf M_p\right),\]

where M_s and M_p are assembled from seven scalar Sommerfeld integrals. The layer stack enters those integrals only through recursive effective mirrors r_minus and r_plus and the source-layer Airy denominator

\[D_\sigma(q) = 1-r_{j-}^{\sigma}(q)r_{j+}^{\sigma}(q) \exp(2i\beta_j d_j), \quad \sigma\in\{s,p\}.\]

Parameters:

• layers – Bottom-to-top layer stack.

• source_layer – Zero-based index of the layer containing source and observer. The layer must have finite thickness unless it is used in a one-sided limit.

• omega – Angular frequency in rad/s.

• qmax – Optional finite upper limit for the Sommerfeld integral.

• epsabs – Absolute quadrature tolerance.

• epsrel – Relative quadrature tolerance.

• limit – Maximum number of quad_vec subintervals.

• split_propagating – If true, split the integration at the source-layer light line sqrt(eps_j) * k0 when it is real.

__init__(layers: Sequence[LayerSpec], source_layer: int, omega: float, qmax: float | None = None, epsabs: float = 1e-09, epsrel: float = 1e-09, limit: int = 400, split_propagating: bool = False, integration_method: str = 'direct', dcim_q_start: float = 0.0, dcim_q_stop: float | None = None, dcim_sample_count: int = 128, dcim_image_count: int = 16, hybrid_direct_q_stop: float | None = None, hybrid_tail_q_stop: float | None = None, hybrid_sample_count: int | None = None, hybrid_image_count: int | None = None, hybrid_validation_rtol: float = 0.05, hybrid_validate_tail: bool = True, hybrid_fallback_to_direct: bool = True, fixed_grid_points: int = 2048)

static _beta_phys(eps: complex, k0: float, q)

Physical z-wavevector branch for one homogeneous layer.

\[\beta(q) = \sqrt{\epsilon k_0^2-q^2+i0^+}.\]

The branch is chosen with Im(beta) >= 0 so evanescent waves decay away from interfaces. For nearly real propagating waves, Re(beta) is kept non-negative.

_effective_reflection(layer: int, step: int, q, polarization: str)

Recursive Airy reflection coefficient for an arbitrary sub-stack.

For a step toward the next layer n = layer + step, the effective reflection is

\[\widetilde R_{l,n}^{\sigma}= \frac{R_{l,n}^{\sigma}+\widetilde R_{n,n+step}^{\sigma} \exp(2i\beta_n d_n)} {1+R_{l,n}^{\sigma}\widetilde R_{n,n+step}^{\sigma} \exp(2i\beta_n d_n)}.\]

The recursion terminates at the semi-infinite outer layer, where the effective reflection is the bare Fresnel coefficient.

_finite_phase(layer: int, q)

Round-trip propagation factor through a finite layer.

\[P_l(q)=\exp(2i\beta_l(q)d_l).\]

This factor appears in the recursive effective-reflection formula.

_fresnel(layer_i: int, layer_j: int, q, polarization: str)

Bare Fresnel reflection coefficient from layer i to layer j.

TE polarization:

\[R_{ij}^{s}=\frac{\beta_i-\beta_j}{\beta_i+\beta_j}.\]

TM polarization:

\[R_{ij}^{p}=\frac{\epsilon_j\beta_i-\epsilon_i\beta_j} {\epsilon_j\beta_i+\epsilon_i\beta_j}.\]

_kz(layer: int, q)

Layer-indexed vertical wavevector.

\[\beta_l(q)=\sqrt{\epsilon_l k_0^2-q^2+i0^+}.\]

Parameters:

• layer – Zero-based layer index in the bottom-to-top stack.

• q – In-plane Sommerfeld wave number.

airy_denominator(q, polarization: str)

Source-layer multiple-reflection denominator.

\[D_\sigma(q)=1-r_{j-}^{\sigma}(q)r_{j+}^{\sigma}(q) e^{2i\beta_j(q)d_j}.\]

Zeros of this denominator are guided-mode or plasmonic pole candidates.

amplitude_coefficients(q, z_observer: float, z_source: float)

Five scalar amplitudes entering the N-layer scattered tensor.

For source and observer in layer j of thickness d_j:

\[C_s = \frac{e^{i\beta_j(z+z_0-d_j)}r_{j-}^s + e^{-i\beta_j(z+z_0-d_j)}r_{j+}^s + 2r_{j+}^sr_{j-}^s\cos[\beta_j(z-z_0)]e^{i\beta_j d_j}} {D_s}.\]

TM amplitudes use the same direct reflection part and opposite signs for the round-trip term:

\[C_p^{\pm} = \frac{B_p \pm 2r_{j+}^pr_{j-}^p \cos[\beta_j(z-z_0)]e^{i\beta_j d_j}}{D_p},\]

\[S_p^{\pm} = \frac{A_p \pm 2i r_{j+}^pr_{j-}^p \sin[\beta_j(z-z_0)]e^{i\beta_j d_j}}{D_p},\]

with B_p = exp(i beta_j(z+z0-d_j)) r_minus_p + exp(-i beta_j(z+z0-d_j)) r_plus_p and A_p using the same two terms with a minus sign between them.

bessel_free_kernels(q: float, z_observer: float, z_source: float) → numpy.ndarray

Seven Bessel-free scalar kernels before angular integration.

The returned vector is

\[\left[\frac{C_s q e^{i\beta_jd_j}}{2\beta_j}, \frac{C_s q e^{i\beta_jd_j}}{2\beta_j}, \frac{C_p^-q\beta_j e^{i\beta_jd_j}}{2k_j^2}, \frac{C_p^-q\beta_j e^{i\beta_jd_j}}{2k_j^2}, \frac{iS_p^+q^2 e^{i\beta_jd_j}}{k_j^2}, \frac{iS_p^-q^2 e^{i\beta_jd_j}}{k_j^2}, \frac{C_p^+q^3 e^{i\beta_jd_j}}{\beta_jk_j^2}\right].\]

Multiplication by [J0,J2,J0,J2,J1,J1,J0] gives the seven Sommerfeld integrands.

branch_points() → numpy.ndarray

Layer light-line branch points in the complex q plane.

\[q_l = \sqrt{\epsilon_l}\,k_0.\]

These branch points mark sheet changes of beta_l(q) and guide the low-q split used by the hybrid DCIM path.

calculate_total_Green_function(x: float, y: float, z_observer: float, z_source: float)

Total source-layer dyadic Green tensor.

\[\mathbf G = \mathbf G_0^{(j)} + \mathbf G_{sc}^{(j)}.\]

complex_quad(func: Callable[[float], complex | numpy.ndarray], a: float, qmax: float | None = None, epsabs: float | None = None, epsrel: float | None = None, limit: int | None = None, split_propagating: bool | None = None) → complex | numpy.ndarray

Vector-valued Sommerfeld quadrature helper.

Computes

\[\int_a^{q_{\max}} f(q)\,dq\]

with scipy.integrate.quad_vec. When requested, the integral is split at the real source-layer light line to separate propagating and evanescent parts.

compute_integrals(rho: float, z_observer: float, z_source: float) → numpy.ndarray

Compute the seven Sommerfeld integrals for the selected method.

Direct mode evaluates

\[I_m(\rho)=\int_0^{q_{\max}}F_m(q)J_{\nu_m}(q\rho)\,dq.\]

dcim performs the legacy whole-domain exponential fit. hybrid_dcim evaluates the low-q interval directly and fits only the evanescent tail.

compute_integrals_dcim(rho: float, z_observer: float, z_source: float) → numpy.ndarray

Legacy whole-domain complex-image approximation.

This path fits

\[F_m(q)\approx\sum_n a_{mn}e^{-s_{mn}q}, \quad 0\le q\le q_f,\]

then applies analytic Laplace-Bessel transforms on [0,infinity). It is retained as a diagnostic because metal multilayer kernels usually require the split-tail treatment implemented by hybrid_dcim.

compute_integrals_hybrid_dcim(rho: float, z_observer: float, z_source: float) → numpy.ndarray

Hybrid direct/GPOF-DCIM evaluation of the Sommerfeld integrals.

The integral is split as

\[I_m = \int_0^{q_c}F_m(q)J_{\nu_m}(q\rho)dq + \int_{q_c}^{\infty}F_m(q)J_{\nu_m}(q\rho)dq.\]

The first term is direct quadrature. On the tail, samples on [q_c,q_f] are fit to shifted images

\[F_m(q)\approx\sum_{n=1}^{N_i}a_{mn}e^{-s_{mn}(q-q_c)}.\]

The fitted tail is accepted only if its finite interval integral over [q_c,q_f] agrees with direct finite-tail quadrature within hybrid_validation_rtol; otherwise the method records diagnostics and falls back to direct quadrature when configured to do so.

direct_integrals_over_range(rho: float, z_observer: float, z_source: float, lower: float, upper: float, epsabs: float | None = None, epsrel: float | None = None, limit: int | None = None) → numpy.ndarray

Directly integrate the seven scalar Sommerfeld integrals on a finite interval.

For Bessel orders nu = [0,2,0,2,1,1,0] and Bessel-free kernels F_m(q), this returns

\[I_m(\rho; q_a,q_b)=\int_{q_a}^{q_b}F_m(q)J_{\nu_m}(q\rho)\,dq.\]

This finite-interval form is used by both the validated direct solver and the hybrid DCIM finite-tail validation.

reflection_coefficient(q, direction: str, polarization: str)

Effective mirror seen from the source layer.

direction='down' returns r_minus and direction='up' returns r_plus for either TE (s) or TM (p) polarization.

scatter_component(x: float, y: float, z_observer: float, z_source: float)

Scattered Green tensor assembled from TE and TM angular spectra.

\[\mathbf G_{sc}^{(j)}=\frac{i}{4\pi} \left[\mathbf M_s(I_1,I_2)+\mathbf M_p(I_3,\ldots,I_6)\right].\]

scattering_p_component(x: float, y: float, integrals: numpy.ndarray)

TM tensor block after angular integration.

\[\begin{split}\mathbf M_p=\begin{pmatrix} -I_3+\cos2\phi I_4 & \sin2\phi I_4 & -\cos\phi I_5^+\\ \sin2\phi I_4 & -I_3-\cos2\phi I_4 & -\sin\phi I_5^+\\ \cos\phi I_5^- & \sin\phi I_5^- & I_6 \end{pmatrix}.\end{split}\]

In an N-layer stack, I5+ and I5- are generally different because the TM amplitudes S_p_plus and S_p_minus are different.

scattering_s_component(x: float, y: float, integrals: numpy.ndarray)

TE tensor block after angular integration.

\[\begin{split}\mathbf M_s=\begin{pmatrix} I_1+\cos2\phi I_2 & \sin2\phi I_2 & 0\\ \sin2\phi I_2 & I_1-\cos2\phi I_2 & 0\\ 0 & 0 & 0 \end{pmatrix}.\end{split}\]

vacuum_component(x: float, y: float, z_observer: float, z_source: float)

Homogeneous Green tensor inside the source layer.

\[\mathbf G_0^{(j)}(\mathbf R)=\frac{e^{ik_jR}}{4\pi R k_j^2} \left[(\mathbf I-\hat{R}\hat{R})k_j^2 +\frac{3\hat{R}\hat{R}-\mathbf I}{R^2} +\frac{i k_j(\mathbf I-3\hat{R}\hat{R})}{R}\right].\]

At coincident points the usual imaginary self term i k_j/(6 pi) I is returned.