Electromagnetic Scattering by Bianisotropic Spheres

Abstract

Electromagnetic fields in bulk bianisotropic media contain plane waves whose k-vectors can be found using the method of the index of refraction’s operator and belong to the Fresnel wave surfaces that fall into one of the five hyperbolic classes of the Durach et al. taxonomy of bianisotropic media. Linear combinations of vector spherical harmonics can be used as a set of solutions of vector Helmholtz equations in gyrotropic media to develop Mie’s theory of scattering by anisotropic spheres as accomplished by Lin et al. and Li et al. In this study, we introduced electromagnetic orbitals for bianisotropic media as linear combinations of vector spherical harmonics, which represent solutions of Maxwell’s equations in bianisotropic media. Using these bianisotropic orbitals, we developed a theory of the scattering of electromagnetic radiation by bianisotropic spheres with arbitrary effective material parameters and sizes. As a by-product, we obtained a simple expression for the expansion of a vector plane wave over vector spherical harmonics in a more compact form than the frequently used by Sarkar et al. We obtained the polarizability expressions in the Rayleigh limit in agreement with the results of the electrostatic approximation of Lakhtahia and Sihvola.

1. Introduction

Electromagnetism occupies the crowning role in physics, science, and modern technology. As in the cases of the Second and the Third Industrial Revolutions, research into electromagnetism is driving the ongoing Fourth Industrial Revolution [1] related to the transition to renewable energy, telecommunications in the 5G and 6G standards [2], advanced micro-/nanofabrication for novel electronic devices [3], bioelectromagnetics [4], information and electronic warfare [5], machine learning, material training [6,7], and in other realms. The slide towards scientific and technological unification of the physical, chemical, biological, and digital worlds brought by the Fourth Industrial Revolution is due to the inalienable electromagnetic nature of these phenomena, based on the rule of the underlying laws of Coulomb, Gauss, Biot-Savart, Ampere, Kirchhoff, Faraday, and Maxwell [8,9,10]. We, the researchers of modern electromagnetism, are devoted to the development of the new electromagnetic materials collectively called composite artificial materials or metamaterials [11,12,13,14,15,16,17]. Metamaterials are made of arrays of subwavelength scatterers designed to exhibit the desired electromagnetic properties. In many important cases, metamaterials can be described as bianisotropic media [13,14,15,16,17]. In bianisotropic media, both the electric and magnetic responses depend on both the electric and magnetic fields of the external radiation [18,19].

The studies of bianisotropic materials are almost as old as electromagnetism itself, persisting through the 19th and 20th centuries in the work of scientists such as Roentgen, Wilson, Landau, Lifshitz, Dzyaloshinskii, Cheng, and Kong [18,19,20,21,22,23,24]. In the 21st century the field of bianisotropic optical materials has received the name of bianisotropics [25,26,27] and is closely related to the research on electromagnetic metamaterials, since, typically, the desired properties of metamaterials depend on them being anisotropic and bianisotropic media [13,14,15,16,17,18,19]. Despite all these efforts and the rich history of research into bianisotropics until recently, very few general properties of bianisotropic media were established due to the complexity of these media [14,18,28,29,30,31,32,33,34,35]. The bianisotropic media are the most general case of local linear media [18,19,25,36,37], with the effective material parameters combined into a 6 × 6 matrix of material parameters $\widehat{M}$, which characterizes the electric displacement field D and magnetic field B in terms of the fields E and H:

$$\left(\begin{array}{c}\mathit{D}\\ \mathit{B}\end{array}\right)=\widehat{M}\left(\begin{array}{c}\mathit{E}\\ \mathit{H}\end{array}\right)=\left(\begin{array}{cc}\widehat{ϵ}& \widehat{X}\\ \widehat{Y}& \widehat{\mu }\end{array}\right)\left(\begin{array}{c}\mathit{E}\\ \mathit{H}\end{array}\right)$$

(1)

The 3 × 3 matrices $\widehat{ϵ}, \widehat{\mu }, \widehat{X}, \widehat{Y}$ are dielectric permittivity, magnetic permeability, and two magnetoelectric coupling matrices respectively. The inverse relationship can also be formulated:

$$\left(\begin{array}{c}\mathit{E}\\ \mathit{H}\end{array}\right)={\left[\begin{array}{cc}\widehat{ϵ}& \widehat{X}\\ \widehat{Y}& \widehat{\mu }\end{array}\right]}^{-1}\left(\begin{array}{c}\mathit{D}\\ \mathit{B}\end{array}\right)=\left[\begin{array}{cc}{\widehat{\alpha }}_{ED}& {\widehat{\alpha }}_{EB}\\ {\widehat{\alpha }}_{HB}& {\widehat{\alpha }}_{HD}\end{array}\right]\left(\begin{array}{c}\mathit{D}\\ \mathit{B}\end{array}\right)$$

(2)

The relationships in Equations (1) and (2) mean that, unlike in many naturally occurring isotropic media, the electrical polarization and magnetization vectors in bianisotropic media are not directed in the same direction as the electric and magnetic fields. Such differences in the isotropic media can be achieved in engineered metamaterial structures [13,14,15,16,17].

One of the jewels in the crown of electromagnetism is Mie’s theory of scattering by spheres. Originally, Mie’s theory was introduced to describe scattering by isotropic spheres [38,39], but later, it was extended to describe scattering by bi-isotropic [40], rotationally symmetric and anisotropic [41], orthorhombic and dielectric–magnetic [42], magneto- and electro-gyrotropic [43,44,45], and spherically-symmetric bianisotropic [46] spheres. Despite this activity, the theory of scattering by a generic bianisotropic sphere has not yet been constituted [47].

The plane waves which can propagate in bianisotropic media belong to Fresnel wave surfaces, which can be characterized using the method of the index of refraction’s operator [31,32,33]. The Fresnel wave surfaces in bianisotropic media follow quartic dispersion equations, and, therefore, can be classified using the taxonomy of Durach et al. [31,32,33], which includes the five hyperbolic classes: non-, mono-, bi-, tri-, and tetra-hyperbolic materials [31,32,33,48,49,50,51,52,53,54,55,56,57,58,59]. The prefix in the name of each topological class indicates the number of double cones that the iso-frequency’s k-surface has in its high-k limit. Note that hyperbolic metamaterials, which are already known for their applications in optical imaging, hyperlensing, and emission rate and directivity control, utilize the diverging optical density of high-k states [51,52,53,54,55,56,57,58,59]. In Figure 1a, we show an example of an iso-frequency Fresnel wave surface for a tetra-hyperbolic bianisotropic medium with the effective material parameter matrix $\widehat{M}$ color-coded in Figure 1b.

