Dressed Global Polarizability Matrix (dGPM)

• dGPM is a reduced-order framework that represents complex meta-atoms via a compact operator learned directly in stratified media.
• It embeds interface-mediated feedback into the operator, enabling accurate reconstruction of near and far-field scattering with efficient numerical dipole representation.
• The method offers significant speed and memory improvements over full-wave simulations, making it ideal for large-scale metasurface design and parameter sweeps.

Searching arXiv for the primary dGPM paper and its GPM precursor to ground the article in the relevant literature. The dressed Global Polarizability Matrix (dGPM) is a reduced-order electromagnetic framework for modeling complex meta-atoms in heterogeneous photonic environments, especially when scatterers touch, straddle, or are embedded within material interfaces. In this formulation, a complex scatterer is replaced by a compact operator acting on a small set of numerical electric and magnetic dipoles, while interface-mediated feedback is absorbed directly into the operator through learning in the stratified environment rather than being treated only during the subsequent many-body solve. The framework was introduced as an extension of the Global Polarizability Matrix (GPM), which had been formulated for dense ensembles of non-spherical particles in stratified media but calibrated in a homogeneous background (Fu et al., 19 Jun 2026, Bertrand et al., 2019).

1. Definition and conceptual scope

The dGPM extends the original GPM to heterogeneous photonic environments in which meta-atoms interact strongly with nearby planar interfaces. In both GPM and dGPM, a complex scatterer is represented by a compact operator that maps a small basis of numerical electric and magnetic dipoles to effective induced moments, from which the scattered field is reconstructed through appropriate Green tensors (Fu et al., 19 Jun 2026). The essential distinction is that the GPM is learned in a homogeneous background and can be freely translated and rotated in three dimensions, whereas the dGPM is learned directly in the stratified environment and therefore incorporates interface-mediated feedback into the operator itself (Fu et al., 19 Jun 2026).

This distinction is operational rather than merely terminological. In the GPM formalism, environmental effects are included during multiple-scattering simulation through the dyadic Green tensor of the stratified medium and by subtracting the background contribution already encoded in the learned operator (Bertrand et al., 2019). In the dGPM formalism, by contrast, direct, reflected, transmitted, and evanescent contributions from nearby planar interfaces are absorbed into the learned operator during training, so that self-interactions mediated by the interface are transferred from the coupling operator into the dressed scatterer description (Fu et al., 19 Jun 2026).

A plausible implication is that the dGPM should be understood as a renormalized scatterer representation specialized to layered environments. The paper makes this intuition explicit by noting that one may think heuristically of a renormalization of the bare polarizability by interface-induced self-energy, although the dressed operator is learned directly rather than assembled from an explicit self-energy term (Fu et al., 19 Jun 2026).

2. Operator structure and electromagnetic representation

In the dGPM, each meta-atom is represented by $N_p$ numerical dipoles, each carrying six moment channels: three electric and three magnetic. The induced-moment vector at numerical dipole $n$ is

$\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$

and the corresponding exciting field vector is

$\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$

(Fu et al., 19 Jun 2026).

The single-particle dGPM is a $6N_p \times 6N_p$ matrix $\tilde{\mathbf{A}}$ organized into $6\times 6$ blocks $\tilde{\mathbf{A}}_{m\to n}$:

$\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$

with

$$\tilde{\mathbf{A}}_{m\to n}\in\mathbb{C}^{6\times 6}$$

(Fu et al., 19 Jun 2026). Unlike the discrete dipole approximation, these numerical dipoles are not local physical dipoles obeying local constitutive relations. The operator is generally fully populated, and its off-diagonal blocks encode nonlocal internal couplings inside the scatterer (Fu et al., 19 Jun 2026). This nonlocality is central to both GPM and dGPM: it is what allows a small number of numerical dipoles to reproduce the external electromagnetic response of complex, high-index, high-aspect-ratio meta-atoms (Fu et al., 19 Jun 2026, Bertrand et al., 2019).

Higher-order multipolar physics is not introduced through an explicit quadrupolar or higher multipole basis. Instead, it is captured implicitly through the nonlocal coupling between numerical sites and through the layered Green tensors (Fu et al., 19 Jun 2026). This suggests that the compactness of the representation does not derive from truncating the physics to dipolar order in the usual sense, but from compressing the scatterer response into an operator basis calibrated against full-wave data.

3. Learning in a stratified environment

