Cartesian and spherical multipole expansions in anisotropic media

The multipole expansion method has long been an essential mathematical tool successfully applied to a wide range of systems described by partial differential equations involving the Laplace operator. Commonly used in electromagnetism^1,2, this method also sees applications from the quantum level³ up to Newtonian gravitational interaction^4,5 and gravitational waves^6,7; indeed, whenever one needs to describe how fields behave outside the source region. However, the multipole expansion in anisotropic media, where wave propagation becomes direction-dependent, has remained relatively ill-studied. This must be distinguished from the multipole expansion associated with an anisotropic source, which results in direction-dependent radiation or scattering⁸. Here, we assume that wave propagation occurs within an anisotropic medium, while the source itself may also exhibit anisotropic characteristics.

For optical or radio-frequency systems, the rise of advanced wave manipulation, such as via metamaterials^9,10,11 or nanophotonics¹², warrants a study of this method in anisotropic media. For instance, meta-atoms mixing electric and magnetic multipole moments give rise to bianisotropy¹³. The multipole expansion is also readily applied to metasurfaces to express the fields, which are then linked through the T-matrix method¹⁴. For example, this approach has been applied to the analysis of substrate-supported dielectric nanoresonator metasurfaces¹⁵. The multipole expansion is also used to obtain a basis of resonant states in the resonant-state expansion method for the analysis of open optical systems¹⁶. High-order multipoles are also gaining interest in connection to the chosen center of coordinates¹⁷. A variation of the multipole expansion can also be used to describe the nonlocal response of certain metasurfaces¹⁸. While the literature commonly uses the multipole expansion to describe and analyze metamaterials with anisotropic wave propagation, the analysis of the multipole expansion of fields embedded in an anisotropic medium is missing. Implications of this analysis could include the effect of a layered substrate with uniaxial wave propagation on the interaction between neighboring particles in a metasurface.

The analogous nature of electromagnetic and sound waves means that anisotropic wave propagation can also occur within acoustic metamaterials¹⁹. In plasma physics as well, the multipole expansion is used to describe the interaction potentials of dust particles²⁰, though the method has yet to see application to anisotropic cases beyond the first order²¹. On the largest scales, anisotropic baryon acoustic oscillations imprint on the clustering of galaxies and can be observed in their power spectra²², while the cosmic neutrino background applies an anisotropic stress that dampens gravitational waves, a potentially observable effect in the multipole expansion of the cosmic microwave background²³.

In the above-mentioned studies, solutions of wave-like equations featuring the Laplace operator are found using the separation of variables technique and representing the field as a linear combination of harmonic polynomials (for the angular part) and spherical Hankel functions of the first kind, second kind, or a combination of both (for the radial part)²⁴, yielding the spherical multipole expansion. Equivalently, the same solutions can be found by using a Taylor series expansion of Green’s function (see, e.g., ref. ²⁵), yielding the Cartesian multipole expansion. The number of Cartesian basis functions (and moments) grows cubically with the series truncation order, while the same growth is only quadratic for the spherical expansion. The equivalence of both for electrostatics is the subject of Problem 4.3 in ref. ²⁴. Indeed, the solution of this exercise shows that any Cartesian moments tensor of order ℓ can be represented by a linear combination of spherical moments tensors of order ℓ, ℓ − 2, etc., down to 0 or 1. For plasma physics and the Vlasov–Fokker–Planck equation, a detailed computation of the conversion between both approaches is presented in ref. ²⁶; the references therein, e.g., ref. ²⁷, are also of interest. Contrary to ref. ²⁶, we are interested in the singular behavior of multipole expansions, i.e., their tendency to explode at the origin. This allows us to consider all Cartesian moments, including, for example, the trace of the second-order moments tensor, which is set to zero by design in other references. While this is necessary for isotropic media to satisfy the wave propagation spherical symmetry, we need to relax this assumption in anisotropic media to obtain a full representation of fields in a general setting. Note that the presented limitation of the spherical formulation is independent of the (non)existence of analytical closed-form expressions of Green’s functions; rather, it is a consequence of the loss of symmetry in the wave propagation. In isotropic media, the spherical symmetry of the Laplace operator implies that irreducible tensors are sufficient to describe scattered or radiated fields. As we will see, in anisotropic media, all tensor elements are necessary to describe the fields fully. As such, we show that the Cartesian multipole expansion is central to obtaining field solutions in anisotropic media.

