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

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.


I. INTRODUCTION
The electromagnetic spectrum of an open optical system is characterized by its resonances, which is evident for optical cavities such as dielectric toroid [1] or microsphere resonators [2]. Resonances are characterized by their spectral positions and linewidths, corresponding to, respectively, the real and imaginary part of the complex eigenfrequencies of the system. Finite linewidths of resonances are typical for open systems and are due to energy leakage from the system to the outside. Objects in close proximity of the cavity modify the electromagnetic susceptibility and perturb the cavity resonances, changing both their position and linewidth, most noticeably for the high-quality (i.e. narrow-linewidth) resonances. This effect is the basis for resonant optical biosensors [3][4][5] in which the changes in the spectral properties of resonators in the presence of perturbations can be used to characterize the size and shape of attached nanoparticles [6]. The whispering gallery mode (WGM) resonances in microdisks and spherical microcavities have been used in sensors for the characterization of nanolayers [7], protein [8] and DNA molecules [9], as well as for single atom [10] and nanoparticle detection [11,12]. Furthermore, the long photon lifetime of WGMs can result in their strong coupling to atoms [13]. Recently, optical resonances have become the core element of a more accurate modeling of multimode and random lasers [14,15] and of light propagation through random media [16]. In nanoplasmonics, the resonances of metal nanoparticles are used to locally enhance the electromagnetic field [17].
Due to the lack of a suited theory, the electromagnetic properties of such open systems were up to now modeled by using finite element method (FEM) and finite differ-ence in time domain (FDTD) solvers. Only recently, approximate approaches using resonance modes have been reported [18][19][20][21][22]. While the eigenmodes of resonators for a few highly symmetric geometries can be calculated exactly, determining the effect of perturbations which break the symmetry presents a significant challenge as the popular computational techniques in electrodynamics, such as the FDTD [23] or FEM [24], need large computational resources [25] to model high quality WGMs.
To treat such perturbations more efficiently, we have developed [26] a rigorous perturbation theory called resonant state expansion (RSE) and applied it to spherical resonators reducible to effective one-dimensional (1D) systems. We have demonstrated on exactly solvable examples in 1D that the RSE is a reliable tool for calculation of wavenumbers and electromagnetic fields of resonant states (RSs) [27], as well as transmission and scattering properties of open optical systems. We have recently developed the RSE also for effectively two-dimensional (2D) systems [28], and planar waveguides [29].
In this paper we extend the RSE formulation to arbitrary three-dimensional (3D) open optical systems, compare its performance with FDTD and FEM, and introduce a local perturbation approach. The paper is organized as follows. In Sec. II we give the general formulation of the RSE for an arbitrary 3D system. In Sec. III we treat the homogeneous dielectric sphere as unperturbed system and introduce the basis for the RSE, which consists of normalized transverse electric (TE) and transverse magnetic (TM) modes and is complemented by longitudinal zero frequency modes. This is followed by examples given in Sec. IV A-C illustrating the method and comparing results with existing analytic solutions, as well as numerical solutions provided by using available commercial software. In Sec. IV D we demonstrate the performance of the RSE as a local perturbation method for a chosen group of modes by introducing a way to select a suitable subset of basis states. Some details of the general formulation of the method including mode normalization and calculation of the matrix elements are given in Appendices A and B.