The Fresnel wave surfaces only include the plane waves with real wave vectors k. Nevertheless, the inhomogeneous plane waves with an imaginary k are also important. In Figure 1c, we plot the eigenvalues of the index of refraction’s operator [31,32,33] in a complex plane for the material color-coded in Figure 1b for plane waves propagating in directions spanning the entire solid angle with steps in the angles Δθ = π/40 and Δϕ = 2π/40. Note that in Figure 1c, only the eigenvalues on the horizontal axis $\mathrm{Re}\left\{k/k_0\right\}$ correspond to the real wave vectors belonging to the Fresnel wave surfaces depicted in Figure 1a.

In unbounded homogeneous bianisotropic media, plane waves represent the set of solutions of Maxwell’s equations which is typically used, but to consider the scattering of electromagnetic waves by bianisotropic spheres, it is more convenient to express the electromagnetic fields inside such spheres in terms of vector spherical harmonics. In this study, we proposed a theory of scattering by bianisotropic spheres with arbitrary effective media parameters $\widehat{M}$. To accomplish this, we introduced bianisotropic orbitals composed of vector spherical harmonics.

2. Results

2.1. Bianisotropic Orbitals

The plane waves whose k-vectors belong to Fresnel’s wave surfaces represent a complete set of solutions of Maxwell’s equations in the bulk bianisotropic media. A different set of solutions can be found as expansion of vector spherical harmonics. Recently, the close connection between the multipole composition of electromagnetic fields and bianisotropy has been revealed [60].

Due to the solenoidal nature of D- and B-fields, namely ∇⋅D = 0 and ∇⋅B = 0, we can express them as an expansion over the vector spherical harmonics ${\mathit{M}}_{lm}^1, {\mathit{N}}_{lm}^1$ found from solutions of the scalar Helmholtz equations [61,62,63]

$${\psi }_{lm}^{\left(j\right)}=-\frac{1}{i\sqrt{l\left(l+1\right)}}{z}_l^{\left(j\right)}\left(kr\right)Y_{lm},\hspace{1em}\hspace{1em}{\psi }_{00}^{\left(j\right)}=i{z}_0^{\left(j\right)}\left(kr\right)Y_{00}$$

(3)

according to

$${\mathit{L}}_{lm}^{\left(j\right)}=\frac{1}{k}\nabla {\psi }_{lm}^{\left(j\right)},\hspace{1em}\hspace{1em}{\mathit{M}}_{lm}^{\left(j\right)}=\nabla \times \left(\mathit{r}{\psi }_{lm}^{\left(j\right)}\right),\hspace{1em}\hspace{1em}{\mathit{N}}_{lm}^{\left(j\right)}=\frac{1}{k}\nabla \times {\mathit{M}}_{lm}^{\left(j\right)}$$

(4)

Detailed definitions of the vector spherical harmonics used in this study are in Appendix A. We represent the D- and B-fields of a bianisotropic orbital with a wave number $k_q$ as

$${\mathit{D}}_q={\displaystyle \sum }_{lm}\left\{f_{qlm}^{DM}{\mathit{M}}_{lm}^1\left(k_q\right)+f_{qlm}^{DN}{\mathit{N}}_{lm}^1\left(k_q\right)\right\}$$

(5)

$${\mathit{B}}_q={\displaystyle \sum }_{lm}\left\{f_{qlm}^{BM}{\mathit{M}}_{lm}^1\left(k_q\right)+f_{qlm}^{BN}{\mathit{N}}_{lm}^1\left(k_q\right)\right\}$$

(6)

The corresponding E- and H-fields can be expressed using Equation (2) as

$$\begin{array}{c}{\mathit{E}}_q={\displaystyle \sum }_{uv}\left({\mu }_{uv}^{eq}{\mathit{M}}_{uv}^{1q}+{\nu }_{uv}^{eq}{\mathit{N}}_{uv}^{1q}+{\lambda }_{uv}^{eq}{\mathit{L}}_{uv}^{1q}\right)={\widehat{\alpha }}_{ED}{\mathit{D}}_q+{\widehat{\alpha }}_{EB}{\mathit{B}}_q\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \sum }_{lm}\left\{f_{qlm}^{DM}\left({\widehat{\alpha }}_{ED}{\mathit{M}}_{lm}^{1q}\right)+f_{qlm}^{DN}\left({\widehat{\alpha }}_{ED}{\mathit{N}}_{lm}^{1q}\right)\right\}+{\displaystyle \sum }_{lm}\left\{f_{qlm}^{BM}\left({\widehat{\alpha }}_{EB}{\mathit{M}}_{lm}^{1q}\right)+f_{qlm}^{BN}\left({\widehat{\alpha }}_{EB}{\mathit{N}}_{lm}^{1q}\right)\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \sum }_{lm,uv}\left({\mathit{M}}_{uv}^{1q},{\mathit{N}}_{uv}^{1q},{\mathit{L}}_{uv}^{1q}\right)\cdot {\left(\begin{array}{cccc}{g}_{MM}^{{\alpha }_{ED}}& g_{NM}^{{\alpha }_{ED}}& g_{MM}^{{\alpha }_{EB}}& g_{NM}^{{\alpha }_{EB}}\\ g_{MN}^{{\alpha }_{ED}}& g_{NN}^{{\alpha }_{ED}}& g_{MN}^{{\alpha }_{EB}}& g_{NN}^{{\alpha }_{EB}}\\ g_{ML}^{{\alpha }_{ED}}& g_{NL}^{{\alpha }_{ED}}& g_{ML}^{{\alpha }_{EB}}& g_{NL}^{{\alpha }_{EB}}\end{array}\right)}_{lm,uv}\cdot \left(\begin{array}{c}{f}_{qlm}^{DM}\\ f_{qlm}^{DN}\\ f_{qlm}^{BM}\\ f_{qlm}^{BN}\end{array}\right)\hfill \end{array}$$

(7)