In this paper, we describe the link between the spherical and Cartesian multipole expansions in isotropic media and demonstrate this equivalence breaks down in anisotropic media. Next, we also discuss the link to the vector spherical harmonics: the extension of their scalar counterpart to vector fields. We then illustrate the lack of equivalence in a radiation inverse problem, including a detailed analysis of a large-anisotropy regime. We show that, in anisotropic media, the use of the Cartesian multipole expansion improves the modeling efficiency and interpretability.

Spherical and Cartesian multipole expansions in isotropic media

We begin by finding equivalence relations between the spherical and Cartesian multipole expansions by analyzing their singular behavior. This will allow us to extrapolate the approach to anisotropic media by considering the well-defined radiation of point sources. To this end, previous work on the pseudopotential method^28,29,30 has sparked interest in the singular behavior of spherical wave functions. This behavior is challenging to describe as the concept of “radial” derivatives of the Dirac δ function is not trivially defined (see refs. ^31,32,33) because of the singular nature of spherical coordinates at the origin. A rigorous approach in ref. ³⁴ allows to write, for the spherical basis function f_ℓm(r) introduced in Equation (8),

where ∇ is the gradient operator, k₀ is the free-space wavenumber, δ(r) is the Dirac δ-function (a distribution), and the solid harmonics ({Y}_{ell m}({{{bf{k}}}})={k}^{ell }{{{{mathcal{Y}}}}}_{ell m}(theta ,phi )) are the homogeneous and harmonic polynomials obtained from the spherical harmonics ({{{{mathcal{Y}}}}}_{ell m}(theta ,phi )). (k, θ, ϕ) are the spherical coordinates of the three-dimensional vector k. We use the convention

where P_ℓm(z) is the associated Legendre polynomial. Y_ℓm(∇) is a partial differential operator obtained by replacing the Cartesian components of k = (k_x, k_y, k_z) in Y_ℓm(k) by

and likewise for ({tilde{Y}}_{ell m}(nabla )). For instance, ({Y}_{1,1}({{{bf{k}}}})=-1/2sqrt{3/(2pi )}({k}_{x}+{{{rm{i}}}}{k}_{y})), hence ({Y}_{1,1}(nabla )delta ({{{bf{r}}}})=-1/2sqrt{3/(2pi )}(partial delta ({{{bf{r}}}})/partial x+{{{rm{i}}}}partial delta ({{{bf{r}}}})/partial y)).

To obtain the singular behavior in Equation (1), we start by noticing that the function r⁻¹ yields a Dirac singularity under the Laplace operator:

vp denotes the principal value. This generalizes to higher orders as^{34,Lemma 1.(c)]}

where (2ℓ − 1)!! = 1 ⋅ 3 ⋅ 5 ⋯ (2ℓ − 1). Moving to waves and the Helmholtz operator, the inverse radius is replaced with the spherical Neumann function y_ℓ(k₀r) which also behaves in ({({k}_{0}r)}^{-ell -1}) as r → 0. A consequence of the same Lemma is that

Radiating waves are described by spherical Hankel functions of the first or second kind ({h}_{ell }^{(1,2)}({k}_{0}r))^{35,Eq. 10.1.1]:}

where j_ℓ(k₀r) is a spherical Bessel function of the first kind. As this function is regular at the origin, the complete radiating wave is

and its singular behavior is given by Equation (1).

On the other hand, the singular behavior of the Cartesian basis functions is simply given, by definition, by (see, e.g., refs. ^36,37)

where we have used the multi-index notation, i.e., (alpha =({alpha }_{x},{alpha }_{y},{alpha }_{z})in {{mathbb{N}}}^{3}), and ({nabla }^{alpha }={(partial /partial x)}^{{alpha }_{x}}{(partial /partial y)}^{{alpha }_{y}}{(partial /partial z)}^{{alpha }_{z}}). The Cartesian solutions f_α(r) can be obtained directly from Green’s function G(r), which satisfies by definition

Indeed, since

we have

The next step is to assume that under appropriate boundary conditions (e.g., radiation conditions for the spherical Hankel function of the first or second kind), there are unique solutions f_ℓm(r) and f_α(r) to the wave equation with sources, Eqs. (1) and (9). Conversely, the Helmholtz operator applied to the solutions f_ℓm(r) and f_α(r) yields singularities that are entirely determined by the polynomials (tilde{Y}_{ell m}({{{bf{k}}}})) and k^α. Hence, from now on, we identify spherical solutions f_ℓm(r) by their corresponding harmonic polynomial ({tilde{Y}}_{ell m}({{{bf{k}}}})) and the Cartesian solutions f_α(r) by the monomials k^α:

