Magnetoquasistatic Resonances of Small Dielectric Objects

TL;DR

This paper explores magnetoquasistatic resonances in small high-permittivity dielectric objects, revealing their physical properties and providing a simplified analysis method for designing high-index resonators.

Contribution

It introduces a new understanding of magnetoquasistatic resonances in arbitrary-shaped dielectrics using Maxwell's equations approximation.

Findings

Resonances occur at specific permittivities and frequencies linked to eigenvalues of a magnetostatic operator.

Resonances involve interplay between dielectric polarization energy and magnetic field energy.

The approach simplifies the analysis and design of high-permittivity dielectric resonators.

Abstract

A small dielectric object with positive permittivity may resonate when the free-space wavelength is large in comparison with the object dimensions if the permittivity is sufficiently high. We show that these resonances are described by the magnetoquasistatic approximation of the Maxwell's equations in which the normal component of the displacement current density field vanishes on the surface of the particle. They are associated to values of permittivities and frequencies for which source-free quasistatic magnetic fields exist, which are connected to the eigenvalues of a magnetostatic integral operator. We present the general physical properties of magnetoquasistatic resonances in dielectrics with arbitrary shape. They arise from the interplay between the polarization energy stored in the dielectric and the energy stored in the magnetic field. Our findings improve the understanding of resonances in high-permittivity dielectric objects and provide a powerful tool that greatly simplifies the analysis and design of high index resonators.

Tables

Table 1. Table 1: Eigenvalues of a sphere y n subscript 𝑦 𝑛 y_{n} compared with the poles of the Mie coefficients r n ​ l subscript 𝑟 𝑛 𝑙 r_{nl} in the limit x → 0 → 𝑥 0 x\rightarrow 0 . Relative error ϵ italic-ϵ \epsilon .

$\#$	1-3	4-11	12-23	24-26	27-32	33-36	37-38	39-42	43-50
$y_n$	3.16	4.53	5.82	6.38	7.07	7.08	7.09	7.10	7.87
$r_{nl}$	3.14	4.49	5.76	6.28	6.99	6.99	6.99	6.99	7.73
$ϵ[\%]$	0.64	0.89	1.04	1.59	1.14	1.29	1.43	1.57	1.81

Table 2. Table 2: Magnetostatic eigenvalues y n subscript 𝑦 𝑛 y_{n} for an isolated finite-size cylinder with height equal to the radius. y ~ n subscript ~ 𝑦 𝑛 \tilde{y}_{n} are the corresponding values evaluated by the empirical approximated formulas Mongia:94 (if available).

$\#$	1	2-3	4-5	6	7-8	9
name	${\text{TE}}_{01\delta }$	${\text{HEM}}_{11\delta }$	${\text{HEM}}_{12\delta }$	${\text{TM}}_{01\delta }$	${\text{HEM}}_{21\delta }$	${\text{TE}}_{011+\delta }$
$y_n$	3.26	4.05	4.52	4.96	5.02	5.30
${\tilde{y}}_n$	3.22	4.11	-	4.95	-	5.28

Table 3. Table 3: Magnetoquasistatic eigenvalues y n subscript 𝑦 𝑛 y_{n} of an isolated triangular prism with edge L 𝐿 L and height H 𝐻 H , with L = 2 ​ H 𝐿 2 𝐻 L=2H . We assumed l c = L subscript 𝑙 𝑐 𝐿 l_{c}=L .

$\#$	1	2-3	4	5-6	7-8	9	10-11
$y_n$	4.52	4.69	5.76	6.21	6.26	6.48	7.41

Table 4. Table 4: Magnetoquasistatic eigenvalues y n subscript 𝑦 𝑛 y_{n} of an ring with minor radius r 𝑟 r and major radius R = 3 ​ r 𝑅 3 𝑟 R=3r , We assumed l c = 3 ​ R subscript 𝑙 𝑐 3 𝑅 l_{c}=3R .

$\#$	1	2-3	4-5	6-7	8-9	10-11	12
$\#$	${\text{T}}_{010}$	${\text{HTP}}_{211}^{\updownarrow }$	${\text{HTP}}_{111}^{\updownarrow }$	${\text{HTP}}_{111}^{\leftrightarrow }$	${\text{HTP}}_{211}^{\leftrightarrow }$	${\text{HTP}}_{311}^{\updownarrow }$	${\text{P}}_{010}$
$y_n$	3.484	7.079	7.095	7.122	7.154	7.292	7.301

Table 5. Table 5: Details on the meshes used in the numerical calculation of the eigenvalues.

object	$N_{node}$	$N_{elem}$	$N_{edge}$
sphere	6527	6048	11665
finite-size cylinder	6060	5148	9433
triangular prism	2520	2025	3592
ring	29056	26112	49409

Equations

$$\nabla \times \mathbf{A}={\mu }_0\mathbf{H},$$

$$\nabla \times \mathbf{H}={\mathbf{J}}^d,$$

$${\mathbf{J}}^d=\frac{{\beta }^2}{{\mu }_0}\mathbf{A}{\Pi }_{\Omega },$$

$$\mathbf{A}\left(\tilde{\mathbf{r}}\right)={\mu }_0{l}_c^2{\mathcal{L}}_m\left\{{\mathbf{J}}^d\right\}\left(\tilde{\mathbf{r}}\right),$$

$${\mathcal{L}}_m\left\{{\mathbf{J}}^d\right\}\left(\tilde{\mathbf{r}}\right)={\iiint }_{\tilde{\Omega }}{\mathbf{J}}^d\left({\tilde{\mathbf{r}}}^{\prime }\right)g_0\left(\tilde{\mathbf{r}}-{\tilde{\mathbf{r}}}^{\prime }\right)\text{d}\tilde{V}\forall \tilde{\mathbf{r}}\in \tilde{\Omega },$$

$${\mathbf{J}}^d\left(\tilde{\mathbf{r}}\right)=y^2{\mathcal{L}}_m\left\{{\mathbf{J}}^d\right\}\left(\tilde{\mathbf{r}}\right)\forall \tilde{\mathbf{r}}\in \tilde{\Omega },$$

$$y_n=\frac{{\Vert {\nabla }_{\tilde{\mathbf{r}}}\times {\mathbf{A}}_n\Vert }_{{\mathbb{R}}^3}}{{\Vert {\mathbf{A}}_n\Vert }_{\tilde{\Omega }}}$$

$${\omega }_n=\frac{c_0}{l_c\sqrt{{\epsilon }_R}}{y}_n.$$

$$y_n\ge \frac{\sqrt[4]{3}}{2\sqrt{\pi }}\frac{1}{\tilde{a}}\approx \frac{0.37}{\tilde{a}}\forall n\in \mathbb{N}$$

$${\int }_{\tilde{\Omega }}\frac{1}{{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}^2}d{\tilde{V}}^{\prime }$$

$${\int }_{\tilde{\Omega }}\left|{\mathbf{J}}^d\right|d{\tilde{V}}^{\prime }$$

$${\epsilon }_{R,n}={\left(\frac{y_n}{x}\right)}^2.$$