$$\begin{array}{c}{\mathit{H}}_q={\displaystyle \sum }_{uv}\left({\mu }_{uv}^{hq}{\mathit{M}}_{uv}^{1q}+{\nu }_{uv}^{hq}{\mathit{N}}_{uv}^{1q}+{\lambda }_{uv}^{hq}{\mathit{L}}_{uv}^{1q}\right)={\widehat{\alpha }}_{HD}{\mathit{D}}_q+{\widehat{\alpha }}_{HB}{\mathit{B}}_q\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \sum }_{lm}\left\{f_{qlm}^{DM}\left({\widehat{\alpha }}_{HD}{\mathit{M}}_{lm}^{1q}\right)+f_{qlm}^{DN}\left({\widehat{\alpha }}_{HD}{\mathit{N}}_{lm}^{1q}\right)\right\}+{\displaystyle \sum }_{lm}\left\{f_{qlm}^{BM}\left({\widehat{\alpha }}_{HB}{\mathit{M}}_{lm}^{1q}\right)+f_{qlm}^{BN}\left({\widehat{\alpha }}_{HB}{\mathit{N}}_{lm}^{1q}\right)\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \sum }_{lm,uv}\left({\mathit{M}}_{uv}^{1q},{\mathit{N}}_{uv}^{1q},{\mathit{L}}_{uv}^{1q}\right)\cdot {\left(\begin{array}{cccc}{g}_{MM}^{{\alpha }_{HD}}& g_{NM}^{{\alpha }_{HD}}& g_{MM}^{{\alpha }_{HB}}& g_{NM}^{{\alpha }_{HB}}\\ g_{MN}^{{\alpha }_{HD}}& g_{NN}^{{\alpha }_{HD}}& g_{MN}^{{\alpha }_{HB}}& g_{NN}^{{\alpha }_{HB}}\\ g_{ML}^{{\alpha }_{HD}}& g_{NL}^{{\alpha }_{HD}}& g_{ML}^{{\alpha }_{HB}}& g_{NL}^{{\alpha }_{HB}}\end{array}\right)}_{lm,uv}\cdot \left(\begin{array}{c}{f}_{qlm}^{DM}\\ f_{qlm}^{DN}\\ f_{qlm}^{BM}\\ f_{qlm}^{BN}\end{array}\right)\hfill \end{array}$$

(8)

where the coefficients $g$ is found using $\left\{\mathit{U}\left|\widehat{\alpha }\right|\mathit{V}\right\}={\int }_0^{\infty }{\int }_0^{2\pi }{\int }_0^{\pi }\mathit{U}\cdot \widehat{\alpha }\cdot \mathit{V}{r}^2\sin \theta dr d\theta d\phi$, $\int Y_{uv}^{*}\left(\widehat{\mathit{r}}\right)Y_{lm}\left(\widehat{\mathit{r}}\right)d\mathsf{\Omega }={\delta }_{ul}{\delta }_{vm}$, and ${\int }_0^{\infty }{j}_u\left(k^{\prime }r\right)j_l\left(kr\right)r^{2}dr=\frac{\pi }{2k^2}\delta \left(k-k^{\prime }\right)$ as:

$$\begin{array}{c}\left(\begin{array}{ccc}\left\{{\mathit{M}}_{uv}^{\star }\left|\widehat{\alpha }\right|{\mathit{M}}_{lm}\right\}& \left\{{\mathit{N}}_{uv}^{\star }\left|\widehat{\alpha }\right|{\mathit{M}}_{lm}\right\}& \left\{{\mathit{L}}_{uv}^{\star }\left|\widehat{\alpha }\right|{\mathit{M}}_{lm}\right\}\\ \left\{{\mathit{M}}_{uv}^{\star }\left|\widehat{\alpha }\right|{\mathit{N}}_{lm}\right\}& \left\{{\mathit{N}}_{uv}^{\star }\left|\widehat{\alpha }\right|{\mathit{N}}_{lm}\right\}& \left\{{\mathit{L}}_{uv}^{\star }\left|\widehat{\alpha }\right|{\mathit{N}}_{lm}\right\}\\ \left\{{\mathit{M}}_{uv}^{\star }\left|\widehat{\alpha }\right|{\mathit{L}}_{lm}\right\}& \left\{{\mathit{N}}_{uv}^{\star }\left|\widehat{\alpha }\right|{\mathit{L}}_{lm}\right\}& \left\{{\mathit{L}}_{uv}^{\star }\left|\widehat{\alpha }\right|{\mathit{L}}_{lm}\right\}\end{array}\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\left(\begin{array}{ccc}{g}_{MM}^{\alpha }& g_{MN}^{\alpha }& g_{ML}^{\alpha }\\ g_{NM}^{\alpha }& g_{NN}^{\alpha }& g_{NL}^{\alpha }\\ g_{LM}^{\alpha }& g_{LN}^{\alpha }& g_{LL}^{\alpha }\end{array}\right)}_{lm,uv}\left(\begin{array}{ccc}\left\{{\mathit{M}}_{uv}^{\star }|{\mathit{M}}_{uv}\right\}=1& 0& 0\\ 0& \left\{{\mathit{N}}_{uv}^{\star }|{\mathit{N}}_{uv}\right\}=1& 0\\ 0& 0& \left\{{\mathit{L}}_{uv}^{\star }|{\mathit{L}}_{uv}\right\}\end{array}\right)\hfill \end{array}$$

(9)

Please note that, while the D- and B-fields are divergence-free, as indicated by Equations (5) and (6), the E- and H-fields include the longitudinal vector spherical harmonic $\mathit{L}$, as was previously shown for anisotropic media as well [43,44,45].

The bianisotropic orbitals with the radial quantum numbers $k_q$ represent the solutions of Maxwell’s equations in homogeneous bianisotropic media if the expansion coefficients $f_{qlm}$ satisfy the following eigenproblem:

$${\displaystyle \sum }_{lm}{\left(\begin{array}{cccc}-g_{MN}^{{\alpha }_{HD}}& -g_{NN}^{{\alpha }_{HD}}& -g_{MN}^{{\alpha }_{HB}}& -g_{NN}^{{\alpha }_{HB}}\\ -g_{MM}^{{\alpha }_{HD}}& -g_{NM}^{{\alpha }_{HD}}& -g_{MM}^{{\alpha }_{HB}}& -g_{NM}^{{\alpha }_{HB}}\\ g_{MN}^{{\alpha }_{ED}}& g_{NN}^{{\alpha }_{ED}}& g_{MN}^{{\alpha }_{EB}}& g_{NN}^{{\alpha }_{EB}}\\ g_{MM}^{{\alpha }_{ED}}& g_{NM}^{{\alpha }_{ED}}& g_{MM}^{{\alpha }_{EB}}& g_{NM}^{{\alpha }_{EB}}\end{array}\right)}_{lm,uv}\left(\begin{array}{c}{f}_{qlm}^{DM}\\ f_{qlm}^{DN}\\ f_{qlm}^{BM}\\ f_{qlm}^{BN}\end{array}\right)=i\left(\frac{k_0}{k_q}\right)\left(\begin{array}{c}{f}_{quv}^{DM}\\ f_{quv}^{DN}\\ f_{quv}^{BM}\\ f_{quv}^{BN}\end{array}\right)$$

(10)