II. RESONANT STATE EXPANSION
Resonant states of an open optical system with a local time-independent dielectric susceptibility tensorε(r) and permeability µ = 1 are defined as the eigensolutions of Maxwell's wave equation, satisfying the outgoing wave boundary conditions. Here, k n is the wave-vector eigenvalue of the RS numbered by the index n, and E n (r) is its electric field eigenfunction in 3D space. The time-dependent part of the RS wave function is given by exp(−iω n t) with the complex eigenfrequency ω n = ck n , where c is the speed of light in vacuum. As follows from Eq. (1) and the divergence theorem, the RSs are orthogonal according to where the first integral in Eq. (2) is taken over an arbitrary simply connected volume V which includes all system inhomogeneities ofε(r) while the second integral is taken over the closed surface S V , the boundary of V , and contains the gradients ∂/∂s normal to this surface. The RSs of an open system form a complete set of functions. This allows us to use RSs for expansion of the Green's function (GF)Ĝ k (r, r ′ ) satisfying the same outgoing wave boundary conditions and Maxwell's wave equation with a delta function source term, where1 is the unit tensor and k = ω/c is the wave vector of the electromagnetic field in vacuum determined by the frequency ω, which is in general complex. The GF expansion in terms of the direct (dyadic) product of the RS vector fields is given by Ref. [28] This expansion requires that the RSs are normalized according to where E(k, r) is an analytic continuation of the RS wave function E n (r) around the point k n in the complex kplane and δ kn,0 is the Kronecker delta accounting for a factor of two in the normalization of k n = 0 modes. For any spherical surface S R of radius R, the limit in Eq. (5) can be taken explicitly leading for k n = 0 modes to where r = |r|, with the origin at the center of the chosen sphere. Static k n = 0 modes, if they exist in the GF spectrum, are normalized according to 2 = drE n ·εE n .
Their wave functions decay at large distances as 1/r 2 or quicker, and the volume of integration in Eq. (5) can be extended to the full space for which the surface integral is vanishing. The proofs of Eqs. (5) and (6) are given in Appendix A. The completeness of RSs allows us to treat exactly a modified (perturbed) problem in which the RS wave vector κ ν and the electric field E ν are modified as compared to k n and E n , respectively, due to a perturbation ∆ε(r) with compact support. We treat this problem by (i) solving Eq. (8) with the help of the GF, (ii) using in Eq. (9) the spectral representation Eq. (4), and (iii) expanding the perturbed wave functions into the unperturbed ones, This is the RSE method. The use of the of the unperturbed GF is an essential element of the RSE as Eq. (9) guarantees that the perturbed wave functions satisfy the outgoing boundary condition. The result of using Eq. (11) in Eq. (10) is a linear matrix eigenvalue problem which is reduced, using a substitution b nν = c nν κ ν /k n , to the matrix equation [26] n ′ This allows us to find the wave vectors κ ν and the expansion coefficients c nν of the perturbed RSs by diagonalizing a complex symmetric matrix. The matrix elements of the perturbation are given by In our previous works on RSE [26,28] we derived the intermediate result Eq. (10) using Dyson's equation for the perturbed GF. The present way to obtain Eq. (10) is equivalent, but is simplifying the treatment by not dealing explicitly with the perturbed GF. We note that in 2D systems the set of RSs of a system is complemented with a continuum of states on the cut of the GF [28]. In this case, all summations in the above equations include states on the cut which are discretized in numerics to produce a limited subset of isolated poles.

