Resonance frequency and radiative Q-factor of plasmonic and dielectric modes of small objects

The electromagnetic scattering resonances of a non-magnetic object much smaller than the incident wavelength in vacuum can be either described by the electroquasistatic approximation of the Maxwell's equations if its permittivity is negative, or by the magnetoquasistatic approximation, if its permittivity is sufficiently high. Nevertheless, these two approximations fail to correctly account for the frequency shift and the radiative broadening of the resonances when the size of the object becomes comparable to the wavelength of operation. In this manuscript, the radiation corrections to the electroquasistatic and magnetoquasistatic resonances are derived, which only depend on the quasistatic current modes. Closed form expressions of the frequency-shift and the radiative Q-factor of either plasmonic or dielectric modes of small objects are introduced, where the dependencies on the material and size of the object are explicitly factorized.


I. INTRODUCTION
There exist two mechanisms through which a nonmagnetic object, assumed small compared to the wavelength, may resonate.
One resonance mechanism occurs in small objects of high and positive permittivity, arising from the interplay between the polarization energy stored in the dielectric and the energy stored in the magnetic field. Manifestation of this kind of resonance can be found at microwaves where low-loss dielectric materials with relative permittivities of the order of 100 are routinely used for various applications including resonators, filters, and antennas [1][2][3][4] . When the object is very small compared to the free-space wavelength, and the permittivity very high, these resonances are well-described by the magnetoquasistatic approximation of the Maxwells equations 5 , where the normal component of the displacement current density field vanishes on the surface of the particle 5,6 . In particular, these resonances are connected to the eigenvalues of the magnetostatic integral operator expressing the vector potential in terms of the displacement current density 5 . Unfortunately, at optical frequencies, the permittivity of the available dielectric materials is only moderately high [7][8][9][10][11] and the size of the object has to be comparable to the incident wavelength to trigger a resonant response. In this scenario the magnetoquasistatic approximation is unable to describe the frequency shift and the broadening of the resonances.
A dual resonant mechanism is the plasmonic resonance: it occurs in small metal nanoparticles with negative permittivity, arising from the interplay between the energy stored in the electric field and the kinetic energy of the free electrons of the metal. When the object is very small compared to the free-space wavelength, these resonances are described by the electroquasistatic approximation of the Maxwells equations [12][13][14][15][16] and associated to the negative values of susceptibilities in correspondence of which source-free solutions exist. However, it is known that as the size of the object becomes comparable to the incident wavelength this approximation is unable to describe the radiative shift and broadening of the resonances.
In light of these observations, to describe the resonances of objects of size comparable to the incident wavelength, one may be tempted to abandon the quasistatic approximation altogether and then turn to the full-Maxwell's equations. There are certainly some advantages in doing so, including the fact that the fullwave formulations would enable, given unlimited computational resources, the treatment of objects of any size. In facts, scattering resonances have been already investigated by considering full-wave eigenvalue problems based on volume (e.g. Refs. 17,18 ), surface [19][20][21][22][23] , and line integral formulations of the Maxwell's equations 24 , differential formulations (e.g. 25,26 ), Mie Theory (e.g. [27][28][29] etc. for a recent review see Ref. 30 ). However, the resulting resonance frequencies, Q-factors, and resonant modes depend on the morphology, the material, and the size of the object. A change of any of these parameters would require an entirely new calculation. These dependencies are buried below the computational layer, and cannot be factorized.
Closed form expressions of the Q-factor and frequency shift of both plasmonic and dielectric modes, where the dependencies on the material and size of the object are factorized, are highly desirable. They would enable the classification of the resonances, and facilitate their engineering, including the coupling with emitters 31-34 , because they could be used as a target for the design 35,36 . Moreover, in many applications the size of metal or dielectric objects does not exceed the free-space wavelength of operation 10 . These are powerful incentives to pursue the extension of the two quasi-static scattering limits to include radiation effect.
In the literature there already exist closed form expressions for the resonance frequency shift and Q-factors in few scenarios. For plasmonic modes, in Ref. 37 Mayergoyz et. al derived the second order correction for electroquasistatic eigenvalues, starting from Maxwell's equations in differential form. In Ref. 15 Wang and Shen derived the expression of the Q-factor of a plasmonic mode, when the non-radiative losses are dominant. To the authors' best knowledge the derivation of a general expression for the Q-factor of plasmonic modes when radiative losses are dominant, sometimes called anomalous light scattering regime 38 , is still missing.
For high-index dielectric resonators, Van Bladel introduced closed form expressions for the Q-factor of the fundamental magnetic dipole, and provided the Q-factor of higher order modes of a rotationally symmetric object 6 , considering an asymptotic expansion of the Maxwell's equation in differential form in terms of the inverse of the index of refraction. Following Van Bladel's work, De Smedt derived the frequency shift and Q-factor of a rotationally symmetric ring resonator 39 . General expressions for both the frequency shift and radiative factor haven't been derived yet.
In this paper, the radiation corrections of both the electroquasistatic and magnetoquasistatic resonances and resonant modes of arbitrarily-shaped non-magnetic homogeneous and isotropic objects are introduced, starting from an integral formulation of the Maxwell's equations. It is demonstrated that, in the scattering from small particles, the frequency shift of any mode is a quadratic function of the size parameter, whose prefactor depends on the ratio between the second order correction and the quasistatic eigenvalue. Furthermore, the radiative Q-factor is an inverse power function of the size parameter ∝ x −ni whose exponent is the order n i of the first non-vanishing imaginary correction, while the prefactor is the ratio between the quasistatic eigenvalue and the n i -th order imaginary correction.
This manuscript is organized as follows. First, the scattering resonances in the two quasistatic regimes are briefly summarized in Sec. II. Then, in Sec. III, the full wave scattering problem is formulated and an eigenvalue problem governing the scattering resonances is introduced. This eigenvalue problem is solved perturbatively in Secs. IV and V starting from the electroquasistatic and magnetoquasistatic limits, treating the size parameter as a small parameter. Collecting same-order terms, closed form radiation corrections are found. In Sec. VI, the frequency shift and the Q-factor are obtained as a function of these radiation corrections. In Sec. VII, the tables of plasmonic and photonic resonances are introduced. They constitute a synthetic classification of the modes of a homogeneous non-magnetic objects, and depend solely on its morphology, but not on its size, material, and frequency of operations. This classification may help the description of the elementary building blocks of the nano-circuitry envisioned by Engheta et al. in Ref. 40 . Eventually, the introduced formalism is validated by investigating the resonance frequency and Qfactors in the scattering response of metal and dielectric objects of size comparable to the incident wavelength.

