Applying the resonant-state expansion to realistic materials with   frequency dispersion

Hame Sehmi; Wolfgang Langbein; Egor Muljarov

arXiv:1906.07007·physics.optics·January 29, 2020

Applying the resonant-state expansion to realistic materials with frequency dispersion

Hame Sehmi, Wolfgang Langbein, Egor Muljarov

PDF

TL;DR

This paper extends the resonant-state expansion method to realistic dispersive materials like metals and semiconductors, enabling detailed analysis of their resonant states, surface modes, and lasing phenomena.

Contribution

It applies the dispersive resonant-state expansion to complex materials using a generalized Drude-Lorentz model, revealing new insights into surface plasmon polaritons and lasing states.

Findings

01

Resonant states are accurately calculated for dispersive materials.

02

Surface plasmon polariton modes evolve with nanoparticle size.

03

Lasing resonant states are characterized near the band gap.

Abstract

The dispersive resonant-state expansion, developed for an accurate calculation of the resonant states in open optical systems with frequency dispersion, is applied here to realistic materials, such as metallic nanoparticles and semiconductor microspheres. The material permittivity is determined by fitting the measured indices of refraction and absorption with a generalized Drude-Lorentz model containing a number of poles in the complex frequency plane. Each Drude or Lorentz pole generates an infinite series of resonant states. Furthermore, for small nanoparticles, each of these poles produces a distinct surface plasmon polariton mode. The evolution of these multiple surface modes with increasing radius traces the transition from the electrostatic limit to significant retardation and radiation. Treating the optical phonon range in a semiconductor microsphere, a reststrahlen band…

Tables1

Table 1. Table 1: Optimized model parameters for Au and GaAs phonon (*) and band-edge (**) ranges, using the fit function with optimization energy ranges corresponding to the data shown in Figs. 1 , 3 , 5 , and 6 . The number of data values N d subscript 𝑁 𝑑 N_{d} and the resulting error S 𝑆 S are also given.