With this formalism in mind, we now strive to write equivalences between the spherical and Cartesian solutions by finding equivalences between their corresponding polynomials.

Spherical to Cartesian multipole expansions

Clearly, any harmonic polynomial can be written as a linear combination of monomials; hence, any spherical solution can be represented as a linear combination of Cartesian solutions. For example, since

we have

Cartesian to spherical multipole expansions

Nevertheless, it is generally not possible to write an arbitrary monomial k^α as a linear combination of harmonic polynomials ({tilde{Y}}_{ell m}({{{bf{k}}}})), as the monomial is not necessarily harmonic. For any monomial, we can nonetheless write by Theorem 5.21 in ref. ³⁸:

where the harmonic polynomials p_{∣α∣−2n}(k) are given by

∣α∣ = α_x + α_y + α_z, ⌊ ⋅ ⌋ denotes the floor function, ({{nabla }_{{{{bf{k}}}}}}^{2}) is the Laplace operator with respect to k, and c_j,n are coefficients specified in the cited theorem, which are omitted here for brevity.

Next, we project the harmonic polynomials p_{∣α∣−2n}(k) onto the solid harmonics Y_{∣α∣−2n,m}(r)

by using the orthonormality of the spherical harmonics ({{{{mathcal{Y}}}}}_{ell ,m}(theta ,phi )):

where ⋅ ^* denotes complex conjugation, and p_{∣α∣−2n}(k)/k^{∣α∣−2n} is the angular part of the polynomial p_{∣α∣−2n}(k) (a function of θ and ϕ). Combining Eqs. (17) and (19), we find

Applying the scaling defined in Eq. (1), we have

where

Now, by using the equivalence between the polynomials and their corresponding multipole expansion (see Eqs. (13) and (14)), we can write any Cartesian solution as

Next, for wave-like equations, like the Helmholtz equation, the Laplacian ∇² of solutions can be re-written from Eq. (8):

Ignoring the origin, we have

Inserting into Eq. (24), we obtain

outside the origin. In other words, any Cartesian solution of order ∣α∣ = ℓ can be written as a linear combination of spherical solutions of order ℓ, ℓ − 2, etc., down to ℓ = 0 or ℓ = 1. For example, the non-harmonic octupole α = (0, 0, 2) (identified by the monomial ({{{{bf{k}}}}}^{alpha }={k}_{z}^{2})) is projected onto solid spherical harmonics of orders zero and two as

which corresponds to

By Eq. (26), multiplying the polynomial ({tilde{Y}}_{0,0}({{{bf{k}}}})) by k² is equivalent to multiplying the solution f_0,0(r) by (-{k}_{0}^{2}), and we obtain

outside the origin. At the origin, derivatives of the Dirac δ function must also be added in compliance with Eq. (25).

The case of anisotropic media

In anisotropic media, the Laplace operator (as in Eq. (25)) is replaced with weighted second-order derivatives, reflecting different propagation characteristics along different directions. For example, electromagnetic fields obey a modified vector wave equation^{39,Equation (5.6)]} whose anisotropy is determined by the constitutive parameters’ tensorial nature. In such cases, as we illustrate in the section “Illustration on an inverse problem,” some sources radiate fields that cannot be fully represented by a spherical multipole expansion. In turn, the more numerous Cartesian moments are necessary to model these fields accurately and efficiently.

Relation to the vector spherical harmonics

Many applications such as optics require the use of vector spherical harmonics. Here, we wish to outline how the link between the spherical and Cartesian approaches can be extended to the vector spherical harmonics and, thus, the vector wave equation.

From vector spherical harmonics to the Cartesian multipole expansion

Fortunately, it is possible to write the components of vector spherical harmonics in terms of scalar spherical harmonics. This allows to generalize the reasoning applied to the scalar wave equation. Using quantum theory formalism and the Condon–Shortley phase, the vector spherical harmonics can be written as^{40,Eq. 7.3.3]}

where j = ℓ, ℓ ± 1 is the total angular momentum quantum number, ℓ the orbital angular momentum quantum number, m = − j, …, j, the sum runs over ({m}^{{prime} }=-ell ,ldots ell), q = ±1, 0, ({C}_{ell m1q}^{jm}) denotes the Clebsch-Gordan coefficient, and the e_q form the covariant spherical basis defined in Cartesian coordinates as