II. RESONANCES IN THE QUASISTATIC REGIME
A homogeneous, isotropic, non-magnetic linear material occupies a volume Ω, of characteristic linear dimension l c , bounded by a closed surface ∂Ω with an outwardpointing normaln. The material has relative permittivity ε R (ω), and it is surrounded by vacuum.
There exist two mechanisms through this object may resonate in the quasistatic regime 5

A. Electroquasistatic resonances
The first resonance mechanism is the electroquasistatic resonance, occurring in metals (more in general, in objects whose dielectric permittivity has a negative real part) where the induced electric charge plays a central role. These resonances are connected to the eigenvalues χ h of the electrostatic integral operator L e that gives the electrostatic field as a function of the surface charge density 41 the expression of L e is 4π|r−r | . In Eq. 2 the spatial coordinates have been normalized by l c , i.e.r = r/l c ,Ω is the corresponding scaled domain, ∂Ω is the boundary ofΩ, and ∇ is the scaled gradient operator. The quasistatic oscillations represented by the electroquasistatic (EQS) current modes j h arise from the interplay between the energy stored in the electric field and the kinetic energy of the free electrons in the metal. The spectrum of the operator L e is discrete 37,41 . The EQS mode j h is characterized by a real and negative eigenvalue χ h , which is size-independent 37 . The modes j h h∈N are longitudinal: they are both curl-free and div-free within the object, but have non-vanishing normal component to the object surface 37,41 . This normal component is related to the induced surface charge density on ∂Ω, and satisfy the charge-neutrality condition ‹ ∂Ω j h (r ) ·ñ (r ) dS = 0.
Moreover, EQS modes associated to non-degenerate eigenvalues are orthonormal accordingly to the scalar product The EQS mode j h may exhibit zero dipole moment ˚Ω j h (r ) dṼ = 0, (6) in this case it is dark, it is bright otherwise.

B. Magnetoquasistatic resonances
The second mechanism is the magnetoquasistatic resonance, occurring in high-permittivity dielectric objects, where the displacement current density field plays a central role. These resonances are connected to the eigenvalues κ ⊥ h of the magnetostatic integral operator L m that gives the vector potential as a function of the current density 5 : with Equation 7 holds in the weak form in the functional space equipped with the inner product 5 constituted by the functions which are div-free withinΩ and having zero normal component to ∂Ω.
The quasistatic oscillations represented by the magnetoquasistatic (MQS) current density modes j ⊥ h arise from the interplay between the polarization energy stored in the dielectric and the energy stored in the magnetic field. The spectrum of the magnetoquasistatic operator 7 is discrete, too 5 . The MQS current mode j ⊥ h is characterized by a real and positive eigenvalue κ ⊥ h , which is size-independent. The current modes j ⊥ h h∈N are transverse modes: they have a non-zero curl within the object, but are div-free and have a vanishing normal component on the object surface 5 . Each current mode j ⊥ h has zero electric dipole moment, namely:Ω Furthermore, the current density modes associated to non-degenerate eigenvalues are orthonormal The current mode j ⊥ h may generate a vector potential with zero normal component to ∂Ω, namely: n (r) ·˚Ω g 0 (r −r ) j ⊥ h (r ) dṼ = 0 ∀r ∈Ω. (12) In this case, the MQS mode j ⊥ h is called A ⊥ -mode.
The longitudinal set of EQS current modes j h h∈N and the transverse set of MQS j ⊥ h h∈N modes are orthogonal accordingly to the scalar product 5, and together are a complete basis of the vector space of square integrable div-free vector fields in Ω.