III. EIGENMODES OF A DIELECTRIC SPHERE AS BASIS FOR THE RSE
To apply the RSE to 3D systems we need a known basis of RSs. We choose here the RSs of a dielectric sphere of radius R and refractive index n R , surrounded by vacuum, since they are analytically known. For any spherically symmetric system, the solutions of Maxwell's equations split into four groups: TE, TM, and longitudinal electric (LE) and longitudinal magnetic (LM) modes [30]. TE (TM) modes have no radial components of the electric (magnetic) field, respectively. Longitudinal modes are curl free static modes satisfying Maxwell's wave equation for k n = 0. Longitudinal magnetic modes have zero electric field, and since we limit ourself in this work to perturbations in the dielectric susceptibility only, they are not mixed by the perturbation to other types of modes and are thus ignored in the following. Furthermore, owing to the spherical symmetry, the azimuthal index m and longitudinal index l are good quantum numbers of the angular momentum operator and take integer values corresponding to the number of field oscillations around the sphere. For each l value there are 2l + 1 degenerate modes with m = −l..l.
Splitting off the time dependence ∝ e −iωt of the electric fields E and D and magnetic field H, the first pair of Maxwell's equations can be written in the form where k = ω/c and D(r) =ε(r)E(r). Combining them leads to Eq. (1) for the RSs and to Eq.
is automatically satisfied, since ∇ × ∇ = 0. However, if k = 0, it is not guaranteed that solutions of Eq. (15) satisfy also Eq. (16). The spectrum of the GF given by LE : where f (r) is a scalar function satisfying the Helmholtz equation with the permeability of the dielectric sphere in vacuum given by Owing to the spherical symmetry of the system, the solution of Eq. (18) splits in spherical coordinates r = (r, θ, ϕ) into the radial and angular components: where Ω = (θ, ϕ) with the angle ranges 0 θ π and 0 ϕ 2π. The angular component is given by the spherical harmonics, which are the eigenfunctions of the angular part of the Laplacian,Λ where P m l (x) are the associated Legendre polynomials. Note that the azimuthal functions are defined here as in order to satisfy the orthogonality condition without using the complex conjugate, as required by Eq. (2). The radial components R l (r, k) satisfy the spherical Bessel equation, (24) and have the following form in which j l (z) and h l (z) ≡ h (1) l (z) are, respectively, the spherical Bessel and Hankel functions of the first kind.
In spherical coordinates, a vector field E(r) can be written as where e r , e θ , and e ϕ are the unit vectors. The electric field of the RSs then has the form for TE modes, for TM modes, and for LE modes. All the wave functions are normalized according to Eqs. (5)- (7), leading to the following normalization constants: The Maxwell boundary conditions following from Eq. (15), namely the continuity of the tangential components of E and H across the spherical dielectric-vacuum interface, lead to the following secular equations determining the RS wavenumbers k n : for TE modes and for TM modes, where z = k n R and j ′ l (z) and h ′ l (z) are the derivatives of j l (z) and h l (z), respectively. While the LE modes are the RSs easiest to calculate due to a simple power-law form of their radial functions, it is convenient to treat them in the RSE as part of the TM family of RSs. Indeed, for r R they coincide with the TM modes taken in the limit k n → 0: Note that k n = 0 is not a solution of the secular equation (31) for TM modes. However, using the analytic dependence of the wave functions of TM modes on k n [see Eqs. (25), (27), and (29)], the limit Eq. (33) can be taken in the calculation of the matrix elements containing LE modes. The same limit k n → 0 has to be carefully approached in the matrix eigenvalue problem Eq. (13) of the RSE, as the matrix elements are divergent, due to the 1/ √ k n factor introduced in the expansion coefficients. We found that adding a finite negative imaginary part to static poles, k n R = −iδ, with δ typically of order 10 −7 (determined by the numerical accuracy) is suited for the numerical results presented in the following section. We have verified this by comparing the results with the ones of the RSE in the form of a generalized linear eigenvalue problem Eq. (12), which has no such divergence, but its numerical solution is a factor of 2-3 slower in the NAG library implementation.

IV. APPLICATION TO 3D SYSTEMS WITH SCALAR DIELECTRIC SUSCEPTIBILITY
In this section we discuss the application of the RSE to 3D systems described by a scalar dielectric function ε(r) + ∆ε(r) =1[ε(r) + ∆ε(r)]. As unperturbed system we use the homogeneous dielectric sphere of radius R with ε(r) given by Eq. (19), having the analytical modes discussed in Sec. III. We use the refractive index n R = 2 of the unperturbed sphere throughout this section and consider several types of perturbations, namely, a homogeneous perturbation of the whole sphere in Sec. IV A, a half-sphere perturbation in Sec. IV B, and a quarter-sphere perturbation in Sec. IV C. We demonstrate in Sec. IV D the performance of the RSE as a local perturbation method for a chosen group of modes by introducing a way to select a suitable subset of basis states. Explicit forms of the matrix elements used in these calculations are given in Appendix B.