The vector spherical harmonics used by Jackson²⁴({{{{bf{{{{mathcal{Y}}}}}}}}}_{ell m}^{(q)}(theta ,phi )), q = 0, ±1, are obtained as linear combinations of ({{{{bf{{{{mathcal{Y}}}}}}}}}_{jm}^{ell }(theta ,phi )) (see supplementary Eqs. (S1–S3)). They correspond to the azimuthal (q = 0), polar (q = 1), and radial (q = − 1) components. In the context of Mie scattering, ({{{{bf{{{{mathcal{Y}}}}}}}}}_{ell m}^{(1)}(theta ,phi )) represents electric multipoles, whereas ({{{{bf{{{{mathcal{Y}}}}}}}}}_{ell m}^{(0)}(theta ,phi )) corresponds to magnetic multipoles. In homogeneous media, the radial component is not allowed for electromagnetic fields. Together with Eq. (31), we show in the Supplementary Note 1 how to express the azimuthal or polar vector spherical harmonics using the Cartesian multipole expansion of an appropriate current source. We show that a current density given by

radiates an electric field whose angular dependence is given by the azimuthal vector spherical harmonic ({{{{bf{{{{mathcal{Y}}}}}}}}}_{ell m}^{(0)}(theta ,phi )). The solid vector harmonics correspond to ({{{{bf{Y}}}}}_{ell m}^{ell }({{{bf{k}}}})={k}^{ell }{{{{bf{{{{mathcal{Y}}}}}}}}}_{ell m}^{ell }({{{bf{k}}}})). On the other hand, the polar vector spherical harmonic is radiated by

where ({{{{bf{Y}}}}}_{ell pm 1,m}^{ell }({{{bf{k}}}})={k}^{ell }{{{{bf{{{{mathcal{Y}}}}}}}}}_{ell pm 1,m}^{ell }({{{bf{k}}}})). The Cartesian components of the current densities ({{{{bf{J}}}}}_{ell m}^{(0,1)}({{{bf{r}}}})) correspond to singularities of the scalar wave equation analyzed in the “Spherical and Cartesian multipole expansions in isotropic media” subsection of the “Methods” section. Hence, any component of the vector spherical harmonics can be represented by an appropriate sum of Cartesian multipoles.

From the Cartesian multipole expansion to the vector spherical harmonics

As expected from the findings for the scalar wave equation, the inverse transformation—i.e., writing an arbitrary vector Cartesian multipole expansion as a linear combination of vector spherical harmonics—is not as straightforward, and a complete treatment goes beyond the scope of this research. To understand the issue at hand, one approach is to count the number of vector spherical harmonics at a given order and to compare with the corresponding number of Cartesian multipole moments. As we will see, as in the scalar case, there is an excess of Cartesian multipole moments.

Let us now count the vector spherical harmonics. The azimuthal vector harmonics correspond to a current density given by the vector spherical harmonics ({{{{bf{{{{mathcal{Y}}}}}}}}}_{ell m}^{ell }(theta ,phi )). In light of Eq. (31), there are 2ℓ + 1 different azimuthal waves of order ℓ > 0 (none for ℓ = 0). By the same argument, there are 2(ℓ ± 1) + 1 polar vector harmonics of order ℓ > 0 corresponding to j = ℓ ± 1. But the polar harmonics are degenerate, in that both values of j yield the same waves, by virtue of the contraction in Eq. (26) (see the Supplementary Notes 1.2.3). Thus, in total, there are

vector spherical harmonics of order ℓ > 0. For ℓ = 0, only the polar solutions exist, and their count is given by 2 ⋅ (0 + 1) + 1 = 3.

Next, we count the moments for the vector Cartesian multipole expansion. Since there are three possible polarizations (x, y, and z) and 2ℓ + 1 spherical harmonics of order ℓ per polarization, one would expect

different vector spherical harmonics (VSH). The missing harmonics are the degenerate polar harmonics with j = ℓ − 1, and their count is given by 2(ℓ − 1) + 1:

This hints towards a similar issue to the scalar case when going to anisotropic media, as the Laplace operator is responsible for the degeneracy. Another viewpoint is that in anisotropic media, the degeneracy disappears, and new degrees of freedom in wave polarization become available.

To note, we circumvent this issue in the section “Illustration on an inverse problem” by opting for the electric Green’s tensor-valued function instead of vector spherical harmonics. This approach enables a direct connection between the electric current and the vector-valued field.