III. ELECTROMAGNETIC MODES
The full-wave scattering problem can be formulated by considering as unknown the current density field J induced in the object. This current density particularizes into conduction current in metals at frequencies below interband transitions, polarization current in dielectrics, sum of conduction and polarization currents in metals in the frequency ranges where interband transitions occur. The object is illuminated by a time harmonic electromagnetic field Re E inc (r) e iωt . In the frequency domain the field J (r) is related to the electric field E (r) by is the electric susceptibility of the object, and ε 0 is the vacuum permittivity. Both the vector field E and J are div-free in the region Ω occupied by the object due to the homogeneity and isotropy of the material. The induced current density is solutions of the full wave volume integral equation [42][43][44] where µ 0 is the vacuum permeability, g (r) = e −ik0r /4πr is the Green function in vacuum, k 0 = ω/c 0 and c 0 = 1/ √ ε 0 µ 0 . The surface and volume integrals represent the contributions to the induced electric field of the scalar and vector potentials, respectively. Then, equation 13 is rewritten as 45 where the spatial coordinates are normalized asr = r/l c , Ω is the scaled domain, ∂Ω is boundary ofΩ,∇ is the scaled gradient operator, x is the size parameter x = 2πl c /λ, and g (r −r , x) is the Green function in vacuum and ∆r = |r −r |.
The spectral properties of the linear operator L play a very important role in any resonant scattering mechanism. Since L is compact its spectrum is countable infinite. The operator L is symmetric but not self-adjoint. For any value of the size parameter x its eigenvalues are complex with positive imaginary part. The eigenfunctions corresponding to two different eigenvalues are not orthogonal in the usual sense: they are bi-orthogonal 27,46 The eigenvalue problem 45 splits into the two eigenvalue problems 1 and 7 (see Sections II A and Section II B) in the quasi-static regime x 1 (small object). This fact was already shown for 2D objects in Ref. 23 and for 3D objects in Ref. 5 . The eigenfunctions of L that in the limit x → 0 tend to the EQS modes are indicated with {u h (r)} and the corresponding eigenvalues are indicated as {χ h }. These eigenfunctions are called plasmonic modes. Dually, the set of eigenfunctions of L that in the limit x → 0 tend to the MQS modes are indicated with {v h (r)} and the corresponding eigenvalues are indicated with κ h /x 2 . Although in the limit x → 0, the eigenvalues κ h /x 2 diverge, the quantities κ h remain constant. These eigenfunctions are called dielectric modes.
Forestiere and Miano et al. in Ref. 27,28 used the adjectives plasmonic and photonic mode instead of plasmonic and dielectric mode to identify the same two sets, while in Ref. 29 the authors called them longitudinal and transverse modes. All these nomenclatures are equivalent. It was shown that this two sets of modes, even if this distinction is made in the long-wavelength regime, remain well distinguishable and have different properties even in the full-wave regime 23,27,28 .
The union of the two sets {u h (r)} and {v h (r)} is a basis for the unknown current density field in equation 13. Its solution is expressed as where both the set of modes {u h (r)} and {v h (r)} are normalized, u * h , u h = 1 and v * h , v h = 1 for any h. This expansion is very useful because it separates the dependence on the material from the dependence on the geometry 23,24,27,45,46 , and has been used in different contexts 28,47,48 .
In the next two sections we develop a perturbation theory to evaluate the plasmonic and dielectric resonances and resonant modes of an object with arbitrary shape and size parameter x 1 by starting from the corresponding modes in the quasistatic regime.

IV. PLASMONIC RESONANCES
To evaluate the plasmonic resonances of small particles, it is convenient to recast the eigenvalue problem 17 as When the free-space wavelength λ = 2πc 0 /ω is large in comparison with the characteristic dimension l c , the size parameter x can be treated as a small parameter and the Green function g (r), the current mode u h , and the eigenvalue χ h can all be expanded in terms of x in the neighborhood of the EQS resonance with eigenvalue χ h and mode j h By using Eqs. 20, 21, and 22, Eq. 19 becomes where u In the Supplemental Material all the details of the deriva-tion of radiation corrections for plasmonic resonances are reported 49 .
Matching the first-order terms in Eq. 23, it is obtained that the first order corrections vanish regardless of the object's shape Collecting the second order terms in Eq. 23, and applying the normal solvability condition of Fredholm integral equations 50,51 , the second order correction χ (2) h is derived and A n (r) = A ·n, B n (r) = B ·n on ∂Ω. Accordingly to Eq. 26, χ h is real. Moreover, the first term in parenthesis in Eq. 26 (starting from the left) is associated to the radiative self-coupling of the surface charge density of the EQS mode j h through by the scalar potential, the second term is instead associated to the quasistatic self-coupling of the EQS mode j h through the vector potential. A second order correction to the EQS modes has been already derived in Ref. 37 by expanding the Maxwell's equation in differential form. It will be demonstrated in Eq. 57 that χ (2) h is connected to the frequency-shift of the plasmonic mode.
The second order correction of the associated plasmonic mode u (2) h has both longitudinal and transverse components, denoted as u where the longitudinal part u is represented in terms of the EQS modes basis j k k∈N , and the transverse part in terms of the MQS modes basis j ⊥ k k∈N . The expansion coefficients are α (2) h,h = 0 and ∀k = h Matching the third order terms in Eq. 23, the third order correction is obtained, which is purely imaginary The first term in parenthesis in Eq. 32, starting from the left, arise from the radiative self-coupling of the surface charges of the EQS mode j h through the scalar potential. The second term is instead associated to the magnitude of the net-dipole moment of the mode j h , and thus it vanishes for dark modes. If χ (3) h = 0, the third order correction determines, as it will be demonstrated in Eq. 58, the radiative Q-factor of h-the plasmonic mode.
However, χ h may vanish due to the symmetry. In this case, to retrieve information about the radiative Q-factor, it is therefore mandatory to consider the fifth order perturbation χ (5) h . When χ where α (2) h,k has been defined in Eq. 30.

V. DIELECTRIC RESONANCES
To evaluate the dielectric resonances beyond the quasistatic regime it is convenient to recast the eigenvalue problem 17 as As for the plasmonic scenario, the mode v and the corresponding eigenvalue κ are expanded at x = 0 in the By substituting Eqs. 35, 36, and 22 in Eq. 34 where In the Supplemental Material all the details on the derivation of radiation corrections for dielectric resonances are reported 49 . Here, only the main results are highlighted.
By matching the terms of corresponding order in Eq. 37, it is possible to demonstrate that first order corrections vanish regardless of the shape of the object: The second order correction κ (2) h is a real quantity 40) The first term in parenthesis in Eq. 40 represents the radiative self-coupling of the MQS mode j ⊥ h through by the vector potential, the second term accounts for the quasistatic coupling of the h-th MQS mode j ⊥ h with the EQS modes. As it will be demonstrated in Eq. 64, κ (2) h is connected to the frequency-shift of dielectric modes.
As shown in Eq. 12, the MQS mode j ⊥ h may be an A ⊥ -mode, generating a vector potential with zero normal component to ∂Ω. In this case, the term j k |∆r −1 |j ⊥ h Ω vanishes ∀k, and j ⊥ h does not couple with any EQS modes j k . For A ⊥ -modes, Eq. 40 further simplifies The second order correction v (2) h to the current density mode has both longitudinal and transverse components, denoted as v For A ⊥ -modes the longitudinal part of v (2) h vanishes. The third order correction κ (3) h is purely imaginary: and originates from the radiative self-coupling of the MQS mode j h through by the vector potential. The mode correction v where β h,h = 0, and As it will be shown in Eq. 45, the correction κ h = 0 is associated to the radiative Q-factor of dielectric modes.
The third order correction κ (3) h may vanish due to the symmetry. In this case, it is mandatory to consider the fifth order correction χ (5) h . If κ h has the following simplified expression: For A ⊥ -modes Eq. 49 further simplifies

