Resonant-state expansion applied to three-dimensional open optical systems: A complete set of static modes

TL;DR

This paper develops and applies a complete set of static modes for the resonant-state expansion in 3D open optical systems, improving accuracy in modeling static electric fields and resonant states.

Contribution

It introduces two alternative complete static mode sets for the RSE in 3D systems and demonstrates their effectiveness in improving perturbation results.

Findings

Convergence of RSE towards exact results for size-reduction perturbations.

Complete static modes improve agreement with finite-element simulations.

Application to dielectric cylinder resonant states confirms method's accuracy.

Abstract

We present two alternative complete sets of static modes of a homogeneous dielectric sphere, for their use in the resonant-state expansion (RSE), a rigorous perturbative method in electrodynamics. Physically, these modes are needed to correctly describe the static electric field of a charge redistribution within the optical system due to a perturbation of the permittivity. We demonstrate the convergence of the RSE towards the exact result for a perturbation describing a size reduction of the basis sphere. We then revisit the quarter-sphere perturbation treated in [Doost {\it et al.}, Phys. Rev. A {\bf 90}, 013834 (2014)], where only a single static mode per each angular momentum was introduced, and show that using a complete set of static modes leads to a small, though non-negligible correction of the RSE result, improving the agreement with finite-element simulations. As another example of applying the RSE with a complete set of static modes, we calculate the resonant states of a dielectric cylinder, also comparing the result with a finite-element simulation.

Tables

Table 1. Table 1: Scaling parameters C 1 subscript 𝐶 1 C_{1} ( C 2 subscript 𝐶 2 C_{2} ) of the number of basis RSs (static modes) for ϵ = 4 italic-ϵ 4 \epsilon=4 and different dimensionalities.

Dim.	1D	2D	3D
$C_1^{\mathrm{TE}}$	0.3125	0.1141	0.04258
$C_1^{\mathrm{TM}}$	0.325	0.1144	0.04287
VC $C_2$	6.712	3.434	1.945
VSC $C_2$	6.196	3.276	1.857

Equations

$$\nabla \times \nabla \times {\text{E}}_n=k_n^2\hat{\mathbf{\epsilon }}(\text{r}){\text{E}}_n$$

$$\nabla \times {\mathbf{E}}_{\lambda }^S=0.$$

$$\nabla \cdot \mathbf{D}=0$$

$$\nabla \cdot \mathbf{E}=-\mathbf{E}\cdot \frac{\nabla \epsilon }{\epsilon }\ne 0$$

$$\nabla \cdot {\text{E}}_n=0$$

$${\text{E}}_{\lambda }^S(\text{r})=-\nabla {\psi }_{\lambda }(\text{r}).$$

$$\int {\text{E}}_{\lambda }^S(\text{r})\cdot {\text{D}}_{\lambda }^S(\text{r})d\text{r}=1,$$

$$I_{\lambda {\lambda }^{\prime }}={\int }_V{\text{E}}_{\lambda }^S\cdot {\text{D}}_{{\lambda }^{\prime }}^{S}d\text{r}=-{\oint }_{S_V}{\psi }_{\lambda }{\text{D}}_{{\lambda }^{\prime }}^S\cdot d\mathbf{S}+{\int }_V{\psi }_{\lambda }\nabla \cdot {\text{D}}_{{\lambda }^{\prime }}^{S}d\text{r},$$

$${{\psi }_{\lambda }|}_{S_V}=0,$$

$$\hat{L}\psi (\text{r})=-\nabla \cdot \hat{\mathbf{\epsilon }}(\text{r})\nabla \psi (\text{r})=\nabla \cdot \text{D}(\text{r})=\rho (\text{r}).$$