Note that the eigenproblem of Equation (10) differs from the eigenproblems formulated for the anisotropic media, which stem from the vector Helmholtz equations [43,44,45]. Such eigenproblems could not be used in the case of the bianisotropic media considered here and we used the pair of Maxwell’s equations composed of Faraday’s and Maxwell-Ampere’s equations instead to obtain Equation (10).

The bianisotropic orbitals provided by the solutions of Equation (10) can be represented as expansions over plane wave solutions of the method of the index of refraction’s operator [31,32,33] with the indexes of refraction $n=k_q/k_0$ (Appendix B). In Figure 2a, we show the inverse eigenvalues $k_q/k_0$ of Equation (10) for the eigenproblems which are truncated and have $l=4$ (black dots), 10 (red), and 40 (green). Note the direct correspondence of the eigenvalues of Equation (10) in Figure 2a with the eigenvalues of the index of refraction’s operator plotted in Figure 1b.

In other words, the inclusion of higher multipoles in the bianisotropic orbitals corresponds to higher angular resolution in the Fourier expansion of the fields in bianisotropic materials with a smaller Δθ and Δϕ with the inclusion of both the real wave vectors belonging to the Fresnel wave surface, as well as the complex eigenvalues $k/k_0$ of the index of refractioná operator.

This correspondence between the plane waves and the bianisotropic orbitals introduced in this study shows the relationship between the solution of the scattering problem presented in this study with the method of plane wave expansion proposed for the scattering by the uniaxial anisotropic spheres [64,65]. In Figure 2b–e, we plot the components ${\left|f_{qlm}\right|}^2$ of the eigenproblem of Equation (10) for a bianisotropic orbital with $k_q/k_0=2.2$ in the angular momentum space $l$-$m$.

2.2. Scattering Cross-Section of Bianisotropic Spheres in a Vacuum

The bianisotropic orbitals can be used to solve a large range of problems with spherically shaped bianisotropic media from spheres and spherical shells, to spherical voids, combinations of such geometries, and so forth. Here, we considered a bianisotropic sphere with the arbitrary effective material parameters $\widehat{M}$ with a radius $R$. We studied the scattering of an electromagnetic plane wave by such a sphere. The field of the plane wave is given by

$${\mathit{E}}_{in}=\mathbf{ϵ}\left(\mathit{k}\right)e^{i\mathit{k}\mathit{r}}={\displaystyle \sum }_{lm}\left(q_{lm}{\mathit{M}}_{lm}^1+p_{lm}{\mathit{N}}_{lm}^1\right)$$

(11)

$${\mathit{H}}_{in}=\mathit{h}\left(\mathit{k}\right)e^{i\mathit{k}\mathit{r}}=\frac{1}{i}{\displaystyle \sum }_{lm}\left(p_{lm}{\mathit{M}}_{lm}^1+q_{lm}{\mathit{N}}_{lm}^1\right),\hspace{1em}\hspace{1em}\mathit{h}=\widehat{\mathit{k}}\times \mathbf{ϵ}$$

(12)

We defined the orientation of the incident electric and magnetic fields with the polarization angle $\alpha$ in the spherical coordinates with respect to the direction of incidence $\widehat{\mathit{k}}=\left(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta \right)$, $\mathbf{ϵ}\left(\mathit{k}\right)=\left(\sin \alpha \widehat{\mathit{\theta }}-\cos \alpha \widehat{\mathit{\varphi }}\right)$, $\mathit{h}\left(\mathit{k}\right)=\left(\cos \alpha \widehat{\mathit{\theta }}+\sin \alpha \widehat{\mathit{\varphi }}\right)$, where $\widehat{\mathit{\theta }}=\left(\cos \theta \cos \varphi ,\cos \theta \sin \varphi ,-\sin \theta \right)$, $\widehat{\mathit{\varphi }}=\left(-\sin \varphi ,\cos \varphi \right)$.

Expansions of the vector plane wave over vector spherical harmonics exist in the literature [43,66]. Nevertheless, we derived a compact expansion of a vector plane wave for our work (Appendix C):

$$\mathit{a}{e}^{i\mathit{k}\mathit{r}}=4\pi {\displaystyle \sum }_{lm}{i}^l\left(-\sqrt{l\left(l+1\right)}\left\{\mathit{a}\cdot {\mathit{Y}}_{lm}^{\left(-1\right)*}\left(\widehat{\mathit{k}}\right)\right\}{\mathit{L}}_{lm}^{\left(1\right)}+\left\{\mathit{a}\cdot {\mathit{Y}}_{lm}^{l*}\left(\widehat{\mathit{k}}\right)\right\}{\mathit{M}}_{lm}^{\left(1\right)}-\left\{\mathit{a}\cdot {\mathit{Y}}_{lm}^{\left(+1\right)*}\left(\widehat{\mathit{k}}\right)\right\}{\mathit{N}}_{lm}^{\left(1\right)}\right)$$

Correspondingly, the coefficients in the expansions Equations (11) and (12) are

$$q_{lm}=4\pi i^l\left\{\mathbf{ϵ}\left(\mathit{k}\right)\cdot {\mathit{Y}}_{lm}^{l*}\left(\widehat{\mathit{k}}\right)\right\},\hspace{1em}\hspace{1em}{p}_{lm}=4\pi i^{l+1}\left\{\mathit{h}\left(\mathit{k}\right)\cdot {\mathit{Y}}_{lm}^{l*}\left(\widehat{\mathit{k}}\right)\right\}$$

(13)

The scattered fields outside of the sphere are given by

$${\mathit{E}}_{sc}={\displaystyle \sum }_{lm}\left(b_{lm}{\mathit{M}}_{lm}^3+a_{lm}{\mathit{N}}_{lm}^3\right)$$

(14)

$${\mathit{H}}_{sc}=\frac{1}{i}{\displaystyle \sum }_{lm}\left(a_{lm}{\mathit{M}}_{lm}^3+b_{lm}{\mathit{N}}_{lm}^3\right)$$

(15)

The scattered fields consist only of the outgoing vector spherical harmonics containing the spherical Hankel functions $h_l^{\left(1\right)}$.

We represent the field inside the bianisotropic sphere as a linear combination of the bianisotropic orbitals described by Equation (10).

$${\mathit{E}}_{sph}={\displaystyle \sum }_q{A}_q{\mathit{E}}_q,\hspace{1em}\hspace{1em}{\mathit{H}}_{sph}={\displaystyle \sum }_q{A}_q{\mathit{H}}_q,$$

(16)