VI. RESONANCE FREQUENCY AND Q-FACTOR
In the previous section, the second order corrections χ (2) h , κ (2) h and non-vanishing imaginary connections χ of the lowest order, i.e. n i , are derived in closed form for both plasmonic and dielectric modes. They depend neither on the size of the object nor on its permittivity, but they solely depend on the morphology of the EQS and MQS modes. In this section, closed form expressions of the resonance frequency and the Q-factors are obtained in terms second order and imaginary corrections of the lowest order in both metal dielectric objects.

A. Plasmonic Resonances
It is now assumed that the object is made of a timedispersive metal described by the Drude model 52,53 where ω p and ν are the plasma and collision angular frequencies, and ν ω p . It is also useful to define the quantity where λ p is the plasma wavelength. The EQS resonance frequency ω h of the h-th mode, is defined as the frequency at which the real part of the metal susceptibility Re {χ (ω)} matches the EQS eigenvalue χ h , i.e.
where x h is the size parameter at EQS resonance.
In the full-wave scenario, by looking at the denominator in the first summation of Eq. 18, the value x h of the size-parameter at the plasmonic resonance is the value of x at which the real part of the metal susceptibility χ (ω) matches the real part of the corresponding eigenvalue and ω h is the corresponding resonance frequency. Eq. 54 is the resonance condition of plasmonic modes. For small particles x p 1, by retaining only the real and imaginary non-zero corrections of the lowest order in Eq. 20, the plasmonic eigenvalue χ h (x) is approximated as where n i is the order of the first non-zero imaginary correction χ (ni) h . By using Eq. 55 in 54, and solving the resulting biquadratic equation (56) In the limit x p 1, the frequency shift ∆ω h = ω h − ω h , and corresponding shift in the size parameter at the resonance ∆x h = x h − x h can be approximated as In conclusion, the relative frequency shift of any plasmonic mode is a quadratic function of x h , whose prefactor is one half the ratio between the second order correction χ (2) h and the EQS eigenvalue χ h . The radiative Q-factor Q r h of the h-th plasmonic mode is obtained by considering the first summation in 18 when the resonance condition 54 is met and assuming negligible non-radiative losses The Q-factor is an inverse power function of the size parameter at the resonance whose exponent is the order n i of the first non-vanishing imaginary correction χ (ni) h , while the prefactor is the ratio between the EQS eigenvalue χ h and χ (ni) h . The regime in which the radiative losses are dominant occurs rarely in metal nano object, and for this reason it is also called anomalous scattering regime 38 .
The opposite regime, dominated by non-radiative losses, is more common. In this case the non-radiative Q-factor Q ⊥nr h is obtained from Eq. 18 Equation 59 is not new, but it was already shown by Wang and Shen in Ref. 15 .
In an intermediate regime, the resulting Q-factor, indicated with Q h , can be obtained as 54

B. Dielectric Resonances
It is now assumed that the object is made of a dielectric material with positive susceptibility χ ≥ 0, with Im {χ} Re {χ}. The size parameter x ⊥ h = 2πω ⊥ h /c 0 at the MQS resonance is defined as the value of x at which the real part of the susceptibility χ matches the eigenvalue κ ⊥ h /x 2 , namely: and ω ⊥ h is the corresponding MQS resonance frequency. In the full-wave regime, by looking at the denominator in the second summation of Eq. 18, the value of size parameter x h = ω h /c 0 l c at the dielectric resonance is the value of x at which the real part of the eigenvalue κ h /x 2 matches the quantity Re {χ} This is the resonance condition for dielectric modes, and ω h is the dielectric resonance frequency. For small particles x 1, by keeping only the real and imaginary non-zero corrections of the lowest order in Eq. 35, the dielectric eigenvalue κ h (x) is approximated as where n i is the order of the first non-zero imaginary correction κ (ni) h . By using Eq. 63 in 62, and solving the resulting quadratic equation For high-index dielectrics χ 1, the relative frequency shift is In conclusion, the relative frequency-shift of any dielectric mode is a quadratic function of x ⊥ h , whose prefactor is approximately the ratio between the second order correction κ (2) h and the quasistatic eigenvalue κ ⊥ h . The radiative Q-factor Q ⊥r h of the h-th dielectric mode is obtained by considering the second summation in 18 when the resonance condition 62 is met, and assuming negligible non-radiative losses The radiative Q-factor is an inverse power function of the size parameter ∝ x ni whose exponent is the order n i of the first non-vanishing imaginary correction κ (ni) h , while the prefactor is the ratio between the quasistatic eigenvalue κ ⊥ h and κ h . In dielectric resonators the opposite regime, dominated by non-radiative losses, is less common. In this case the non-radiative Q-factor Q ⊥nr In an intermediate regime, the Q-factor, indicated with Q ⊥ h , is obtained as 54 (68)