A. Homogeneous sphere perturbation
The perturbation we consider here is a homogeneous change of ε over the whole sphere, given by where Θ is the Heaviside function, with the strength ∆ǫ = 5 used in the numerical calculation. so that the perturbed modes obey the same secular equations Eq. (30) and Eq. (31) with the refractive index n R of the sphere changed to n 2 R + ∆ǫ, and the perturbed wavenumbers κ ν calculated using the RSE can be compared with the exact values κ (exact) ν obtained from the secular equations.
We choose the basis of RSs for the RSE in such a way that for a given orbital number l and m we select all RSs with |k n | < k max (N ) using a maximum wave vector k max (N ) chosen to result in N RSs. We find that as we increase N , the relative error κ ν /κ (exact) ν − 1 decreases as N −3 . Following the procedure described in Ref. 27 we can extrapolate the perturbed wavenumbers. The resulting perturbed wavenumbers for N = 1000 (corresponding to k max R = 800) are shown in Fig. 1 for the TM RSs and Fig. 2 for the TE RSs. The perturbation is strong, leading to WGMs with up to 2 orders of magnitude narrower linewidths. The RSE reproduces the wavenumbers of about 100 RSs to a relative error in the 10 −7 range, which is improving further by one to two orders of magnitude after extrapolation. The homogeneous perturbation does not couple LE modes to TE modes as LE modes have the symmetry of TM modes [see Eq. (33)] leading to vanishing overlap integrals with TE RSs. The contribution of the LE-mode RS in the TM polarization is significant, as is shown in Fig. 1 by the large decrease of the relative error by up to 8 orders of magnitude when adding them to the basis. This validates the analytical treatment of the LE-mode RSs in the RSE developed in this work. We have verified that taking a finite imaginary value of δ = 10 −7 in Eq. (13) for the LE-modes instead of using strict k n = 0 poles in Eq. (12), as done throughout this work, changes the relative error of the TM mode calculation by less than 10% and within the range of 10 −9 only. For practical applications, this limitation should not be relevant as the error in the measured geometry will typically be significantly larger.

B. Hemisphere Perturbation
We consider here a hemisphere perturbation as sketched in Fig. 3 which mixes TE, TM, and LE modes with different l, while conserving m. The perturbation is given by and increases ε in the northern hemisphere by ∆ǫ, while leaving the southern hemisphere unchanged. In our numerical simulation, we use ∆ǫ = 0.2. The calculation of the matrix elements is done using Eqs. (B7)-(B12) of Appendix B which require numerical integration. Owing to the symmetry of the perturbation, matrix elements between TM and TE RSs can only be non-zero when the RSs have m of opposite sign and equal magnitude, i.e. they are are sine and cosine states of equal |m|. Similarly, matrix elements between two TE RSs or two TM RSs can only be non-zero if both states have the same m. We can therefore restrict the basis to m = 3 TM states and m = −3 TE states for the numerical calculations of this section. We treat the LE RSs as TM modes with k n R = −i10 −7 and a normalization factor modified according to Eq. (33). The resulting RS wavenumbers are shown in Fig. 3. Due to the smaller perturbation compared to that considered in Sec. IV A, the mode positions in the spectrum do not change as much. The imaginary part of most of the WGMs decreases due to the higher dielectric constant in the perturbed hemisphere. However, some of the modes also have an increased imaginary part due to the scattering at the edge of the perturbation.
To the best of our knowledge, an analytic solution for this perturbation is not available and thus we cannot calculate the relative error of the RSE result with respect to the exact solution. However, we can investigate the convergence of the method in order to demonstrate how the RSE works in this case, for the perturbation not reducible to an effective one-dimensional problem. We accordingly show in Fig. 3(a) the perturbed modes for two different values of basis size N and in Fig. 3 are the RS wavenumbers calculated for basis sizes of N 1 ≈ N/2, We see that the perturbed resonances are converging with increasing basis size, approximately following a power law with an exponent between −2 and −3.