Please note that for reciprocal media, the eigenvalues of Equation (10) come in pairs with the radial wavenumber $\pm k_q$, and the corresponding modes in Equation (16) are equivalent. Therefore, among the orbitals with a real $k_q$, we only included the modes with a positive $k_q$. For the orbitals with a complex $k_q$, we only included the modes with $\mathrm{Im} k_q>0$. Note, however, that in non-reciprocal media, the symmetry of the reciprocity between the modes is broken, and all the bianisotropic orbitals of Equation (10) should be included into Equation (16). Correspondingly, additional boundary conditions (ABC) are needed, as described in [34].

The continuity of the tangential components of the E- and H-fields at the surface of the sphere lead to the following boundary conditions, expressed in terms of the Riccati–Bessel functions ${\psi }_l\left(x\right)=x{j}_l\left(x\right)$ and ${\xi }_l\left(x\right)=x{h}_l^{\left(1\right)}\left(x\right)$, where $j_l\left(x\right)$ and $h_l^{\left(1\right)}\left(x\right)$ are spherical Bessel functions of the first and third kinds, and the parameters $x=k_{0}R$ and $x_q=k_{q}R$:

$$E_{in\theta }+E_{sc\theta }=E_{sph\theta },\hspace{1em}{q}_{lm}+b_{lm}\left(\frac{{\xi }_l\left(x\right)}{{\psi }_l\left(x\right)}\right)={\displaystyle \sum }_q{A}_q\left(\frac{x}{x_q}\right){\mu }_{lm}^{eq}\left(\frac{{\psi }_l\left(x_q\right)}{{\psi }_l\left(x\right)}\right)$$

(17)

$$H_{in\theta }+H_{sc\theta }=H_{sph\theta },\hspace{1em}{p}_{lm}+a_{lm}\left(\frac{{\xi }_l\left(x\right)}{{\psi }_l\left(x\right)}\right)=i{\displaystyle \sum }_q{A}_q\left(\frac{x}{x_q}\right){\mu }_{lm}^{hq}\left(\frac{{\psi }_l\left(x_q\right)}{{\psi }_l\left(x\right)}\right)$$

(18)

$$E_{in\varphi }+E_{sc\varphi }=E_{sph\varphi },\hspace{1em}\hspace{1em}{p}_{lm}+a_{lm}\left(\frac{{\xi }_l^{\prime }\left(x\right)}{{\psi }_l^{\prime }\left(x\right)}\right)={\displaystyle \sum }_q{A}_q\left(\frac{x}{x_q}\right)\left\{{\nu }_{lm}^{eq}\left(\frac{{\psi }_l^{\prime }\left(x_q\right)}{{\psi }_l^{\prime }\left(x\right)}\right)+{\lambda }_{lm}^{eq}\left(\frac{j_l\left(x_q\right)}{{\psi }_l^{\prime }\left(x\right)}\right)\right\}$$

(19)

$$H_{in\varphi }+H_{sc\varphi }=H_{sph\varphi },\hspace{1em}\hspace{1em}{q}_{lm}+b_{lm}\left(\frac{{\xi }_l^{\prime }\left(x\right)}{{\psi }_l^{\prime }\left(x\right)}\right)=i{\displaystyle \sum }_q{A}_q\left(\frac{x}{x_q}\right)\left\{{\nu }_{lm}^{hq}\left(\frac{{\psi }_l^{\prime }\left(x_q\right)}{{\psi }_l^{\prime }\left(x\right)}\right)+{\lambda }_{lm}^{hq}\left(\frac{j_l\left(x_q\right)}{{\psi }_l^{\prime }\left(x\right)}\right)\right\}$$

(20)

The boundary Equations (17)–(20) can be represented in the matrix form as

$${\widehat{{\rm M}}}_{\psi }\overrightarrow{A}-\widehat{\Psi }\left(ab\right)=\left(pq\right)$$

(21)

$${\widehat{{\rm N}}}_{\psi }\overrightarrow{A}-\widehat{\Phi }\left(ab\right)=\left(pq\right)$$

(22)

where $\overrightarrow{A}$ is the column of the amplitudes of the bianisotropic orbitals $A_q$ inside the sphere, and $\left(ab\right)$ and $\left(pq\right)$ include the amplitudes of the scattered vector spherical harmonics and the known incident amplitudes of Equation (13).

The matrices $\widehat{\Psi }, \widehat{\Phi },$ are diagonal and contain the functions $\left(\frac{{\xi }_l\left(x\right)}{{\psi }_l\left(x\right)}\right)$ and $\left(\frac{{\xi }_l^{\prime }\left(x\right)}{{\psi }_l^{\prime }\left(x\right)}\right)$, while the matrices ${\widehat{{\rm M}}}_{\psi }$, ${\widehat{{\rm N}}}_{\psi }$ represent the coupling between different multipoles due to the bianisotropy in the bianisotropic orbitals, and correspond to the sums on the right-hand sides of Equations (17)–(20).

Excluding the coefficients $\left(pq\right)$, we obtained the relationship between $\left(ab\right)$ and $\overrightarrow{A}$

$$\left(ab\right)=\widehat{\Omega }\overrightarrow{A},\hspace{1em}\widehat{\Omega }={\left\{\widehat{\Psi }-\widehat{\Phi }\right\}}^{-1}\left\{{\widehat{{\rm M}}}_{\psi }-{\widehat{{\rm N}}}_{\psi }\right\}$$

(23)

The amplitudes of the bianisotropic orbitals $A_q$ were found from Equations (21) and (23) as follows

$$\overrightarrow{A}=\widehat{\Xi }\left(pq\right)={\left\{{\widehat{{\rm M}}}_{\psi }-\widehat{\Psi }\widehat{\Omega }\right\}}^{-1}\left(pq\right)$$

(24)

Substituting the amplitudes $A_q$ from Equation (24) into Equation (23), we obtained the T-matrix for a bianisotropic sphere with arbitrary effective medium parameters and the scattering amplitudes $\left(ab\right)$ in terms of the parameters of the incident wave $\left(pq\right)$ as

$$\left(ab\right)=\widehat{T}\left(pq\right),\hspace{1em}\widehat{T}=\widehat{\Omega }\widehat{\Xi }$$

(25)

The scattering cross-section $Q_s$ can be found from the scattering amplitudes as

$$Q_s=\frac{1}{k_0^2}{\displaystyle \sum }_{lm}\left({\left|a_{lm}\right|}^2+{\left|b_{lm}\right|}^2\right)$$

(26)