Here, we show that the loss of equivalence between spherical and Cartesian multipole expansions has measurable and significant effects on an electromagnetic radiation inverse problem. We consider an electric uniaxial medium characterized by the relative permittivity tensor

On top of appearing in natural materials (e.g., calcite⁴¹), electromagnetic anisotropic media commonly describe metamaterials (e.g., through effective medium theory⁴²), which are seeing an increase of adoption from optical⁴³ to radio frequencies⁴⁴. High anisotropy can be achieved through 3D printing of carbon fiber composite filament⁴⁵.

To illustrate, we consider a radiation inverse problem—i.e., we try to predict simulation results (the ground-truth data) by performing a multipole expansion in spherical and Cartesian coordinates. For ε_zz≤ 10, the ground truth data consists of the electric field on a sphere surrounding the source

at a distance of one vacuum wavelength λ₀ simulated using Meep⁴⁶ in the frequency domain. The simulation domain consists of a cube of side 2λ₀ surrounded by a perfectly matched layer of thickness λ₀/5. The source consists, first, of a pure dipole (α = (0, 0, 0)) and then of a non-harmonic octupole (α = (0, 0, 2)). The latter is approximated by four dipole sources of alternating signs separated by a simulation voxel along the z axis.

To avoid the “inverse crime,”⁴⁷ the fitting data (i.e., the prediction through a multipole expansion) is obtained from the analytical Green’s function⁴⁸. The Cartesian multipole expansion basis functions are obtained by differentiating Green’s function up to a fixed order ({ell }_{max }), i.e.,

where G(r) is Green’s tensor-valued function. G is defined as the solution of

where I is the identity matrix. In turn,

and μ is the medium’s magnetic permeability. We use the following convention for the Fourier transform f(ω) of a function f(t):

For ε_zz > 10, the ground truth data is also obtained from the analytical Green’s function. Doing so, while the “inverse crime” is committed, the data range can be extended to high degrees of anisotropy (up to ε_zz = 10³) while avoiding prohibitively large numerical problems.

In parallel, the spherical multipole expansion basis functions are derived as linear combinations of the Cartesian basis functions, as in Eq. (16).

The inverse problem consists in finding the spherical (resp. Cartesian) moments C_ℓm (resp. C_α) that minimize the squared L²-norm of the error between the ground truth and the multipole expansion, normalized to the ground truth squared L²-norm (energy). The minimization is carried out using the Nelder-Mead solver. After convergence, we evaluate the proportion of the ground truth energy explained by the multipole expansion. A schematic illustration of the problem is shown in Fig. 1.

The results in Fig. 2 for an anisotropic medium with ε_zz = 4 show that the dipole radiation (i.e., α = (0, 0, 0)) is perfectly explained by both the spherical and Cartesian multipole expansions. This is expected, as the corresponding solutions f_0,0 and f_(0, 0, 0) are proportional to each other and the monomial given by a constant is harmonic. However, the spherical multipole expansion fails to explain the non-harmonic octupole (α = (0, 0, 2)) radiation, even with increased truncation orders ({ell }_{max }). Indeed, Fig. 3 shows that, while the original Cartesian multipole moment is accurately recovered by the inverse problem, the spherical multipole expansion needs four dominating moments ((ℓ, m) = (0, 0), (2, − 2), (2, 0), (2, 2)) to explain 75.4% of the ground truth energy. Moreover, Fig. 4 presents the effect of the degree of anisotropy on the representation power of the spherical multipole expansion. This figure shows that this lack of representation power worsens with increasing anisotropy up to ε_zz = 8 for ({ell }_{max }=2). Indeed, in this regime, we observe the superposition of ordinary and extraordinary waves⁴⁹ given by the dispersion relations

where k is the wavenumber vector, k^⊤ its transpose, and k₀ the free-space wavenumber. The second relation does not feature the Laplacian  − k^⊤k, but an anisotropic propagation (-{{{{bf{k}}}}}^{top }overline{overline{varepsilon }}{{{bf{k}}}}) incompatible with the contraction in Equation (24).

However, surprisingly, the representation power of the spherical multipole expansion for the non-harmonic octupole α = (0, 0, 2) increases with increasing anisotropy ε_zz above ε_zz = 8. Indeed, Fig. 4 shows asymptotes to a perfect representation for all orders ({ell }_{max }) as ε_zz grows to 10³. The spherical expansion involves the moments ℓ = 0, 2, m = 0 to represent the fields in this large-anisotropy regime.