VII. RESULTS AND DISCUSSION
Given the sole shape of an homogeneous object, the table of plasmonic modes and the table of dielectric modes can be introduced. The two tables of a sphere are shown in Figs. 1 and 2; the tables of a finite-cylinder with height equal to the radius are shown in Figs. 5 and 6. The tables illustrate the EQS and MQS current modes, ordered according to their real quasistatic eigenvalue χ h and κ ⊥ h , respectively. Besides the quasistatic eigenvalue, each resonance is also characterized by the second order correction κ (2) h and χ (2) h , and by the lowest-order (nonvanishing) imaginary correction κ h , where n i is odd and n ≥ 3. In both tables, the value of n i is also highlighted, enclosed in a circle on the top-right of each box. Plasmonic modes are also labeled with a "D" if dark, while dielectric resonances are labeled with "A ⊥ " if the vector potential generated by the current mode has zero normal component to ∂Ω. The information contained in these two tables depends neither on the permittivity of the object nor on its size nor on the frequency of operation; it depends solely on the morphology of the quasistatic modes.
The tables of plasmonic and dielectric resonances contain essential information to characterize and engineer the electromagnetic scattering of small objects. Specifically, the frequency shift of both plasmonic and dielectric resonances is a quadratic function of the size parameter at the quasistatic resonance, whose prefactor is −χ (2) h /2 χ h for plasmonic modes and κ (2) h /κ ⊥ h for dielectric resonances, respectively. Furthermore, the radiative Q-factor is an inverse power function of the size parameter 1/ ∝ x n i at the resonance whose exponent is exactly the order n i , while the prefactor is the ratio h . In this manuscript, the electrostatic eigenvalue problem 1 is solved by the surface integral method outlined in Refs. 12,37 . The magnetoquasistatic eigenvalue problem 1 is solved by the numerical method outlined in Refs. 5 by using loop basis functions. Then, the radiation correction for both plasmonic (Eqs. 26, 32, and 33) and dielectric eigenvalues (Eqs. 40, 45, and 48) are computed using standard quadrature formulas, and, if singular, using the formulas provided by Graglia 55,56 .

A. Sphere
The plasmonic and dielectric resonances of a sphere of radius R are now investigated. It is assumed that l c = R. The sphere is the ideal shape to numerically validate the radiation corrections, because they also have analytic expression.
In particular, the EQS eigenvalues and corresponding corrections of a sphere are: where n = 1, 2, 3, . . . is the multipolar order. The MQS modes are divided in two sets. The first set is constituted by the A ⊥ -modes. Their eigenvalues and the corresponding corrections are the remaining modes have a non-vanishing normal component of the vector potential on ∂Ω, where r n,l denotes the l-th zero of the spherical Bessel function j n . The formulas in Eqs. 69, 70, and 71 can be also extrapolated by the perturbing the denominators of the Mie coefficients, and they also coincide with Pad expansion of the Mie coefficients found by D. C. Tzarouchis and A. Sihvola in Refs. 57,58 In Table I, the radiative corrections of EQS eigenvalues numerically obtained by Eqs. 26, 32, and 33 are compared against their analytic value in Eq. 69, with a resulting error below 3%. An analogous comparison is also done for in Table II for MQS eigenvalues, where the numerical values obtained by Eqs. 40, 45, and 48 are compared against their analytic counterpart Eqs. 70-71, obtaining an error below 5%.

Table of plasmonic resonances
The table of plasmonic resonances of a sphere is shown in Fig. 1. The first three EQS modes j 1−3 of a sphere are degenerate, in Fig. 1 (a) an exemplification of them is shown. They are bright and represent three electric dipoles oriented along mutually orthogonal directions. The second order correction χ  the sum of two contributions associated to the first and second term (from left to right) of Eq. 26, whose values are 1.6 and −3.5. The latter, prevailing, contribution is associated to the quasistatic self-coupling of j 1−3 , while the former, weaker, term is associated to the radiative self-coupling of the surface charge densities. The third order correction χ 1−3 = +2.92 i , arises from the destructive interference of the two terms in Eq. 32. The first one from the left represents the radiative self-coupling of the surface charges through by the scalar potential (−1.0), while the second, prevailing, term is proportional to the net-dipole moment of j 1−3 (+3.0).
The next five EQS modes j 4−8 are degenerate. They have quadrupolar character, thus they are dark. One exemplification is shown in Fig. 1 (b). As for the dipole modes j 1−3 , the second order correction χ Therefore, to get information on the radiative Q-factor of j 4−8 , it is mandatory to calculate the next imaginary correction, namely χ (5) h through Eq. 33. This expression is greatly simplified because the terms [[j n,h |∆r 2 |j n,h ]] ∂Ω and j k |j h Ω vanish (72) the two terms in Eq. 72 have values of opposite sign, namely −0.055 i and +0.137 i, respectively.

Table of dielectric resonances
The table of dielectric resonances of a sphere is shown in Fig. 2. The fundamental MQS modes j ⊥ 1−3 of the sphere, shown in Fig. 2 (a), are three degenerate magnetic dipoles oriented along three orthogonal axis. The current modes j ⊥ 1−3 are A ⊥ -modes since they generate a vector potential with zero normal component to ∂Ω. Thus, the second order correction κ (2) 1−3 has the simplified expression given by Eq. 41, because the term j k |∆r −1 |j ⊥ 1−3 Ω vanishes ∀k, and j ⊥ 1−3 do not couple with any EQS modes j k . The second order correction depends on the radiative self-coupling of j ⊥ 1−3 through the vector potential. The third order correction κ Then, the three degenerate modes j ⊥ 4−6 are shown in Fig. 2 (b), are referred to as electric toroidal dipoles in the nano-optics community 59 . Contrarily to the previous modes, each of the modes j ⊥ 4−6 generates a vector potential with a non-vanishing normal component on ∂Ω, enabling the coupling between MQS and EQS modes, namely j ⊥ k |∆r −1 |j 4−6 Ω = 0. In particular, j ⊥ 4−6 couple with the three EQS dipole modes j 1−3 , shown in Fig. 1  (a), whereas the coupling with the remaining EQS modes is zero. Both the contributions to κ  It arises from four contribution the third one from the left is negligible, the remaining contributions contribute with the same sign.
Next, the degenerate modes are j ⊥ 7−11 are shown in Fig. 2 (c). As already pointed out in Ref. 5 the modes j ⊥ 7−11 have the same MQS eigenvalue of electric toroidal dipoles 5 j ⊥ 4−6 . Nevertheless, differently from the electric toroidal dipoles they are A ⊥ -modes. For this reason the second order correction κ (2) 7−11 has the simplified expression shown in Eq. 41. The symmetry of the current mode implies j ⊥ h |∆r 2 |j ⊥ h Ω = 0 and thus a vanishing third order correction. The fifth order correction κ (5) 7−11 can be easily calculated by 50, where, the second term from the left, which accounts for the coupling among different MQS modes is also negligible.
Similar considerations apply also to the degenerate modes j ⊥ 12−16 , shown in Fig. 2 (d).