$$g_S\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)=\frac{1}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}-\frac{{\epsilon }_S-1}{{\epsilon }_S+1}\frac{1}{\left|\mathbf{r}-{\mathbf{r}}^{\prime \prime }\right|}$$

$$\begin{gathered} \frac{\mathbf{J}\left({\mathbf{r}}^{\prime }\right)}{i\omega {\epsilon }_0\chi }=-\frac{1}{i\omega {\epsilon }_0}\nabla \iint g\left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)\mathbf{J}\left({\mathbf{r}}^{\prime }\right)\cdot \hat{\mathbf{n}}\left({\mathbf{r}}^{\prime }\right)dS \\ -i\omega {\mu }_0{\iiint }_{V}g\left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)\mathbf{J}\left({\mathbf{r}}^{\prime }\right)dV+{\mathbf{E}}_{inc}\left(\mathbf{r}\right)\forall \mathbf{r}\in \Omega , \end{gathered}$$

$$\frac{\mathbf{J}\left(\tilde{\mathbf{r}}\right)}{\chi }-\mathcal{L}\left\{\mathbf{J}\right\}\left(\tilde{\mathbf{r}}\right)=i\omega {\epsilon }_0{\mathbf{E}}_{inc}\left(\tilde{\mathbf{r}}\right)\forall \tilde{\mathbf{r}}\in \tilde{\Omega }$$

$$\begin{gathered} \mathcal{L}\left\{\mathbf{W}\right\}\left(\tilde{\mathbf{r}}\right)=-\tilde{\nabla }\iint {\gamma }^{\perp }{W}_n\left({\tilde{\mathbf{r}}}^{\prime }\right)d{\tilde{S}}^{\prime } \\ +x^2{\iiint }_{\tilde{\Omega }}{\gamma }^{\perp }\mathbf{W}\left({\tilde{\mathbf{r}}}^{\prime }\right)d{\tilde{V}}^{\prime } \end{gathered}$$

$$\begin{gathered} \mathcal{L}\left\{\mathbf{W}\right\}\left(\tilde{\mathbf{r}}\right)=-\tilde{\nabla }\iint g_0\left(\tilde{\mathbf{r}}-{\tilde{\mathbf{r}}}^{\prime }\right)\left(1-ix\left|\tilde{\mathbf{r}}-{\tilde{\mathbf{r}}}^{\prime }\right|+\frac{x^2}{2}{\left|\tilde{\mathbf{r}}-{\tilde{\mathbf{r}}}^{\prime }\right|}^2\right)W_n\left({\tilde{\mathbf{r}}}^{\prime }\right)d{\tilde{S}}^{\prime } \\ +x^2{\iiint }_{\tilde{\Omega }}{g}_0\left(\tilde{\mathbf{r}}-{\tilde{\mathbf{r}}}^{\prime }\right)\mathbf{W}\left({\tilde{\mathbf{r}}}^{\prime }\right)d{\tilde{V}}^{\prime }+O(x^3) \end{gathered}$$

$$g_0\left(\mathbf{r}\right)=\frac{1}{4\pi r}$$

$${\mathcal{L}}_e\left\{{\mathbf{W}}_h^{\parallel }\right\}=\frac{1}{{\gamma }_h^{\parallel }}{\mathbf{W}}_h^{\parallel },$$

$${\mathcal{L}}_e\left\{\mathbf{W}\right\}=-\tilde{\nabla }\iint g_0\left(\tilde{\mathbf{r}}-{\tilde{\mathbf{r}}}^{\prime }\right)W_n\left({\tilde{\mathbf{r}}}^{\prime }\right)d{\tilde{S}}^{\prime }.$$

$$⟨\mathbf{A},\mathbf{B}⟩={\iiint }_{\Omega }\mathbf{A}\cdot \mathbf{B}dV.$$

$${\mathcal{L}}_m\left\{{\mathbf{W}}_h^{\perp }\right\}=\frac{1}{{\gamma }_h^{\parallel }}{\mathbf{W}}_h^{\perp }$$

$${\mathcal{L}}_m\left\{{\mathbf{W}}_h^{\perp }\right\}\left(\tilde{\mathbf{r}}\right)={\iiint }_{\tilde{\Omega }}{g}_0\left(\tilde{\mathbf{r}}-{\tilde{\mathbf{r}}}^{\prime }\right){\mathbf{W}}_h^{\perp }\left({\tilde{\mathbf{r}}}^{\prime }\right)d{\tilde{V}}^{\prime }.$$

$$\mathbf{J}\left(\tilde{\mathbf{r}}\right)=\sum _h\left(I_h^{\parallel }{\mathbf{W}}_h^{\parallel }\left(\tilde{\mathbf{r}}\right)+I_h^{\perp }{\mathbf{W}}_h^{\perp }\left(\tilde{\mathbf{r}}\right)\right).$$

$$I_h^{\parallel }$$

$$I_h^{\perp }$$

$$\tilde{\nabla }\times \tilde{\nabla }\times {\mathbf{A}}^{\perp }={\gamma }^{\perp }{\mathbf{A}}^{\perp }\tilde{\Pi }\left({\tilde{\mathbf{r}}}^{\prime }\right)$$

$$\frac{{\omega }_h^2{l}_c^2}{c_0^2}\chi ={\gamma }_h^{\perp }.$$

$${\gamma }_h^{\perp }=\frac{{\Vert \tilde{\nabla }\times {\mathbf{A}}_h^{\perp }\Vert }_{{\mathbb{R}}^3}^2}{{\Vert {\mathbf{A}}_h^{\perp }\Vert }_{\tilde{\Omega }}^2}$$

$$\frac{1}{2{\mu }_0}{\Vert \tilde{\nabla }\times {\mathbf{A}}_h^{\perp }\Vert }_{{\mathbb{R}}^3}^2=\frac{{\epsilon }_0\chi }{2}{\Vert {\omega }_h{\mathbf{A}}_h^{\perp }\Vert }_{\Omega }^2.$$

$${\mathcal{J}}_L=\left\{\mathbf{w}\in H\left(\mathrm{div},\Omega \right)|\nabla \cdot \mathbf{w}=0\text{ in }\Omega ,\mathbf{w}\cdot \hat{\mathbf{n}}=0\text{ on }\partial \Omega \right\}$$

$${\mathbf{J}}^d=\sum _{k=1}^{N_L}{I}_k^L{\mathbf{w}}_k^{\text{L}}$$

Full text

It is well established that small metal objects with negative permittivity may resonate when the free-space wavelength is large in comparison with their dimensions Mayergoyz:03 ; mayergoyz2005electrostatic . These resonances can be predicted by the electroquasistatic approximation of the Maxwell-equations, and they are associated to the values of permittivity for which source-free electrostatic fields exist Mayergoyz:03 ; mayergoyz2005electrostatic .

Small dielectric objects with positive permittivity may also resonate when the free-space wavelength is large in comparison with their dimensions, providing that their permittivity is sufficiently high VanBladel:75 ; hulst1981light ; Videen . At microwaves, low-loss dielectric materials with relative permittivities $\varepsilon_{R}$ of the order of $\approx 100$ are routinely used for various applications including resonators and filters cava2001dielectric ; mongia1994low ; Reaney06 , while the relative permittivity of certain titanates schlicke1953quasi ; PhysRev.135.A748 ; sethares1966design can reach values higher than 1000.