The dressed operator is retrieved from full-wave data in the stratified medium using source and learning surfaces. The construction employs a source vector $n$0, a stratified Green matrix $n$1 mapping the source surface to the numerical dipoles, and a stratified Green matrix $n$2 mapping the induced moments at the numerical dipoles to the learning surface (Fu et al., 19 Jun 2026). The dGPM-predicted scattered field on the learning surface is

$n$3

and the dressed operator is retrieved through the regularized least-squares problem

$n$4

(Fu et al., 19 Jun 2026).

The precomputation pipeline reported for a single meta-atom consists of three steps: generating the reference dataset with COMSOL using random electric and magnetic dipolar sources on a source surface, building the stratified Green matrices $n$5 and $n$6, and solving the inverse problem for $n$7 (Fu et al., 19 Jun 2026). Typical runtime for $n$8 is reported as approximately $n$9–$\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$0 minutes with fine meshes ensuring negligible discretization error (Fu et al., 19 Jun 2026).

For the ZnS-capped GaAs pillar at $\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$1, the reported training configuration used $\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$2, $\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$3 sources with six polarizations each, and approximately $\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$4 learning points. This yielded mean relative error $\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$5 with standard deviation $\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$6 on training and $\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$7 with standard deviation $\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$8 on generalization (Fu et al., 19 Jun 2026). The error metric is given explicitly by

$\mathbf{P}_n = [p_x,p_y,p_z,m_x,m_y,m_z]^\top,$9

with variance

$\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$0

(Fu et al., 19 Jun 2026).

The 2019 GPM precursor used an inverse procedure with two-stage least-squares inversion based on SVD pseudoinverses, learning the operator from scattered fields sampled on a learning surface in a homogeneous background (Bertrand et al., 2019). That earlier work did not explicitly coin the term dGPM, but it established the compact nonlocal operator paradigm that the 2026 work extends to interface-coupled training (Bertrand et al., 2019).

4. Multiple scattering, Green operators, and formal dressing

For a collection of meta-atoms, the many-body system is assembled through a block-diagonal operator

$\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$1

where each block may be either a dGPM for an interface-coupled scatterer or a GPM for a free-space scatterer (Fu et al., 19 Jun 2026). The self-consistent multiple-scattering equation is

$\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$2

or equivalently

$\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$3

(Fu et al., 19 Jun 2026).

The inter-particle coupling operator is decomposed as

$\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$4

where $\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$5 contains the direct contributions within a homogeneous layer and $\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$6 contains reflected, transmitted, and evanescent couplings (Fu et al., 19 Jun 2026). For planar interfaces, the interface-mediated Green tensor may be written in spectral form using Sommerfeld integrals over in-plane wavevectors and an $\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$7 (TE/TM) decomposition (Fu et al., 19 Jun 2026). In practice, the source-to-dipole, dipole-to-learning-surface, and many-body Green operators are evaluated with the stratified-medium formalism of Ref. [24] in the 2026 paper, so that all direct, reflected, transmitted, and evanescent components are included (Fu et al., 19 Jun 2026).

The conceptual dressing relation is written as

$\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$8

where $\mathbf{F}_n = [E_x,E_y,E_z,H_x,H_y,H_z]^\top$9 denotes the interface-mediated self-interaction (Fu et al., 19 Jun 2026). The paper emphasizes, however, that the dressed operator is not computed from an explicit $6N_p \times 6N_p$0; instead, the effect of $6N_p \times 6N_p$1 is encoded implicitly through the stratified Green operators used during training (Fu et al., 19 Jun 2026).

This point connects directly to the 2019 GPM formalism. There, for a single particle in a stratified medium, the environmental contribution appears through

$6N_p \times 6N_p$2

with the equivalent dressed form

$6N_p \times 6N_p$3

derived as an operational consequence of the self-consistent equation (Bertrand et al., 2019). The 2026 dGPM may therefore be viewed as an explicit learning-based realization of dressing that had been implicit in the earlier GPM-plus-environment formulation.

5. Observables, benchmarks, and computational performance

Once the induced moments are known, scattered near fields at a point $6N_p \times 6N_p$4 are reconstructed by applying the stratified dyadic Green tensors:

$6N_p \times 6N_p$5

$6N_p \times 6N_p$6

(Fu et al., 19 Jun 2026). The framework reconstructs fields anywhere outside the learning surfaces; inside the learning surface, point-dipole singularities appear near numerical dipoles, but this does not affect reconstruction outside (Fu et al., 19 Jun 2026).