Material	Au	GaAs*	GaAs**
$ ε_{\infty}$	$0.5$	$11.0$	$8.6013$
$γ (eV)$	$- 0.065748$	$-$	$-$
$σ (eV)$	$1133.0$	$-$	$-$
$ Ω_{1}^{'} (eV)$	$2.5936$	$0.033314$	$1.497$
$Ω_{1}^{''} (eV)$	$- 0.41875$	$1.4904 \times 10^{- 4}$	$- 0.03665$
$σ_{1}^{'} (eV)$	$1.4029$	$0.0$	$0.0$
$σ_{1}^{''} (eV)$	$0.76857$	$0.033262$	$0.01224$
$ Ω_{2}^{'} (eV)$	$3.8192$	$-$	$1.5612$
$Ω_{2}^{''} (eV)$	$- 1.3246$	$-$	$- 0.05643$
$σ_{2}^{'} (eV)$	$0.41939$	$-$	$0.0$
$σ_{2}^{''} (eV)$	$4.5468$	$-$	$0.02432$
$ Ω_{3}^{'} (eV)$	$9.6899$	$-$	$1.6463$
$Ω_{3}^{''} (eV)$	$- 4.2933$	$-$	$- 0.0457$
$σ_{3}^{'} (eV)$	$0.012244$	$-$	$0.0$
$σ_{3}^{''} (eV)$	$14.817$	$-$	$0.02404$
$ Ω_{4}^{'} (eV)$	$-$	$-$	$2.2853$
$Ω_{4}^{''} (eV)$	$-$	$-$	$- 0.00778$
$σ_{4}^{'} (eV)$	$-$	$-$	$0.0$
$σ_{4}^{''} (eV)$	$-$	$-$	$2.9302$
$ ℏ ω_{1} (eV)$	$0.64$	$0.031$	$1.3$
$ℏ ω_{N} (eV)$	$6.6$	$0.0372$	$1.65$
$N_{d}$	$49$	$47$	$8$
$S$	$1.4795$	$0.0372$	$0.0055$

Equations107

\nabla \times \nabla \times E_{n} (r) = \frac{ω _{n}^{2}}{c ^{2}} \hat{ε} (r, ω_{n}) E_{n} (r)

\nabla \times \nabla \times E_{n} (r) = \frac{ω _{n}^{2}}{c ^{2}} \hat{ε} (r, ω_{n}) E_{n} (r)

\hat{ε} (r, ω) = \hat{ε}_{\infty} (r) + j \sum \frac{i σ ^ _{j} ( r )}{ω - Ω _{j}},

\hat{ε} (r, ω) = \hat{ε}_{\infty} (r) + j \sum \frac{i σ ^ _{j} ( r )}{ω - Ω _{j}},

\nabla \times \nabla \times E (r) = \frac{ω ^{2}}{c ^{2}} [\hat{ε} (r, ω) + Δ \hat{ε} (r, ω)] E (r)

\nabla \times \nabla \times E (r) = \frac{ω ^{2}}{c ^{2}} [\hat{ε} (r, ω) + Δ \hat{ε} (r, ω)] E (r)

Δ \hat{ε} (r, ω) = Δ \hat{ε}_{\infty} (r) + j \sum \frac{i Δ σ ^ _{j} ( r )}{ω - Ω _{j}}

Δ \hat{ε} (r, ω) = Δ \hat{ε}_{\infty} (r) + j \sum \frac{i Δ σ ^ _{j} ( r )}{ω - Ω _{j}}

E (r) = - \frac{ω ^{2}}{c ^{2}} \int \hat{G}_{ω} (r, r^{'}) Δ \hat{ε} (r^{'}, ω) E (r^{'}) d r^{'} .

E (r) = - \frac{ω ^{2}}{c ^{2}} \int \hat{G}_{ω} (r, r^{'}) Δ \hat{ε} (r^{'}, ω) E (r^{'}) d r^{'} .

E (r) = n \sum c_{n} E_{n} (r),

E (r) = n \sum c_{n} E_{n} (r),

\hat{G}_{ω} (r, r^{'}) = c^{2} n \sum \frac{E _{n} ( r ) \otimes E _{n} ( r ^{'} )}{ω _{n} ( ω - ω _{n} )} .

\hat{G}_{ω} (r, r^{'}) = c^{2} n \sum \frac{E _{n} ( r ) \otimes E _{n} ( r ^{'} )}{ω _{n} ( ω - ω _{n} )} .

1

1

ω_{n} m \sum (δ_{nm} - U_{nm}) c_{m} = ω m \sum (δ_{nm} + V_{nm}) c_{m},

ω_{n} m \sum (δ_{nm} - U_{nm}) c_{m} = ω m \sum (δ_{nm} + V_{nm}) c_{m},

V_{nm} = \int E_{n} (r) \cdot Δ \hat{ε}_{\infty} (r) E_{m} (r) d r

V_{nm} = \int E_{n} (r) \cdot Δ \hat{ε}_{\infty} (r) E_{m} (r) d r

U_{nm}

U_{nm}

= j \sum \frac{i}{ω _{n} - Ω _{j}} \int E_{n} (r) \cdot Δ \hat{σ}_{j} (r) E_{m} (r) d r,

ω_{n} = Ω_{j} + ξ q_{n},

ω_{n} = Ω_{j} + ξ q_{n},

\hat{ε} (r, ω_{n}) = \hat{ε}_{\infty} (r) + j^{'} \neq = j \sum \frac{i σ ^ _{j^{'}} ( r )}{Ω _{j} - Ω _{j^{'}}} + \frac{s ^ _{j} ( r )}{q _{n}},

\hat{ε} (r, ω_{n}) = \hat{ε}_{\infty} (r) + j^{'} \neq = j \sum \frac{i σ ^ _{j^{'}} ( r )}{Ω _{j} - Ω _{j^{'}}} + \frac{s ^ _{j} ( r )}{q _{n}},

\frac{\partial ε ^ ( r , ω )}{\partial ( ω )}_{ω_{n}} = - \frac{1}{ξ} \frac{s ^ _{j} ( r )}{q _{n}^{2}},

\frac{\partial ε ^ ( r , ω )}{\partial ( ω )}_{ω_{n}} = - \frac{1}{ξ} \frac{s ^ _{j} ( r )}{q _{n}^{2}},

ξ = - \frac{Ω _{j}}{q _{n}^{2}} \int_{V} E_{n} (r) \cdot \hat{s}_{j} (r) E_{n} (r) d r,

ξ = - \frac{Ω _{j}}{q _{n}^{2}} \int_{V} E_{n} (r) \cdot \hat{s}_{j} (r) E_{n} (r) d r,

\alpha_{n}=\left\{\begin{array}[]{cl}1&\text{for non-pRSs}\\ &\\ \displaystyle\sqrt{\frac{\omega_{n}-\Omega_{j}}{\Omega_{j}}}&\text{for pRSs of $\Omega_{j}$}\\ \end{array}\right.

\alpha_{n}=\left\{\begin{array}[]{cl}1&\text{for non-pRSs}\\ &\\ \displaystyle\sqrt{\frac{\omega_{n}-\Omega_{j}}{\Omega_{j}}}&\text{for pRSs of $\Omega_{j}$}\\ \end{array}\right.

E_{n} (r) = α_{n} \tilde{E}_{n} (r),

E_{n} (r) = α_{n} \tilde{E}_{n} (r),

V_{nm}

V_{nm}

ω_{n} m \sum (δ_{nm} - S_{nm}) b_{m} = ω m \sum (δ_{nm} + α_{n}^{2} Q_{nm}) b_{m},

ω_{n} m \sum (δ_{nm} - S_{nm}) b_{m} = ω m \sum (δ_{nm} + α_{n}^{2} Q_{nm}) b_{m},

Q_{nm} = \int \tilde{E}_{n} (r) \cdot Δ \hat{ε}_{\infty} (r) \tilde{E}_{m} (r) d r,

Q_{nm} = \int \tilde{E}_{n} (r) \cdot Δ \hat{ε}_{\infty} (r) \tilde{E}_{m} (r) d r,

S_{nm} = \int \tilde{E}_{n} (r) \cdot [Δ \hat{ε} (r, ω_{n}) - Δ \hat{ε}_{\infty} (r)] \tilde{E}_{m} (r) d r,

S_{nm} = \int \tilde{E}_{n} (r) \cdot [Δ \hat{ε} (r, ω_{n}) - Δ \hat{ε}_{\infty} (r)] \tilde{E}_{m} (r) d r,

S_{nm} = \frac{i}{Ω _{j}} \int \tilde{E}_{n} (r) \cdot Δ \hat{σ}_{j} (r) \tilde{E}_{m} (r) d r,

S_{nm} = \frac{i}{Ω _{j}} \int \tilde{E}_{n} (r) \cdot Δ \hat{σ}_{j} (r) \tilde{E}_{m} (r) d r,

E (r) = n \sum b_{n} \tilde{E}_{n} (r)

E (r) = n \sum b_{n} \tilde{E}_{n} (r)

ω = ω_{n} \frac{1 - U _{nn}}{1 + V _{nn}}

ω = ω_{n} \frac{1 - U _{nn}}{1 + V _{nn}}

U_{nm} = j = \pm 1 \sum [\frac{i}{ω _{n} - Ω _{j}} - \frac{i}{ω _{n} - Ω _{j}}] \int E_{n} (r) \cdot \hat{σ}_{j} (r) E_{m} (r) d r,

U_{nm} = j = \pm 1 \sum [\frac{i}{ω _{n} - Ω _{j}} - \frac{i}{ω _{n} - Ω _{j}}] \int E_{n} (r) \cdot \hat{σ}_{j} (r) E_{m} (r) d r,

ε (ω) = ε_{\infty} + j = 1 \sum L (\frac{i σ _{j}}{ω - Ω _{j}} + \frac{i σ _{j}^{*}}{ω + Ω _{j}^{*}}),

ε (ω) = ε_{\infty} + j = 1 \sum L (\frac{i σ _{j}}{ω - Ω _{j}} + \frac{i σ _{j}^{*}}{ω + Ω _{j}^{*}}),

ε (ω) = ε_{\infty} - \frac{γ σ}{ω ( ω + iγ )} + j = 1 \sum L (\frac{i σ _{j}}{ω - Ω _{j}} + \frac{i σ _{j}^{*}}{ω + Ω _{j}^{*}}),

ε (ω) = ε_{\infty} - \frac{γ σ}{ω ( ω + iγ )} + j = 1 \sum L (\frac{i σ _{j}}{ω - Ω _{j}} + \frac{i σ _{j}^{*}}{ω + Ω _{j}^{*}}),

\frac{iσ}{ω} - \frac{iσ}{ω + iγ} = - \frac{γ σ}{ω ( ω + iγ )} .

\frac{iσ}{ω} - \frac{iσ}{ω + iγ} = - \frac{γ σ}{ω ( ω + iγ )} .

\mathbf{E}(\mathbf{r})=\dfrac{A_{nl}}{n_{r}k_{n}r}\left(\begin{array}[]{ccc}l(l+1)\,\psi_{nl}(r)Y_{lm}(\theta,\varphi)\\[5.0pt] \dfrac{\partial}{\partial r}r\,\psi_{nl}(r)\dfrac{\partial}{\partial\theta}Y_{lm}(\theta,\varphi)\\ \dfrac{\partial}{\partial r}r\,\psi_{nl}(r)\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\varphi}Y_{lm}(\theta,\varphi)\\ \end{array}\right),

\mathbf{E}(\mathbf{r})=\dfrac{A_{nl}}{n_{r}k_{n}r}\left(\begin{array}[]{ccc}l(l+1)\,\psi_{nl}(r)Y_{lm}(\theta,\varphi)\\[5.0pt] \dfrac{\partial}{\partial r}r\,\psi_{nl}(r)\dfrac{\partial}{\partial\theta}Y_{lm}(\theta,\varphi)\\ \dfrac{\partial}{\partial r}r\,\psi_{nl}(r)\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\varphi}Y_{lm}(\theta,\varphi)\\ \end{array}\right),

ψ_{n l} (r) = \frac{j _{l} ( n _{r} k _{n} r )}{j _{l} ( n _{r} k _{n} R )} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Applying the resonant-state expansion to realistic materials with frequency dispersion

H. S. Sehmi

W. Langbein

E. A. Muljarov

School of Physics and Astronomy, Cardiff University, Cardiff CF24 3AA, United Kingdom

Abstract

The dispersive resonant-state expansion, developed for an accurate calculation of the resonant states in open optical systems with frequency dispersion, is applied here to realistic materials, such as metallic nanoparticles and semiconductor microspheres. The material permittivity is determined by fitting the measured indices of refraction and absorption with a generalized Drude-Lorentz model containing a number of poles in the complex frequency plane. Each Drude or Lorentz pole generates an infinite series of resonant states. Furthermore, for small nanoparticles, each of these poles produces a distinct surface plasmon polariton mode. The evolution of these multiple surface modes with increasing radius traces the transition from the electrostatic limit to significant retardation and radiation. Treating the optical phonon range in a semiconductor microsphere, a reststrahlen band separating the resonant states is found. Considering a small energy range around the semiconductor band gap, the transition from absorption to gain is described by inverting the Lorentz pole weight, which results in the formation of lasing resonant states. Interestingly, the series of resonant states converging towards the absorption pole from the lower frequency side reshapes for a gain pole into a clockwise loop approaching the pole from the higher frequency side, being separated from a series spanning from low to high frequencies and containing the lasing modes.

I Introduction

Applying the concept of resonant states (RSs) to optical systems is a mathematically rigorous way of treating their electromagnetic resonances. RSs are the eigenmodes of an optical system which can be found by solving Maxwell’s wave equation with outgoing wave boundary conditions. In open systems, the RS eigenfrequencies $\omega_{n}$ , numbered by the index $n$ , are generally complex, which physically reflects the fact that the energy leaks out of the system. For an isolated resonance in an optical spectrum of the system, the real part of the corresponding $\omega_{n}$ determines the center frequency of the resonance and the imaginary part its half width at half maximum, with the quality factor of the resonance given by half of the ratio of real to imaginary part.

Based on the concept of RSs, a rigorous perturbative method of calculation of the RSs of an arbitrary finite optical system (called perturbed system) has been developed Muljarov et al. (2010). This method, called the resonant-state expansion (RSE), converts the Maxwell wave equation into a linear matrix eigenvalue problem using the RSs of a simpler and usually analytically solvable system (called unperturbed system) as a basis for expansion. In this way, the RSE is not limited to small perturbations and finds numerically exactly all the eigenmodes of the perturbed system. The RSE has been recently generalized Muljarov and Langbein (2016a) to optical systems with frequency dispersion, such as metal or semiconductor nanoparticles. This dispersive RSE (dRSE) is suited for materials with dispersion of the permittivity described by a generalized Drude-Lorentz model Sehmi et al. (2017a). In practice, one can find the parameters of this model – generalized conductivities and complex poles of the permittivity – by fitting the measured indices of refraction and absorption, as we have recently demonstrated for metals Sehmi et al. (2017a) and semiconductors Sehmi et al. (2017b).

The RSE can treat perturbations of arbitrary strength Muljarov et al. (2010), and can be superior to existing computational methods in electrodynamics, such as finite difference in time domain (FDTD) and finite element method (FEM), in terms of accuracy and efficiency Doost et al. (2014); Lobanov et al. (2017). One of the main advantages of the RSE method is that it calculates the full spectrum of RSs in a wide frequency range striving towards completeness, and it does not depend on further approximations, i.e. it is numerically exact. Indeed, the only parameter of the RSE is the number of the RSs included in the basis, which can in principle be made arbitrarily large. While the formalism of the RSE is analytically advanced, its technical implementation is rather straightforward, as it reduces solving Maxwell’s wave equation to a diagonalization of a matrix containing a perturbation. For the latter, one needs a complete set of RSs of an unperturbed system, ideally chosen in such a way that the matrix elements of the perturbation can be efficiently computed.

In the present work we apply the dRSE to metal and semiconductor nanoparticles, using as basis system a non-dispersive dielectric sphere. This approach has been used in Ref. [Muljarov and Langbein, 2016a] for a dispersion described by the Drude model. Here, we extend the method to systems with a dispersion of the permittivity described by up to four pairs of Lorentz poles. Including these poles in the dispersion leads to new physical insights into the properties of metal and semiconductor particles, such as a coexistence of multiple surface plasmon polariton modes, which is of significance for the field of plasmonics Shahbazyan and Stockman (2014). We also explore how the dRSE can describe the onset of lasing when material absorption changes to gain in semiconductor microstructures, important for the field of optical microcavities Vahala (2004); Kavokin et al. (2007) and semiconductor lasers Chow et al. (2012).

The paper is organized as follows. In Sec. II we outline the formalism of the dRSE, providing general expressions for the normalization of dispersive RSs and for linear matrix eigenvalue problem determining the perturbed resonant states. The dRSE formalism is based on the condition that the dispersion of the unperturbed and perturbed systems has the same pole structure. In order to allow the unperturbed system to have no dispersion, or the perturbation to introduce additional poles in the dispersion, we present in Sec. III a variation of the dRSE called infinitesimal dispersive RSE (idRSE), which extends the basis of the unperturbed RSs to include a subset of RSs, called pole RSs, created by the additional poles of the permittivity in the limit of infinitesimal conductivities. We then apply in Sec. IV the dRSE to metal and semiconductor spherically symmetric systems. Using the idRSE, we find in Sec. IV.1 the RSs of a gold sphere described by the Drude-Lorentz model with 3 pairs of Lorentz poles. This is done using the RSs of a dielectric sphere as basis extended to include the pole RSs. We concentrate our study on the surface plasmon polariton RSs and investigate their dependence on the sphere radius changing from 1 to 100 nm. In Sec. IV.2 we consider the optical-phonon range of the permittivity of GaAs, which can be accurately fitted by 1 pair of Lorentz poles, and find the RSs of a 50 $\mu$ m-radius GaAs microsphere. We study a small shift of these Lorentz poles in Sec. IV.3 by using a diagonal approximation of the dRSE which is corresponding to a 1st-order perturbation theory. Finally, in Sec. IV.4 we calculate the lasing modes of GaAs microspheres with 10 $\mu$ m and 1 $\mu$ m radius, having the permittivity fitted close to the GaAs band edge by 4 pairs of classical Lorentz poles. The optical gain is introduced by inverting the weight of the pole pair at the band edge from absorption to gain. Details of the permittivity models and the calculation of basis RSs and their matrix elements are given in Appendices A–C.

II Formalism of the dispersive RSE

The basis for the dRSE is formed from the RSs of an unperturbed system which are the eigensolutions of the Maxwell wave equation

[TABLE]

with outgoing wave boundary conditions (here $c$ is the speed of light in vacuum). The permittivity tensor $\hat{{\boldsymbol{\varepsilon}}}(\mathbf{r},\omega)$ of the unperturbed system is taken as analytic function in the complex frequency plane, with a countable number of simple poles. According to the Mittag-Leffler theorem Arfken and Weber (2013), the frequency dispersion of the permittivity tensor can be expressed as

[TABLE]

where $\hat{{\boldsymbol{\varepsilon}}}_{\infty}(\mathbf{r})$ is the high-frequency value of the permittivity and $\Omega_{j}$ are the resonance frequencies, at which the permittivity has simple poles, with the residues at the poles given by generalized conductivity tensors $\hat{\boldsymbol{\sigma}}_{j}(\mathbf{r})$ . Equation (2) is known as generalized Drude-Lorentz model Sehmi et al. (2017a); more details about the model and its special cases (Drude and Lorentz models) are provided in Appendix A.

The dRSE solves the Maxwell wave equation

[TABLE]

for the perturbed system (the system of interest) described by a modified permittivity tensor $\hat{{\boldsymbol{\varepsilon}}}(\mathbf{r},\omega)+\Delta\hat{{\boldsymbol{\varepsilon}}}(\mathbf{r},\omega)$ , with the perturbation

[TABLE]

having the same pole structure as $\hat{{\boldsymbol{\varepsilon}}}(\mathbf{r},\omega)$ , i.e. with the poles at the same frequencies $\Omega_{j}$ and frequency-independent tensors $\Delta\hat{{\boldsymbol{\varepsilon}}}_{\infty}(\mathbf{r})$ and $\Delta\hat{\boldsymbol{\sigma}}_{j}(\mathbf{r})$ . The electric field $\mathbf{E}(\mathbf{r})$ and the frequency $\omega$ in Eq. (3) are, respectively, the wave function and the eigenfrequency of a RS of the perturbed system. The perturbed Maxwell’s wave equation (3) can be solved with the help of the Green’s function (GF) $\hat{\mathbf{G}}_{\omega}(\mathbf{r},\mathbf{r}^{\prime})$ of the unperturbed Eq. (1) in the following way:

[TABLE]

The completeness of the basis of the unperturbed RSs allows us to convert this integral equation to a matrix eigenvalue problem. This can be done by expanding the perturbed RS field into the unperturbed ones,

[TABLE]

and by using the Mittag-Leffler theorem again, this time expanding the Green’s dyadic,

[TABLE]

Here $\otimes$ denotes the dyadic product of the vector fields, and the RS wave functions $\mathbf{E}_{n}(\mathbf{r})$ are properly normalized Doost et al. (2014); Muljarov and Langbein (2016a, b); Muljarov and Langbein (2017). For non-static RSs, i.e. those having $\omega_{n}\neq 0$ , the correct normalization linked to Eq. (7) is given by Muljarov and Langbein (2016a, b); Muljarov and Weiss (2018)

[TABLE]

where $\mathbf{F}_{n}=(\mathbf{r}\cdot\nabla)\mathbf{E}_{n}$ , $V$ is an arbitrary volume containing the system, surface $S_{V}$ is its boundary, and $\partial/\partial s$ is the spatial derivative along the outer surface normal. Note that the electric field normalized according to Eq. (8) is a factor $\sqrt{2}$ smaller than the one used in the previous works Muljarov et al. (2010); Muljarov and Langbein (2016a). This has been introduced in order to symmetrize the normalization with respect to the electric and magnetic field (for more details see Ref. [Muljarov and Weiss, 2018]).

The sum in Eq. (7) is taken over all RSs; their eigenfrequencies $\omega_{n}$ play the role of simple poles of the GF. Owing to the frequency dispersion Eq. (2), the GF has alternative spectral representations linked to Eq. (7) via sum rules which the RSs of the dispersive system obey Muljarov and Langbein (2016a); Muljarov and Weiss (2018). Using in Eq. (5) a proper combination of the different spectral representations of the GF and equating the coefficients at the basis functions $\mathbf{E}_{n}(\mathbf{r})$ , this integral equation reduces to the following matrix eigenvalue problem, linear in frequency $\omega$ :

[TABLE]

which presents the key equation of the dRSE. Here, $V_{nm}$ and and $U_{nm}$ are the matrix elements of the non-dispersive and dispersive parts of the perturbation, given by

[TABLE]

and

[TABLE]

respectively. For a detailed derivation of the dRSE equation (9), see Refs. Muljarov and Langbein (2016a); Muljarov and Weiss (2018).

III Infinitesimal-dispersive basis

The dRSE formally requires the permittivity of the basis system and the perturbed system to have the same poles. While this is looking like a limitation of the dRSE, in reality, the perturbation $\Delta\hat{{\boldsymbol{\varepsilon}}}(\mathbf{r},\omega)$ can introduce new poles in the permittivity. An extreme case of this situation (illustrated in Sec. IV) is having an unperturbed system without frequency dispersion, such as a dielectric sphere in vacuum, but its RSs are used as a basis for the dRSE treating a system with dispersion, such as a metal or a semiconductor nanoparticle. However, in the functional space of dispersive solutions, the basis of non-dispersive RSs is incomplete and needs to be complemented (for completeness) with so-called pole RSs (pRSs).

The pRSs originate from a group of RSs of a dispersive system with a permittivity pole at $\omega=\Omega_{j}$ having an infinitesimal conductivity, $\hat{\boldsymbol{\sigma}}_{j}\to 0$ . All the pRSs due to this pole then become degenerate, $\omega_{n}\to\Omega_{j}$ , while having different values of the effective refractive index, and hence different wave functions. These wave functions, however, have vanishing normalization constants, $A_{n}\to 0$ , so that without perturbing the pole to a finite conductivity, the pRSs do not contribute to the modes of the system, as both $V_{nm}$ and $U_{nm}$ are vanishing for $n$ or $m$ referring to pRS.

To treat the pRSs of the $j$ -th pole of the permittivity, we introduce $i\hat{\boldsymbol{\sigma}}_{j}(\mathbf{r})=\xi\hat{\boldsymbol{s}}_{j}(\mathbf{r})$ with $\xi\to 0$ . $\hat{\boldsymbol{s}}_{j}(\mathbf{r})$ is an arbitrary finite tensor, which can be chosen, without loss of generality, as $\hat{\boldsymbol{s}}_{j}(\mathbf{r})=\hat{\mathbf{1}}$ within the volume $V$ of the system and zero outside it. In the limit $\xi\to 0$ , the eigenfrequencies of the pRSs become

[TABLE]

where $q_{n}$ are finite eigenvalues determined by the geometry of the system and its material parameters. For example, for transverse-magnetic (TM) modes of a sphere in vacuum they are found from solving a secular equation (37), see Appendix B. Note that for brevity of notations, here and below we have omitted index $j$ which would appear in all the eigenvalues and matrix elements due to the pRSs of pole $j$ .

According to Eqs. (2) and (13), the permittivity calculated at the pRS frequency $\omega_{n}$ , contributing both to Eqs.(8) and (11), is given by

[TABLE]

which is a finite quantity, determining the effective refractive index and hence the wavelength of the pRS. However, the frequency derivative of the permittivity, contributing to the RS normalization Eq. (8), is divergent,

[TABLE]

which leads to a vanishing normalization coefficient of the pRS field, as follows from Eq. (8):

[TABLE]

since $\xi\to 0$ but ${\Omega_{j}}/{q_{n}^{2}}$ is finite.

For any finite perturbation of the conductivity $\Delta\hat{\boldsymbol{\sigma}}_{j}\neq 0$ , the vanishing normalization constants of pole RSs are compensated in the perturbation matrix by divergencies due to the dispersion, $i\Delta\hat{\boldsymbol{\sigma}}_{j}/(\omega_{n}-\Omega_{j})\propto 1/\xi$ , thus making $U_{nm}$ finite. To take this limit properly, it is convenient to introduce the factors

[TABLE]

and to re-define the RS wave functions as

[TABLE]

where all $\tilde{\mathbf{E}}_{n}(\mathbf{r})$ have finite normalization. Expressing the perturbation matrices and eigenvectors as

[TABLE]

the dRSE matrix equation (9) takes the form:

[TABLE]

in which the elements $Q_{nm}$ and $S_{nm}$ are finite. Explicitly, the matrix $Q_{nm}$ is defined according to

[TABLE]

while the matrix $S_{nm}$ is given by

[TABLE]

if $n$ refers to a non-pRS, and by

[TABLE]

if $n$ refers to a pRS with $\omega_{n}=\Omega_{j}$ . Accordingly, the perturbed field is given by

[TABLE]

in terms of the modified basis RS wave functions $\tilde{\mathbf{E}}_{n}(\mathbf{r})$ and the eigenvector components $b_{n}$ determined by Eq. (20). Clearly, if $\hat{\boldsymbol{\sigma}}_{j}(\mathbf{r})=\Delta\hat{\boldsymbol{\sigma}}_{j}(\mathbf{r})=0$ , none of the corresponding pRSs is expected to contribute to the perturbed field. Indeed, in this case $S_{nm}=0$ for any $j$ -th pRS $n$ , according to Eq. (23), while $\alpha_{n}=0$ is leading to $\alpha_{n}^{2}Q_{nm}=0$ in Eq. (20), even though $V_{nm}\neq 0$ . Hence $b_{n}=0$ , showing that this pRS does not contribute to the perturbed RSs.

The method presented in this section is called infinitesimal-dispersive RSE (idRSE). The described procedure is suited for any number of new poles introduced in the perturbation. It has been already applied to systems with dispersion described by the Drude model Muljarov and Langbein (2016a). In Sec. IV below we illustrate the idRSE on several examples, for the dispersion of the permittivity including both Drude and Lorentz poles, or Lorentz poles only.

IV Results

In this section we present a selection of applications of the dRSE. To provide a clearer illustration of the physical effects treated, and to compare results of the dRSE with available exact solutions, we concentrate here on spherically symmetric systems, but note that the RSE can also be applied to non-spherical systems Doost et al. (2014). We consider TM polarization, and treat $l=1$ modes of metallic particles, showing localized surface plasmons (SP), where $l$ is the orbital quantum number. For semiconductor microspheres instead, we use $l\gg 1$ in order to obtain sharp whispering-gallery (WG) modes. The wave functions, normalization coefficients, and the secular equation determining the RS eigenvalues of a uniform sphere in vacuum are provided in Appendix B, along with the matrix elements for uniform perturbation. Approximate recursive relations determining the start values for solving the secular equation are given in Appendix C. The basis for the dRSE is truncated in such a way that all unperturbed RSs satisfying the condition $|n_{r}(\omega_{n})\omega_{n}|<\Omega_{\rm c}$ are included in the basis, where $n_{r}(\omega)=\sqrt{\varepsilon(\omega)}=n_{r}^{\prime}+in_{r}^{\prime\prime}$ is the complex refractive index, and $\Omega_{\rm c}$ is the cut-off frequency determining the basis size $N$ . For each $l$ , a longitudinal static mode Doost et al. (2014) with $\omega_{n}=0$ has been included in the basis; this is required for completeness of a non-dispersive basis or a basis without Ohm’s law pole in the dispersion.

To apply the dRSE to spherical particles made of realistic materials, such as gold and GaAs considered below, we have generated the corresponding parameters of the Drude-Lorentz model using the earlier developed fit program Sehmi et al. (2017a). This fit program is using the experimentally measured indices of refraction and absorption, $n_{r}^{\prime}$ and $n_{r}^{\prime\prime}$ , respectively, in order to find the optimal values of the pole frequencies $\Omega_{j}$ and their associated conductivities $\sigma_{j}$ . The value of the high-frequency permittivity $\varepsilon_{\infty}$ can be used as an adjustable parameter which can be chosen to generate a set of the basis RSs providing quick convergence with $N$ . We have found in particular that the value of $\varepsilon_{\infty}=1$ used for the basis system corresponds to the slowest convergence as in this case the refractive-index contrast between the system and the environment (vacuum) is minimized for high-frequency modes. In the example provided in Sec. IV.1 we transform a dielectric (silica) sphere with $\varepsilon_{\infty}=2.1272$ to a gold sphere with $\varepsilon_{\infty}+\Delta\varepsilon_{\infty}=0.5$ , while in Secs. IV.2–IV.4 we use $\varepsilon_{\infty}=11.0$ or 8.6013 for calculating the RSs of a GaAs microsphere, while keeping $\Delta\varepsilon_{\infty}=0$ . Table 1 of Appendix A shows all the fit parameters of gold and GaAs used in the dRSE calculations presented in this section below.

IV.1 Infinitesimal-dispersive RSE: Plasmons

To construct a basis for the idRSE treating a gold nanosphere, we add to the set of the RSs of a non-dispersive dielectric nanosphere with $n_{r}=1.4585$ the pRSs of a Drude pole $\Omega_{0}=-i\gamma$ and of 3 pairs of Lorentz poles $\Omega_{j}$ ( $j=\pm 1$ , $\pm 2$ and $\pm 3$ ) with infinitesimal conductivities. The poles are shown in Fig. 1(b), lined up with the corresponding fit of $n_{r}^{\prime}$ and $n_{r}^{\prime\prime}$ of gold Johnson and Christy (1972), which is illustrated in Fig. 1(a).

The spectrum of the perturbed system – a gold sphere of radius $R=10$ nm – shown in Fig. 1(b) reveals RSs to the high energy side of each pole in the dispersion. These are known as surface plasmon (SP) polariton resonances. In the electrostatic limit, they occur at $\varepsilon(\omega)=-2$ for $l=1$ , as given by Eq. (38) of Appendix B, and are well known in the literature Bohren and Huffman (1998). This equation provides one SP state for each pole in $\varepsilon$ . Accordingly, using the Drude model, only one such RS is present Muljarov and Langbein (2016a). Notably, the RSE calculates these poles beyond the electrostatic limit – the resulting evolution of the SPs versus radius is discussed below. The validity of the response predicted by the SPs at complex frequencies depends on the reliability of the DL model of $\varepsilon(\omega)$ at these frequencies, keeping in mind that the parameters of the model are fitted to describe the permittivity at real frequencies. At small imaginary parts of the SP frequencies, such as for SP0, the validity is expected to be better than at larger imaginary parts, such as for SP1 and SP2. Importantly, for correct calculation of the response of the system at real frequencies, such as the scattering and absorption spectra, the complete set of its RSs is required Lobanov et al. (2018), which necessarily includes all the SP modes generated by the model of $\varepsilon(\omega)$ used. The absorption lineshape is therefore not just a single resonance, but is more complex, as also observed in the experiment.

In addition to the SP mode around each pole of the permittivity, the full solution of the secular equation (see Eq. (37) of Appendix B) contains a countable infinite number of modes in the vicinity of the Drude and each Lorentz pole, originating from the pRSs. The inset in Fig. 1(b) shows these RSs approaching the first Lorentz pole $\Omega_{1}$ . We can also see from the figure that the idRSE matches the exact perturbed solution; however, in the presence of Lorentz poles, it is difficult to find all modes within a frequency range reliably from the secular equation (37). The idRSE can find these solutions due to the completeness of RSs and is expected to work equally well for systems without known analytic solutions Doost et al. (2014). The relative error of the idRSE, given in Fig. 1(c) for the SPs, shows a $N^{-3}$ dependence on the number of states in the basis, as observed previously for spherically symmetric systems Muljarov et al. (2010); Muljarov and Langbein (2016a).

The evolution of the SP modes with nanosphere radius $R$ changing from 1 nm to 100 nm is given in Fig. 2. We find that for $R<15$ nm, the SP modes keep close to the electrostatic limit, as expected for sizes much smaller than the light wavelength, $R\omega/c\ll 1$ , noting that the SP wavelengths are in the 250 to 500 nm range. Increasing $R$ , different SP modes change differently. Generally, they shift to longer wavelength, as expected from modes determined by size rather than dispersion. SP0, related to the Drude pole and thus representing the classical SP of gold, is not changing much in the full range of the radii shown, which is due to the dominant role of the Drude part in the dispersion of gold and a rather small value of $\gamma$ for the Drude pole, see Table 1. For small $R$ , SP0 has an energy around 2.4 eV and a linewidth of about 0.3 eV, comparable to the features seen in the measured absorption Bohren and Huffman (1998). Increasing $R$ , SP0 moves along a small arc, first shifting to the red, then increasing the magnitude of the imaginary part of the wave number due to the increasing role of the radiative broadening, and then shifting back to the blue, possibly due to a repulsion from SP1, the next closest surface plasmon mode originating from the Lorentz pole at $\Omega_{1}$ . The next SP mode, SP1, has an energy of 2.9 eV in the static regime and a linewidth of 1.5 eV; this mostly accounts for the experimentally observed intraband absorption features. Note however that in calculating the absorption and scattering properties, the SPs will interfere Lobanov et al. (2018), so that the lineshape is expected to be non-Lorentzian, again in agreement with experiment. Calculating these properties using the dRSE will be the topic of a future work. With increasing $R$ , SP1 is approximately keeping its linewidth, but is red-shifting, passing SP0 at around 80 nm. This leads to a redshift and broadening of the observed SP resonance in absorption, in qualitative agreement with experiment Myroshnychenko et al. (2008). SP2 shows a similar behavior as SP1, but with a twice broader linewidth, therefore describing a broad absorption in the interband absorption range of gold.

IV.2 Infinitesimal-dispersive RSE: Optical phonons

The optical phonons in GaAs provide a sharp resonance in the permittivity over the energy range of $28-40$ meV. This is well described by a single Lorentz pole with a purely imaginary positive conductivity, corresponding to a classical Lorentz oscillator having a purely real negative weight of the permittivity pole. The refractive and absorption index spectra and their fit are shown in Fig. 3(a). The fit parameters are given in Table 1 of Appendix A. We note that from about 33 to 36 meV, the absorption index supersedes the refractive index, which is the so-called reststrahlen region where a negative real part of the permittivity is causing an increased reflectivity of the material.

We use the idRSE with a basis system given by a dielectric microsphere with the same $\varepsilon_{\infty}$ as used in the above fit. We choose $l=15$ and $R=50\,\mu$ m in order for the Lorentz pole $\Omega_{1}$ to be close to WG mode frequencies. The RS frequencies of the basis system and of the perturbed system calculated via idRSE and exactly via Eq. (37), are shown in Fig. 3(c). Due to the large $l$ , the basis is clearly separated into WG, Fabry-Perot (FP), and leaky RSs. In the new system, we find the WG RSs strongly modified. There is a series of RSs approaching $\Omega_{1}$ from the lower frequency side, see also a zoom in Fig. 3(b). This series is countable infinite, limited in the numerical calculation by the finite basis size used. It resembles the series of basis RSs with increasing frequency and accordingly reducing wavelength. In fact, approaching from the lower frequency side the pole at $\Omega_{1}$ , which is also the transversal optical phonon frequency, the real part of the permittivity increases without bound, and so does the refractive index $n_{r}^{\prime}$ , in this way reducing the light wavelength within the sphere towards zero. The imaginary part of the perturbed RSs is, however, different from that of the WG modes, as it is dominated by the material absorption of the perturbed system. On the right side of $\Omega_{1}$ instead, the refractive index is reduced when approaching the longitudinal optical phonon frequency, separated from $\Omega_{1}$ by the reststrahlen band. Consequently, the light wavelength within the sphere is increased. This produces a set of RSs which can be labeled with a radial quantum number reducing towards the frequency at which the refractive index is zero. The leaky RSs are, in turn, nearly unaffected by the perturbation within the sphere, which is understandable considering their little spatial overlap with the sphere.

We again note that it is difficult to find all the exact solutions corresponding to the series of RSs near the pole of the permittivity, which is seen as the absence of some of the exact RS frequencies in Fig. 3(b). However, we found that using the RS frequencies determined by the idRSE as starting points for solving numerically the secular equation (via e.g. a gradient decent method) enabled us to overcome this difficulty. We also note that it is important for the accuracy and convergence of the dRSE to the exact result with increasing basis size, that no RSs are missed in the basis, which is selected using the high spatial frequency cutoff condition $|n_{r}(\omega_{n})\omega_{n}|<\Omega_{\rm c}$ . The relative error for $N=319$ , shown in Fig. 3(c), is in the $10^{-8}$ range. This error is significantly smaller than the ones observed for perturbations of $\varepsilon_{\infty}$ for a similar basis size Doost et al. (2014); Muljarov et al. (2010), a finding which we attribute the perturbation being localized in frequency. Once the frequency cutoff of the basis is much larger than the frequency of the perturbed pole (here $\Omega_{\rm c}\gg\Omega_{1}$ ), the effect of the perturbation on the RSs not included in the basis is small, as suggested by the denominator $\omega_{n}-\Omega_{j}$ in the corresponding matrix elements given in Eq. (12).

IV.3 Pole-shift: diagonal dRSE

To treat a shift in the pole positions of the permittivity, the idRSE requires the shifted pRSs to be included in the basis, along with the RSs of the unshifted poles. This can add a significant numerical effort, as the matrix diagonalization required to solve Eq. (9) has the computational complexity ${\cal O}(N^{3})$ . To avoid this, one can use a diagonal version of the dRSE, only retaining the diagonal elements in Eq. (9). We call this method ddRSE. It is linked to the first-order perturbation theory Doost et al. (2014); Weiss et al. (2016), which is also neglecting all off-diagonal elements. Accordingly, we expect this approach to be suited to treat small shifts in resonance frequency of the permittivity, such as the change of optical phonon frequencies with pressure or temperature. Neglecting the off-diagonal elements of Eq. (9) yields

[TABLE]

Considering only the resonance shift of one pole $\Omega_{1}$ , without changing conductivities or $\varepsilon_{\infty}$ , all $V_{nm}$ vanish, and

[TABLE]

where $\widetilde{\Omega}_{1}$ is the perturbed pole position.

To evaluate the suitability of this approach, we use the optical phonon-induced dispersion in GaAs discussed in Sec. IV.2. We consider a shift of the phonon frequency by $-1$ %, which is applied to the real part of $\Omega_{1}$ , in order to obtain $\widetilde{\Omega}_{1}$ . We note that the GaAs phonon frequency under isotropic pressure shifts by 0.49 meV/GPa Adachi (1994), so that a 1% change would corresponds to exerting about 0.7 GPa pressure.

The resulting RSs close to $\Omega_{1}$ are shown in Fig. 4(a). The exact perturbed RSs converge towards the shifted pole $\widetilde{\Omega}_{1}$ , as expected. The diagonal (first order) approach is reproducing the exact RSs well when they are sufficiently far from the pole. Close to the pole instead, large errors are encountered. Notably, using all off-diagonal elements defined by Eq. (26) and solving the eigenvalue problem Eq. (9) does not significantly reduce the errors. This indicates that the effect of the basis incompleteness, missing the pRSs of $\widetilde{\Omega}_{1}$ , is dominating the remaining error.

The relative error of the ddRSE results compared to the exact solutions is shown in Fig. 4(b). We find that the error for the WG modes is decreasing with increasing separation from $\Omega_{1}$ . The leaky modes have a small error since they have little spatial overlap with the sphere, and are spectrally broad, making them insensitive to a shift of $\Omega_{1}$ . We can thus conclude that the diagonal approximation, equivalent to a first-order perturbation theory, is a valid approach for small pole shifts, and separations of the RSs of interest from the poles much larger than these shifts.

The physical properties in scattering and absorption are expected to be dominated by the RSs away from the poles of the permittivity, as close to the poles the permittivity is very high, screening the outside field. Furthermore, for the response at real frequencies, the specific position of the RSs close to the poles is having a minor influence, since the poles typically have a significant imaginary part, and are therefore separated from the real axis.

IV.4 Dispersive RSE: Absorption to gain

Another interesting feature one can address using the dRSE is the transition from absorption to gain, leading to lasing. In terms of the RSs, lasing occurs once the imaginary part of the RS frequency becomes positive, indicating an exponential increase of the field with time. This increase of the field comes with a corresponding increase of stimulated transitions, which in realistic systems results in gain saturation, reducing the imaginary part of the RS frequency to zero in a steady-state state regime of lasing Türeci et al. (2008); Chow et al. (2012).

In the classical model of Lorentz oscillators, poles representing absorption have positive imaginary conductivity Sehmi et al. (2017a). Continuously changing this conductivity to negative imaginary values transforms the absorption to gain. When representing the absorption due to interband transitions in a solid, a prefactor $1-f_{\mathrm{e}}-f_{\mathrm{h}}$ is applied to the conductivity, where $f_{\mathrm{e}}$ and $f_{\mathrm{h}}$ are the occupation functions of, respectively, the related electron state in the conduction band and hole state in the valence band, which in thermal equilibrium are given by the Fermi distribution function Chow et al. (2012). Creating the population inversion changes the sign of the conductivity by changing the occupation factor $1-f_{\mathrm{e}}-f_{\mathrm{h}}$ from 1 for no occupation ( $f_{\mathrm{e}}=f_{\mathrm{h}}=0$ ) to $-1$ for full occupation ( $f_{\mathrm{e}}=f_{\mathrm{h}}=1$ ).

As an example, we use GaAs, which is a direct band-gap semiconductor widely used in semiconductor lasers. To model the transition from absorption to gain, we concentrate on the energy region close to the band gap, where a population inversion is created by optical or electrical pumping, and fit the permittivity data Sop in the range 1.3 to 1.65 eV with four pairs of Lorentz poles having purely imaginary conductivities. The measured absorption index $n_{r}^{\prime\prime}$ and the resulting fit of the Lorentz model are shown in Fig. 5(a). We see that the fit is not as good as for the optical-phonon range shown in Fig. 3(a), and allowing the conductivities to have non-zero real parts could lead to a better agreement, as in the case of gold. However, we intentionally fit the permittivity with a model of classical Lorentz oscillators, in order to enable a simulation of the population inversion. The fit function contains three poles close to the band-gap having a small linewidth to describe the sharp increase of the absorption at the band-edge, and the fourth pole at higher energy outside the fit range, providing a background. When changing the sign of the first-pole conductivity $\sigma_{1}$ , in order to simulate the change from no to full occupation, the resulting absorption index shows a negative region around 1.49 eV, indicating the region where the material exhibits gain.

We first consider a rather large microsphere of radius $R=10\,\mu$ m, and choose $l=10$ . The wavelength in the medium resonant to the gain is about $\lambda=0.24\,\mu$ m, much smaller than the radius, and thus the RSs close to the gain pole are of FP type for the chosen $l$ . To provide a round-trip amplification, the losses due to the transmission across the sphere surface must be compensated by the gain through the sphere. The amplitude reflection coefficient for normal incidence is $r=(n_{r}^{\prime}-1)/(n_{r}^{\prime}+1)\approx 0.55$ , using the refractive index of $n_{r}^{\prime}=3.5$ . The amplification of the wave amplitude propagating one diameter is $g=\exp(-2Rn_{r}^{\prime\prime}\omega/c)$ . The condition for round-trip amplification is then $gr>1$ , which requires $n_{r}^{\prime\prime}<-0.004$ . We see in Fig. 5(a) that this condition is fulfilled around 1.49 eV.

The RSs under this gain, calculated using the dRSE with the dispersive basis of RSs in the absorption regime (using a basis size $N=676$ corresponding to $\Omega_{c}=26.78$ eV), are shown in Fig. 5(b) over the same energy range as used in Fig. 5(a). Away from the poles, the basis shows nearly equidistant FP modes, which have an imaginary part following the losses due to the absorption index $n_{r}^{\prime\prime}$ . In the absorption regime, close to the poles, the RSs converge from the lower frequency side towards the poles, which are at the transversal optical resonance frequencies $\omega_{\rm T}$ , and show gaps at higher frequencies up to the longitudinal optical frequency $\omega_{\rm L}$ , where $n_{r}=0$ , similar to the case of the optical phonons in Fig. 3(b). Note that in the present case, the set of RSs approaching the longitudinal optical frequencies is rather large, in the order of 100. This is due to the large radial quantum number of the basis FP modes in the region of the pole frequency. This number can be estimated as $4R/\lambda\approx 166$ , using $\lambda=0.24\,\mu$ m.

After switching to gain, some RSs in the gain spectral region show a positive imaginary part, as expected from the above discussion, and are thus lasing RSs. Furthermore, a curious arrangement of RSs is observed close to the gain pole, which occurs due to the inverted influence of the pole – the refractive index now increases approaching the pole from higher frequencies, so that the infinite series of RSs now converges towards the pole from the higher frequency side. These RSs, however, form a loop, moving clockwise over the pole, starting at a value around the corresponding longitudinal optical frequency $\omega_{\rm L}$ , which is now on the lower frequency side of the pole, as shown in the inset of Fig. 5(b). The original, unperturbed series of FP modes, present without the pole, is instead repelled from the pole. This is qualitatively different from the case of absorption, where the RS series is attracted to the pole, splitting the RSs into an infinite set of all radial quantum numbers below the pole frequency $\omega_{\rm T}=\Omega_{1}$ , and a second infinite set of all radial quantum numbers above $\omega_{\rm L}$ .

To investigate the regime of WG mode lasing at $l=10$ , we consider in Fig. 6 a much smaller sphere of radius $R=0.94\,\mu$ m. In this case, the gain region lines up with a WG mode in the basis (WG3), moving it to a positive imaginary energy. The evolution of the imaginary part of WG3 with the change from absorption to gain is shown in the inset. We observe that with the build-up of the gain, given by the decreasing imaginary part of the conductivity, the imaginary part of the WG3 energy moves continuously towards positive values, eventually reaching the lasing regime. The other RSs only change slightly since they are not close to $\Omega_{1}$ , so that WG3 remains the only lasing mode of the system.

V Conclusions

We have shown a set of applications of the dispersive RSE, a method developed for accurate and efficient calculation of the resonant states of optical systems exhibiting frequency dispersion. We have also presented the infinitesimal-dispersive RSE, which uses the resonant states of a non-dispersive system as a basis for treating an optical system with dispersion, adding new poles in the frequency dispersion in the limit of zero conductivity. This requires an extension of the RSE basis, in order to include the so-called pole resonant states which are all at the pole frequency.

We apply the method to calculate the resonant states of spheres made of materials described by a frequency-dependent permittivity in the form of a generalized Drude-Lorentz model. Using our fit program Sehmi et al. (2017a) accurately determining the parameters of the model, we fit gold with one Drude and three pairs of Lorentz poles. Using the infinitesimal-dispersive RSE, we perturb a 10 nm-radius nanosphere of sand into gold, with the error in the resonant-state frequencies scaling as the inverse cube of the basis size. We also analyse the evolution of the surface plasmon polaritons with the size of the gold nanosphere, showing the effects of radiation and retardation with increasing size.

We accurately describe the optical phonon range in GaAs with a single pair of Lorentz poles and perturb a 50 $\mu$ m-radius dielectric sphere to GaAs showing the reststrahlen band near the pole and a series of resonant states approaching the pole. Fitting the band gap range 1.3–1.65 eV of GaAs with four pairs of classical Lorentz model poles, we perturb the lowest energy pole from absorption to gain and identify lasing resonant states. In general, each Lorentz pole of the permittivity creates an additional infinite series of resonant states with the full set of radial quantum numbers corresponding to frequencies starting from a value around the longitudinal frequency at which the refractive index vanishes and extending towards infinity or converging to the next Lorentz pole at which the index diverges. Adding an absorption pole attracts the non-dispersive series, splitting it into two infinite series, above and below the pole. Adding a gain pole instead repels the non-dispersive series, and creates an additional series forming a clockwise loop in complex frequency, starting from the longitudinal frequency left of the pole and converging towards the pole from the right.

Finally, we introduce a diagonal dispersive RSE, a first-order perturbation method for treating small changes in the pole position, which is suited for the resonant states sufficiently separated from the pole of the dispersion and which has a much lower numerical complexity than the full dispersive RSE.

We believe that the dispersive RSE, once implemented in a accessible software package, including arbitrary 3D open systems Doost et al. (2014), and scattering and absorption with external excitation Lobanov et al. (2018), will be a powerful and widely used tool for treating electromagnetic problems.

Acknowledgements.

This work was supported by the Sêr Cymru National Research Network in Advanced Engineering and Materials and by the EPSRC under grant EP/M020479/1. E.A.M. acknowledges support by RBRF (Grant No. 16-29-03282).

Appendix A Drude-Lorentz model of the permittivity and fit for Au and GaAs

Implementing causality, the dispersive permittivity obeys the relation $\hat{{\boldsymbol{\varepsilon}}}^{\ast}(\mathbf{r},\omega)=\hat{{\boldsymbol{\varepsilon}}}(\mathbf{r},-\omega^{\ast})$ . This implies that poles appear in pairs in the permittivity Eq. (2). In particular, each pole at $\Omega_{j}$ has a partner at $\Omega_{-j}=-\Omega^{\ast}_{j}$ with $\hat{\boldsymbol{\sigma}}_{-j}=\hat{\boldsymbol{\sigma}}^{\ast}_{j}$ . A pole with Re $\Omega_{j}=0$ , such as that of the Drude model, can be treated as a pair of poles merged to produce a single pole on the imaginary $\omega$ -axis with a purely real conductivity. Furthermore, Lorentz reciprocity requires that all tensors in Eq. (2) are symmetric. For isotropic materials, such as gold and GaAs, these tensors can be replaced by scalars.

In the calculations presented in this paper, we use for GaAs a Lorentz model of the permittivity consisting of $L$ pairs of poles,

[TABLE]

and for gold a Drude-Lorentz model,

[TABLE]

in which the Drude term of the dispersion has been added, consisting of Ohm’s law pole at $\omega=0$ and a Drude pole at $\omega=\Omega_{0}=-i\gamma$ :

[TABLE]

Here $\varepsilon_{\infty}$ is the permittivity at high frequencies and $\sigma$ is the real DC conductivity. The generalized conductivities $\sigma_{j}=\sigma^{\prime}_{j}+i\sigma^{\prime\prime}_{j}$ at the Lorentz poles $\Omega_{j}=\Omega_{j}^{\prime}+i\Omega_{j}^{\prime\prime}$ are complex.

Table 1 contains the parameters of the fits shown in Figs. 1, 3, and 5, using the Drude-Lorentz and the Lorentz models with $L=3$ , 1, and 4, respectively. The resulting root-mean square error $S$ of the fit is also shown. Its definition, with and without experimental errors taken into account in the optimization, is provided in Ref. [Sehmi et al., 2017a]. This difference in the definition leads, in particular, to a much higher value of $S$ for gold, in which case the fit takes the measurement errors Johnson and Christy (1972) into account.

Appendix B TM modes of a sphere, their normalization and matrix elements

The electric field of a TM mode of a sphere of radius $R$ in vacuum has the form (in spherical coordinates)

[TABLE]

where $Y_{lm}(\theta,\varphi)$ are real-valued spherical harmonics Doost et al. (2014), and the radial wave functions within the sphere, $r\leqslant R$ , are given by

[TABLE]

Here $k_{n}=\omega_{n}/c$ is the RS wave number and $n_{r}$ is the refractive index, which is given by

[TABLE]

with the permittivity $\varepsilon(\omega)$ . The normalization constant $A_{nl}$ is found from the general RS normalization condition with dispersion, Eq. (8), and has the form Muljarov and Langbein (2016b)

[TABLE]

with

[TABLE]

and

[TABLE]

where $j_{l}(x)$ are the spherical Bessel functions, $x=n_{r}\omega_{n}R/c$ ,

[TABLE]

and $n_{r}$ is given by Eq. (32).

The RS eigenfrequencies $\omega_{n}$ are solutions of the secular equation for the TM polarization

[TABLE]

in which $h_{l}(z)\equiv h^{(1)}_{l}(z)$ is the spherical Hankel function of the first kind, with $z=\omega_{n}R/c$ , and $n_{r}(\omega_{n})$ is defined by Eq. (32).

For small values of $R$ , corresponding to the electrostatic limit, the spherical Bessel and Hankel functions can be expanded into power series of their arguments and keeping only the leading terms, which is valid for $|\omega|R/c\ll 1$ and $|n_{r}(\omega)\omega|R/c\ll 1$ , one can obtain from Eq. (37) a well-known equation

[TABLE]

determining for each $l$ the eigenfrequency $\omega$ of the surface plasmon in the low-frequency (electrostatic) limit.

For the pRSs, due to a pole of the dispersion at $\Omega_{j}$ , all the eigenfrequencies are the same, $\omega_{n}=\Omega_{j}$ , but the eigenvalues of the refractive index $n_{r}(\omega_{n})$ are different. They are again determined by the secular equation (37) which is now solved for $n_{r}$ while keeping $z=\Omega_{j}R/c$ fixed. These eigenvalues $n_{r}$ are linked to $q_{n}$ , as defined by Eq. (14),

[TABLE]

and to the values $\eta_{n}$ in the normalization constants Eq. (36):

[TABLE]

in accordance with Eqs. (15) and (17), using $\xi=\alpha_{n}^{2}\Omega_{j}/q_{n}$ . Clearly, $\xi\to 0$ makes the normalization constants $A_{nl}$ vanishing. Re-defining these constants as

[TABLE]

in accordance with Eq. (18), we find the new (finite) normalization constants of the pRSs:

[TABLE]

For any homogeneous perturbation of the sphere, all the matrix elements, contributing to Eq. (9) or Eq. (20) are given by the integrals

[TABLE]

where $V_{R}$ is the volume of the sphere. These integrals have the following analytic form:

[TABLE]

where $x=n_{r}(\omega_{n})\omega_{n}R/c$ and $y=n_{r}(\omega_{m})\omega_{m}R/c$ , and both RSs $n$ and $m$ belong to the same orbital quantum number $l$ . The function $F_{l}(x,y)$ is defined as

[TABLE]

which for coinciding arguments becomes

[TABLE]

Note that using the normalization coefficients Eq. (42) for the pRSs, the diagonal elements Eq. (44) become $W_{nn}=-q_{n}$ , in accordance with their normalization condition

[TABLE]

see Eq. (16).

Appendix C Calculating the basis RSs

To generate a basis of unperturbed RSs, the secular equation (37) is solved numerically using the Newton-Raphson algorithm. This relies on a suited choice of guess values, which for $|n_{r}z|\lesssim 1$ are found on a two-dimensional grid in the complex $z$ plane. For large values of $n_{r}z$ , however, we use approximate recursive relations given below, which provide reliable guess values based on the eigenvalues found for the neighboring RSs. To derive these relations, it is convenient to modify Eq. (37) to the following form:

[TABLE]

where $j^{\prime}_{l}(x)$ and $h^{\prime}_{l}(x)$ are the derivatives of the spherical Bessel and Hankel functions with respect to their arguments. Noting that

[TABLE]

at large $x$ , we find

[TABLE]

which provides an estimate for the left-hand side of Eq. (48). Its right-hand side is in turn a slowly changing function. Therefore, the neighboring RSs are separated by $\Delta x\approx\pi$ , which is the period of the tangent function. This separation can be further evaluated as

[TABLE]

where $n_{r}^{\prime}=dn_{r}/dz$ . Then we find approximate relations for the neighboring roots of Eq. (37) in the following three cases:

(i) No dispersion of permittivity or RS eigenfrequencies are far from any poles of the dispersion, i.e. $n_{r}$ is nearly constant:

[TABLE]

(ii) Strong dispersion, such as that for the RS near the pole of the dispersion:

[TABLE]

(iii) Pole RSs – which are required for constructing an infinitesimal-dispersive basis. In this case Eq. (37) is solved for $n_{r}$ with $z$ fixed. For the neighboring pRSs, we find

[TABLE]

Bibliography23

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Muljarov et al. (2010) E. A. Muljarov, W. Langbein, and R. Zimmermann, Europhys. Lett. 92 , 50010 (2010).
2Muljarov and Langbein (2016 a) E. A. Muljarov and W. Langbein, Phys. Rev. B 93 , 075417 (2016 a).
3Sehmi et al. (2017 a) H. S. Sehmi, W. Langbein, and E. A. Muljarov, Phys. Rev. B 95 , 115444 (2017 a).
4Sehmi et al. (2017 b) H. S. Sehmi, W. Langbein, and E. A. Muljarov, in 2017 Progress In Electromagnetics Research Symposium - Spring (PIERS) (IEEE, United States, 2017 b), pp. 994–1000.
5Doost et al. (2014) M. B. Doost, W. Langbein, and E. A. Muljarov, Phys. Rev. A 90 , 013834 (2014).
6Lobanov et al. (2017) S. V. Lobanov, G. Zoriniants, W. Langbein, and E. A. Muljarov, Phys. Rev. A 95 , 053848 (2017).
7Shahbazyan and Stockman (2014) T. Shahbazyan and M. Stockman, Plasmonics: Theory and Applications (Springer, Netherlands, 2014), ISBN 9789400778054.
8Vahala (2004) K. Vahala, Optical Microcavities (World Scientific, Singapore, 2004).