In the visible and NIR spectral ranges, resonances in nanoscale high-permittivity objects, such as AlGaAs, Si, Ge nanoparticles, have been observed experimentally, e.g. evlyukhin2012demonstration ; kuznetsov2012magnetic ; kuznetsov2016optically , and have been exploited for several applications, e.g. Smirnova:16 ; Schmidt:12 ; PhysRevB.85.245432 ; Feng:16 ; doi:10.1021/acsphotonics.7b00509 ; baranov2017modifying ; staude2013tailoring . In the nano-optics community, these resonances are known as “Mie resonances” (e.g. kuznetsov2016optically ) and they are described in framework of the full-Maxwell equations (Mie theory Bohren1998 , quasi-normal modes Lalanne:18 , characteristic modes garbacz1965modal , material-independent-modes Forestiere:16 ; Pascale2019 , etc.).

However, resonances in dielectric objects have a long history which begins at the dawn of the twentieth century with the work on Debye on natural resonant frequencies of free dielectric spheres debye1909lichtdruck . In 1939, R. D. Richtmyer showed that suitably shaped objects made of a dielectric material can function as electrical resonators at high frequency richtmyer1939dielectric . Since the 1960s, they have been used as high-Q elements for microwave filters and oscillators designs kajfez1998dielectric , and following the work of Long et. al Long:83 also as antennas. Dielectric resonators has been traditionally analyzed by using perfect magnetic wall boundary conditions (PMW) schlicke1953quasi ; sethares1966design ; Mongia:94 . However, because electromagnetic fields do exist beyond the geometrical boundary of the cavity this condition is unable to accurately predict resonances yee1965natural ; Guillon:77 ; kajfez1998dielectric . Many ad-hoc corrections to the PMW conditions have been proposed, including the Cohn model cohn1968microwave , where an idealized waveguide with PMW walls is considered, and the Itoh-Rudokas model Itoh:77 . Instead, Van Bladel investigated these resonances without ad-hoc assumptions using an asymptotic expansion of the Maxwell’s equation in differential form in terms of the inverse of the index of refraction VanBladel:75 . Glisson et al. Glisson:83 ; Kajfez:84 obtained the resonant frequencies of rotationally symmetric dielectric bodies by searching the frequencies at which the determinant of impedance matrix is zero; they assembled the impedance matrix by discretizing a surface integral formulation of the full-Maxwell equation.

In this manuscript, we show by using an integral formulation that the resonances in high-permittivity small dielectric objects can be predicted by the magnetoquasistatic approximation haus1989electromagnetic of the Maxwell’s equations, in which the normal component of the displacement current density field vanishes on the surface of the object. These resonances are associated to values of permittivities and frequencies for which source-free quasistatic magnetic fields exist, and they are in one-to-one correspondence to the eigenvalues of the magnetostatic integral operator relating the vector potential inside the object and the displacement current density induced in the object itself. By studying this operator, we derive the general physical properties of magnetoquasistatic resonances in small high-permittivity objects of arbitrary shape. They arise from the interplay between the polarization energy stored in the dielectric and the energy stored in the magnetic field. Their resonance frequencies constitute an infinite countable set accumulating at $+\infty$, and, for any given shape of the object, they are inversely proportional to its linear dimension and to the square root of its relative permittivity. The eigenmodes corresponding to different eigenvalues are orthogonal in the usual sense. We also introduce an a-priori lower bound for the minimum resonance frequency. From a numerical perspective, this approximated integral formulation gives, in its range of validity, a great simplification with respect to full-wave approaches and with respect to asymptotic differential formulations.

I Resonances in a Magnetoquasistatic Approximation

Let us consider a homogeneous and isotropic dielectric object with a bounded arbitrary shape $\Omega$ and relative permittivity $\varepsilon_{R}$. We define its characteristic size $l_{c}$ as a chosen linear dimension of the object, and its size parameter $x=2\pi l_{c}/\lambda$, where $\lambda$ is the vacuum wavelength. We look for source-free solutions of the Maxwell’s equations in high-permittivity dielectric objects, in the limit $x\rightarrow 0$ (small object). Under these conditions, the electromagnetic field is primarily determined by the displacement current density field ${\bf J}^{d}$ induced inside the object itself VanBladel:75 . In the appendix A, we show that this assumption is well-founded. Therefore, we look for the values of the parameter $\beta=\left(\omega/c_{0}\right)\sqrt{\varepsilon_{R}}$ for which there exists a non-trivial solution of the source-free magnetoquasistatic problem haus1989electromagnetic

[TABLE]

with the constitutive relation

[TABLE]

where $\Pi_{\Omega}$ is the characteristic function on the set $\Omega$, i.e. $\Pi_{\Omega}\left(\mathbf{r}\right)=1$ for $\mathbf{r}\in\Omega$, [math] otherwise, $\omega$ is the angular frequency, and $c_{0}=1/\sqrt{\varepsilon_{0}\mu_{0}}$ is the light velocity in vacuum. The magnetoquasistatic vector potential $\mathbf{A}$ satisfies the Coulomb gauge in $\Omega$ and $\mathbb{R}^{3}\backslash\Omega$; $\mathbf{A}$ and the quasistatic magnetic field $\mathbf{H}$ are regular at infinity. Eq. 2 disregards the effects of the displacement current density field in vacuum. The continuity of the tangential components of $\mathbf{A}$ and $\mathbf{H}$ imply, respectively, the continuity of the normal components of $\mathbf{H}$ and ${\bf J}^{d}$ across the boundary $\partial\Omega$ of $\Omega$. Since the normal component of the current density field ${\bf J}^{d}$ at the boundary $\partial\Omega$ is equal to zero, the current density field ${\bf J}^{d}$ is div-free everywhere in $\mathbb{R}^{3}$; instead, the normal component of the vector potential at $\partial\Omega$ may be discontinuous. The fact that ${\bf J}^{d}$ has a vanishing normal component on $\partial\Omega$ implies that also the normal component of the polarization current density field is zero.

It is convenient to scale the spatial coordinates by the characteristic size of the object $l_{c}$, $\mathbf{r}=l_{c}\tilde{\mathbf{r}}$. Thus, we denote with $\tilde{\Omega}$ the scaled domain $\Omega$. Then, problem 1a,1b,2 is solved by expressing the vector potential $\mathbf{A}$ in terms of the current density ${\bf J}^{d}$ as:

[TABLE]

where we have introduced the magnetostatic integral operator

[TABLE]

and $g_{0}\left(\mathbf{r}\right)=1/(4\pi r)$ is the static Green function in vacuum. In 4 there is the static Green function because we are neglecting the displacement current density in vacuum. By combining Eqs. 2 and 3, we obtain the linear eigenvalue problem

[TABLE]