To validate our formulas and codes and check the accuracy of the numerical results obtained, in Figure 3, we compare our results with the published results of Lin et al. and Li et al. [43,45]. In Figure 3a, we show the scattering cross-section $Q_s/\left(\pi R^2\right)$ as a function of the angles $\theta$ for the anisotropic sphere with $\widehat{ϵ}=1$ and $\widehat{\mu }=\widehat{1}+\left({\mu }_s-1\right)\widehat{\mathit{z}}\widehat{\mathit{z}}$ for various values of ${\mu }_s$ in linear polarization for $\alpha =0$. In Figure 3b, we show the scattering cross-section $Q_s/\left(\pi R^2\right)$ as a function of the angles $\theta$ for the gyromagnetic sphere with $\widehat{ϵ}=1$ and $\widehat{\mu }=\widehat{1}-i{\mu }_g\left(\widehat{\mathit{x}}\widehat{\mathit{y}}-\widehat{\mathit{y}}\widehat{\mathit{x}}\right)$ for various values of ${\mu }_g$ in left- and right-handed circular polarizations. Figure 3 is an exact match with the results obtained in Figure 2 of Ref. [43] and Figure 2 and Figure 3 of Ref. [45].

In Figure 4a,b, we plot the dependence of the reduced scattering cross-section $Q_s/\left(\pi R^2\right)$ for all incidence angles $\theta$ and $\varphi$ for a bianisotropic sphere with the effective material parameters color-coded in Figure 1b and with $x=k_{0}R=1$. Please note that to emphasize that our method is applicable to arbitrary bianisotropic materials, the effective parameters matrix $\widehat{M}$ color-coded in the bottom of panel Figure 1b corresponds to a reciprocal material with effective parameters randomly generated in the range between −5 and 5. This material features anisotropic dielectric permittivity, magnetic permeability, and chirality tensors, which can be engineered by combining split-ring, helix, omega, fishnet, parallel-plate, and wire metaatoms [16,17].

Figure 4a corresponds to the incidence polarization angle $\alpha =0$, while in Figure 4b, the polarization angle is $\alpha =\pi /2$.

In Figure 4, one can see the strong dependence of scattering by the bianisotropic spheres in the incidence direction and polarization.

2.3. Expressions of Polarizability in the Rayleigh Limit

In the Rayleigh limit of electromagnetically small spheres $x\ll 1$, only the dipole terms of Mie’s theory are retained. For $l, u=1$, $g_{MN}^{\alpha }=g_{ML}^{\alpha }=g_{NM}^{\alpha }=0$ and the E- and H-fields, Equations (7) and (8) can be written as

$$\mathit{E}={\displaystyle \sum }_{mv}{\mathit{M}}_{1v}^1{\left(g_{MM}^{{\alpha }_{ED}},g_{MM}^{{\alpha }_{EB}}\right)}_{m,v}\left(\begin{array}{c}{f}_{1m}^{DM}\\ f_{1m}^{BM}\end{array}\right)+{\displaystyle \sum }_{mv}{\mathit{N}}_{1v}^1{\left(g_{NN}^{{\alpha }_{ED}},g_{NN}^{{\alpha }_{EB}}\right)}_{m,v}\left(\begin{array}{c}{f}_{1m}^{DN}\\ f_{1m}^{BN}\end{array}\right)+{\displaystyle \sum }_{mv}{\mathit{L}}_{1v}^1{\left(g_{NL}^{{\alpha }_{ED}},g_{NL}^{{\alpha }_{EB}}\right)}_{m,v}\left(\begin{array}{c}{f}_{1m}^{DN}\\ f_{1m}^{BN}\end{array}\right)$$

(27)

$$\mathit{H}={\displaystyle \sum }_{mv}{\mathit{M}}_{1v}^1{\left(g_{MM}^{{\alpha }_{HD}},g_{MM}^{{\alpha }_{HB}}\right)}_{m,v}\left(\begin{array}{c}{f}_{1m}^{DM}\\ f_{1m}^{BM}\end{array}\right)+{\displaystyle \sum }_{mv}{\mathit{N}}_{1v}^1{\left(g_{NN}^{{\alpha }_{HD}},g_{NN}^{{\alpha }_{HB}}\right)}_{m,v}\left(\begin{array}{c}{f}_{1m}^{DN}\\ f_{1m}^{BN}\end{array}\right)+{\displaystyle \sum }_{mv}{\mathit{L}}_{1v}^1{\left(g_{NL}^{{\alpha }_{HD}},g_{NL}^{{\alpha }_{HB}}\right)}_{m,v}\left(\begin{array}{c}{f}_{1m}^{DN}\\ f_{1m}^{BN}\end{array}\right)$$

(28)

Applying Maxwell’s equation to the E- and H-fields of Equations (20) and (21), we can obtain

$$-\mathit{i}({\mathit{k}}_0/\mathit{k})\mathit{D}={\displaystyle \sum }_{\mathbf{1}\mathit{m},\mathbf{1}\mathit{v}}{\mathit{N}}_{\mathbf{1}\mathit{v}}^{\mathbf{1}}\cdot {\left(g_{MM}^{{\alpha }_{HD}},g_{MM}^{{\alpha }_{HB}}\right)}_{m,v}\cdot \left(\begin{array}{c}{\mathit{f}}_{\mathbf{1}\mathit{m}}^{\mathit{D}\mathit{M}}\\ {\mathit{f}}_{\mathbf{1}\mathit{m}}^{\mathit{B}\mathit{M}}\end{array}\right)+{\displaystyle \sum }_{\mathbf{1}\mathit{m},\mathbf{1}\mathit{v}}{\mathit{M}}_{\mathbf{1}\mathit{v}}^{\mathbf{1}}\cdot {\left(g_{NN}^{{\alpha }_{HD}},g_{NN}^{{\alpha }_{HB}}\right)}_{m,v}\cdot \left(\begin{array}{c}{\mathit{f}}_{\mathbf{1}\mathit{m}}^{\mathit{D}\mathit{N}}\\ {\mathit{f}}_{\mathbf{1}\mathit{m}}^{\mathit{B}\mathit{N}}\end{array}\right)$$

$$\mathit{i}({\mathit{k}}_0/\mathit{k})\mathit{B}={\displaystyle \sum }_{\mathbf{1}\mathit{m},\mathbf{1}\mathit{v}}{\mathit{N}}_{\mathbf{1}\mathit{v}}^{\mathbf{1}}\cdot {\left(g_{MM}^{{\alpha }_{ED}},g_{MM}^{{\alpha }_{EB}}\right)}_{m,v}\cdot \left(\begin{array}{c}{\mathit{f}}_{\mathbf{1}\mathit{m}}^{\mathit{D}\mathit{M}}\\ {\mathit{f}}_{\mathbf{1}\mathit{m}}^{\mathit{B}\mathit{M}}\end{array}\right)+{\displaystyle \sum }_{\mathbf{1}\mathit{m},\mathbf{1}\mathit{v}}{\mathit{M}}_{\mathbf{1}\mathit{v}}^{\mathbf{1}}\cdot {\left(g_{NN}^{{\alpha }_{ED}},g_{NN}^{{\alpha }_{EB}}\right)}_{m,v}\cdot \left(\begin{array}{c}{\mathit{f}}_{\mathbf{1}\mathit{m}}^{\mathit{D}\mathit{N}}\\ {\mathit{f}}_{\mathbf{1}\mathit{m}}^{\mathit{B}\mathit{N}}\end{array}\right)$$