$$\hat{L}\psi_{\lambda}(\textbf{r})=\left\{\begin{array}[]{cl}\Lambda\psi_{\lambda}&{\rm within}\ V_{0}\\ 0&{\rm otherwise}\,,\end{array}\right.$$

$$I_{\lambda {\lambda }^{\prime }}$$

$${\int }_V{\text{E}}_{\lambda }^S\cdot {\text{D}}_{n}d\text{r}=-{\oint }_{S_V}{\psi }_{\lambda }{\text{D}}_n\cdot d\mathbf{S}+{\int }_V{\psi }_{\lambda }\nabla \cdot {\text{D}}_{n}d\text{r}=0,$$

$${\rho }_{\lambda }(\text{r})=\Lambda {\psi }_{\lambda }(\text{r}),$$

$${\sigma }_{\lambda }=-{{\text{D}}_{\lambda }^S|}_{S_V^{-}}$$

$$\epsilon (r)=(ϵ-1)\Theta (R-r)+1$$

$${\psi }_{\lambda }(\text{r})=f_{\lambda }(r)Y_{lm}(\Omega ),$$

$$({\nabla }_r^2+{\lambda }^2)f_{\lambda }(r)=0\mathrm{for}r<R$$

$${\nabla }_r^2{f}_{\lambda }(r)=0\mathrm{for}r>R,$$

$${\nabla }_r^2=\frac{d^2}{d{r}^2}+\frac{2}{r}\frac{d}{dr}-\frac{l(l+1)}{r^2}$$

$$f_{\lambda}(r)=A_{\lambda}\times\left\{\begin{array}[]{ll}j_{l}(\lambda r)&r<R\\ j_{l}(\lambda R)(R/r)^{l+1}&r>R\,,\end{array}\right.$$

$${\text{E}}_{\lambda }^S(\text{r})$$

$$I_{\lambda {\lambda }^{\prime }}$$

$$F_{\lambda }$$

$$G_{\lambda }$$

$$\alpha F_{\lambda }+\beta G_{\lambda }=0$$

$$\alpha {\lambda }^2{f}_{\lambda }(R)-\beta R{\left[\epsilon (r)\frac{d{f}_{\lambda }(r)}{dr}\right]}_{R_{-}}^{R_{+}}=0,$$

$$I_{\lambda \lambda }=\frac{R}{{\lambda }^2}{F}_{\lambda }{G}_{\lambda }+ϵ{\lambda }^2{\int }_0^R{f}_{\lambda }^2(r)r^{2}dr=1,$$

$$G_{\lambda }=0$$

$${\rho }_{\lambda }(\text{r})=ϵ{\lambda }^2{\psi }_{\lambda }(\text{r})\Theta (R-r),$$

$$A_{\lambda }^2=\frac{2}{ϵ{\lambda }^2{R}^3}{\left[j_l^2(\lambda R)-j_{l-1}(\lambda R)j_{l+1}(\lambda R)\right]}^{-1}.$$

Full text

I Introduction

The resonant-state expansion (RSE) is a powerful theoretical method, recently developed in electrodynamics MuljarovEPL10 for accurate calculation of resonant states (RSs) of an arbitrary open optical system. The concept of RSs presents a mathematically rigorous way of describing physical resonances of an open system, seen in its optical spectra, for example in the scattering cross-section. Using the Mittag-Leffler theorem, one can determine these spectra by expanding them into the RSs of the system LobanovPRA18 . The RSs thus contain the full information about the system and ideally present a complete set of functions suited for expansion of any vector field within the volume of the optical system. This completeness of the RSs as well as the Mittag-Leffler expansion of the Green’s diadic are at the heart of the RSE which allows one to accurately calculate the RSs of an optical system, using as basis the RSs of an unperturbed system which differs from the system of interest by a perturbation of e.g. its permittivity DoostPRA14 . To determine a complete set of RSs, it is advantageous to choose the basis system to be exactly solvable, such as a homogeneous dielectric sphere in vacuum, for which the analytical solutions in the form of Mie resonances are well known in the literature LaiPRA90 ; DoostPRA14 .

When applying the RSE to three-dimensional (3D) open optical systems DoostPRA14 , it has been recognized that in addition to the RSs $\textbf{E}_{n}(\textbf{r})$ describing physical resonances and satisfying Maxwell’s wave equations

[TABLE]

with outgoing boundary conditions, one needs to include in the basis for the RSE also the zero-frequency modes $\mathbf{E}^{S}_{\lambda}(\textbf{r})$, satisfying the static Maxwell’s equation

[TABLE]

Here, index $n$ is used to label RSs, $\lambda$ to label static modes, $k_{n}=\omega_{n}/c$ is the RS wave number in vacuum, $\omega_{n}$ is the RS eigen frequency, $\hat{\boldsymbol{\varepsilon}}(\textbf{r})$ is the permittivity tensor, and the permeability is assumed to be $\mu=1$ for simplicity. The electric fields $\mathbf{E}^{S}_{\lambda}$ of the static modes do not satisfy the full set of Maxwell’s equations, and therefore do not represent any physical states of the system. However, they are required for the Mittag-Leffler expansion of the Green’s function, playing the role similar to that of the cut poles included for completeness into the RSE of 2D systems, for which the Green’s function has cuts in the complex frequency plane DoostPRA13 . In the recent work DoostPRA14 , a single static longitudinal electric mode was added for each angular momentum. The added modes satisfy Maxwell’s equation

[TABLE]

both within and outside a dielectric sphere but violate Maxwell’s boundary condition of the continuity of the normal component of the displacement $\mathbf{D}$ across the sphere boundary. Adding only these modes to the basis for the RSE was suited to describe a homogeneous perturbation of the sphere permittivity, and only required for transversal magnetic modes which have a finite electric field normal to the interface. However, we found that adding only these static modes was insufficient to treat another simple perturbation of the dielectric sphere – reducing its size.

In general, the RSs of any optical system satisfy Eq. (3) at any point in space. This implies, in particular, that their electric fields have a non-zero divergence

[TABLE]

in the regions of inhomogeneity of the permittivity $\varepsilon(\textbf{r})$. Obviously, this non-zero divergence cannot be reproduced by the RSE if all basis RSs respect

[TABLE]

inside the system. This is the case, however, for the RSs of a homogeneous sphere, which are known in the literature as transverse-electric (TE) and transverse-magnetic (TM) modes responsible for the resonant behavior in Mie scattering Bohren1998 . Indeed, due to Maxwell’s equation $\nabla\times\textbf{H}_{n}=-ik_{n}\varepsilon\textbf{E}_{n}$, Eq. (5) holds everywhere except the sphere boundary. The static modes introduced in Ref. DoostPRA14, also have this property. One therefore needs additional modes in the basis with $\nabla\cdot\mathbf{E}^{S}_{\lambda}\neq 0$, which can only be static modes, having $k^{S}_{\lambda}=0$.

In this paper, we present in Sec. II a general formulation of the problem of static modes of an open optical system and then introduce in Sec. III two sets of longitudinal electric static modes of a dielectric sphere in vacuum, both suited for a correct treatment of an arbitrary perturbation of its permittivity. One such set of modes is complementing the static modes already introduced in Ref. DoostPRA14, . Physically, these modes carry both volume charges within the sphere and surface charges on its boundary, so we call them volume-surface charge (VSC) modes. The other set introduces only volume charges, as they satisfy Maxwell’s boundary condition for D on the sphere boundary, and are therefore called volume-charge (VC) static modes. Note that for treating perturbations of the permeability $\mu$, similar sets of longitudinal static modes can be introduced for the magnetic field. In Secs. IV.1 and IV.2, we test both sets on the exactly solvable cases of, respectively, strength and size perturbations of a homogeneous sphere, also studying the convergence of the RSE towards the exact solution. We then revisit in Sec. IV.3 the quarter-sphere perturbation and compare results with the previous calculation DoostPRA14 and with finite-element simulations. Finally, in Sec. IV.4 we use the RSE with static modes for calculating the RSs of a dielectric cylinder, also comparing results with a finite-element simulation.

II Static modes of an arbitrary open optical system

We first consider an arbitrary finite open optical system in vacuum or in a homogeneous medium. Static modes $\textbf{E}_{\lambda}^{S}$ of this system satisfy Eq. (2), and therefore can be written as

[TABLE]

Considering only square integrable solutions, the proper normalization of static modes, as it was derived in Ref. DoostPRA14, , takes the form

[TABLE]

where $\textbf{D}_{\lambda}^{S}(\textbf{r})=\hat{\boldsymbol{\varepsilon}}(\textbf{r})\textbf{E}_{\lambda}^{S}(\textbf{r})$, and the integration is performed over the entire 3D space. Note that the factor of 2 difference compared to Eq. (7) of Ref. DoostPRA14 was introduced in a generalized formulation of the RSE MuljarovOL18 , which we use in the present work. Here we assume a frequency-independent permittivity; a generalization of Eq. (7) for systems with frequency dispersion can be found in Ref. MuljarovOL18, . We note that a zero frequency pole in the permittivity of a conductive material results in a vanishing amplitude of the normalized static modes and thus in their vanishing contribution.

Let us now consider an expression similar to Eq. (7) for two static modes, $\lambda$ and $\lambda^{\prime}$, including the case of the same mode $\lambda=\lambda^{\prime}$,

[TABLE]

where $V$ is an arbitrary (finite or infinite) volume which includes all the inhomogeneities of the optical system, and $S_{V}$ is its boundary. The right-hand side of Eq. (8) is obtained by using Eq. (6) and the divergence theorem. If we require, as boundary condition for the static modes, that $\psi_{\lambda}$ is vanishing on the surface $S_{V}$,

[TABLE]

and outside it, for any state $\lambda$, then the surface term on the RHS of Eq. (8) vanishes, and the volume term can be seen as a scalar product which introduces a linear operator $\hat{L}$ such that

[TABLE]

The above quantity has the physical meaning of the free-charge density $\rho(\textbf{r})$ described by the electrostatic potential $\psi(\textbf{r})$. Assuming these free charges can exist only within the system volume $V_{0}$ (included in $V$), results in the eigenvalue equation

[TABLE]

where $\Lambda$ is the eigenvalue of the operator $\hat{L}$ corresponding to the eigenstate $\psi_{\lambda}$ . Equation (11) together with the boundary condition Eq. (9) thus form a generalized Sturm-Liouville problem. We note that by its physical meaning, the electrostatic potential $\psi_{\lambda}$ must be continuous across the system boundary. However, its spatial derivative does usually have a break across the system boundary which may lead to a presence of surface charges on the boundary, as discussed in detail in Sec. III below.

From Eq. (11) follows the orthonormality of static modes,

[TABLE]

where $\delta_{\lambda\lambda^{\prime}}$ is the Kronecker delta, since $\Lambda\neq\Lambda^{\prime}$ for different modes, or in case of degeneracy the modes can be made orthogonal by symmetry. Furthermore, taking the same scalar product with the wave functions of the RSs, we see that static modes are orthogonal to all of the RSs of the optical system:

[TABLE]

due to the boundary condition Eq. (9) for static modes and $\nabla\cdot\textbf{D}_{n}=0$ for the RSs. Note that this argument is valid for a finite volume $V$ only. For $V\to\infty$, the integral Eq. (13) may diverge. However, a proper orthogonality condition for this case, involving surface integrals, is provided in Ref. DoostPRA14, .

Finally, we note that in addition to the volume charge density of a static mode $\lambda$, given by

[TABLE]

which is present only within the system volume $V_{0}$, there is a surface charge density

[TABLE]

on the inner side of the boundary $S_{V}$, which is due to the fact that the normal component of the displacement $\textbf{D}_{\lambda}$ is discontinuous across the boundary, owing to the boundary condition Eq. (9). The volume $V$ can be any, and if it coincides with the system volume $V_{0}$, there is a surface charge density on the system boundary, associated with each static mode. If instead $V\to\infty$, the surface charge density $\sigma_{\lambda}\to 0$, due to $\textbf{D}_{\lambda}^{S}$ vanishing quick enough at $\textbf{r}\to\infty$. Below we consider these situations in detail for the analytically solvable case of a dielectric sphere.

III Static modes of a dielectric sphere in vacuum

Consider a dielectric sphere in vacuum, having radius $R$ and permittivity $\epsilon$. The system is described by the permittivity tensor $\hat{\boldsymbol{\varepsilon}}(\textbf{r})=\varepsilon(r)\hat{\bm{1}}$, where

[TABLE]

and $\Theta(x)$ is the Heaviside step function. Let us initially take the volume $V$ in Eq. (8) to be the infinite volume of the full 3D space. Owing to the spherical symmetry of the system, we can make the ansatz

[TABLE]

where $Y_{lm}(\Omega)$ are spherical harmonics (for definition see Appendix A). For a given fixed angular momentum $l$, the eigenvalue problem Eq. (11) then takes the form

[TABLE]

and

[TABLE]

where

[TABLE]

and $l$ is the angular momentum. Here, we have redefined for convenience the eigenvalue according to $\Lambda=\epsilon\lambda^{2}$.

Using the boundary conditions that the electrostatic potential $\psi_{\lambda}(\textbf{r})$ is continuous everywhere, finite at $\textbf{r}\to 0$, and vanishing at $\textbf{r}\to\infty$, we find the solution of the radial equations (18) and (19) with the permittivity Eq. (16) in the following form:

[TABLE]

where $j_{l}(x)$ is the spherical Bessel function of order $l$ and $A_{\lambda}$ is the normalization constant.

To find the secular equation determining the eigenvalues $\lambda$ we require that the scalar product Eq. (8) between different static modes $\lambda\neq\lambda^{\prime}$ is vanishing: $I_{\lambda\lambda^{\prime}}=0$. For the potentials Eq. (17) we obtain

[TABLE]

where $\textbf{Y}_{ilm}(\Omega)$ are the vector spherical harmonics (for their definition and properties, see Ref. LobanovPRA18, and Appendix A). Using their orthonormality, we obtain

[TABLE]

where

[TABLE]

and $j^{\prime}_{l}(x)$ is the derivative of the spherical Bessel function.

The secular equation (23) determines the eigenvalues $\lambda$ of the static modes. Obviously, it is fulfilled if $F_{\lambda}=0$ or $G_{\lambda}=0$, which are the two special cases determining, respectively, the VSC and VC static modes, considered below in detail. A more general solution of Eq. (23) is given by

[TABLE]

with arbitrary $\lambda$-independent constants $\alpha$ and $\beta$. The condition Eq. (26) can also be written as

[TABLE]

where $R_{\pm}=R\pm 0_{+}$ with a positive infinitesimal $0_{+}$.

Finally, the normalization constants $A_{\lambda}$ are found from the diagonal elements $I_{\lambda\lambda}$, defined by Eq. (8), which have the following explicit form

[TABLE]

with an analytical integral in the second term.

III.1 Volume-charge (VC) static modes

The set of VC static modes is generated by the condition

[TABLE]

with $G_{\lambda}$ given by Eq. (25), which determines the eigenvalues $\lambda$. This condition physically implies that Maxwell’s boundary condition of the continuity of the normal component of the displacement $\textbf{D}_{\lambda}$ across the sphere boundary is fulfilled, meaning that there are no free surface charges carried by the static modes. This can be easily seen from Eq. (9) used for $\alpha=0$ and $\beta=1$. The (volume) charge density is given by Eq. (14):

[TABLE]

where $\psi_{\lambda}(\textbf{r})$ is provided by Eqs. (17) and (21). We therefore call this set of static modes the volume-charge (VC) basis. The normalization condition Eq. (28) reduces to the second term only, and the normalization constants are given by

[TABLE]

III.2 Volume-surface charge (VSC) static modes

The set of VSC static modes is generated by the condition

[TABLE]

with $F_{\lambda}$ given by Eq. (24), determining another set of eigenvalues $\lambda$. This condition, in turn, splits into two cases:

[TABLE]

The first one determines an infinite set of modes with $\lambda\neq 0$ found from zeros of the Bessel function $j_{l}(x)$, and the wave functions are given by

[TABLE]

The second one is nothing else than the single static mode introduced in Ref. DoostPRA14, for each $l$. Its wave function can be obtained by taking the limit $\lambda\to 0$ in Eq. (21), which gives

[TABLE]

Both type of modes, $\lambda=0$ and $\lambda\neq 0$, violate Maxwell’s boundary condition for $\textbf{D}_{\lambda}$ at $r=R$ which implies the presence of a free surface charge for each static mode at the sphere boundary. The electric charge density in this case is given by

[TABLE]

where the first term, having the same form as in Eq. (30), corresponds to the volume charge, while the second term describes the surface charge. We therefore call this set of modes the volume-surface charge (VSC) basis.

For $\lambda\neq 0$ modes, the normalization is given by Eq. (31), which can be simplified, using Eq. (33), to

[TABLE]

The $\lambda=0$ mode normalization is instead given by

[TABLE]

which is produced by the surface term in Eq. (28), evaluated at $\lambda\to 0$. Using these normalization constants allows us to obtain explicit expressions for the surface charge density

[TABLE]

which appears in Eq. (36):

[TABLE]

We see that all $\lambda\neq 0$ modes are zero outside the volume $V_{0}$ of the dielectric sphere and produce volume and surface charges, while the $\lambda=0$ mode is non-zero in all space and produces only a surface charge on the sphere surface $S_{V_{0}}$.

As it was shown in Ref. DoostPRA14, , to use merely the $\lambda=0$ modes, having only the surface free charge as it is clear from Eq. (36), is sufficient for treating a homogeneous perturbation across the sphere with a step at $r=R$, which creates effective charges only at the sphere surface. Since this step is present in most applications using the homogeneous sphere as basis, it is expected to be more computationally efficient to use the VSC basis of static modes, as they provide a surface charge density on the surface of the basis sphere, which in the VC basis is harder to reproduce.

IV Application to systems with scalar dielectric susceptibility

In this section we consider the application of the RSE with both VS and VSC sets of static modes to various perturbations of a dielectric sphere. The perturbed system is described by a scalar permittivity, $\hat{\boldsymbol{\varepsilon}}(\mathbf{r})+\Delta\hat{\boldsymbol{\varepsilon}}(\mathbf{r})=\hat{\mathbf{1}}[\varepsilon(\mathbf{r})+\Delta\varepsilon(\mathbf{r})]$, with $\Delta\varepsilon(\mathbf{r})$ being the perturbation of the dielectric constant. As basis system we use a homogeneous dielectric sphere of radius $R$ and $\epsilon=4$. We consider several types of perturbations, namely, a homogeneous increase of the permittivity of the sphere in Sec. IV.1, a size reduction of the sphere in Sec. IV.2, a quarter-sphere perturbation in Sec. IV.3, and a deformation of the sphere into a cylinder in Sec. IV.4.

We use the standard formalism of the RSE, as described in MuljarovEPL10 ; DoostPRA12 ; DoostPRA13 ; DoostPRA14 ; MuljarovOL18 ; LobanovPRA18 . Some details of the technical implementation and optimization of the inclusion of a large number of static modes are given in Appendix B.

IV.1 Homogeneous sphere perturbation

The perturbation we consider here is a homogeneous change of the permittivity over the whole sphere, given by

[TABLE]

with the strength of $\Delta\epsilon=5$ used in the numerical calculation. This perturbation is spherically symmetric, so that RSs of different angular quantum numbers $(l,m)$, and different transverse polarization (TE or TM) do not mix, and are degenerate in $m$, making the RSE problem effectively 1D. We show in Fig. 1, for illustration, RSs with $l=5$, as in Ref. DoostPRA14, . Spherically symmetric perturbations do not couple static modes to TE RSs, since for the latter only the $\textbf{Y}_{1lm}$ component is non-vanishing in the basis of the vector-spherical harmonics (see LobanovPRA18 for details), while the static modes have the $\textbf{Y}_{1lm}$ component vanishing, see Eq. (LABEL:Elam). We therefore show here only the TM RSs, which couple to the static modes by a spherically symmetric perturbation, as both types of fields have a non-vanishing radial component. For the VSC set, only the $\lambda=0$ mode couples to the TM RSs by the perturbation Eq. (40). Indeed, the matrix elements between TM RSs and all other statics modes are proportional to the volume integral in Eq. (13) and thus vanish. For the VC set instead, all static modes are coupled to the TM RSs. The perturbed TM RSs obey the same secular equation as the basis system (see Eq.(31) in Ref. DoostPRA14, ) with the new permittivity of the sphere $\epsilon+\Delta\epsilon=9$, so that the new wave numbers $k$ calculated using the RSE can be compared with the exact values $k^{\mathrm{(exact)}}$ obtained from that secular equation.

Without any static modes, the error in the wave numbers $k$ of the perturbed RSs remains large, in the $10^{-1}$ to $10^{-2}$ range, see Fig. 1(b), consistent with the results shown in Ref. DoostPRA14, . Using the VSC set, the error of $k$ is reduced to values below $10^{-7}$ for $Rk_{\mathrm{max}}=4096$. We find that as we increase $N$, the error \bigl{|}k-k^{\mathrm{(exact)}}\bigr{|} scales as $k_{\mathrm{max}}^{-3}$. Figure 1(a) shows the resulting new wave numbers $k$ and Fig. 1(b) shows their errors for $Rk_{\mathrm{max}}=$ 64, 256, 1025, 4096 (corresponding to $N_{1}=$ 40, 164, 652, 2608, respectively) and $Rk_{\mathrm{max}}^{\mathrm{S}}=$ 397, 1586, 6344, 25377 (corresponding to $N_{2}=$124, 503, 2017, 8076, respectively), see Appendix B for details. Note, however, that the only static mode relevant for this perturbation is $\lambda=0$, as explained above (but technically we include here for consistency the full VSC set). Following the procedure described in Ref. DoostPRA12, , we extrapolate the new wave numbers to infinite $k_{\mathrm{max}}$ yielding $k^{(\infty)}$. We find that this extrapolation provides 1 to 2 orders of magnitude further reduction of the error, see Fig. 1(b).

Using instead the VC static modes (with $Rk_{\mathrm{max}}^{\mathrm{S}}$=430, 1718, 6873, 27492, corresponding to $N_{2}=$134, 544, 2185, 8748, respectively), the error of $k$, especially for RSs with small wave numbers (see Fig. 1(b)) is reduced to only $10^{-4}$. Again, the extrapolation to $k^{(\infty)}$ provides 1 to 2 orders of magnitude reduction of the error. We find that the error for small $|Rk|$, which is dominated by the static modes, scales as $k_{\mathrm{max}}^{-1}$, much slower than with the VSC basis. This is a consequence of the slow convergence in description of the effect of a surface charge induced by the perturbation (represented by a delta function $\delta(R-r)$) when no surface charges are present in the basis, which is the case of the VC set.

IV.2 Size perturbation of a sphere

We now consider a perturbation amounting to a size reduction of the sphere by 20%, given by

[TABLE]

This perturbation is also spherically symmetric, and we again show $l=5$ TM RSs in Fig. 2. The perturbed RSs obey the same secular equation for a sphere (Eq.(31) in Ref. DoostPRA14, ), with the radius reduced to $0.8R$, so that the new wave numbers $k$ calculated using the RSE can be again compared with the exact values $k^{\mathrm{(exact)}}$. The RSE wave numbers for $Rk_{\mathrm{max}}=4096$ are shown in Fig. 2(a). Without static modes, the error stays large, above $10^{-2}$, similar to the homogeneous perturbation. Adding either the VSC or the VC set, the error of $Rk$ is reduced to below $10^{-4}$, and the extrapolation to $k^{(\infty)}$ provides 1 to 2 orders of magnitude further reduction of the error. We find that for large $k_{\mathrm{max}}$, the relative error scales as $(k_{\mathrm{max}}^{\mathrm{S}})^{-1}$ in both cases. This is due to the step of the perturbation at $r=0.8R$, which leads to a surface charge inside the basis sphere which is not provided by a single state in both sets. Notably, the error using the VC set shows oscillations versus Re$(Rk)$, which are absent for the VSC set. This is attributed to the contributions of the two surface charges at $r=R$ and $r=0.8R$, which lead to an interference in the matrix elements in the VC set. For the VSC set instead, the surface charge at $r=R$ is coupling only to the $\lambda=0$ mode, and the surface charge at $r=0.8R$ is coupling only to the $\lambda\neq 0$ modes, so that no interference is present.

IV.3 Quarter-sphere perturbation

We now revisit the example of the quarter-sphere perturbation treated in Ref. DoostPRA14, . This perturbation breaks the continuous rotational symmetry for both spherical angles $\theta$ and $\varphi$ and is thus not reducible to an effective one- or two-dimensional system. The perturbation is given by

[TABLE]

and corresponds physically to a uniform increase of the dielectric constant in a quarter-sphere volume, as sketched in Fig. 3. For the results shown here, we use $\Delta\epsilon=1$, as in Ref. DoostPRA14, . The calculation of the matrix elements requires numerical integration, as detailed in Ref. DoostPRA14, . The perturbation mixes RSs of different $l$, $m$, and polarization. The remaining mirror symmetry $\varphi\rightarrow-\varphi$ of the system decouples $m\geqslant 0$ TE and $m<0$ TM RSs, having fields of odd parity ($-$), from $m<0$ TE and $m\geqslant 0$ TM RSs, having fields of even parity ($+$). For this calculation we use the VSC basis, split according to the parity of the fields, with the same selection rules as TM RSs.

The lifting of the $m$-degeneracy of the new RSs provides a splitting of resonances in Fig. 3(a) and (b). An analytic solution for this perturbation is not available, so that we estimate the error using as exact solution the extrapolated values of the largest $k_{\mathrm{max}}$. A convergence with a power law around $(k_{\mathrm{max}})^{-3}$ is observed, resulting in relative errors in the $10^{-4}$ to $10^{-5}$ range. The irregular arrangement of the perturbed whispering-gallery modes (WGMs) in the right-bottom part of Fig. 3(a) is due to this remaining small error which manifests itself in fluctuating $\mathrm{Im}\,k$ well seen in the log scale used.

Exactly the same system was treated in Ref. DoostPRA14, , using only the $\lambda=0$ static modes, and was compared with FEM calculations using the commercial solver COMSOL, see Fig. 4 of Ref. DoostPRA14, . These COMSOL results are also shown in the present Fig. 3 with $N_{\rm G}=25$k, 50k, 100k and 200k finite elements. With increasing $N_{\rm G}$, the COMSOL wave numbers converge, with an error scaling approximately as $N_{\rm G}^{-1}$. However, it was found in Ref. DoostPRA14, that increasing the precision of the RSE and COMSOL, a relative difference remained in the few times $10^{-4}$ range. One possible origin of this systematic deviation could be due to the incomplete static-mode basis used in Ref. DoostPRA14, . We have therefore repeated the RSE calculation using the VSC basis, as shown here. We have extrapolated both the COMSOL and the RSE results, as shown in the insets of Fig. 3(b). We find that the remaining relative difference between the extrapolated COMSOL and RSE results is in the $10^{-4}$ range, which is still somewhat larger that the estimated error of the extrapolations shown in Fig. 3(c). Nevertheless, the deviation between FEM and RSE results reported in Ref. DoostPRA14, was partly due to the incomplete basis used in that work. To explicitly identify the role of the static modes, we show in the insets of Fig. 3(b) also the RSE results using no static mode, or using only the $\lambda=0$ static modes as in Ref. DoostPRA14, . We see that the full effect of the static modes on the perturbed RSs shown is about $5\times 10^{-4}$ relative change.

IV.4 Sphere to cylinder perturbation

We consider here a perturbation which transforms a sphere into a cylinder of a height $2h$ equal to its diameter $2a$. It retains axial symmetry and inversion symmetry, and thus represents an effective 2D system. The perturbation is given by

[TABLE]

with $R=\sqrt{h^{2}+a^{2}}$ being the radius of the basis sphere. Owing to the axial symmetry and the mirror symmetry, the RSs of the cylinder have well defined $m$, and parity ($+,-$) of the field under the $z\rightarrow-z$ mirror imaging. The full perturbation matrix of the RSE is thus separated into blocks corresponding to these quantum numbers, significantly reducing the size of the eigenvalue problem. Also note that the RSs of $m$ and $-m$ are degenerate due to the mirror symmetry $\varphi\rightarrow-\varphi$.

The basis sphere of radius $R$ has the same permittivity and just encloses the cylinder. To calculate the values of the symmetry-allowed matrix elements of the perturbation, we discretize the integrals into shell-segments. Shells of radial thickness $R/s$ are used, where $s$ is the number of shells. The polar angular ranges for each shell are then determined using the intersection of the center radius of the shell with the cylinder surface. The number of shells is chosen using the Nyquist criterion of sampling, with 2 shells per period of $k_{\mathrm{max}}$, $s=\sqrt{\epsilon}k_{\mathrm{max}}R/\pi$. This choice ensures that with increasing $k_{\mathrm{max}}$, the spatial sampling of the structure is refined according to the spatial resolution of the basis RSs.

The resulting RSs with $m=0,1$, and 7, and parity $+$ and $-$ are displayed in Fig. 4(a) for the largest basis size considered, $Rk_{\mathrm{max}}=75$ and $Rk_{\mathrm{max}}^{\mathrm{S}}=258$. The RSs with large losses, $-\mathrm{Im}(Rk)>1.5$, are called leaky modes. For the other RSs, we see that for small $m$-values ($m=0$ and 1), the losses $-\mathrm{Im}(Rk)$ are in the 0.1 range, increasing with $k$ their spread towards smaller losses and reaching Q-factors above 1000, where $Q=-\mathrm{Re}\,k/(2\mathrm{Im}\,k)$. The low-loss RSs for these small values of $m$ can be formed at large $k$ by a field localized away from the edges and experiencing total internal reflection at the surfaces (center top and bottom, and center of the cylinder side wall). For the larger $m=7$ instead, we find that the lowest-frequency RSs have the character of WGMs, with $-\mathrm{Im}(Rk)<10^{-3}$ and Q-factors around $10^{4}$. The error of the wave numbers, defined as the absolute difference to the values obtained for $Rk_{\mathrm{max}}=75$, are given in Fig. 4(b), not including the leaky modes. We find that RSs with $|Rk|\lesssim 5$ have the errors in the $10^{-2}$ range for $Rk_{\mathrm{max}}=20$, and in the $10^{-3}$ range for $Rk_{\mathrm{max}}=41$, thus scaling approximately as $k_{\mathrm{max}}^{-3}$, similar to the example with the quarter-sphere perturbation discussed above. Notably, as known about the RSE MuljarovEPL10 , for a given $k_{\mathrm{max}}$, the error is increasing with increasing wave number, limiting the range of small error to RSs with $|k|\lesssim k_{\mathrm{max}}/4$. We can see that within this range, the RSE is determining a large number of RS simultaneously, a few hundred for $Rk_{\mathrm{max}}=75$.

To verify the wave numbers obtained for the cylinder using the RSE with static modes, we concentrate on the RS with $m=1$ and odd parity which has been calculated in Ref. MuljarovPRB16, using COMSOL, with a resonance frequency of $Rk=4.16275-0.24382i$ and an estimated error of $0.002$ in real and imaginary part. We show in Fig. 5 the wave number of this RS, and compare it with the results of the RSE using the VSC or the VC basis, as functions of $k_{\mathrm{max}}$, as well as with their extrapolated values including errors (see Ref. DoostPRA12, for the extrapolation procedure). We find that the convergence is somewhat different for the VC and VSC static basis, but the extrapolated values are equal within the estimated errors. The COMSOL result is in agreement with the extrapolated RSE within the estimated error, indicating that including the used sets of static modes indeed provides a complete basis for general 3D geometries.

V Summary

In summary, we have introduced a general procedure for determining a full set of static modes for an arbitrary open optical system, supplementing the set of RSs, in order to form a complete basis for the RSE. Including this set is required to achieve high accuracy in the RSE in most cases, and specifically for general three-dimensional perturbations. We have shown in particular that static modes are required to treat perturbations which have a gradient of the permittivity in the direction of the RS electric field, thus creating effective free charges which are incompatible with Maxwell’s boundary conditions of the basis system, and therefore not described by the RSs alone. Two alternative sets of static modes have been introduced for a dielectric sphere in vacuum, with one of them including free charges on the surface of the sphere, and thus being somewhat advantageous when treating perturbations with a permittivity step at the sphere surface. Using this complete basis in the RSE was then shown to result in convergence towards the analytical RSs of a reduced-size sphere, and towards the numerically determined RSs of a cylinder and a sphere with a quarter-sphere perturbation. Using the RS basis extended by including either of the two static mode sets, the RSE is expected to be numerically exact for general three-dimensional confined geometries.

Acknowledgements.

This work was supported by the EPSRC under grant EP/M020479/1. S. L. and E. A. M. acknowledge support from RFBR (Grant No. 16-29-03283).

Appendix A Scalar and vector spherical harmonics

Following DoostPRA14 , we define the scalar spherical harmonics $Y_{lm}(\Omega)$ as real functions:

[TABLE]

where $\Omega=(\theta,\varphi)$ in the spherical coordinate system, $P^{m}_{l}(x)$ are the associated Legendre polynomials, $l$ and $m$ are the spherical quantum numbers, and

[TABLE]

This is done in order to satisfy the orthonormality condition without using the complex conjugate:

[TABLE]

where $d\Omega=\sin\theta d\theta d\varphi$.

The vector spherical harmonics are defined following LobanovPRA18 as

[TABLE]

They satisfy the following orthonormality condition:

[TABLE]

Appendix B RSE with static modes

In Ref. DoostPRA14, , the RSE equation given by Eq.(13) therein was used to take into account the LE ($\lambda=0$) static mode, using a non-zero wave vector $Rk=10^{-7}$, since it was numerically more efficient than solving the generalized eigenvalue problem given by Eq.(12) therein. As we are now dealing with a large number of static RS, we instead follow a different approach, exploiting the zero wave vector of all static modes. We start with the generalized eigenvalue problem, given by Eq.(12) of Ref. DoostPRA14, , written in matrix form as LobanovPRA18

[TABLE]

where $\mathbbm{k}$ is a diagonal matrix containing the wave numbers of the basis RSs and static modes, $\mathbbm{M}=\mathbbm{1}+\mathbbm{V}$ with $\mathbbm{1}$ being the unit matrix and $\mathbbm{V}$ the perturbation matrix. $\varkappa$ and $\mathbf{c}$ are, respectively, the eigenvalue and the eigenvector of a perturbed RS. The matrix elements of the perturbation matrix $\mathbbm{V}$ are given by

[TABLE]

with $i$ and $j$ labeling both the RSs and the static modes of the unperturbed system. The RSs are normalized in accordance with MuljarovOL18 , with the already mentioned factor of 2 difference compared to the earlier used normalization.

We now separate the notation explicitly into the $N_{1}$ RSs and $N_{2}$ static modes, so that the matrices $\mathbbm{k}$ and $\mathbbm{M}$ split into four sub-matrices, with $N_{1}$-dimensional square top-left sub-matrices $\mathbbm{k}_{1}$ and $\mathbbm{M}_{11}$ corresponding to the RSs. The eigenvector $\mathbf{c}$ also splits into two sub-vectors, with an $N_{1}$-dimensional top sub-vector $\mathbf{c}_{1}$ corresponding to the RSs and $N_{2}$-dimensional bottom sub-vector $\mathbf{c}_{2}$ corresponding to the static modes. The matrix equation (51) then reads

[TABLE]

where $\mathbbm{M}_{12}$ is the transpose of $\mathbbm{M}_{21}$. The above matrix equation then splits into a pair of matrix equations,

[TABLE]

which can be written as an $N_{1}\times N_{1}$ eigenvalue problem for the RSs,

[TABLE]

with the effective perturbation matrix

[TABLE]

and an auxiliary equation for the static mode amplitudes

[TABLE]

The reduced generalized eigenvalue problem Eq. (56) for the RSs can be further modified to the standard eigenvalue problem solved by matrix diagonalization, as it was done in MuljarovEPL10 ; DoostPRA14 . Calculating the effective perturbation matrix $\widetilde{\mathbbm{M}}_{11}$ requires inversion of a $N_{2}\times N_{2}$ matrix $\mathbbm{M}_{22}$ with the compute time of about $T^{\mathrm{inv}}=T^{\mathrm{inv}}_{0}(N_{1}N_{2})^{3/2}$, and we find using Matlab 2017a on 16 CPU cores (dual Intel E5-2640 v3) the value $T^{\mathrm{inv}}_{0}=7.5$ ps. Solving Eq. (56) involves a matrix diagonalization which requires a compute time of about $T^{\mathrm{dia}}=T^{\mathrm{dia}}_{0}N_{1}^{3}$, with $T^{\mathrm{dia}}_{0}=115$ ps. For an equal number of RSs and static modes, the treatment of the static modes is thus about 15 times faster.

Now consider that the influence of static modes on a given $\varkappa$ is decreasing with increasing $\lambda$ only due to a reduction of the overlap matrix elements in $M$, owing to the increasing spatial frequency of the static RSs, see Eq. (52). For the RSs instead, there is additionally an increase of the wave number difference $|k_{n}-\varkappa|$, helping to reduce their influence with increasing $k_{n}$, see Eq. (38) of Ref. DoostPRA14, . It is therefore beneficial for the accuracy of the RSE at a given compute time to use the lower numerical complexity of treating static modes to increase their number. We balance the numerical complexity of static modes and RSs in the RSE by choosing $T^{\mathrm{dia}}\approx T^{\mathrm{inv}}$, adjusting $N_{1}$ and $N_{2}$ appropriately.

To choose the basis RSs we use a cut-off in their spatial frequency inside the system, $|\sqrt{\epsilon}k_{n}|<k_{\mathrm{max}}$, with the maximum wave number $k_{\mathrm{max}}$ in vacuum, and the refractive index $\sqrt{\epsilon}$ of the basis sphere. The number of RSs is then approximately given by $N_{1}\approx C_{1}^{(p)}\left(\sqrt{\epsilon}Rk_{\mathrm{max}}\right)^{d}$, with the dimensionality $d\in\{1,2,3\}$, and the polarization $p\in\{{\rm TE,TM}\}$. The values of $C_{1}^{(p)}$ determined for $\epsilon=4$ and $Rk_{\mathrm{max}}\gg 1$ are given in Table 1.

For the static modes, we use $\lambda<k_{\mathrm{max}}^{\mathrm{S}}$, with a separate maximum wave number $k_{\mathrm{max}}^{\mathrm{S}}$. The $\lambda=0$ static modes are always included in the VSC set. The resulting number of static modes is $N_{2}\approx C_{2}(Rk_{\mathrm{max}}^{\mathrm{S}})^{d}$, with the values of $C_{2}$ given in Table 1, determined for $Rk_{\mathrm{max}}^{\mathrm{S}}\gg 1$. We note that the scaling constants are not significantly different between VSC and VC sets. Requiring $T^{\mathrm{dia}}=T^{\mathrm{inv}}$, we then find

[TABLE]

where $C_{1}=C_{1}^{\rm TM}$ for the 1D TM case, $C_{1}=C_{1}^{\rm TE}$ for the 1D TE case, and $C_{1}=C_{1}^{\rm TM}+C_{1}^{\rm TE}$ for the 2D and 3D case. We use this relation between $k_{\mathrm{max}}^{\mathrm{S}}$ and $k_{\mathrm{max}}$ for all numerical results reported in the present work.