where $y=l_{c}\left(\omega/c_{0}\right)\sqrt{\varepsilon_{R}}=x\sqrt{\varepsilon_{R}}$. Equation 5 holds in the weak form in the functional space constituted by the functions which are div-free within $\tilde{\Omega}$ and having zero normal component on $\partial\tilde{\Omega}$, and equipped with the inner product $\langle\mathbf{w},\mathbf{v}\rangle_{\tilde{V}}=\iiint_{\tilde{V}}\mathbf{w}^{*}\cdot\mathbf{v}\,\mbox{d}\tilde{V}.$

The operator $\mathcal{L}_{m}$ is compact, positive-definite, and self-adjoint. As a consequence: Equation 5 admits a countable set of eigenvalues $\left\{y_{n}^{2}\right\}_{n\in\mathbb{N}}$ and eigenmodes ${\bf J}^{d}_{n}$; the eigenvalues $y_{n}^{2}$ are real and positive and accumulate at infinity. The eigenmodes corresponding to different eigenvalues are orthogonal in the usual sense; and they constitute a complete basis of the considered functional space. Furthermore, the eigenvalue $y_{n}$ is proportional to the magnetic energy of the corresponding eigenmode:

[TABLE]

$\mathbf{A}_{n}$ is the magnetic vector potential generated by ${\bf J}^{d}_{n}$ in the whole space, and $\left\|{\bf v}\right\|_{V}^{2}=\langle{\bf v},{\bf v}\rangle_{V}$.

The resonance angular frequencies $\omega_{n}$ are given by

[TABLE]

The mathematical structure of the integral operator 4 does not depend on the linear characteristic dimension of the dielectric object $l_{c}$, namely it is scale invariant. This fact combined with Eq. 7 leads to the important property of the magnetoquasistatic resonances: for any given shape of the object the resonance frequencies are always inversely proportional to both $l_{c}$ and $\sqrt{\varepsilon_{R}}$. Furthermore, the resonance frequencies accumulate at infinity. By contrast, we recall that in small metal particles the electrostatic resonance frequencies accumulate at finite frequencies Mayergoyz05 (for arbitrarily shaped particles with a simple Drude dispersion they accumulate at $\omega_{p}/\sqrt{2}$ where $\omega_{p}$ is the plasma frequency of the metal.)

Moreover, the following bound on the eigenvalues hold:

[TABLE]

where $\tilde{a}$ is the radius of a sphere $\tilde{B}_{a}$ having the same volume of $\tilde{\Omega}$. The inequality 8 is obtained by multiplying both members of Eq. 5 by ${\bf J}^{d}$, integrating over $\tilde{\Omega}$ and using the Cauchy-Schwarz inequality and the inequality jeffreys1999methods

[TABLE]

where $\tilde{B}_{a}$ is centered in ${\bf r}$. Let us apply the bound 8 to the case of an oblate spheroid: by keeping the major semi-axis fixed, and decreasing the minor one, the volume of the particle decreases and so does $\tilde{a}$: this leads to an increase of lower bound for the eigenvalues (and for the magnetostatic energies) of all the modes. This is in contrast with what it is observed for a metal spheroid, where the same decrease of the minor semi-axis implies a decrease of the electrostatic energy of the fundamental electrostatic (plasmon) mode maier2007plasmonics .

In a “material resonance picture” Forestiere:16 ; Forestiere:19 , once the operating frequency is assigned, the resonance permittivities are given by:

[TABLE]

Since $\mathcal{L}_{m}$ is positive-definite, source-free magnetoquasistatic field may exist only for positive permittivities. Moreover, for any given shape, they $\varepsilon_{R,n}$ are inversely proportional to $x^{2}$, they also constitute accumulate at $\infty$. It is now apparent the difference between the magnetoquasistatic resonances and the plasmon resonances, which only exist for negative permittivities, are size-independent, and accumulate at $-1$ mayergoyz2005electrostatic .

In the Appendix A, we show that in the limit $x\rightarrow 0$ the Maxwell equations admit two orthogonal sets of current modes. The modes of the first set are div-free and curl-free within the object, but have a non-vanishing normal component on the object surface. The modes of the second set are div-free within the object, have a vanishing normal component on the object surface and non-zero curl. The latter set is solution of the eigenvalue problem 5 and the corresponding resonance frequencies $\omega_{n}$ are given by 7.

In a magnetoquasistatic resonance the energy oscillates back and forth between the polarization energy stored in the dielectric and the energy stored in the magnetic field, as shown in Appendix B.

The eigenvalue problem 5 can be numerically solved through a finite element approach, briefly outlined in Appendix C, where, unlike differential formulations, only the spatial domain $\Omega$ is discretized and the radiation condition are automatically satisfied. This approach only requires the calculation of the eigenvalues of a real symmetric matrices, for which efficient routines exists van1983matrix ; laug .

The magnetostatic formulation can be easily extended to the scenario where the object is standing on a substrate with relative permittivity $\varepsilon_{S}$, which is very relevant for the applications. Since the thickness of the substrate is typically much larger than the dimensions of the object, we can safely assume it to be semi-infinite. Therefore, by using the method of images jackson1999classical , the resonance are still associated to the eigenvalues of the operator 4 provided that the Green function $g_{0}$ in Eq. 4 is replaced with

[TABLE]

where $\mathbf{r}^{\prime\prime}$ is the mirror image of $\mathbf{r}^{\prime}$ with respect to the substrate plane.

II Results and Discussions

Sphere. To benchmark the magnetoquasistatic approximation, we now consider a dielectric sphere for which there exists a full-wave analytic solution, i.e. the Mie theory Bohren1998 . Van De Hulst gives hulst1981light the resonant conditions of a high-permittivity small sphere by finding the poles of the Mie coefficients in the limit $x\rightarrow 0$. The resonances occur at $y_{n}=r_{n,l}$ for the TM multipoles and at $y_{n}=r_{n-1,l}$ for the TE multipoles, for any $n,l\in\mathbb{N}$, where $r_{n,l}$ denotes the $l$-th zero of the spherical Bessel function $j_{n}$. In Fig. 1 (a) we compare the first 100 magnetoquasistatic eigenvalues of the operator $\mathcal{L}_{m}$ with the poles of the Mie coefficients, while in Table 1 we show the corresponding values and the numerical error. The details on the hexahedral mesh used for the computation are given in Appendix C. In particular, we note that for the first 50 eigenvalues the numerical error is below $2\%$. The magnetoquasistatic model correctly predicts the peculiar degeneracies of the TE and TM resonances: the $l^{th}$ resonance of the $n^{th}$ TE multipole is degenerate with the $l^{th}$ resonance of the $\left(n-1\right)^{th}$ TM multipole: namely the eigenvalue $r_{0,l}$ has a degeneracy 3, while $r_{n,l}$ with $n\geq 1$ has a degeneracy $4(n+1)$.