C. Quarter-Sphere Perturbation
We consider here a perturbation which breaks both continuous rotation symmetries of the sphere and is thus is not reducible to an effective one or two-dimensional system. The perturbation is given by and corresponds physically to a uniform increase of the dielectric constant in a quarter-sphere area, as sketched in Fig. 4. In our numerical simulation, we take ∆ǫ = 1. Again, the calculation of the matrix elements requires numerical integration. Owing to the reduced symmetry of the perturbation as compared to that treated in the previous section we now have modes of different l, m, and polarization mixing, although TE sine (TM cosine) and TE cosine (TM sine) modes are decoupled, owing to the mirror symmetry of the system. This allows us to split the simulation of all modes into two separate simulations called A and B, respectively, each of size N . The lifting of the m-degeneracy of the unperturbed modes can be seen as splitting off resonances in Fig. 4(a) and (b). In most cases the splitting in the real part of the resonant wavenumber is greater than the linewidth of the modes. The convergence of the RSE is well seen in Fig. 4(a) and (b) showing the perturbed RS wavenumbers for two different basis sizes N . An analytic solution for this perturbation is not available, so that we use the method described in Sec. IV B to estimate the error, and show in Fig. 4(c) the resulting absolute errors M ν for several values of N . A convergence with a power law exponent between −2 and −3 is again observed, resulting in relative errors in the 10 −4 to 10 −5 range for N = 8000.
To verify the RSE results, we have simulated the system using the commercial solver ComSol (http://www.comsol.com) which uses the finite element method and Galerkin's method, approximating the openness of the system with an absorbing perfectly matched layer (PML). We have surrounded the sphere with a vacuum shell followed by a PML shell of equal thickness D. The results are shown in Fig. 4(b) using D = R/2, and a "physics controlled" mesh with N G = 25k, 50k, 100k and 200k finite elements. We used the nearest unperturbed RS wave vector as linearization point (i.e. the input value) for the ComSol solver, and requested the determination of 40 eigenfrequencies, which we found to be the minimum number reliably returning all 15 nondegenerate modes deriving from the l = 7 unperturbed fundamental WGM. With increasing N G , the ComSol RS wavenumbers tend towards the RSE poles, with an error scaling approximately as N −1 G . This is verifying the validity of the RSE results.
To make a comparison between the RSE and Com-Sol in terms of numerical complexity we use the poles computed by an N = 16000 RSE simulation as "exact solution" to calculate the average relative errors of the poles shown in Fig. 4(b) versus effective processing time on an Intel E8500 CPU. The result is shown in Fig. 5, including ComSol data for different shell thicknesses D of R/2, R/4, and R/8, revealing that D = R/4 provides the best performance. This comparison shows that the RSE is 2-3 orders of magnitude faster than ComSol for  Fig. 4 the present example, and at the same time determines significantly more RSs. The RSE computing time includes the calculation of the matrix elements which were done evaluating the 1dimensional integrals (see Appendix B 2) using 10000 equidistant grid points. The computing time of the matrix elements is significant only for N 2000, while for larger N the matrix diagonalization time, scaling as N 3 , is dominating. We have verified that the accuracy of the matrix element calculation is sufficient to not influence the relative errors shown.
We also include in Fig. 5 the performance of FDTD calculations using the commercial software Lumerical (http://www.lumerical.com). They were undertaken using a simulation cube size from 2.5R to 4R, exploiting the reflection symmetry, and for grid steps between R/8 and R/80, with a sub-sampling of 32. The simulation area was surrounded by a PML of a size chosen automatically by the software. The excitation pulse had a center wavenumber of kR = 5.1 and a relative bandwidth of 10% to excite the relevant modes, and the simulation was run for 360 oscillation periods. The calculated timedependent electric field after the excitation pulse was transformed into a spectrum and the peaks were fitted with a Lorentzian to determine the real and imaginary part of the mode. The parameters used were chosen to optimize the performance, and in the plot the results with the shortest computation time for a given relative error are given.
We can conclude that the RSE is about two orders of magnitude faster than both FEM and FDTD for this specific problem, showing its potential to supersede presently used methods. A general analysis of the performance of RSE relative to FEM and FDTD is beyond the scope of this work and will be presented elsewhere.
To illustrate how a particular perturbed RS is created as a superposition of unperturbed RSs, we show in Fig. 6 the contributions of the unperturbed RSs to the perturbed WGM indicated by the arrow in Fig. 4(b) with index ν and wavenumber κ ν , given by the open star in Fig. 6. The contribution of the basis states to this mode are visualized by circles of a radius proportional to 6 |c nν | 2 , where the sum is taken over the 2l + 1 degenerate basis RSs of a given eigenfrequency, centered at the positions of the RS wavenumbers in the complex k-plane. The expansion coefficients c nν decrease quickly with the distance between the unperturbed and perturbed RS wavenumbers, with the dominant contribution coming from the nearest unperturbed RS, a typical feature of perturbation theory in closed systems. The unperturbed RS nearest to the perturbed one in Fig. 6 has the largest contribution, and is a l = 7 TE WGM with the lowest radial quantum number. Other WGMs giving significant contributions have the same radial quantum number and the angular quantum numbers ranging between l = 6 and l = 9, see the small stars in Fig. 6 corresponding to l = 7 basis states. This is a manifestation of a quasi-conservation of the angular momentum l for bulky perturbations like the quarter-sphere perturbation considered here.
Generally, we see that a significant number of unperturbed RSs are contributing to the perturbed RS, which is indicating that previous perturbation theories for open systems would yield large errors for the strong perturbations treated in this work since they are limited to low orders [31,32] or to degenerate modes only [33].

D. Local Perturbation
The weights of the RSs shown in Fig. 6 indicate that a perturbed mode can be approximately described by a subset of the unperturbed modes, which typically have wavenumbers in close proximity to that of the perturbed mode. It is therefore expected that a local perturbation approach based on the RSE is possible. We formulate here such an approach.
We commence with a small subset S of modes of the unperturbed system which are of particular interest, for example because they are used for sensing. To calculate the perturbation of these modes approximately, we consider a global basis B as used in the previous sections, with a size N providing a sufficiently small relative error. We then choose a subset S + ⊂ B with N ′ < N elements containing S, i.e. S ⊂ S + , and solve the RSE Eq. (13) restricted to S + . The important step in this approach is to find a numerically efficient method to choose the additional modes in S + which provide the smallest relative error of the perturbed states deriving from S for a given N ′ . Specifically, the method should be significantly faster than the matrix diagonalization Eq. (13).
To develop such a method, we consider here the Rayleigh-Schrödinger perturbation theory based on the RSE and expand the RS wave vector κ up to second order, (38) as directly follows from Eq. (13). Note that the secondorder result in Eq. (38) is different from that given in Ref. 31.
We expect that the second-order correction given by Eq. (38) is a suited candidate to estimate the importance of modes. We therefore sort the modes in B according to the weight W n given by where D is the set of modes degenerate with the mode n in B. The summation over all degenerate modes is motivated by their comparable contribution to the perturbed mode, as known from degenerate perturbation theory. We add modes of B to S + in decreasing W n order. Groups of degenerate modes D are added in one step as they have equal W n . A special case are the LE modes in the basis of the dielectric sphere, which are all degenerate having k n = 0. They are added in groups of equal l in the order of reducing weight.
To exemplify the local perturbation method, we use the quarter sphere perturbation with two different perturbations strengths ∆ǫ = 1 and ∆ǫ = 0.2, and choose the degenerate l = 7 modes shown in Fig. 4(b) as S. The perturbed RSs deriving from S are shown in Fig. 7(a) and (b), as calculated by RSE using either a global basis B with N = 16000, or a minimum local basis S + = S with N ′ ∼ 10, or a larger S + with N ′ ∼ 100. As in the previous section we show the results separately for each class of RSs (A and B) decoupled by symmetry. We find that for ∆ǫ = 0.2 (∆ǫ = 1) the perturbation lifts the degeneracy of S by a relative wavenumber change of about 1% (5%), and that the minimum local basis S + = S of only degenerate modes reproduces the wavenumbers with a relative error of about 10 −4 (10 −3 ), i.e. the perturbation effect is reproduced with an error of a few %. Increasing the local basis size to N ′ ∼ 100 the error reduces by a factor of three, by similar absolute amounts in the real and the imaginary part of the wavenumber [see insets of Fig. 7(a) and (b)].
The relative error of the local-basis RSE is generally decreasing with increasing basis size, as shown in Fig. 7(c). It can however be non-monotonous on the scale of individual sets of degenerate modes. This is clearly seen for for ∆ǫ = 0.2 and small N ′ , where adding the second group increases the error, which is reverted when the third group is added. These groups are the l = 6 and l = 8 fundamental WGMs as expected from Fig. 6(b), which are on opposite sides of S (l = 7 WGMs) in the complex frequency plane. Adding only one of them therefore imbalances the result, leading to an increase of the relative error. Comparing results in Fig. 7(c) for two different values of ∆ǫ, we see that the second-order correction dominates the relative error, as in the wide range of N ′ the error scales approximately like a square of the perturbation strength. The global-basis RSE, also shown in Fig. 7(b), has for a given basis size significantly larger errors. Fur-thermore, a minimum basis size is required for the basis to actually contain S, in the present case N ≈ 500. The local basis thus provides a method to calculate the perturbation of arbitrary modes with a small basis size.
The local perturbation method described in this section enables the calculation of high frequency perturbed modes which have previously been numerically inaccessible to FDTD and FEM due to the necessity of the corresponding high number of elements needed to resolve the short wavelengths involved and inaccessible to the RSE with a global basis due to the prohibitively large N required. The example we used for the illustration shows that a basis of ∼ 100 RSs in the local RSE can be sufficient to achieve the same accuracy as provided by FDTD and FEM in a reasonable computational time [see Figs. 5 and 7(c)]. For this basis size, solving the RSE Eq. (13) is 6 orders of magnitude faster than FDTD and FEM, and the computational time in our numerical implementation is dominated by the matrix element calculation which can be further optimized. A detailed evaluation of the performance of the local basis RSE and a comparison of selection criteria different from Eq. (39) will be given in a forthcoming work.

V. SUMMARY
We have applied the resonant state expansion (RSE) to general three-dimensional (3D) open optical systems. This required including in the basis both types of transversal polarization states, TE and TM modes, as well as longitudinal electric field modes at zero frequency. Furthermore, a general proof of the mode normalization used in the RSE is given. Using the analytically known basis of resonant states (RSs) of a dielectric spherea complete set of eigenmodes satisfying outgoing wave boundary conditions -we have applied the RSE to perturbations of full-, half-and quarter-sphere shapes. The latter does not have any rotational or translational symmetry and is thus not reducible to lower dimensions, so that their treatment demonstrates the applicability of the RSE to general 3D perturbations.
We have compared the performance of the RSE with commercially available solvers, using both the finite element method (FEM) and finite difference in time domain (FDTD), and showed that for the geometries considered here, the RSE is several orders of magnitude more computationally efficient, showing its potential to supersede presently used computational methods in electrodynamics. We have furthermore introduced a local perturbation method for the RSE, which is restricting the basis in order to treat a small subset of modes of interest. This further reduces computational efforts and improves on previous local perturbation methods.
We prove in this section that the spectral representation Eq. (4) leads to the RS normalization condition Eq. (5) and further to Eq. (6). To do so, we consider an analytic continuation E(k, r) of the wave function E n (r) around the point k = k n in the complex k-plane (k n is the wavenumber of the given RS). We choose the analytic continuation such that it satisfies the outgoing wave boundary condition and Maxwell's wave equation n )σ(r) (A1) with an arbitrary source term corresponding to the current density j(r) = σ(r)ic(k 2 − k 2 n )/(4πk). The source σ(r) has to be zero outside the volume V of the inhomogeneity ofε(r) for the electric field E(k, r) to satisfy the outgoing wave boundary condition. It also has to be non-zero somewhere inside V , as otherwise E(k, r) would be identical to E n (r). We further require that σ(r) is normalized according to V E n (r) · σ(r) dr = 1 + δ kn,0 , with the Kronecker delta δ kn,0 = 1 for k n = 0 and δ kn,0 = 0 for k n = 0. This ensures that the analytic continuation reproduces E n (r) in the limit k → k n . Indeed, solving Eq. (A1) with the help of the GF and using the GF spectral representation Eq. (10), we find: and using Eq. (A2) obtain lim k→kn E(k, r) = E n (r) .
We now consider the integral (A4) and evaluate it by using Maxwell's wave Eqs. (1) and (A1) for E n and E, respectively, and the source term normalization Eq. (A2): On the other hand, rearranging the integrand in Eq. (A4) and using the divergence theorem, we obtain with S V being the the boundary of V . Here, we used that for two arbitrary vector fields, a(r) and b(r), we can write The divergence theorem therefore allows us to convert all volume integrals in Eq. (A4) into surface integrals over the closed surface S V , the boundary of V , taken with an infinitesimal extension to the outside area whereε(r) is homogeneous, so that both ∇ · E and ∇ · E n vanish on that surface leaving only the integral shown in Eq. (A6). Finally, using Eq. (A5) in Eq. (A6) and taking the limit k → k n we obtain the normalization condition Eq. (5). The limit in Eq. (5) can be taken explicitly for any spherical surface [26]. In fact, outside the system, wherê ε(r) =1 (or a constant) the wave function of any k n = 0 mode is given by E n (r) = F n (k n r), where F n (q) is a vector function satisfying the equation and the proper boundary conditions at system interfaces and at q → ∞. The analytic continuation of E n (r) can be therefore be taken in the form E(k, r) = F n (kr) .
We use a Taylor expansion at k = k n to obtain E(k, r) ≈ F n (k n r) + (k − k n )r ∂F n (kr) ∂(kr) k=kn = E n (r) + k − k n k n r ∂E n (r) ∂r (A9) and ∂E(k, r) ∂r ≈ ∂E n (r) ∂r + k − k n k n ∂ ∂r r ∂E n (r) ∂r , where r = |r| is the radius in the spherical coordinates.
Choosing the origin to coincide with the center of the sphere of integration S V = S R we note that ∂/∂s = ∂/∂r in Eq. (5). Substituting Eqs. (A9) and (A10) into Eq. (5) and taking the limit k → k n obtain Eq. (6).
the form of a piece of a homogeneous spherical shell layer. The latter is suitable for treating an arbitrary symmetric or asymmetric perturbation of the sphere and is used in particular for half-and quarter-sphere perturbations considered in Sec. IV B and IV C, respectively.

Homogeneous sphere perturbation
The homogeneous perturbation Eq. (34) does not mix different m or l values, nor does it mix TE modes with TM or LE modes. Using the definition Eq. (14) we calculate the matrix elements between TE RSs performing the angular integration which leads to the lm-orthogonality: The radial integration can also be done analytically, so that the matrix elements take the form for identical basis states n = n ′ and for different basis states n = n ′ , where x = n R k n R and y = n R k n ′ R . Similarly, for TM RSs we find × R 0 l(l + 1)R l (r, k n )R l (r, k n ′ ) + ∂[rR l (r, k n )] ∂r ∂[rR l (r, k n ′ )] ∂r dr , and after analytic integration we obtain j l (x) (B3) for identical basis states n = n ′ and for different basis states n = n ′ , where with x = n R k n R and y = n R k n ′ R . Note that LE and TM modes are mixed by the perturbation, and nonvanishing matrix elements between them are calculated using Eqs. (B3) and (B4), treating the LE modes as TM modes with k n = 0 and the normalization constants multiplied by l(n 2 R − 1), in agreement with Eq. (33).
Factorizing the radial and angular integrals and using the fact that χ ′ m (ϕ) = mχ −m (ϕ), the matrix elements of the perturbation Eq. (B6) become between TE modes,