Far fields are obtained through a near-to-far transformation for layered media, with RETOP given as an example (Fu et al., 19 Jun 2026). For the ZnS-capped GaAs pillar at $6N_p \times 6N_p$7, angular radiation diagrams match COMSOL to within $6N_p \times 6N_p$8 integrated relative error and reproduce features including critical-angle cutoffs and evanescent-to-propagating conversion (Fu et al., 19 Jun 2026). Reflection and transmission cross sections $6N_p \times 6N_p$9 are computed by normalizing radiated powers by the incident flux, while absorption for lossy meta-atoms is evaluated by energy balance from $\tilde{\mathbf{A}}$0, thereby avoiding volumetric integration (Fu et al., 19 Jun 2026).

The published benchmarks span single particles, dimers, and hybrid ensembles. For a single interface-straddling pillar, the reported training error is $\tilde{\mathbf{A}}$1 with standard deviation $\tilde{\mathbf{A}}$2, and generalization error is $\tilde{\mathbf{A}}$3 with standard deviation $\tilde{\mathbf{A}}$4; near-field component errors are reported as low as $\tilde{\mathbf{A}}$5 across all six components, and the far-field integrated error is $\tilde{\mathbf{A}}$6 (Fu et al., 19 Jun 2026). For a dimer at separation $\tilde{\mathbf{A}}$7 and incidence angle $\tilde{\mathbf{A}}$8, the mean relative error in near fields is reported as $\tilde{\mathbf{A}}$9 for dominant $6\times 6$0 and $6\times 6$1 across six components, with $6\times 6$2 far-field error (Fu et al., 19 Jun 2026). Dimer cross-section benchmarks show mean relative errors below $6\times 6$3 for separations $6\times 6$4, rising to $6\times 6$5 at $6\times 6$6 and $6\times 6$7 at $6\times 6$8 (Fu et al., 19 Jun 2026).

The reported speed and memory advantages over full-wave simulation are substantial. For the dimer benchmark, COMSOL requires approximately $6\times 6$9 minutes and $\tilde{\mathbf{A}}_{m\to n}$0 MB per angle, while dGPM requires approximately $\tilde{\mathbf{A}}_{m\to n}$1 s per angle after approximately $\tilde{\mathbf{A}}_{m\to n}$2 s initialization with approximately $\tilde{\mathbf{A}}_{m\to n}$3 MB, corresponding to approximately $\tilde{\mathbf{A}}_{m\to n}$4 speed and approximately $\tilde{\mathbf{A}}_{m\to n}$5 memory reduction (Fu et al., 19 Jun 2026). In the precursor GPM paper, a dense ensemble of $\tilde{\mathbf{A}}_{m\to n}$6 cylinders in a metallo-dielectric stack required approximately $\tilde{\mathbf{A}}_{m\to n}$7 GB for the COMSOL reference and approximately $\tilde{\mathbf{A}}_{m\to n}$8 MB for the MATLAB GPM solution (Bertrand et al., 2019). These figures establish the method as a reduced-order surrogate for repeated scattering calculations rather than a direct replacement for full-wave training.

6. Comparison with T-matrix and other scattering formalisms

A central argument for dGPM is that conventional T-matrix methods rely on vector spherical wave expansions outside circumscribing spheres. For high-aspect-ratio meta-atoms such as micropillars, the circumscribing spheres can be much larger than the physical scatterers, and at practical separations the spheres may overlap, causing loss of convergence and preventing near-gap field resolution (Fu et al., 19 Jun 2026). Standard T-matrix plus S-matrix combinations also cannot handle scatterers straddling interfaces (Fu et al., 19 Jun 2026).

The dGPM addresses these issues in three ways stated explicitly in the 2026 paper: it uses conformal learning surfaces, allowing close packing without circumscribing-sphere overlap; it learns directly in layered media, embedding reflection, transmission, and evanescent feedback in the operator; and it preserves reduced-order compactness and portability for translations parallel to the interface and rotations around the interface normal while remaining compatible with many-body Green-tensor multiple scattering (Fu et al., 19 Jun 2026).

