Multipole-Frequency Decomposition

• Multipole-frequency decomposition is a frequency-domain representation that isolates complex oscillatory behavior into a few multipolar channels, poles, or rational functions.
• It underpins metallic GW methods by reducing dense frequency sampling to a limited number of complex poles, achieving quasiparticle energy accuracy within a few meV compared to full-frequency techniques.
• The approach accelerates boundary-element and nano-optics simulations by decoupling phase factors, enabling efficient extraction of mode densities, scattering spectra, and optimal illumination conditions.

Multipole-frequency decomposition refers, in the literature summarized here, to a set of frequency-domain constructions in which a response function, operator family, or scattered field is represented through multipolar channels, complex poles, or analytically extracted oscillatory factors. In metallic many-body perturbation theory it appears as the multipole approximation for dielectric screening in full-frequency $GW$; in boundary-element discretizations of the Helmholtz equation it appears as frequency extraction through a Hadamard-factorized phase; and in nano-optics and photonics it appears as exact field-, source-, and current-based multipole expansions that resolve frequency-dependent resonances, mode densities, and illumination dependence (Leon et al., 2023, Dirckx et al., 2020, Majorel et al., 2022, Loulas et al., 2024, Shevchenko et al., 22 Aug 2025, Hernández-Sarria et al., 23 Oct 2025).

1. Recurring mathematical structure

Across these domains, the common operation is not a single universal formula but a related pattern: the explicitly frequency-dependent part of a problem is isolated into a representation that is cheaper to sample, easier to interpret, or directly tied to modal content. In some cases the reduced representation is a pole sum, in some a phase-extracted matrix family, and in some a multipole basis expansion of radiated or induced fields.

Setting	Core representation	Stated role
Metallic $GW$	$$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$ or $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$	Accurate and efficient alternative to current full-frequency methods
Helmholtz BEM	$B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$ with $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$ in the far-field	Data-sparse representation valid for a range of frequencies
Generalized polarizabilities	$$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$	Exact multipole decomposition for any illumination
2-D photonics	$$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$	Full frequency-dependent multipole decomposition and scattering spectrum
Exact current multipoles	$$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$	Exact and general current-moment description beyond the point-multipole approximation

This comparison suggests an umbrella usage of the term: a multipole-frequency decomposition is any formulation that rewrites a frequency-dependent problem into a small set of poles, harmonics, rational functions, or multipolar amplitudes whose variation with $\omega$ or $GW$0 can be sampled or interpreted directly (Leon et al., 2023, Dirckx et al., 2020, Majorel et al., 2022, Loulas et al., 2024, Shevchenko et al., 22 Aug 2025).

2. Full-frequency $GW$1 for metals: multipole approximation of screening

In the $GW$2 approximation, the self-energy is

$GW$3

with the screened Coulomb interaction obtained from the Dyson equation for the polarizability. In reciprocal space one computes the irreducible RPA polarizability $GW$4, solves $GW$5, and forms $GW$6. The metallic multipole approximation replaces dense real-frequency sampling by a small number $GW$7 of complex poles, either in a Lehmann-like form for $GW$8 or in a linear-pole form for $GW$9, and the two forms are stated to be fully equivalent for a given $$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$0 (Leon et al., 2023).

For metals, the numerical procedure uses a “double-parallel” sampling of $$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$1 complex frequencies,

$$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$2

with typical $$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$3, $$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$4, and a nonuniform “power-$$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$5” grid

$$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$6

for metals. The nonlinear fit for $$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$7 is solved by linearization, with an iterative “secant” or Newton-Raphson scheme in $$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$8 and a linear least-squares solution for $$\Pi_{GG'}(q,\omega)\approx \sum_{p=1}^{N_p}\dfrac{2\,R_p^{GG'}(q)\,\Omega_p(q)}{\omega^2-[\Omega_p(q)]^2}$$9 at each step; the stated convergence criterion is a residual $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$0. The metallic case also requires a treatment of the $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$1 intraband contribution. Two ab-initio routes are described: an $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$2-sum-rule extrapolation that reconstructs a Drude term in the head of $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$3, and a constant-dielectric approximation that simply replaces $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$4 (Leon et al., 2023).