To understand this regime, we start from the anisotropic vector wave equation for Green’s function G(k) in wavenumber space, which reads (omitting the frequency dependence)

with I the identity matrix. Solving for Green’s function and taking the limit as ε_zz → ∞, we obtain

where ({{{{bf{I}}}}}_{perp }={{{{bf{e}}}}}_{x}{{{{bf{e}}}}}_{x}^{top }+{{{{bf{e}}}}}_{y}{{{{bf{e}}}}}_{y}^{top }) and k_⊥ = k_xe_x + k_ye_y. As the propagation dynamics are determined by the denominator, we apply the contraction in Eq. (25) to g(k), obtaining

Multiplying by the right-hand side denominator and transforming back to spatial variables, the remainder function h must satisfy

Solving this ordinary differential equation, we obtain (see Supplementary Note 2)

where c_± are constants. Thus, the contraction reads, by combining Eqs. (48) and (49) in the spatial domain,

where ({nabla }_{perp }={{{{bf{e}}}}}_{x}frac{partial }{partial x}+{{{{bf{e}}}}}_{y}frac{partial }{partial y}). In other words, in the large anisotropy regime, we observe the same structure as in isotropic media. However, the singular region is no longer at the origin but is now defined by the function h. Arguably, the non-decaying traveling plane waves in h can be discarded by setting c_± = 0 to satisfy a radiation-like condition. In that case, the remaining support of h shrinks to the z axis because of the Dirac function δ(x, y). In turn, we find the anisotropic contraction outside this axis:

which explains the representation power of the spherical multipole expansion in the large anisotropic regime.

Finally, one might argue that it is always possible to represent fields on a sphere by a projection onto spherical harmonics. While this is mathematically true, it is usually not helpful in physical problems. Indeed, the moments recovered by the projection do not correspond to those of the source responsible for the radiation. To illustrate, we project the xx component of Green’s function on the sphere described at the beginning of this section in isotropic (ε_zz = 1) and anisotropic (ε_zz = 4) media. Figure 5 shows the change in explained ground truth energy with increasing truncation orders ({ell }_{max }). As expected, in isotropic media, a finite number of moments (up to second order) are needed to represent the field fully. However, moments of order up to eight are necessary in the considered anisotropic medium, even though the source is the same—only the field propagation changed.

It would also be possible to repeat this comparison by projecting the fields onto the restriction on the sphere of an arbitrary basis of polynomials of three variables. Since the spherical harmonics are formed by a restriction on the unit sphere of harmonic and homogeneous polynomials of three variables, this would consist of a Cartesian counterpart to the spherical harmonics. However, such a Cartesian projection would still be blind to the propagation anisotropy. Moreover, the added degrees of freedom would also be redundant, as it has been long shown that the spherical harmonics span all square-integrable functions on the sphere³⁸.

By analyzing the singular behavior of spherical and Cartesian multipole expansions, we applied harmonic function theory to find explicit links between both approaches in isotropic media while uncovering a lack of equivalence in anisotropic media. In those, the Cartesian approach is necessary, as illustrated in a radiation problem for electromagnetic waves. We showed that a spherical multipole expansion cannot represent the fields radiated by some sources, while a Cartesian multipole expansion can. Although a radical order increase improves the representation of the spherical expansion, the Cartesian expansion is still more efficient outside the large anisotropy regime.

Furthermore, by approaching the equivalence problem from the singular behavior of the solutions, we have shown how the Cartesian multipole expansion—obtained straightforwardly from Green’s function—can be used to obtain the spherical one, which is not readily available in anisotropic media. This matters for cases where the dominant source feature is harmonic and, therefore, well-represented by the spherical approach. Also, the Laplacian-involving mapping from Cartesian to spherical approaches presented in Eq. (24) can be applied to many wave-like systems where the multipole expansion is already applied. Anytime the propagation is anisotropic, the Cartesian approach is expected to offer a better representation power and, ultimately, an improved understanding of the underlying physics.

Note that the discrepancy between both approaches only appears for sources of order above one; this, together with the relatively recent rise of manufactured anisotropic media such as metamaterials, may explain the lack of prior research. In addition to the electromagnetic application to optical and radio-frequency systems, our results are relevant to describing anisotropic interactions in plasma physics as discussed in the introduction. Similarly, the implications for cosmological studies on anisotropic baryon acoustic oscillations and the cosmic neutrino background also highlight key examples where our findings can have an impact.