The magnetoquasistatic modes are ordered according to their magnetic energy, which does not necessarily follow their multipolar order. The first eigenvalue $y=r_{01}$ is associated to three degenerate TE displacement current density modes, $\mathbf{M}_{pm1}^{\left(1\right)}\left(r_{01}\tilde{\mathbf{r}}\right)$ with $p=e,o$ and $m=0,1$, which correspond to magnetic dipoles oriented along the three coordinate axis. As an example, we show $\mathbf{M}_{o11}^{\left(1\right)}\left(r_{01}\tilde{\mathbf{r}}\right)$ in Fig. 2 (a). The fact that the lowest-energy magnetoquasistatic mode of a high-permittivity small sphere is the magnetic dipole is in agreement with experimental studies kuznetsov2012magnetic . The second eigenvalue has an eight fold degeneracy: it is associated to five TE modes $\mathbf{M}_{pm2}^{\left(1\right)}\left(r_{11}\tilde{\mathbf{r}}\right)$ with $p=e,o$ and $m=0-2$ (in Fig. 2 (b’) $\mathbf{M}_{o12}^{\left(1\right)}\left(r_{11}\tilde{\mathbf{r}}\right)$ is shown) and three TM modes $\mathbf{N}_{pm1}^{\left(1\right)}\left(r_{11}\tilde{\mathbf{r}}\right)$ with $m=0-1$ and $p=e,o$ (in Fig. 2 (b”) $\mathbf{N}_{e11}^{\left(1\right)}\left(r_{11}\tilde{\mathbf{r}}\right)$ is shown). The third eigenvalue is associated to seven TE modes $\mathbf{M}_{pm3}^{\left(1\right)}\left(r_{21}\tilde{\mathbf{r}}\right)$ with $p=e,o$ and $m=0-3$ (in Fig. 2 (c’) $\mathbf{M}_{o13}^{\left(1\right)}\left(r_{21}\tilde{\mathbf{r}}\right)$ is shown) and five TM modes $\mathbf{N}_{pm2}^{\left(1\right)}\left(r_{21}\tilde{\mathbf{r}}\right)$ with $p=e,0$ and $m=0-2$ (in Fig. 2 (c”) $\mathbf{N}_{e12}^{\left(1\right)}\left(r_{21}\tilde{\mathbf{r}}\right)$ is shown). The fourth eigenvalue is associated to three TE modes $\mathbf{M}_{pm1}^{\left(1\right)}\left(r_{02}\tilde{\mathbf{r}}\right)$ with $m=0-1$ and $p=e,o$, associated to two counter-rotating current loops (we show $\mathbf{M}_{o11}^{\left(1\right)}\left(r_{02}\tilde{\mathbf{r}}\right)$ in Fig. 2 (d)). The fifth eigenvalue is associated to nine TE modes $\mathbf{M}_{pm4}\left(r_{31}\mathbf{r}^{\prime}\right)$ with $m=0-4$ and $p=e,o$ (in Fig. 2 (e’) $\mathbf{M}_{o14}^{\left(1\right)}\left(r_{31}\tilde{\mathbf{r}}\right)$ is shown) and seven TM modes $\mathbf{N}_{pm3}^{\left(1\right)}\left(r_{31}\tilde{\mathbf{r}}\right)$ with $m=0-3$ and $p=e,0$ (in Fig. 2 (e”) $\mathbf{N}_{e13}^{\left(1\right)}\left(r_{31}\tilde{\mathbf{r}}\right)$ is shown). The sixth eigenvalue is associated to five TE modes $\mathbf{M}_{pm2}^{\left(1\right)}\left(r_{12}\tilde{\mathbf{r}}\right)$ with $m=0-2$ and $p=e,o$ (in Fig. 2 (f’) $\mathbf{M}_{o12}^{\left(1\right)}\left(r_{12}\tilde{\mathbf{r}}\right)$ is shown) and three TM modes $\mathbf{N}_{pm1}^{\left(1\right)}\left(r_{22}\tilde{\mathbf{r}}\right)$ with $m=0-3$ and $p=e,0$ (in Fig. 2 (f”) $\mathbf{N}_{e13}^{\left(1\right)}\left(r_{31}\tilde{\mathbf{r}}\right)$ is shown). All the magnetoquasistatic modes have vanishing normal component and non-zero tangential component on the sphere boundary.

We now show that the magnetoquasistatic eigenvalues predict the occurrence of the resonance peaks in the scattering response of a high-permittivity small sphere. In Fig. 3 we show with a continuous line the power absorbed by the sphere normalized by $\pi R^{2}$ as a function of $y=x\sqrt{\varepsilon_{R}}$ for different values of $\varepsilon_{R}$, calculated by the Mie theory Bohren1998 , and with four vertical dashed lines the first non-degenerate eigenvalues $y_{n}$ of the magnetostatic volume integral operator listed in Tab 1. Thus, the first vertical line represents the position of the mode $\mathbf{M}_{o11}^{\left(1\right)}\left(r_{01}\tilde{\mathbf{r}}\right)$ , the second line is associated to one TE mode $\mathbf{M}_{o12}^{\left(1\right)}\left(r_{11}\tilde{\mathbf{r}}\right)$ and one TM mode $\mathbf{N}_{e11}^{\left(1\right)}\left(r_{11}\tilde{\mathbf{r}}\right)$, the third is associate to $\mathbf{M}_{o13}^{\left(1\right)}\left(r_{21}\tilde{\mathbf{r}}\right)$ and $\mathbf{N}_{e12}^{\left(1\right)}\left(r_{21}\tilde{\mathbf{r}}\right)$, the forth line to $\mathbf{M}_{o11}^{\left(1\right)}\left(r_{02}\tilde{\mathbf{r}}\right)$. They are shown in Fig. 2 (a-d). In Fig. 3 (a) we consider $x\in\left[0.025,0.065\right]\ll 1$ and $\varepsilon_{R}=10^{4}+i10^{2}$. We find very good agreement between the eigenvalues $y_{n}$ and the absorption peaks positions. This is expected because the investigated values of $x$ are much less than one. Then in Fig. 3 (b) we consider $x\in\left[0.25,0.65\right]<1$ and $\varepsilon_{R}=10^{2}+i$. We note a red-shift of the peaks with respect to the magnetoquasistatic prediction, because the values of $x$ starts to be comparable to one. Eventually, in Fig. 3 (c) we investigate a Silicon sphere with $x\in\left[0.64,1.6\right]$ and $\varepsilon_{R}=15.45+0.1456i$. Although $x$ is now comparable to one, there is still a correlation between the peaks and the magnetoquasistatic predictions. However, the peaks now show a broadening and a red-shift. The third and fourth peaks appearing in Fig. 3 (c) arise from the splitting of the third peak of panel (b) because the degeneracy of the modes $\mathbf{M}_{o13}^{\left(1\right)}\left(r_{21}\tilde{\mathbf{r}}\right)$ and $\mathbf{N}_{e12}^{\left(1\right)}\left(r_{21}\tilde{\mathbf{r}}\right)$ is lifted by the retardation.

Finite-Size Cylinder. We now investigate a finite-size cylinder of radius $R$ and height $h=R$, which is very common among nanofabricated structures, because it is compatible with planar nanofabrication processes. Although no analytic solution exists in this case, semi-empirical formulas have been proposed for low-index resonances obtained by brute-force numerical calculations Mongia:94 . We assume a characteristic size $l_{c}=R$.

