Cartesian and spherical multipole expansions in anisotropic media

Introduction

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 electromagnetism1,2, this method also sees applications from the quantum level3 up to Newtonian gravitational interaction4,5 and gravitational waves6,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 scattering8. 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 metamaterials9,10,11 or nanophotonics12, warrants a study of this method in anisotropic media. For instance, meta-atoms mixing electric and magnetic multipole moments give rise to bianisotropy13. The multipole expansion is also readily applied to metasurfaces to express the fields, which are then linked through the T-matrix method14. For example, this approach has been applied to the analysis of substrate-supported dielectric nanoresonator metasurfaces15. 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 systems16. High-order multipoles are also gaining interest in connection to the chosen center of coordinates17. A variation of the multipole expansion can also be used to describe the nonlocal response of certain metasurfaces18. 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 metamaterials19. In plasma physics as well, the multipole expansion is used to describe the interaction potentials of dust particles20, though the method has yet to see application to anisotropic cases beyond the first order21. On the largest scales, anisotropic baryon acoustic oscillations imprint on the clustering of galaxies and can be observed in their power spectra22, 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 background23.

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)24, 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. 25), 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. 24. 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. 26; the references therein, e.g., ref. 27, are also of interest. Contrary to ref. 26, 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.

Methods

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 method28,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. 34 allows to write, for the spherical basis function fm(r) introduced in Equation (8),

$$({nabla }^{2}+{k}_{0}^{2}){f}_{ell m}({{{bf{r}}}})=frac{{(-1)}^{ell }4pi {{{rm{i}}}}}{{k}_{0}^{ell +1}}{Y}_{ell m}(nabla )delta ({{{bf{r}}}})=:{tilde{Y}}_{ell m}(nabla )delta ({{{bf{r}}}})$$
(1)

where is the gradient operator, k0 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

$${{{{mathcal{Y}}}}}_{ell m}(theta ,phi )=sqrt{frac{2ell +1}{4pi }}sqrt{frac{(l-m)!}{(l+m)!}}{P}_{ell m}left[cos (theta )right]{{{{rm{e}}}}}^{{{{rm{i}}}}mphi }$$
(2)

where Pm(z) is the associated Legendre polynomial. Ym() is a partial differential operator obtained by replacing the Cartesian components of k = (kxkykz) in Ym(k) by

$${k}_{x}to partial /partial x,,{k}_{y}to partial /partial y,,{k}_{z}to partial /partial z$$
(3)

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−1 yields a Dirac singularity under the Laplace operator:

$${nabla }^{2}{{{rm{vp}}}}({r}^{-1})=-4pi delta ({{{bf{r}}}})$$
(4)

vp denotes the principal value. This generalizes to higher orders as34,Lemma 1.(c)]

$${nabla }^{2}{{{rm{vp}}}}left[{r}^{-ell -1}{{{{mathcal{Y}}}}}_{ell m}(theta ,phi )right]=frac{{(-1)}^{ell -1}4pi }{(2ell -1)!!}{Y}_{ell m}(nabla )delta ({{{bf{r}}}})$$
(5)

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(k0r) which also behaves in ({({k}_{0}r)}^{-ell -1}) as r → 0. A consequence of the same Lemma is that

$$({nabla }^{2}+{k}_{0}^{2}){{{rm{vp}}}}left[{y}_{ell }({k}_{0}r){{{{mathcal{Y}}}}}_{ell m}(theta ,phi )right]=frac{{(-1)}^{ell }4pi }{{k}_{0}^{ell +1}}{Y}_{ell m}(nabla )delta ({{{bf{r}}}})$$
(6)

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]:

$${h}_{ell }^{(1)}({k}_{0}r)={j}_{ell }({k}_{0}r)+{{{rm{i}}}}{y}_{ell }({k}_{0}r)$$
(7)

where j(k0r) is a spherical Bessel function of the first kind. As this function is regular at the origin, the complete radiating wave is

$${f}_{ell m}({{{bf{r}}}}):={{{rm{vp}}}}left[{h}_{ell }^{(1)}({k}_{0}r){{{{mathcal{Y}}}}}_{ell m}(theta ,phi )right]$$
(8)

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)

$$({nabla }^{2}+{k}_{0}^{2}){f}_{alpha }({{{bf{r}}}})={nabla }^{alpha }delta ({{{bf{r}}}})$$
(9)

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

$$({nabla }^{2}+{k}_{0}^{2})G({{{bf{r}}}})=delta ({{{bf{r}}}})$$
(10)