Dipole Excitation
The sphere is now excited by an electric dipole. Specifically, the sphere of radius R is centered in the origin, while the exciting dipole is oriented alongx and it is positioned in r d = (0, 0, 1.5R), namely Under the same excitation conditions, two different scenarios are investigated. In the first one, shown in Fig. 3, the sphere is made of a Drude metal, in the second one, shown in Fig. 4, of a high/moderate-index dielectric.
In both cases, low-losses are assumed: this hypothesis is essential for a quantitative comparison between the estimated Q-factor and the peaks broadening, which would have been otherwise dominated by non-radiative losses. The absorbed power spectrum W abs is calculated by the Mie analytic solution 60 , and the translation-addition theorem for vector spherical wave functions (VSWFs) 61,62 is used to translate the exciting dipole of Eq. 73 into the corresponding VSWFs set centered in the sphere's center. The reason behind the choice of the dipole as excitation and of the absorbed power as physical observable is that without these two hypothesis, some modes would not be excited or probed. The radiative shift of the peaks of W abs and their Qfactors are then investigated. In particular, the plasmonic and dielectric resonance frequencies are compared against the frequenciesω h , associated to the curve peaks. Similarly, Q-factors of plasmonic and dielectric modes are validated against the corresponding heuristic Q-factors, given by the ratio of the resonance frequencyω h to the width ∆ω FWHM,h of the resonance curve between two points, at the either side of the resonance, where the ordinate is the half of the maximum absorbed power, namely the full width at half maximum (FWHM) 63 In the W abs spectra of Fig. 3,4 a segment joining the two ordinates at half maximum is also shown. Metal Sphere. A metal sphere is investigated in Fig.  3, assuming a low-loss Drude metal with ν = 10 −4 ω p . The absorbed power spectra W abs are evaluated as a function of x/x p = ω/ω p for three different values of x p ,   1 (c). The normalization adopted in Fig. 3 makes the corresponding curves independent of the particular value of ω p . Nevertheless, it is useful to contextualize the chosen values of x p to actual materials: for a gold sphere 64 with ω p ≈ 6.79 T rad/s, they correspond to R = 4.5nm (a), R = 22nm (b), and R = 45nm (c). The positions predicted by Eq. 56 for first three sets of degenerate EQS modes, u 1−3 , u 4−8 , and u 9−15 , are shown with vertical dashed lines.
Based on the information contained in Fig. 1, the ratio between the second order correction and the EQS eigenvalue is significantly larger for the u 1−3 than for u 4−8 or u 9−15 , thus the modes u 1−3 experiences the largest shift as x p increases. Moreover, the lowest non-vanishing imaginary correction n i = 3 for u 1−3 , 5 for u 4−8 , and 7 for u 9−15 ; thus the modes u 1−3 also undergoes the largest radiative broadening.
In Fig. 3 (a) the radius R is small compared to the plasma wavelength λ p , and the EQS approximation works well: Eq. 51 exactly predicts the occurrence of the W abs peaks.
In Fig. 3 (b) x p is increased to 0.5, and the W abs peaks undergo a broadening and deviate from their quasistatic position, especially the peaks associated to u 1−3 . Nevertheless, the plasmonic resonances obtained through Eq. 56, which incorporates the radiation corrections, accurately predict the occurrence of the W abs peaks.  In Table III, the resonance frequencies ω h are compared against the corresponding peak positionsω h and the Q-factors against their heuristic counterparts. As expected, for x p = 0.5 the Q-factor of the dipole modes u 1−3 is the lowest one, being, unlike the others, dominated by radiative losses. For the investigated peaks Eqs. 58, 59, and 60 correctly predict the Q-factors.
Eventually, in Fig. 3, R is comparable to the plasma wavelength λ p , namely x p = 1. The peaks experience a further shift, nevertheless thanks to the radiation corrections, Eq. 56 is still able to accurately locate the resonances with an error < 1.5%, as also quantified by Tab. III. Also the Q-factors of the modes u 1−3 and u 4−8 are predicted with good accuracy. Now, they are both dominated by radiative losses, Q with color matching the box of the corresponding modes of the table of dielectric resonances in Fig. 2.
Based on the information contained in Fig. 2, the ratio between the second order correction and the MQS eigenvalue is larger for the v 1−3 than for v 4−6 or v 7−11 , thus the modes v 1−3 experiences the largest shift as the x p increases. Moreover, the lowest non-vanishing imaginary correction n i = 3 for v 1−3 , n i = 5 for v 4−7 and for v 7−11 thus the magnetic dipole modes v 1−3 undergoes the largest radiative broadening.
In Fig. 4 (a), the size parameter is small, i.e. x ∈ [0.025, 0.065] 1, and the MQS approximation works well: Eq. 61 exactly predicts the occurrence of the W abs peaks.
Next, in Fig. 4 (b) x ∈ [0.25, 0.65] < 1 and χ = 99 + 0.01i: the peaks position red-shifts against their MQS position, but Eq. 64, which includes the radiation corrections, predicts their occurrence. In Tab. IV the expected resonance frequencies and Q-factors are compared against the peak positions and the corresponding heuristic Q-factors. The radiative damping determines the broadening of the magnetic dipole modes v 1−3 , while non-radiative mechanisms dominate the remaining peaks. The agreement of the overall Q-factors Q ⊥ h with its heuristic counterpartsQ h is good.
Eventually, in Fig. 4 (c), a silicon sphere with χ = 14.45−0.1456i is investigated in the range x ∈ [0.64, 1.6]. The shift of the peaks against the MQS position is significant, but including the radiation corrections, Eq. 61 predicts their occurrence with an error ≤ 2%. In Tab. IV shows that the radiative damping causes of the broadening of the first and third peaks, and the values of the Q r h calculated by Eq. 58 are very close to their heuristic counterparts. The second peak is not visible anymore. The table of plasmonic resonances of a cylinder with H = R is presented in Fig. 5. It is assumed that l c = R. In Fig. 5 (a) the first set of two degenerate EQS current density modes j 1−2 is described. They represent two electric dipoles oriented along mutually orthogonal directions, thus they are bright. The second order correction χ (2) 1−2 = −3.94 results from the competition of the two terms in Eq. 26, which account for the quasistatic selfcoupling of j 1−2 through the vector potential (−5.68), and for the radiative self-coupling of the surface charge densities through the scalar potential (1.74). Also the third order correction χ 1−2 = +2.92 i arises from the interplay of two competing terms in Eq. 32, accounting for the radiative self-coupling of the surface charge density (−1.47i), and the magnitude of the net-dipole moment of j 1−2 (+4.39i).
The next EQS modes are two degenerate couples, namely j 3−4 , and j 5−6 shown in Fig. 5 (b)-(c). They have quadrupolar character with zero dipole moment, thus they are dark. In both cases, the third order corrections in Eq. 32 vanishes: in particular j h |j h Ω = 0 because j 3−4 and j 5−6 are dark, while [[j n,h |∆r 2 |j n,h ]] ∂Ω = 0 because of their symmetry. Therefore, to obtain information on the radiative Q-factor, it is mandatory to calculate the fifth order correction χ (5) h . It can be calculated by Eq. 72 because [[j n,h |∆r 2 |j n,h ]] ∂Ω = 0 and j k |j h Ω = 0. In both cases, the two terms in Eq. 72 contribute with opposite signs, and the latter one, associated to the radiative self-coupling of the EQS current modes, is dominant.
In Fig. 5 (d) the EQS mode j 7 is shown. This mode is associated to a vertical electric dipole, thus it is bright. The second order correction χ