In Fig. 1 (b) we show the first $100$ magnetoquasistatic eigenvalues of the investigated cylinder. In Fig. 4 we plot the displacement current density modes corresponding to the first seven distinct eigenvalues. A similar catalog has been produced by Kajfez et al. Kajfez:84 by using a surface integral equation formulation of the full-Maxwell equations for bodies of revolution. Here, we follow the classification introduced of Glisson et al. Glisson:83 , which is in turn borrowed from the literature on cylindrical waveguides snitzer1961cylindrical ; Standards . Specifically, the modes which are symmetric along the azimuthal direction are denoted either as $\text{TM}_{0n\delta}$ or as $\text{TM}_{0n\delta}$. All the remaining modes are known as hybrid modes with respect to the axis of rotation, and they are denoted by $\text{HEM}_{mnp}$, where the subscripts $m$, $n$, $p$ are associated to the number of oscillation of the mode along the azimuthal, radial, and axial directions. It is also worth to note the third subscript is denoted as $\delta$ if smaller than unity. The fundamental mode of the finite size cylinder shown in Fig. 4 (a) is $\text{TE}_{01\delta}$, i.e. a magnetic dipole oriented along the vertical axis. The next two degenerate modes, e.g. $\text{HEM}^{\prime}_{11\delta}$- $\text{HEM}^{\prime\prime}_{11\delta}$, shown in Fig. 4 (b’)-(b”), radiate like magnetic dipoles oriented along two orthogonal horizontal directions. We note that while in a high-permittivity small sphere the three magnetic dipoles are degenerate, in the investigated cylinder they are split up by the symmetry breaking, with $\text{TE}_{01\delta}$ having lower magnetoquasistatic energy with respect to $\text{HEM}^{\prime}_{11\delta}$- $\text{HEM}^{\prime\prime}_{11\delta}$. Moving further in the direction of increasing magnetostatic energy, we encounter the hybrid modes $\text{HEM}^{\prime}_{12\delta}$,$\text{HEM}^{\prime\prime}_{12\delta}$ followed by the azimuthally symmetric mode $\text{TM}_{01\delta}$. Then, we encounter the hybrid modes $\text{HEM}_{21\delta}^{\prime}$, $\text{HEM}_{21\delta}^{\prime\prime}$, then $\text{TE}_{011+\delta}$, followed by the modes $\text{HEM}_{22\delta}^{\prime}$,$\text{HEM}_{22\delta}^{\prime\prime}$.

In Table 2 we compare the values of the magnetoquasistatic eigenvalues $y_{n}$ against semi-empirical formulas Mongia:94 (when available) obtained by fitting the results of different numerical methods for the modes $\text{TE}_{01\delta}$ DeSmedt:84 , $\text{HEM}_{11\delta}$ kishk1993computed , $\text{TM}_{01\delta}$ tsuji1984analytical , $\text{TE}_{011+\delta}$ DeSmedt:84 .

As for the sphere, we show in Fig. 5 that from the knowledge of the magnetoquasistatic eigenvalues we can predict the occurrence of the resonance peaks in the absorbed power spectra. The absorbed power has been calculated by an independent full-wave numerical method, i.e. the PMWCHT approach harrington1993field ; ForestiereSIE . We show with vertical dashed lines the first six non-degenerate eigenvalues $y_{n}$ associated to the modes $\text{TE}_{01\delta}$, $\text{HEM}_{11\delta}$, $\text{HEM}_{12\delta}$, $\text{TM}_{01\delta}$, $\text{HEM}_{21\delta}$, $\text{TE}_{011+\delta}$, whose values are listed in Tab. 2. In Fig. 5 (a) we consider $\varepsilon_{R}=10^{4}+10^{2}i$ and $x\in\left[0.025,0.065\right]\ll 1$: The magnetoquasistatic eigenvalues exactly predict the absorption peaks because the hypotheses of the magnetoquasistatic model are verified. Next, in Fig. 5 (b) we consider $x\in\left[0.25,0.65\right]<1$ and $\varepsilon_{R}=10^{2}+10i$: we note a red-shift and a broadening of the peaks because the values of $x$ are approaching one. Eventually, in Fig 5 (c), we consider a silicon cylinder with $x\in\left[0.64,1.6\right]\approx 1$ and $\varepsilon_{R}=15.45+0.1456i$. Even if $x\approx 1$, there is still a correlation between the peaks and the magnetoquasistatic predictions. However, the shift and the broadening of the peaks are now significant.

Right Triangular Prism. We now present the magnetoquasistatic analysis of a right triangular prism, which is not axisymmetric as the objects investigated so far, but belongs to the symmetry group $D_{3h}$. The triangular prism has height $H$, while its basis is an equilater triangle of edge $L=2H$. In Fig. 1 (c) we show the first 100 magnetostatic eigenvalues, while in Fig. 6 we show the corresponding modes. It is apparent from Tab. 3 that the eigenvalues have multiplicity one or two, which is consistent with the fact that the prism is invariant under transformation of the group $D_{3h}$ lyubarskii2013application . The fundamental mode of a triangular prism is a magnetic dipole oriented along the vertical axis, shown in Fig. 6 (a). The next two degenerate modes are magnetic dipoles lying on the horizontal plane and oriented along one height and the corresponding orthogonal edge, as shown in Fig. 6 (a)-(b). By moving further in the direction of increasing magnetostatic energy we note that the modes in Fig. 6 (c) and (d)-(d’) resemble the modes $\text{TM}_{01\delta}$ and $\text{HEM}_{12\delta}$ of the cylinder but the order of the corresponding magnetostatic energy is inverted. This fact suggests that the order of the magnetoquasistatic resonances can be tailored by a convenient design of the geometry of the object.

Ring. Eventually, we investigate a dielectric ring (solid torus) with minor radius $r$ and major radius $R=3r$, whose boundary has a different genus with respect to the object that we have investigated so far. Similar rod-ring shaped dielectric resonators have been investigated in the context of antennas, e.g. ong2004rod . In Fig. 1 (d) we show the first 100 magnetostatic eigenvalues, while in Tab. 4 we list the values only of the low-order ones.

We catalog these modes using a magnetic coordinate system wesson2011tokamaks where $\bf\hat{\phi}$, $\bf\hat{r}$, $\bf\hat{\theta}$ are the toroidal, radial and poloidal directions, respectively. We describe the number of oscillation of the displacement current density modes along the $\bf\hat{\phi}$, $\bf\hat{r}$, $\bf\hat{\theta}$ direction by the toroidal $t$, radial $r$, and poloidal $p$ numbers. The modes which are invariant along the toroidal direction ($t=0$) are denoted either as toroidal $\text{T}_{0rp}$ if the displacement current is directed along $\bf\hat{\phi}$ or poloidal $\text{P}_{0rp}$ if the displacement current is directed along $\bf\hat{\theta}$. We denote all the remaining modes as hybrid toroidal-poloidal modes $\text{HTP}_{mrp}$.