Indeed, since

$${nabla }^{alpha }left[({nabla }^{2}+{k}_{0}^{2})G({{{bf{r}}}})right]=({nabla }^{2}+{k}_{0}^{2}){nabla }^{alpha }G({{{bf{r}}}})={nabla }^{alpha }delta ({{{bf{r}}}})$$
(11)

we have

$${f}_{alpha }({{{bf{r}}}})={nabla }^{alpha }G({{{bf{r}}}}).$$
(12)

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 fm(r) and fα(r) to the wave equation with sources, Eqs. (1) and (9). Conversely, the Helmholtz operator applied to the solutions fm(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 fm(r) by their corresponding harmonic polynomial ({tilde{Y}}_{ell m}({{{bf{k}}}})) and the Cartesian solutions fα(r) by the monomials kα:

$${f}_{ell m}({{{bf{r}}}})={{{rm{vp}}}}left[{h}_{ell }^{(1)}({k}_{0}r){{{{mathcal{Y}}}}}_{ell m}(theta ,phi )right] iff tilde{Y}_{ell m}({{{bf{k}}}})$$
(13)
$${f}_{alpha }({{{bf{r}}}})={nabla }^{alpha }G({{{bf{r}}}}) iff {{{{bf{k}}}}}^{alpha }$$
(14)

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

$${tilde{Y}}_{2,0}({{{bf{k}}}})={{{rm{i}}}}sqrt{5pi }/{k}_{0}^{3}(2{k}_{z}^{2}-{k}_{x}^{2}-{k}_{y}^{2})$$
(15)

we have

$${f}_{2,0}({{{bf{r}}}})={{{rm{i}}}}sqrt{5pi }/{k}_{0}^{3}left[2{f}_{(0,0,2)}({{{bf{r}}}})-{f}_{(2,0,0)}({{{bf{r}}}})-{f}_{(0,2,0)}({{{bf{r}}}})right].$$
(16)

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. 38:

$${{{{bf{k}}}}}^{alpha }={sum }_{n=0}^{lfloor | alpha | /2rfloor }{({k}^{2})}^{n}{p}_{| alpha | -2n}({{{bf{k}}}})$$
(17)

where the harmonic polynomials pα−2n(k) are given by

$${p}_{| alpha | -2n}({{{bf{k}}}})={sum }_{j=n}^{lfloor | alpha | /2rfloor }{c}_{j,n}{k}^{2(j-n)}{({{nabla }_{{{{bf{k}}}}}}^{2})}^{j}{{{{bf{k}}}}}^{alpha }$$
(18)

α = αx + αy + αz, denotes the floor function, ({{nabla }_{{{{bf{k}}}}}}^{2}) is the Laplace operator with respect to k, and cj,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)

$${p}_{| alpha | -2n}({{{bf{k}}}})={sum }_{m=-| alpha | +2n}^{| alpha | -2n}{A}_{| alpha | -2n,m}{Y}_{| alpha | -2n,m}({{{bf{k}}}})$$
(19)

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

$${A}_{| alpha | -2n,m}=int_{0}^{pi }{{{rm{d}}}}theta int_{0}^{2pi }{{{rm{d}}}}phi sin (theta ){p}_{| alpha | -2n}({{{bf{k}}}})/{k}^{| alpha | -2n}{{{{mathcal{Y}}}}}_{| alpha | -2n,m}^{* }(theta ,phi )$$
(20)

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

$${{{{bf{k}}}}}^{alpha }={sum }_{n=0}^{lfloor | alpha | /2rfloor }{({k}^{2})}^{n}{sum }_{m=-| alpha | +2n}^{| alpha | -2n}{A}_{| alpha | -2n,m}{Y}_{| alpha | -2n,m}({{{bf{k}}}})$$
(21)

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

$${{{{bf{k}}}}}^{alpha }={sum }_{n=0}^{lfloor | alpha | /2rfloor }{({k}^{2})}^{n}{sum }_{m=-| alpha | +2n}^{| alpha | -2n}{tilde{A}}_{| alpha | -2n,m}{tilde{Y}}_{| alpha | -2n,m}({{{bf{k}}}})$$
(22)

where

$${tilde{A}}_{| alpha | -2n,m}=frac{{k}_{0}^{|alpha|-2n +1}}{{(-1)}^{|alpha|-2n }4pi {{{rm{i}}}}}{A}_{| alpha | -2n,m}$$
(23)

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