The reported validation is unusually explicit. For simple metals such as Al and Na, using $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$5 poles and the constant-dielectric intraband correction on a $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$6 $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$7-mesh, MPA quasiparticle energies differ by $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$8 from a converged full-frequency contour-deformation approach. For states far from $$W_{GG'}(q,\omega)=v_{GG'}(q)+\sum_{p=1}^{N_p}\dfrac{B_p^{GG'}(q)}{\omega-\omega_p(q)}$$9, plasmon-pole errors are $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$0, while MPA reduces them to few meV. In Cu, single-pole PPA fails, with $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$1 errors in quasiparticle energies together with wrong $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$2 shape and spurious satellites; with $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$3 poles and quadratic grid $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$4 up to $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$5, MPA reproduces full-frequency $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$6 and $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$7 while using only $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$8 frequency points versus $B(\kappa)=H(\kappa)\circ \widehat B(\kappa)$9 in full frequency, corresponding to a $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$0 speed-up. Converged Cu quasiparticle energies agree with experiment to within $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$1. The physical interpretation is also made explicit: in Cu the fitted MPA poles of $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$2 coincide within $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$3 with peaks in the measured electron-energy-loss spectrum, and the four strongest loss peaks at $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$4, $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$5, $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$6, and $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$7 match the real parts of the first four MPA poles at $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$8, $H_{ij}(\kappa)=e^{\imath\kappa d_{ij}}$9, $$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$0, and $$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$1 (Leon et al., 2023).

3. Frequency extraction for Helmholtz BEM matrices

For the 3-D scalar Helmholtz equation, the discretized Galerkin BEM matrix $$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$2 is written entrywise as

$$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$3

where the far-field phase factor is

$$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$4

and $$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$5 is the de-oscillated remainder. The phase is exactly the phase of the Helmholtz Green’s function evaluated at the cluster-center distance. The paper states that rigorous estimates show that every far-field block of $$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$6 has numerical rank growing only like

$$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$7

when $$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$8 increases, with mesh size $$M^K(\omega)=\sum_{j=1}^{N_c}\alpha_j^K(\omega)\cdot E_{0j}(\omega)$$9 (Dirckx et al., 2020).

The construction proceeds in two layers. First, $$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$0 at any fixed $$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$1 is approximated by a standard $$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$2-matrix using geometric cluster trees, an admissibility criterion, and Adaptive Cross Approximation on admissible blocks,

$$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$3

Because $$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$4 is far less oscillatory, the numerical ranks $$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$5 remain essentially constant in $$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$6, so the same block skeletons can be reused for all $$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$7. Second, each admissible block is treated as a matrix-valued function of $$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$8 and approximated by a low-tensor-rank expansion

$$\mathbf E_{\mathrm{scat}}(\rho,\varphi;\omega)=\sum_m\big[-iZ_0B_m(\omega)\mathbf M_m^{(4)}+A_m(\omega)\mathbf N_m^{(4)}\big]$$9

where the $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$0 are rational functions without poles on the interval $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$1. The algorithm is an alternating ACA-AAA procedure: a pivot $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$2 and frequency $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$3 are chosen, the scalar function $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$4 is passed to AAA, a block sample at $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$5 is obtained by ACA, and orthogonalization plus a residual-based stopping criterion in Frobenius norm determine termination (Dirckx et al., 2020).

The computational consequences are reported quantitatively. Sampling a new frequency requires evaluating the scalar weights $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$6 in $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$7, forming $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$8 in $$\mathbf M^{(l)}_{\rm exact}=\dfrac{i}{\omega}\dfrac{(2l-1)!!}{(l-1)!}\int_{\mathbb R^3}\big[\mathbf J(\mathbf r)\otimes \mathbf r^{\otimes(l-1)}\big]\dfrac{j_{l-1}(kr)}{(kr)^{l-1}}\,d^3r$$9, reconstructing each far block in $\omega$0, and then multiplying entrywise by $\omega$1. Four 3-D meshes—sphere, Netgen crankshaft, Stanford bunny, and NASA almond—of size $\omega$2 were tested for $\omega$3. Over that range, the classical $\omega$4-matrix memory grows roughly $\omega$5, whereas the extracted $\omega$6-matrix stays essentially flat; build time is stated to be $\omega$7 faster at high $\omega$8. A compact rational model on $\omega$9, using $GW$00 Chebyshev samples, degree-$GW$01 rationals, and tolerance $GW$02, took about $GW$03 times the cost of one $GW$04-matrix build, but thereafter any new frequency could be recovered in a few seconds, often yielding a $GW$05 speedup compared to rebuilding the $GW$06-matrix. Asymptotic tests on refined spheres with constant $GW$07 confirm $GW$08 scaling of memory and time for the extracted method, with nearly constant ACA ranks (Dirckx et al., 2020).

4. Illumination-independent multipole spectra in nanophotonics

In the generalized polarizabilities formalism, the starting point is the exact size-dependent multipole expansion together with a generalized field propagator derived from the Lippmann-Schwinger integral equation. After discretization into $GW$09 cubical cells of volume $GW$10, the total polarization obeys

$GW$11

where $GW$12 and $GW$13 is the generalized propagator. Inserting the exact multipole kernels into the discretized multipole integrals yields generalized polarizability tensors $GW$14 such that, for each multipole $GW$15,

$GW$16

The electric dipole is split into a “primary” part $GW$17 and a toroidal part $GW$18, and the magnetic dipole $GW$19 is given in exact, non-long-wavelength form through spherical-Bessel-weighted volume integrals over the induced polarization $GW$20 (Majorel et al., 2022).

A central consequence of this construction is the definition of a multipole-mode density that does not depend on a chosen illumination: $GW$21 The paper identifies $GW$22 with the total density of mode $GW$23 at frequency $GW$24, and states that it is proportional to the maximum energy that can be radiated by mode $GW$25, irrespective of illumination. Once $GW$26 are known, the multipole moment for any illumination is obtained by an $GW$27 tensor-vector contraction. Optimal illumination for a given multipole under the constraint $GW$28 follows from the leading eigenvector of

$GW$29

with the largest eigenvalue giving the maximum achievable $GW$30 (Majorel et al., 2022).

The formalism is presented as a remedy for a specific limitation of the conventional multipole expansion: the latter gives access only to the modes excited by a specific illumination, and the study of various illuminations requires multiple, costly numerical simulations. After a single pre-computation of the generalized propagator, the generalized polarizabilities allow one to obtain the exact multipole decomposition for any illumination, the total density of available multipole modes, and the optimum illumination field distributions that maximally couple to specific multipole modes. The reported validation includes scattering spectra and radiation patterns of a $GW$31 cuboid, for which full-field direct simulation and the $GW$32-tensor reconstruction show perfect agreement up to quadrupolar order, as well as mode-density spectra $GW$33, $GW$34, $GW$35, and $GW$36 for a $GW$37 spheroid in air (Majorel et al., 2022).

5. Electromagnetic multipole decompositions in three and two dimensions

For arbitrarily shaped magneto-dielectric scatterers, the scattered electromagnetic field can be expanded in vector spherical harmonics with multipole amplitudes $GW$38 and $GW$39. The formulation explicitly includes both polarization currents $GW$40 and magnetization currents $GW$41, yielding exact volume-integral expressions

$GW$42

and from these one forms multipolar scattering efficiencies $GW$43 and $GW$44. In Faraday geometry, magnetic circular dichroism is decomposed into multipolar contributions by writing $GW$45, $GW$46, and $GW$47 as sums over $GW$48, with separate dielectric, magnetic, and interference terms. The paper states that multipole resonances associated with magnetization currents can be even stronger than multipole contributions from conventional dielectric currents, and that the analytical results were verified through comparison with finite-element simulations in COMSOL Multiphysics (Hernández-Sarria et al., 23 Oct 2025).

For two-dimensional photonics, the full-wave theory is based on an expansion of the outgoing scattered electric field in fourth-kind cylindrical vector wave functions,

$GW$49

The coefficients $GW$50 and $GW$51 are obtained either by contour integrals over an enclosing circle or by explicit source-current volume integrals after inserting an electric-field volume integral equation with equivalent electric and magnetic currents. Once the coefficients are known, the partial scattering cross section is

$GW$52

and the total two-dimensional scattering cross section is

$GW$53

The stated motivation is a “critical gap” in the literature, namely that existing multipole theories are primarily for three-dimensional particles of isotropic materials, whereas the new formulation addresses isolated, arbitrarily shaped, inhomogeneous, anisotropic cylinders or collections thereof (Loulas et al., 2024).

A third line of development replaces field amplitudes by exact current multipoles. The exact Cartesian current-multipole tensor of order $GW$54 is

$GW$55

which reduces to the familiar point-multipole moment in the long-wavelength limit. Closed-form linear combinations of the tensor elements reproduce the classical electric and magnetic multipole coefficients $GW$56 and $GW$57, allowing scattering and extinction cross sections to be decomposed by multipole order. The paper emphasizes that current multipoles include nonradiating current configurations absent in the classical field-based expansion; certain tensor combinations appear in $GW$58 or $GW$59, while others do not appear in any classical coefficient and are therefore strictly dark. This gives, in the authors’ formulation, a unique, unambiguous decomposition of $GW$60 at all frequencies (Shevchenko et al., 22 Aug 2025).

6. Interpretation, validation, and recurrent misconceptions

Several recurrent misconceptions are explicitly addressed in this literature. In metallic $GW$61, the complex poles of the multipole approximation are not presented as arbitrary fitting parameters: in Cu, comparison with electron-energy-loss spectra is used to argue that each MPA pole corresponds to a bona-fide collective excitation rather than a mere numerical artifact, with the imaginary parts reproducing finite widths without ad-hoc damping parameters (Leon et al., 2023).

In nano-optics, the conventional multipole expansion is not an illumination-independent modal inventory. The generalized polarizability formalism states that the conventional expansion only gives access to modes excited by a specific illumination, whereas the generalized tensors $GW$62 allow the total density of available multipole modes to be computed regardless of illumination and directly provide the optimum illumination for specific modes (Majorel et al., 2022).

In current-based electromagnetic theory, field multipoles are not a complete description of internal current structure. The exact current-multipole expansion states that dark or nonradiating modes, including the anapole as a special combination of dipole and octupole currents, are absent from the classical field-based expansion and can store energy while influencing the frequency dependence and linewidth of radiative multipole resonances (Shevchenko et al., 22 Aug 2025).

In boundary-element frequency extraction, the method is not merely a restatement of low-rank compression for oscillatory matrices. The essential step is the extraction of the known phase of the Green’s function, after which the de-oscillated matrix $GW$63 is approximated by a frequency-independent $GW$64-matrix skeleton combined with adaptive rational approximation in the continuous frequency dimension (Dirckx et al., 2020).

Taken together, these works indicate a shared computational strategy. This suggests that, despite the different operators and physical observables involved, multipole-frequency decomposition is most useful when it relocates the hard part of frequency dependence into a representation with low pole count, low tensor rank, or directly interpretable multipolar amplitudes, while preserving the quantities of interest—quasiparticle energies, spectral functions, BEM matrices, scattering spectra, mode densities, or dichroic observables.