The GPM literature had already positioned the operator approach relative to discrete dipole approximations, T-matrix methods, and BEM/FEM. The GPM uses a small number of non-locally coupled numerical dipoles with $\tilde{\mathbf{A}}_{m\to n}$9 degrees of freedom rather than local dipoles with $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$0 degrees of freedom, and thereby can achieve accurate scattering predictions for high-index, resonant, non-spherical particles with dramatically fewer unknowns than DDA (Bertrand et al., 2019). In the 2019 benchmarks, a resonant Si cylinder represented by $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$1 numerical dipoles achieved approximately $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$2 error in scattering cross section with $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$3 s runtime and $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$4 MB, whereas DDSCAT with $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$5 dipoles achieved $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$6 error with $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$7 s runtime and $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$8 MB (Bertrand et al., 2019). The 2019 paper also reports that T-matrix methods fail badly for two cylinders with strong near-field coupling and touching learning surfaces (Bertrand et al., 2019).

A plausible implication is that dGPM is best situated between full-wave solvers and traditional multipole formalisms: more specialized than general FEM or BEM, but more robust than circumscribing-sphere-based T-matrix schemes for interface-coupled, high-aspect-ratio, and densely packed structures.

7. Applications, assumptions, and extensions

The framework is intended for rapid parameter sweeps and large-scale design of ordered and disordered metasurfaces. The 2026 paper states that one may precompute operators once, possibly in shared libraries, assemble heterogeneous scenes mixing buried, suspended, and interfacial meta-atoms, perform angle and frequency sweeps, and evaluate scattering, reflection, transmission, absorption, and radiation patterns efficiently (Fu et al., 19 Jun 2026). It further states that the operator-based solver can be integrated with automatic differentiation, with TorchGDM given as an example, to support inverse-design pipelines and GPU-accelerated optimization (Fu et al., 19 Jun 2026).

The formal assumptions are linear, time-harmonic electromagnetism, layered planar environments, and material dispersion as represented in the full-wave training stage (Fu et al., 19 Jun 2026). Formal validity additionally requires non-overlap of learning surfaces across different meta-atoms (Fu et al., 19 Jun 2026). Several limitations are identified explicitly. Extreme near-touching configurations with significant learning-surface overlap degrade accuracy; very high-order multipolar content may require increasing $\mathbf{P}=\tilde{\mathbf{A}}\,\mathbf{F},$9 or refining dipole placement; strongly nonlocal or spatially dispersive materials are not explicitly modeled; and symmetry portability is reduced relative to GPM, since only translations parallel to the interface and rotations about the interface normal are retained (Fu et al., 19 Jun 2026).

The implementation scaling is also spelled out. Assembly of the coupling operator scales approximately as $$\tilde{\mathbf{A}}_{m\to n}\in\mathbb{C}^{6\times 6}$$0 in the number of meta-atoms because of pairwise coupling, dense direct solves scale approximately as $$\tilde{\mathbf{A}}_{m\to n}\in\mathbb{C}^{6\times 6}$$1, and memory scales approximately as $$\tilde{\mathbf{A}}_{m\to n}\in\mathbb{C}^{6\times 6}$$2 (Fu et al., 19 Jun 2026). In ordered arrays, identical Green blocks can be reused, making initialization cheaper and leaving the dense solve dominant from about $$\tilde{\mathbf{A}}_{m\to n}\in\mathbb{C}^{6\times 6}$$3; in disordered arrays, initialization dominates until about $$\tilde{\mathbf{A}}_{m\to n}\in\mathbb{C}^{6\times 6}$$4 (Fu et al., 19 Jun 2026). In the reported MATLAB implementation without GPU acceleration, systems with hundreds to thousands of meta-atoms are tractable (Fu et al., 19 Jun 2026).

The extensions proposed in the 2026 paper include nonlinearities via augmented learned operators, time modulation, anisotropic substrates, multilayer stacks already supported by stratified Green tensors, and non-planar interfaces provided suitable Green operators are available (Fu et al., 19 Jun 2026). Straightforward directions for improved scalability are also enumerated: iterative solvers, preconditioning, GPU implementations, and hierarchical low-rank compression of the coupling matrix $$\tilde{\mathbf{A}}_{m\to n}\in\mathbb{C}^{6\times 6}$$5 (Fu et al., 19 Jun 2026).

Taken together, these features define the dGPM as a compact operator formalism in which interface-mediated self-interactions are learned rather than appended, enabling accurate and efficient simulation of large collections of interface-coupled meta-atoms in layered media (Fu et al., 19 Jun 2026). The broader significance lies in extending operator-based electromagnetic modeling from homogeneous backgrounds to heterogeneous photonic environments while preserving the compact scatterer-based representation established by the GPM methodology (Fu et al., 19 Jun 2026, Bertrand et al., 2019).