In Fig. 7 we show the displacement current modes of the first seven non-degenerate eigenvalues. The first mode, shown in Fig 7 (a), is the fundamental toroidal mode $\text{T}_{010}$: its magnetostatic energy is significantly lower with respect to the remaining modes. The fundamental poloidal mode $\text{P}_{010}$ is the seventh non-generate mode, shown in Fig. 7 (g). The remaining modes in Fig. 7 are hybrid. Specifically, the modes in Fig. 7 (c)-(c’) and (d)-(d’) are characterized by the same numbers $\left(t,r,p\right)=\left(1,1,1\right)$: we distinguish them as up-down $\text{HTP}_{111}^{\updownarrow}$ (c)-(c’), and inboard-outboard $\text{HTP}_{111}^{\leftrightarrow}$ (d)-(d’). Similarly, the modes in Fig. 7 (b)-(b’) and (e)-(e’) are denoted as $\text{HTP}_{211}^{\updownarrow}$ and $\text{HTP}_{211}^{\leftrightarrow}$. Next, we show in Fig. 7 (f) $\text{HTP}_{311}^{\updownarrow}$.

III Conclusions

There exist two mechanisms through which a small non-magnetic homogeneous object may resonate SupMat . The first is the electroquasistatic resonance Mayergoyz:03 where the induced electric charge plays a central role. These resonances are connected to the eigenvalues of the electrostatic integral operator that gives the electrostatic field as a function of the charge density. They physically arise from the interplay between the energy stored in the electric field and the kinetic energy of the electrons. Each resonance is characterized by a negative eigenpermittivity, which is size-independent. The eigenpermittivities constitute a countable infinite set, and accumulate at the point $-1$. The induced current density fields are both curl-free and div-free within the particle, but have non-vanishing (and discontinuous) normal components to the particle surface.

The second mechanism is the magnetoquasistatic resonance, described in this paper, where the displacement current density field is the main player. These resonances are connected to the eigenvalues of a magnetostatic integral operator that gives the vector potential as a function of the current density. They arise from the interplay between the polarization energy stored in the dielectric and the energy stored in the magnetic field. These resonances are only possible for positive permittivity. For any given shape, the resonance frequencies are inversely proportional to the characteristic size of the object, and inversely proportional to the square root of the permittivity. They are an infinite but countable set accumulating at $+\infty$. The induced current density fields have a non-zero curl within the particle, but are div-free and have a vanishing normal component on the particle surface.

The applicability of the magnetoquasistatic approximation can be extended to accurately describe the radiation damping of the modes and the frequency-shift due to the finite particle particle size by using perturbation techniques.

Appendix A Solution of the electromagnetic scattering in the quasi-static limit

Let us consider an isotropic and homogeneous material occupying a volume $\Omega$, which is bounded by a closed surface $\partial\Omega$ with outward-pointing normal $\hat{\mathbf{n}}$ . The medium is nonmagnetic with relative permittivity $\varepsilon_{R}$, surrounded by vacuum. The object is illuminated by a time harmonic electromagnetic field incoming from infinity $\text{Re}\left\{\mathbf{E}_{inc}\left(\mathbf{r}\right)e^{i\omega t}\right\}$, where $\omega$ is the angular frequency. Here, we derive the behavior of the electromagnetic scattering in the quasi-static limit starting from the full wave model.

The scattering problem is formulated by considering as unknown the effective current density field $\mathbf{J}=\mathbf{J}\left(\mathbf{r}\right)$ induced in the body (which particularizes into conduction current in metals at frequencies below interband transitions, polarization current in dielectrics, sum of conduction and polarization currents in metals in frequency ranges where interband transitions occur). We have $\mathbf{J}=i\omega\varepsilon_{0}\chi\mathbf{E}$ where $\mathbf{E}=\mathbf{E}\left(\mathbf{r}\right)$ is the total electric field (induced + incident) and $\chi=\left(\varepsilon_{R}-1\right)$ is the electric susceptibility. Both the vector fields $\mathbf{E}$ and $\mathbf{J}$ are div-free in $\Omega$ due to the homogeneity and isotropy of the material. The current density $\mathbf{J}$ is governed by the full wave volume integral equation van2007electromagnetic ; Schaubert:84 ; jin2011theory ; rubinacci2006broadband ; miano2010numerical :

[TABLE]

where $\varepsilon_{0}$ is the vacuum permittivity, $\mu_{0}$ is the vacuum permeability, $g\left(\mathbf{r}\right)=e^{-ik_{0}r}/4\pi r$ is the Green function in vacuum, $k_{0}=\omega/c_{0}$ and $c_{0}=1/\sqrt{\varepsilon_{0}\mu_{0}}$. The surface integral represents the contribution of the scalar potential to the induced electric field and the volume integral represents the contribution of the vector potential. We introduce the dimensionless size parameter $x=\omega l_{c}/c_{0}$ where $l_{c}$ is a characteristic linear length of the region $\Omega$. Then, equation 12 is rewritten as

[TABLE]

where $\tilde{\mathbf{r}}=\mathbf{r}/l_{c}$,

[TABLE]

$\tilde{\Omega}$ is the scaled domain, $\tilde{\nabla}$ is the scaled gradient operator and $W_{n}={\bf W}\cdot\hat{\mathbf{n}}$. In the quasi-static limit ($x\rightarrow 0$) the operator $\mathcal{L}$ has the following expression up to the third order in x:

[TABLE]

$\text{as}\,x\rightarrow 0$, where

[TABLE]

is the static Green function for the vacuum.

We now study the solution of equation 13 in the quasi-static limit $x\rightarrow 0$. To achieve this purpose we introduce a complete basis for representing the unknown, which is obtained by joining two orthogonal sets. The first set $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ is given by the solution of the eigenvalue problem

[TABLE]

where $\mathcal{L}_{e}$ is the electrostatic integral operator

[TABLE]

The eigenfunctions $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ are both div-free and curl-free in the limit $x\rightarrow 0$ in $\Omega$, but they have non-vanishing normal components to $\partial\Omega$. Since $\mathcal{L}_{e}$ is Hermitian and definite negative the eigenvalues $\gamma_{h}^{\parallel}$ are real and negative, and the eigenfunctions are orthonormal according to the scalar product

[TABLE]

Both the eigenfunctions $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ and the eigenvalues $\left\{\gamma_{h}^{\parallel}\right\}$ do not depend on the size of object, but only on its shape. However, the set $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ is not sufficient to represent the vector space of square integrable div-free functions in $\Omega$. To complete the basis it is sufficient to add to $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ the set of solenoidal vector fields $\left\{\mathbf{W}_{h}^{\perp}\right\}$ with vanishing normal components to $\partial\Omega$ and that are solution of the eigenvalue problem in weak form

[TABLE]

where $\mathcal{L}_{m}$ is the magnetostatic integral operator

[TABLE]