This translates into a system of equations

$${\widehat{G}}_M{F}_{Mq}=i(k_0/k_q)\widehat{U}{F}_{Nq}$$

(29)

$${\widehat{G}}_N{F}_{Nq}=i(k_0/k_q)\widehat{V}{F}_{Mq}$$

(30)

where $F_{Mq}={\left(f_{1-1q}^{DM},f_{1-1q}^{BM},f_{10q}^{DM},f_{10q}^{BM},f_{11q}^{DM},f_{11q}^{BM}\right)}^T, F_{Mq}={\left(f_{1-1q}^{DN},f_{1-1q}^{BN},f_{10q}^{DN},f_{10q}^{BN},f_{11q}^{DN},f_{11q}^{BN}\right)}^T$, and the matrices ${\widehat{G}}_M$ and ${\widehat{G}}_N$ are composed of the coefficients $g_{MM}$ and $g_{NN}$

$${\widehat{G}}_M=\left(\begin{array}{cccccc}{g}_{MM--}^{{\alpha }_{HD}}& g_{MM--}^{{\alpha }_{HB}}& g_{MM0-}^{{\alpha }_{HD}}& g_{MM0-}^{{\alpha }_{HB}}& g_{MM+-}^{{\alpha }_{HD}}& g_{MM+-}^{{\alpha }_{HB}}\\ g_{MM--}^{{\alpha }_{ED}}& g_{MM--}^{{\alpha }_{EB}}& g_{MM0-}^{{\alpha }_{ED}}& g_{MM0-}^{{\alpha }_{EB}}& g_{MM+-}^{{\alpha }_{ED}}& g_{MM+-}^{{\alpha }_{EB}}\\ g_{MM-0}^{{\alpha }_{HD}}& g_{MM-0}^{{\alpha }_{HB}}& g_{MM00}^{{\alpha }_{HD}}& g_{MM00}^{{\alpha }_{HB}}& g_{MM+0}^{{\alpha }_{HD}}& g_{MM+0}^{{\alpha }_{HB}}\\ g_{MM-0}^{{\alpha }_{ED}}& g_{MM-0}^{{\alpha }_{EB}}& g_{MM00}^{{\alpha }_{ED}}& g_{MM00}^{{\alpha }_{EB}}& g_{MM+0}^{{\alpha }_{ED}}& g_{MM+0}^{{\alpha }_{EB}}\\ g_{MM-+}^{{\alpha }_{HD}}& g_{MM-+}^{{\alpha }_{HB}}& g_{MM0+}^{{\alpha }_{HD}}& g_{MM0+}^{{\alpha }_{HB}}& g_{MM++}^{{\alpha }_{HD}}& g_{MM++}^{{\alpha }_{HB}}\\ g_{MM-+}^{{\alpha }_{ED}}& g_{MM-+}^{{\alpha }_{EB}}& g_{MM0+}^{{\alpha }_{ED}}& g_{MM0+}^{{\alpha }_{EB}}& g_{MM++}^{{\alpha }_{ED}}& g_{MM++}^{{\alpha }_{EB}}\end{array}\right)$$

$${\widehat{G}}_N=\left(\begin{array}{cccccc}{g}_{NN--}^{{\alpha }_{ED}}& g_{NN--}^{{\alpha }_{EB}}& g_{NN0-}^{{\alpha }_{ED}}& g_{NN0-}^{{\alpha }_{EB}}& g_{NN+-}^{{\alpha }_{ED}}& g_{NN+-}^{{\alpha }_{EB}}\\ g_{NN--}^{{\alpha }_{HD}}& g_{NN--}^{{\alpha }_{HB}}& g_{NN0-}^{{\alpha }_{HD}}& g_{NN0-}^{{\alpha }_{HB}}& g_{NN+-}^{{\alpha }_{HD}}& g_{NN+-}^{{\alpha }_{HB}}\\ g_{NN-0}^{{\alpha }_{ED}}& g_{NN-0}^{{\alpha }_{EB}}& g_{NN00}^{{\alpha }_{ED}}& g_{NN00}^{{\alpha }_{EB}}& g_{NN+0}^{{\alpha }_{ED}}& g_{NN+0}^{{\alpha }_{EB}}\\ g_{NN-0}^{{\alpha }_{HD}}& g_{NN-0}^{{\alpha }_{HB}}& g_{NN00}^{{\alpha }_{HD}}& g_{NN00}^{{\alpha }_{HB}}& g_{NN+0}^{{\alpha }_{HD}}& g_{NN+0}^{{\alpha }_{HB}}\\ g_{NN-+}^{{\alpha }_{ED}}& g_{NN-+}^{{\alpha }_{EB}}& g_{NN0+}^{{\alpha }_{ED}}& g_{NN0+}^{{\alpha }_{EB}}& g_{NN++}^{{\alpha }_{ED}}& g_{NN++}^{{\alpha }_{EB}}\\ g_{NN-+}^{{\alpha }_{HD}}& g_{NN-+}^{{\alpha }_{HB}}& g_{NN0+}^{{\alpha }_{HD}}& g_{NN0+}^{{\alpha }_{HB}}& g_{NN++}^{{\alpha }_{HD}}& g_{NN++}^{{\alpha }_{HB}}\end{array}\right)$$

while $\widehat{U}=\left(\begin{array}{cccccc}-1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& -1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)$, $\widehat{V}=\left(\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ -1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& -1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& -1& 0\end{array}\right)$.

Excluding $F_{Mq}=i(k_q/k_0)\widehat{V}{\widehat{G}}_N{F}_{Nq}$ from the system Equations (29) and (30), we obtained an eigenproblem for $F_{Nq}$, which corresponded to six bianisotropic orbitals in the Rayleigh limit

$$\left(\widehat{U}{\widehat{G}}_M\widehat{V}{\widehat{G}}_N\right)F_{Nq}={\left(k_0/k_q\right)}^2{F}_{Nq}$$

(31)

The boundary conditions of Equations (17)–(20) turn into