$${f}_{alpha }({{{bf{r}}}})={sum }_{n=0}^{lfloor | alpha | /2rfloor }{({nabla }^{2})}^{n}{sum }_{m=-| alpha | +2n}^{| alpha | -2n}{tilde{A}}_{| alpha | -2n,m}{f}_{| alpha | -2n,m}({{{bf{r}}}})$$
(24)

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

$${nabla }^{2}{f}_{ell m}({{{bf{r}}}})=-{k}_{0}^{2}{f}_{ell m}({{{bf{r}}}})+{tilde{Y}}_{ell m}(nabla )delta ({{{bf{r}}}})$$
(25)

Ignoring the origin, we have

$${nabla }^{2}{f}_{ell m}({{{bf{r}}}})=-{k}_{0}^{2}{f}_{ell m}({{{bf{r}}}})$$
(26)

Inserting into Eq. (24), we obtain

$${f}_{alpha }({{{bf{r}}}})={sum }_{n=0}^{lfloor | alpha | /2rfloor }{sum }_{m=-| alpha | +2n}^{| alpha | -2n}{tilde{A}}_{| alpha | -2n,m}{(-{k}_{0}^{2})}^{n}{f}_{| alpha | -2n,m}({{{bf{r}}}})$$
(27)

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

$${k}_{z}^{2}=2/3left[sqrt{pi }{k}^{2}{Y}_{0,0}({{{bf{k}}}})+2sqrt{pi /5}{Y}_{2,0}({{{bf{k}}}})right]$$
(28)

which corresponds to

$${k}_{z}^{2}=-{{{rm{i}}}}{k}_{0}/6left[1/sqrt{pi }{k}^{2}{tilde{Y}}_{0,0}({{{bf{k}}}})+2{k}_{0}^{2}/sqrt{5pi }{tilde{Y}}_{2,0}({{{bf{k}}}})right]$$
(29)

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

$${f}_{(0,0,2)}({{{bf{r}}}})={{{rm{i}}}}{k}_{0}^{3}/6left[1/sqrt{pi }{f}_{0,0}({{{bf{r}}}})-2/sqrt{5pi }{f}_{2,0}({{{bf{r}}}})right]$$
(30)

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 equation39,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 as40,Eq. 7.3.3]

$${{{{bf{{{{mathcal{Y}}}}}}}}}_{jm}^{ell }(theta ,phi )={sum}_{{m}^{{prime} },q}{C}_{ell {m}^{{prime} }1q}^{jm}{{{{mathcal{Y}}}}}_{ell {m}^{{prime} }}(theta ,phi ){{{{bf{e}}}}}_{q}$$
(31)

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 eq form the covariant spherical basis defined in Cartesian coordinates as

$${{{{bf{e}}}}}_{pm 1}=mp frac{1}{sqrt{2}}left({{{{bf{e}}}}}_{x}pm {{{rm{i}}}}{{{{bf{e}}}}}_{y}right)$$
(32)
$${{{{bf{e}}}}}_{0}={{{{bf{e}}}}}_{z}$$
(33)