The eigenfunctions $\left\{\mathbf{W}_{h}^{\perp}\right\}$ are div-free in $\Omega$ and have vanishing normal components to $\partial\Omega$, but their curl in $\Omega$ is different from zero. Since $\mathcal{L}_{m}$ is Hermitian and definite positive the eigenvalues $\gamma_{h}^{\perp}$ are real and positive, and the eigenfunctions are orthonormal according to the scalar product 19. Both the eigenfunctions $\left\{\mathbf{W}_{h}^{\perp}\right\}$ and the eigenvalues do not depend on the size of object, but only on its shape. Furthermore, the set of eigenfunctions $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ is orthogonal to set of eigenfunctions $\left\{\mathbf{W}_{h}^{\perp}\right\}$ . The union of the two sets $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ and $\left\{\mathbf{W}_{h}^{\perp}\right\}$ is a complete basis for the vector space of square integrable div-free vector fields in $\Omega$. Therefore, we represent the unknown current density $\mathbf{J}$ as

[TABLE]

We substitute expression 22 into equation 13 and we obtain to the leading order in $x$ the following expressions for the expansion coefficients $I_{h}^{\perp}$ and $I_{h}^{\perp}$:

[TABLE]

Expression 22 with 23 is the solution of the scattering problem 13 in the quasi-static limit $x\rightarrow 0$. The eigenfunction $\mathbf{W}_{h}^{\parallel}$ may be resonantly exited when $\text{Re}\left\{\chi\right\}=\gamma_{h}^{\perp}$ and the eigenfunction $\mathbf{W}_{h}^{\perp}$ may be resonantly exited when $\text{Re}\left\{\chi\right\}=\gamma_{h}^{\parallel}/x^{2}$. For these reasons the sets $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ and $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ can be interpreted as the current modes of the body in the quasi-static limit $x\rightarrow 0$, and $\left\{\gamma_{h}^{\parallel}\right\}$ and $\left\{\gamma_{h}^{\perp}\right\}$ can be interpreted as the corresponding eigen-susceptibilities. Thus, we call the eigenfunctions $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ electroquasistatic current modes and the eigenfunctions $\left\{\mathbf{W}_{h}^{\parallel}\right\}$ magnetoquasistatic current modes. Since the eigen-susceptibilities $\left\{\gamma_{h}^{\parallel}\right\}$ are all negative the electroquasistatic modes can be resonantly excited only in metals (surface plasmons). On the contrary, the eigen-susceptibilities $\left\{\gamma_{h}^{\parallel}/x^{2}\right\}$ are all positive, therefore the magnetoquasistatic modes can be resonantly excited only in dielectrics. The scalar products $\langle\mathbf{W}_{h}^{\parallel},\mathbf{E}_{inc}\left(\mathbf{r}\right)\rangle$ and $\langle\mathbf{W}_{h}^{\perp},\mathbf{E}_{inc}\left(\mathbf{r}\right)\rangle$ in Eq. 23 describe the coupling of the electroquasistatic and magnetoquasistatic current modes with the external excitation.

The magnetoquasistatic eigenvalue problem 20 is the problem 5 of the main manuscript where $\gamma^{\perp}$ is replaced by $y^{2}$. Since $\displaystyle\min_{h\in\mathbb{N}}{\gamma_{h}^{\parallel}}\approx 1$ and $x\ll 1$ the magnetoquasistatic resonances occur in dielectrics with high permittivity, hence the resonance frequency of the $h^{th}$ magnetoquasistatic mode is $\omega_{h}\approx\left(\frac{c_{0}}{l_{c}\sqrt{\varepsilon_{r}}}\right)\sqrt{\gamma_{h}^{\perp}}$: this is the resonant condition 7 of the main manuscript (where $\gamma^{\perp}$ is replaced by $y^{2}$ ). Moreover, the problem 20 can be rewritten in differential form as in the following:

[TABLE]

where $\tilde{\nabla}\times$ is the scaled curl, $\mathbf{A}^{\perp}\left(\mathbf{r}\right)=\mathcal{L}_{m}\left\{\mathbf{W}^{\perp}\right\}\left({\bf r}\right)$ for $\tilde{\mathbf{r}}\in\tilde{\Omega}$ and $\tilde{\Pi}\left(\tilde{\mathbf{r}}\right)=1$ for $\tilde{\mathbf{r}}\in\tilde{\Omega}$ and [math] for $\tilde{\mathbf{r}}\notin\tilde{\Omega}$. Note that at resonance we have $\gamma_{h}^{\perp}=l_{c}^{2}\beta^{2}$. This validates the magnetoquasistatic model defined by Eqs. 1a, 1b, and 2 of the main manuscript.

Appendix B Energy balance in the magnetoquasistatic resonances

In the electroquasistatic resonances of metals below the interband transitions the energy oscillates back and forth between the kinetic energy of the conduction electrons and the Coulomb potential energy arising from the surface charges on the surface of the metal. In the magnetoquasistatic resonances in dielectrics the energy oscillates back and forth between the polarization energy of the dielectric and the magnetic energy. Indeed we now show that in these resonances the energy stored in the magnetic field is balanced by the energy stored in the dielectric in the form of polarization energy. By assuming that the dielectric susceptibility is real, the resonance frequency $\omega_{h}$ for the magnetoquasistatic current mode $\mathbf{W}_{h}^{\perp}$ is given by the condition

[TABLE]

The eigenvalue $\gamma_{h}$ is related to the time average of the magnetic energy. Indeed, we have

[TABLE]

where $\mathbf{A}_{h}^{\perp}$ is the magnetic vector potential generated by the current mode $\mathbf{W}_{h}^{\perp}$ and $\tilde{\nabla}\times$ is the scaled curl. By combining Eqs. 25 and 26 we obtain

[TABLE]

The term on the left hand side is the energy stored in the magnetic field associated to the current mode $\mathbf{W}_{h}^{\perp}$ while term on the right hand side is the energy stored in the dielectric, in the form of polarization energy, at the resonance frequency $\omega_{h}$.

Appendix C Numerical Model

Equation 5 can be discretized by drawing on the standard repertoire of computational electromagnetics for Volume Integral Equations Schaubert:84 ; jin2011theory ; van2007electromagnetic . The unknown of the magnetoquasistatic problem, i.e. the displacement current density field ${\bf J}^{d}$, belongs to the functional space $\mathcal{J}_{L}$ bossavit1998computational :

[TABLE]

We obtain the discretization of Eq. 5 by representing ${\bf J}^{d}$ in terms of $N_{L}$ loop shape functions ${\bf w}^{L}_{k}$

[TABLE]

Each function $\mathbf{w}_{k}^{\text{L}}$ is associated to the $k$-th edge of the finite element discretization of the volume $\Omega$. It is defined as the curl of the $k$-th edge-element shape functions bossavit1998computational $\mathbf{N}_{k}$:

[TABLE]

The discrete generalized eigenvalue problem is obtained by substituting the representation (28) into Eq. (5) and applying the Galerkin method, projecting along the loop shape functions:

[TABLE]

The matrices $\mathbf{R}$ and $\mathbf{L}$ are associated to the l.h.s. and r.h.s. of Eq. 5, respectively. The generic occurrences of these matrices are:

[TABLE]

Eventually, the problem 29 is reduced to a standard symmetric eigenvalue problem by exploiting the LAPACK laug routine DSYGST, then all eigenvalues and eigenvectors of the resulting real symmetric matrix are computed through the routine DSYEV. In Tab. 5 we provide details about the number of nodes, elements, and edges of the hexahedral meshes used in the calculation.