Table of Dielectric Resonances
The table of dielectric resonances of the finite size cylinder is presented in Fig. 6. The fundamental MQS mode j ⊥ 1 , shown in Fig. 6 (a), is a magnetic dipole oriented along the vertical axis. The current mode j ⊥ 1 generates a vector potential with zero normal component to ∂Ω, thus it is an A ⊥ -mode. As a consequence the term j ⊥ 1 does not couple with any EQS mode, i.e. j k |∆r −1 |j ⊥ 1 Ω = 0 ∀k, and the second order correction κ (2) 1 is given by Eq. 41. The third order correction is given by Eq. 45 and it is connected to the radiative selfcoupling of j ⊥ 1 through by the vector potential. The next couple of degenerate modes j ⊥ 2 , j ⊥ 3 , shown in Fig. 6 (b), are magnetic dipoles oriented along two orthogonal horizontal directions. Analogously to j ⊥ 1 , they are A ⊥ -modes, and similar considerations apply.
The next degenerate modes j ⊥ 4−5 are shown in Fig.  6 (c). They are known as electrical toroidal dipoles within the nano-optics community, and as HEM 12δ within the antenna community 4 . Contrarily to the previous modes, the currents j ⊥ 4−5 generate vector potentials with a non-vanishing normal component on ∂Ω, enabling the coupling between MQS and EQS modes, namely j ⊥ k |∆r −1 |j 4−5 Ω = 0. In particular, the MQS modes j ⊥ 4−5 couple with the two horizontal EQS dipoles j 1 , j 2 , shown in Fig. 5 (a), whereas the coupling with the remaining EQS modes is negligible. Both the contributions to κ (2) 4−5 shown in Eq. 40 are thus different from zero and have the same sign. Due to the symmetry, the term j ⊥ h |∆r 2 |j ⊥ h Ω is zero and the third order correction vanishes. Thus, the calculation of κ (5) 4−5 by Eq. 48 is required to retrieve information about the radiative Q-factor. Among the four contributions in 48 the third from the left is negligible, the remaining contributions contribute with the same sign, and the last one is the dominant one.
The next mode is j ⊥ 6 , namely a vertical electrical toroidal dipole, or TM 01δ , is shown in Fig. 6 (d). It generates electric fields with a non-vanishing normal component on ∂Ω, and exhibits a non-negligible coupling with the vertical EQS dipole j 7 , j ⊥ 7 |∆r −1 |j 4−5 Ω = 0. Its properties are analogous to those of j ⊥ 4−5 , and the values of the second and fifth order correction are reported in Fig. 6.
The next couple of degenerate modes j ⊥ 7−8 , shown in Fig. 6 (e), are A ⊥ -modes. Therefore, the second order correction is given by Eq. 41. The symmetry of the current mode implies j ⊥ h |∆r 2 |j ⊥ h Ω = 0 and thus a vanishing third order correction. The fifth order correction κ (5) 4−5 can be calculated by Eq. 50, where, the second term from the left, which accounts for the coupling among different MQS modes, is negligible.