$$j_1\left(k_{0}R\right) \left(pq\right)+h_1\left(k_{0}R\right)\left(ab\right)={\displaystyle \sum }_q{A}_q{j}_1\left(k_{q}R\right)diag\left\{i,1,i,1,i,1\right\}{\widehat{G}}_M{F}_{Mq}$$

(32)

$$\begin{array}{c}\frac{\partial }{\partial r}{\left[r j_1\left(k_{0}r\right)\right]}_{r=R} \left(pq\right)+\frac{\partial }{\partial r}{\left[r h_1\left(k_{0}r\right)\right]}_{r=R}\left(ab\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \sum }_q{A}_q\left(\frac{k_0}{k_q}\right)diag\left\{1,i,1,i,1,i\right\}\times \left(\frac{\partial }{\partial r}{\left[r j_1\left(k_{q}r\right)\right]}_{r=R}{\widehat{G}}_N{F}_{Nq}+j_1\left(k_{q}R\right){\widehat{G}}_L{F}_{Nq}\right)\hfill \end{array}$$

(33)

where ${\widehat{G}}_L$ is composed of the coefficients $g_{NL}$, similar to the matrix ${\widehat{G}}_N$ provided above, $\left(pq\right)={\left(p_{1-1},q_{1-1},p_{10},q_{10},p_{11},q_{11}\right)}^T$, and $\left(ab\right)={\left(a_{1-1},b_{1-1},a_{10},b_{10},a_{11},b_{11}\right)}^T$.

Taking the limit of $k_{q}R,k_{0}R\to 0$ in the spherical Bessel functions, and excluding the scattered amplitudes $\left(ab\right)$, and substituting Equation (29) for ${\widehat{G}}_M{F}_{Mq}$, we obtained

$$3\left(pq\right)=diag\left\{1,i,1,i,1,i\right\}\left(\widehat{1}+2{\widehat{G}}_N+{\widehat{G}}_L\right){\widehat{F}}_N\overrightarrow{A}$$

(34)

Equation (34) can be re-written using $\left(\begin{array}{c}\mathit{D}\\ \mathit{B}\end{array}\right)=\frac{1}{\sqrt{3\pi }}{\widehat{T}}^{-1}{\widehat{F}}_N\overrightarrow{A}$, $\left(pq\right)=\sqrt{3\pi }\widehat{T}\left(\begin{array}{c}\mathbf{ϵ}\\ \mathit{h}\end{array}\right)$, and the matrix $\widehat{T}$ given by $\widehat{T}=\left(\begin{array}{cccccc}-i& +1& 0& 0& 0& 0\\ 0& 0& 0& 1& +i& 0\\ 0& 0& -i\sqrt{2}& 0& 0& 0\\ 0& 0& 0& 0& 0& \sqrt{2}\\ i& +1& 0& 0& 0& 0\\ 0& 0& 0& -1& +i& 0\end{array}\right)$.

In accordance with Equation (34), $\left(\widehat{1}+2{\widehat{G}}_N+{\widehat{G}}_L\right)=\widehat{T}\left(\widehat{1}+2{\widehat{M}}^{-1}\right){\widehat{T}}^{-1}$, or $3{\widehat{T}}^{-1}{\left(\widehat{1}+2{\widehat{G}}_N+{\widehat{G}}_L\right)}^{-1}\widehat{T}=3\widehat{M}{\left(\widehat{M}+2\widehat{1}\right)}^{-1}$. This resulted in

$$\left(\begin{array}{c}\mathit{D}\\ \mathit{B}\end{array}\right)=3\widehat{M}{\left(\widehat{M}+2\widehat{1}\right)}^{-1}\left(\begin{array}{c}\mathbf{ϵ}\\ \mathit{h}\end{array}\right)$$

(35)

Considering that $\left(\begin{array}{c}\mathit{D}\\ \mathit{B}\end{array}\right)=\left(\begin{array}{c}\mathit{E}\\ \mathit{H}\end{array}\right)+4\pi \left(\begin{array}{c}\mathit{P}\\ \mathit{M}\end{array}\right)={\widehat{M}}^{-1}\left(\begin{array}{c}\mathit{D}\\ \mathit{B}\end{array}\right)+4\pi \left(\begin{array}{c}\mathit{P}\\ \mathit{M}\end{array}\right)$, we found $\left(\begin{array}{c}\mathit{D}\\ \mathit{B}\end{array}\right)=4\pi {\left(\widehat{1}-{\widehat{M}}^{-1}\right)}^{-1}\left(\begin{array}{c}\mathit{P}\\ \mathit{M}\end{array}\right)$,

$$\left(\begin{array}{c}\mathit{P}\\ \mathit{M}\end{array}\right)=\frac{3}{4\pi }\left(\widehat{M}-\widehat{1}\right){\left(\widehat{M}+2\widehat{1}\right)}^{-1}\left(\begin{array}{c}\mathbf{ϵ}\\ \mathit{h}\end{array}\right)$$

(36)

From this, we obtained the polarizability $\widehat{\alpha }$ of the sphere of volume $V$ in the Rayleigh approximation

$$\left(\begin{array}{c}\mathit{p}\\ \mathit{m}\end{array}\right)=V\left(\begin{array}{c}\mathit{P}\\ \mathit{M}\end{array}\right)=\widehat{\alpha }\left(\begin{array}{c}\mathbf{ϵ}\\ \mathit{h}\end{array}\right)=\frac{3V}{4\pi }\left(\widehat{M}-\widehat{1}\right){\left(\widehat{M}+2\widehat{1}\right)}^{-1}\left(\begin{array}{c}\mathbf{ϵ}\\ \mathit{h}\end{array}\right)$$

This agrees with the polarizability of the bianisotropic spheres obtained previously in electrostatic approximations [67,68,69].

3. Discussion and Conclusions

Mie’s theory of electromagnetic scattering by spheres is a very important part of electromagnetism. In the existing literature, it has been extended from the original results of Mie to bi-isotropic spheres by Bohren [40], and, recently, to anisotropic spheres by Lin et al. and Li et al. [43,44,45]. Nevertheless, there is no existing theory of scattering by bianisotropic spheres, and the methods used to find the scattering for isotropic or anisotropic spheres are not applicable to bianisotropic spheres.

The extension of our approach to other geometries could be very promising. For example, recently much attention has been paid to the Mie resonances in nanocylinder systems [70,71,72], and we believe that this opens a broad avenue for application.

To conclude, in this study, we introduced the bianisotropic orbitals and presented a theory of the scattering by bianisotropic spheres with arbitrary effective media parameters and sizes. In the Rayleigh limit, we obtained the results known from the electrostatic approximation approach.