The vector spherical harmonics used by Jackson24({{{{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

$${{{{bf{J}}}}}_{ell m}^{(0)}({{{bf{r}}}})propto {{{{bf{Y}}}}}_{ell m}^{ell }(nabla )delta ({{{bf{r}}}})$$
(34)

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

$${{{{bf{J}}}}}_{ell m}^{(1)}({{{bf{r}}}})propto {{{{bf{Y}}}}}_{ell pm 1,m}^{ell }(nabla )delta ({{{bf{r}}}})$$
(35)

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

$$[2ell +1]+[2(ell +1)+1]=4ell +4$$
(36)

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

$$3(2ell +1)=6ell +3$$
(37)

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:

$$underbrace{{{4ell +4}}}_{{{{{rm{VSH}}}}}} + underbrace{{{2(ell -1)+1}}}_{{{{rm{degenerate}}}}}=underbrace{{{6ell +3}}}_{{{{rm{Cartesian}}}}}$$
(38)

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.

Illustration on an inverse problem

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

$$overline{overline{varepsilon }}={{{{bf{e}}}}}_{x}{{{{bf{e}}}}}_{x}^{top }+{{{{bf{e}}}}}_{y}{{{{bf{e}}}}}_{y}^{top }+{varepsilon }_{zz}{{{{bf{e}}}}}_{z}{{{{bf{e}}}}}_{z}^{top }$$
(39)

On top of appearing in natural materials (e.g., calcite41), electromagnetic anisotropic media commonly describe metamaterials (e.g., through effective medium theory42), which are seeing an increase of adoption from optical43 to radio frequencies44. High anisotropy can be achieved through 3D printing of carbon fiber composite filament45.

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

$${{{bf{J}}}}({{{bf{r}}}})={nabla }^{alpha }delta ({{{bf{r}}}}){{{{bf{e}}}}}_{x}$$
(40)

at a distance of one vacuum wavelength λ0 simulated using Meep46 in the frequency domain. The simulation domain consists of a cube of side 2λ0 surrounded by a perfectly matched layer of thickness λ0/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,”47 the fitting data (i.e., the prediction through a multipole expansion) is obtained from the analytical Green’s function48. The Cartesian multipole expansion basis functions are obtained by differentiating Green’s function up to a fixed order ({ell }_{max }), i.e.,

$${{{{bf{f}}}}}_{alpha }({{{bf{r}}}})={nabla }^{alpha }{{{bf{G}}}}({{{bf{r}}}}){{{{bf{e}}}}}_{x}$$
(41)

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

$$(nabla {nabla }^{top }-{nabla }^{2}{{{bf{I}}}}-{k}_{0}^{2}overline{overline{varepsilon }}){{{bf{G}}}}({{{bf{r}}}})={{{bf{I}}}}delta ({{{bf{r}}}})$$
(42)

where I is the identity matrix. In turn,

$${{{bf{E}}}}({{{bf{r}}}})=-{{{rm{i}}}}omega mu int !!!iint ,{{{{rm{d}}}}}^{3}{{{{bf{r}}}}}^{{prime} }{{{bf{G}}}}({{{bf{r}}}}-{{{{bf{r}}}}}^{{prime} }){{{bf{J}}}}({{{{bf{r}}}}}^{{prime} })$$
(43)

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

$$f(omega )=int,{{{rm{d}}}}tf(t){{{{rm{e}}}}}^{-{{{rm{i}}}}omega t}$$
(44)

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 = 103) 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 Cm (resp. Cα) that minimize the squared L2-norm of the error between the ground truth and the multipole expansion, normalized to the ground truth squared L2-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.

Fig. 1: Flowchart of the inverse problem for εzz = 4 and α = (0, 0, 2).
figure 1

a The ground truth data consists of all components (pictured: norm of the x component, max-normalized; warm colors indicate high values) of the electric field radiated by the source J(r) = αδ(r)ex simulated using Meep and sampled on a sphere surrounding the source. The explaining data consists of the analytical Cartesian (b) and spherical (c) multipole expansion basis functions. These basis functions are multiplied by coefficients (d) determined by an optimization algorithm (e), whose goal is to minimize the L2-norm of the error (f) between the ground truth and the multipole expansion (g). This error is normalized to the ground truth L2-norm. Its complement after convergence represents the proportion of ground truth energy explained by the Cartesian (h) and spherical (i) multipole expansions. The fields obtained after convergence of the minimizer are shown in (j) for the Cartesian and (k) for the spherical multipole expansion.

Full size image

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 f0,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 waves49 given by the dispersion relations

$${k}_{0}^{2}-{{{{bf{k}}}}}^{top }{{{bf{k}}}}=0$$
(45)
$${varepsilon }_{zz}{k}_{0}^{2}-{{{{bf{k}}}}}^{top }overline{overline{varepsilon }}{{{bf{k}}}}=0$$
(46)

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

Fig. 2: Effect of the truncation order ({ell }_{max }) on the representation power of both multipole expansions (spherical: solid and solid-dotted; Cartesian: dashed) for both source types (α = (0, 0, 0): no markers; α = (0, 0, 2): markers).
figure 2

The number of degrees of freedom (i.e., the number of parameters determined by the inverse problem) is also given on the right-hand side axis. The results circled in red are illustrated in Fig. 1.

Full size image

Fig. 3: Moments recovered by the inverse problem for εzz = 4 and ({ell }_{max }=4).
figure 3

Spherical multipole expansions for (a) α = (0, 0, 0) and (c) α = (0, 0, 2). Cartesian multipole expansions for (b) α = (0, 0, 0) and (d) α = (0, 0, 2). The dot size is proportional to the maximum field amplitude of the corresponding moment. Warm colors indicate higher orders.

Full size image
Fig. 4: Effect of the degree of anisotropy εzz on the representation power of the spherical multipole expansion for the source J(r) = (0, 0, 2)δ(r)ex = ∂2δ(r)/∂z2ex.
figure 4

Solid yellow line: ({ell }_{max }=2); solid-dashed orange line: ({ell }_{max }=4); dashed red line: ({ell }_{max }=6). For εzz > 10, the analytical Green’s function was used as ground truth data.

Full size image

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 103. 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)

$$left({{{bf{k}}}}{{{{bf{k}}}}}^{top }-{{{{bf{k}}}}}^{top }{{{bf{k}}}}{{{bf{I}}}}+{k}_{0}^{2}overline{overline{varepsilon }}right){{{bf{G}}}}({{{bf{k}}}})=-{{{bf{I}}}}$$
(47)

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

$${{{mathbf{G}}}}({{{mathbf{k}}}})to {{{mathbf{G}}}}_infty({{{mathbf{k}}}})=-underbrace{frac{1}{k_0^2-k_z^2}frac{1}{k_0^2-{{{mathbf{k}}}}^{top}{{{mathbf{k}}}}}}_{g({{{mathbf{k}}}})}left[(k_0^2-k_z^2){{{mathbf{I}}}}_{perp}-{{{mathbf{k}}}}_{perp}{{{mathbf{k}}}}_{perp}^{top}right]$$
(48)

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

$$({k}_{0}^{2}-{{{{bf{k}}}}}^{top }{{{bf{k}}}})g({{{bf{k}}}}):=h({{{bf{k}}}})=frac{1}{{k}_{0}^{2}-{k}_{z}^{2}}$$
(49)

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

$${k}_{0}^{2}h({{{bf{r}}}})+frac{{partial }^{2}h({{{bf{r}}}})}{partial {z}^{2}}=delta (x,y)delta (z)$$
(50)

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

$$h({{{bf{r}}}})={c}_{+}{{{{rm{e}}}}}^{{{{rm{i}}}}{k}_{0}z}+{c}_{-}{{{{rm{e}}}}}^{-{{{rm{i}}}}{k}_{0}z}+frac{1}{2{k}_{0}}delta (x,y){{{rm{sign}}}}(z)sin ({k}_{0}z)$$
(51)

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

$${nabla }^{2}{{{{bf{G}}}}}_{infty }({{{bf{r}}}})=-{k}_{0}^{2}{{{{bf{G}}}}}_{infty }({{{bf{r}}}})-left[{{{{bf{I}}}}}_{perp }left({k}_{0}^{2}+frac{{partial }^{2}}{partial {z}^{2}}right)+{nabla }_{perp }{nabla }_{perp }^{top }right]h({{{bf{r}}}})$$
(52)

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 δ(xy). In turn, we find the anisotropic contraction outside this axis:

$${nabla }^{2}{{{{bf{G}}}}}_{infty }({{{bf{r}}}})=-{k}_{0}^{2}{{{{bf{G}}}}}_{infty }({{{bf{r}}}})$$
(53)

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.

Fig. 5
figure 5

Representation power of the spherical harmonics projected onto the xx component of Green’s function in isotropic (εzz = 1, solid blue line) and anisotropic media (εzz = 4, dashed red line).

Full size image

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 sphere38.

Conclusion

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.

Related Articles

Optical sorting: past, present and future

Optical sorting combines optical tweezers with diverse techniques, including optical spectrum, artificial intelligence (AI) and immunoassay, to endow unprecedented capabilities in particle sorting. In comparison to other methods such as microfluidics, acoustics and electrophoresis, optical sorting offers appreciable advantages in nanoscale precision, high resolution, non-invasiveness, and is becoming increasingly indispensable in fields of biophysics, chemistry, and materials science. This review aims to offer a comprehensive overview of the history, development, and perspectives of various optical sorting techniques, categorised as passive and active sorting methods. To begin, we elucidate the fundamental physics and attributes of both conventional and exotic optical forces. We then explore sorting capabilities of active optical sorting, which fuses optical tweezers with a diversity of techniques, including Raman spectroscopy and machine learning. Afterwards, we reveal the essential roles played by deterministic light fields, configured with lens systems or metasurfaces, in the passive sorting of particles based on their varying sizes and shapes, sorting resolutions and speeds. We conclude with our vision of the most promising and futuristic directions, including AI-facilitated ultrafast and bio-morphology-selective sorting. It can be envisioned that optical sorting will inevitably become a revolutionary tool in scientific research and practical biomedical applications.

Light-matter coupling via quantum pathways for spontaneous symmetry breaking in van der Waals antiferromagnetic semiconductors

Light-matter interaction simultaneously alters both the original material and incident light. Light not only reveals material details but also activates coupling mechanisms. The coupling has been demonstrated mechanically, for instance, through the patterning of metallic antennas, resulting in the emergence of plasmonic quasiparticles and enabling wavefront engineering of light via the generalized Snell’s law. However, quantum-mechanical light-matter interaction, wherein photons coherently excite distinct quantum pathways, remains poorly understood. Here, we report on quantum interference between light-induced quantum pathways through the orbital quantum levels and spin continuum. The quantum interference immediately breaks the symmetry of the hexagonal antiferromagnetic semiconductor FePS3. Below the Néel temperature, we observe the emergence of birefringence and linear dichroism, namely, quantum anisotropy due to quantum interference, which is further enhanced by the thickness effect. We explain the direct relevance of the quantum anisotropy to a quantum phase transition by spontaneous symmetry breaking in Mexican hat potential. Our findings suggest material modulation via selective quantum pathways through quantum light-matter interaction.

Dynamic thermalization on noisy quantum hardware

Emulating thermal observables on a digital quantum computer is essential for quantum simulation of many-body physics. However, thermalization typically requires a large system size due to incorporating a thermal bath, whilst limited resources of near-term digital quantum processors allow for simulating relatively small systems. We show that thermal observables and fluctuations may be obtained for a small closed system without a thermal bath. Thermal observables occur upon classically averaging quantum mechanical observables over randomized variants of their time evolution that run independently on a digital quantum processor. Using an IBM quantum computer, we experimentally find thermal occupation probabilities with finite positive and negative temperatures defined by the initial state’s energy. Averaging over random evolutions facilitates error mitigation, with the noise contributing to the temperature in the simulated observables. This result fosters probing the dynamical emergence of equilibrium properties of matter at finite temperatures on noisy intermediate-scale quantum hardware.

An observational study of pleiotropy and penetrance of amyotrophic lateral sclerosis associated with CAG-repeat expansion of ATXN2

Spinocerebellar ataxia type 2 (SCA2) and amyotrophic lateral sclerosis (ALS) are both associated with a CAG-repeat expansion in ATXN2 and with TDP-43-positive neuronal cytoplasmic inclusions. The two disorders have been viewed as distinct entities, where an intermediate length expansion of 31-33 CAG-repeats is associated with sporadic ALS and a full length expansion of ≥34 CAG-repeats is associated with SCA2. We report the clinical phenotype of ATXN2-positive patients and their relatives, identified in three specialist ALS clinics, which force a reconsideration of this dichotomy. We also report the frequency of ATXN2 expansions in two large cohorts of ALS patients and in a population-matched cohort of controls. We report ten cases of familial ALS in which disease is associated with either an intermediate or a full-length ATXN2 CAG-repeat expansion. Pedigrees and patients feature additional phenotypes including parkinsonism, dementia and essential tremor (ET). We conclude that CAG-repeat expansions in ATXN2 exhibit pleiotropy and are associated with a disease spectrum that includes ALS, SCA2, and parkinsonism; to recognise this complexity we propose the new term ‘ATXN2-related neurodegeneration’. We also observed sporadic ALS associated with full-length expansions. We conclude that ATXN2 CAG-repeat expansions, irrespective of length, should be considered a risk factor for ALS. Interrupted CAG-repeats were associated with an ALS phenotype in our data but we also identified ALS cases with uninterrupted expansions. Our findings have relevance for researchers, patients and families linked to CAG-repeat expansions in ATXN2.

Collective quantum enhancement in critical quantum sensing

Critical systems represent a valuable resource in quantum sensing and metrology. Critical quantum sensing (CQS) protocols can be realized using finite-component phase transitions, where criticality arises from the rescaling of system parameters rather than the thermodynamic limit. Here, we show that a collective quantum advantage can be achieved in a multipartite CQS protocol using a chain of parametrically coupled critical resonators in the weak-nonlinearity limit. We derive analytical solutions for the low-energy spectrum of this unconventional quantum many-body system, which is composed of locally critical elements. We then assess the scaling of the quantum Fisher information with respect to fundamental resources. We demonstrate that the coupled chain outperforms an equivalent ensemble of independent critical sensors, achieving quadratic scaling in the number of resonators. Finally, we show that even with finite Kerr nonlinearity or Markovian dissipation, the critical chain retains its advantage, making it relevant for implementing quantum sensors with current microwave superconducting technologies.

Responses

Your email address will not be published. Required fields are marked *