Dipole Excitation
The cylinder is excited by an electric dipole. Specifically, the cylinder of radius R and height H = R is centered in the origin and oriented alongẑ, while the exciting dipole is oriented alongŷ and it is positioned in r d = (R, 0, 1.5R), namely Assuming the same geometry and excitation conditions, two different scenarios are investigated. In the first one, the cylinder is assumed to be made of a Drude metal, in the second of a high/moderate-index dielectric. The absorbed power W abs spectrum is calculated by an in-house full-Maxwell numerical method based on the Surface Integral Equation Method 65 . The frequency shift and the broadening of the resonances are then investigated. Metal cylinder In Fig. 7 the first scenario is investigated, assuming a low loss Drude metal with ν = 10 −3 ω p . The absorbed power W abs is shown as a function of x/x p = ω/ω p for three different values of x p , namely (a) 0.1, (b) 0.5, and (c) 1. For gold cylinders, by assuming, as in Ref. 64 , ω p ≈ 6.79 Trad/s, the three values of x p correspond to R = 4.5nm (a), R = 22nm (b), and R = 45nm (c). The expected resonance positions (obtained by Eq. 56) of the first three plasmonic modes, namely u 1−2 , u 3−4 , and u 5−6 , are shown with vertical dashed lines.
Based on the table of plasmonic modes, the ratio between the second order correction and the EQS eigenvalue is significantly larger for the u 1−2 than for u 3−4 or u 5−6 , thus the modes u 1−2 experiences the largest shift as x p increases. Moreover, the lowest non-vanishing imaginary correction n i = 3 for u 1−2 , 5 for u 3−4 and u 5−6 ; thus the modes u 1−3 also undergoes the largest radiative broadening.
In Fig. 7 (a) R λ p , the EQS resonances obtained by Eq. 51 accurately predict the occurrence of the W abs peaks. The non-radiative losses cause the broadening of all peaks.
In Fig. 7 (b), by increasing the size parameter to x p = 0.5, the peaks of W abs start to shift, especially the one associated to the horizontal electric dipoles u 1−2 . Nevertheless, including the radiation correction, Eq. 56 predicts their occurrence with a relative error less than 0.5%, as shown in Table V. In this table, the Q-factors obtained by Eqs. 58, 59, 60 are also compared against their heuristic counterparts, defined in Eq. 74. The Qfactor of the the modes u 1−2 is limited by their radiative losses, while the Q-factors of the modes u 3−4 and u 5−6 are limited by non-radiative damping mechanisms. This is because, accordingly to the In Fig. 7 (c), the cylinder has radius comparable to the plasma wavelength. Despite the further shift of the peaks, the radiation corrections are still able to correctly locate the resonances, with an error less than 0.5% (see Tab. V). Moreover, Tab. V shows that while the broadening of the first peak arises from the radiative damping, both radiative and non-radiative damping contribute to the broadening of the second and third peak. The total Q-factor is also in good agreement with the heuristic Q-factor in the investigated cases.
High/moderate index cylinder In Fig. 8 the power absorbed by a cylinder constituted by a non-dispersive high/moderate-index dielectric with low-losses is investigated. The absorbed power W abs is shown as a func- In Fig. 8 (a), the size parameter of the cylinder is small, i.e. x ∈ [0.0275, 0.055] 1, and the MQS approximation accurately predicts the occurrence of the peaks through Eq. 61. Next, in Fig. 8 (b), it is assumed that χ = 99 − 0.01i and x ∈ [0.0275, 0.55] < 1: the peaks undergo a shift from their MQS positions, but, including the radiation corrections, Eq. 64 predicts their occurrence with an error < 0.2%. The resonances and Q-factors of v 1 , v 2−3 , and v 7−8 are compared in Tab. VI against the actual peak positions and heuristic Q-factors. Good agreement is found. Moreover, in this case, the Q-factors of v 1 , v 2−3 are dominated by radiative losses, while for v 7−8 non-radiative losses play also a role.
Eventually, in Fig. 8 (c), a silicon cylinder with x ∈ [0.72, 1.44] and χ = 14.45 − 0.1456i is considered. The size parameter is now of the order of one and the peaks undergo a significant shift and broadening. Nevertheless, Eq. 64 is able to predict the peaks occurrence with an error lower than 2.5%, as shown in Tab. VI. In this same table, the corresponding Q-factors, all dominated by radiative losses, are also listed and are very close to their heuristic counterpart.

VIII. CONCLUSIONS
Maxwell's equations provide an exhaustive description of classical electromagnetic phenomena, from the simplest to the most sophisticated ones. However, in many applications wave phenomena occurring at short time scales are of no practical concern, and the fields may be described by either their magnetoquasistatic or the electroquasistatic approximation, i.e. the approximations behind the description of capacitors and inductors. This is also the case for light scattering: resonances in metal or high-index objects, assumed much smaller than the vacuum wavelength, may be respectively described by the electroquasistatic or the magnetoquasistatic approximation. Unfortunately, both approximations are unable to predict the frequency-shift and radiative Q-factors, which arise from the coupling with the radiation.
In this paper, closed form expressions for the radiation corrections to the real and imaginary part of both electroquasistatic and magnetoquasistatic eigenvalues and of the corresponding modes are derived. These corrections only depend on the quasistatic current mode distributions.
The expression of the radiation corrections may be greatly simplified if the electroquasistatic mode is dark, or if the magnetoquasistatic mode generates a vector potential with vanishing normal component to the surface of the object (A ⊥ -modes).
The frequency shift of any mode is a quadratic function of the size parameter, whose prefactor depends on the ratio between the second order correction and the quasistatic eigenvalue. The radiative Q-factor is an inverse power function of the size parameter whose exponent is the order n i of the first non-vanishing imaginary correction, while the prefactor is the ratio between the the static eigenvalue and its n i -th order imaginary correction.
The resonances of a small objects can be then naturally classified in two tables, denoted as table of plasmonic resonances and table of dielectric resonances, containing the essential information to analyze and engineer the electromagnetic scattering from small objects. In these tables the resonances are sorted accordingly to their real quasistatic eigenvalue, and characterized by the second order correction and by the non-vanishing imaginary correction of the lowest order, i.e. n i . All the quantities contained in these tables do not depend neither on the size nor on the permittivity, but only on the quasistatic mode morphology.
The introduced expressions for the resonance frequency and Q-factor are successfully validated by predicting the resonance peaks and their broadening in the absorption spectra of a sphere and a finite-size cylinder. Both Drude-metals and high/moderate index dielectric are considered.