1. Introduction

One of the most useful optical geometries for generating and controlling electromagnetic radiation is the optical resonator. In general it is a photonic structures that allows a portion of the field to circulate in a closed path within the device. The resonance frequencies, resulting in a constructive interference of the field amplitude, depend on features such as size, materials and geometry. The resonator is a key component in lasers, waveguides, sensors and filtering devices [1–3]. The optical resonator can be configured into a number of equivalent geometries, from two plane mirrors in air [4], macro- and micro-optic ring resonators [1, 2], disk resonators [5], cylindrical resonators [6, 7], photonic crystal cavity resonators [8–10] and the focus of this presentation the spherical resonator [11].

The spherical resonators provides an environment for controlling and confining light in three dimensions and forms the backdrop for a number of more advanced resonator applications; the change of spontaneous emission for atoms placed in spherical dielectric cavities [9]; the second harmonic generation [12]; whispering-gallery mode based structures [13–15]; comb resonators [16]; plasmonics [17]. In addition the spherical shape can be used as an approximate geometry for a number naturally occurring objects and biological specimens; liquid droplet [18, 19], air bubble [20], viruses [21], biological cells [22] and bacteria [23], to which the resonator analysis presented here can be applied.

The computation techniques used to calculate the resonator properties in the spherical structures are highly dependent on the physical size compared to the wavelength of the electromagnetic fields involved. For large lossless objects, much larger than the wavelength of light considered, a ray optics approach can lead to important insight on the modal properties and is usually a simple computation to perform [15]. When the sphere’s size is comparable to the wavelength, or the structure contains internal features comparable to the wavelength, the full vector wave analysis is usually required for accurate results. For spherical structures displaying variations with respect to the radial coordinate only, exact solutions to Maxwell’s wave equation results in a characteristic equation from which the roots provide access to the resonator wavelengths and field profiles [24, 25]. For single and multiple spherical layer structures, the scattering matrix method has been used to obtain the resonator properties [26, 27]. When the spherical resonator includes material variations with respect to the angular coordinate directions, numerical and approximate methods are frequently considered for determining optical resonator properties. The most commonly encountered are the plane-wave-method [28, 29], the finite-difference-time-domain method [30, 31] and the finite-element-method [32]. The general solvers are somewhat inefficient as they do not exploit the spherical space properties of the optical structure being considered. Recently, an efficient numerical mode solver for cylindrically symmetric resonators configurations has been developed from Maxwell’s wave equation cast into an eigenvalue problem [33, 34]. That technique is reconfigured here such that it is applicable to spherical objects of arbitrary internal structure. The relative dielectric profile and the electromagnetic fields are represented by convergent series indexed over a set of basis functions linked to the spherical coordinate system. When introduced into Maxwell’s vector wave equation a matrix is generated which yields resonator state frequencies (wavelengths) and field profiles. The approach is general and as such can be used for any dielectric structure uniform or containing internal features as well as lossy material properties. This paper is organized as follows. The details related to the development of the numerical solver are introduced in sections 2. Section 3 presents applications of the technique to different spherical and commonly encountered structures and results are compared to published values. In section 4, the technique is applied to determining the resonator state properties of a lossless and lossy Bragg shell configurations. Section 5 closes the presentation. The analyses of structures containing dielectric variations in the polar and azimuthal angles are left for another presentation. These can be analysed using the technique presented here without any additional modifications to the expressions provided.

2. Development and application procedure

The numerical technique is developed from the vector wave equation for either the electric or magnetic field in a charge free, current free, isotropic, nonmagnetic and linear medium. Throughout the computation process, the time dependence is taken of the form $e^{-i\omega t}$ with $\omega ={\omega }_r+j{\omega }_i$ the complex angular frequency. The general starting equations, after performing time derivatives, are given in (1) where the different dielectric regions can be represented by the spatial variations in the relative dielectric constant, ${\epsilon }_r$, and product of the free space constitutive parameters replaced by the speed of light in vacuum squared, $c^2$.

The operational forms of the two equations in (1) differ only in the placement of the relative dielectric constant, which is in general a function of the coordinates. The operator form of the E field equation can be obtained from that of the H field when derivatives of the dielectric constant are set to zero. For this reason it is sufficient to provide the mathematical steps for the H equation, appendix A. The E field expressions are a subset of the H expressions as outlined in the appendix.

The spherical symmetric space of the computational domain is employed to define the form of the basis functions for expressing the field and inverse dielectric profiles. The radial coordinate variation is represented using the lowest order spherical Bessel function. Restricting to the lowest order is a departure from the standard approach of employing all orders of the spherical Bessel functions. In addition, the lowest order Bessel function is employed to represent the decaying field regions as well. The benefit to using only the lowest order spherical Bessel function are four fold in developing the matrix. The first is that any radial function can be series expanded in lowest order spherical Bessel functions, thus any higher order spherical Bessel (or decaying function) can be accommodated without loss of accuracy. The second is that the radial coordinate derivatives will only involve the zero and first order spherical Bessel and thus only two variants of the spherical Bessel are needed in populating the matrix elements. The third involves the computational domain integral of the product of three spherical Bessel functions (one form the field, one from the dielectric and one from the orthogonal field). Using only the lowest order spherical Bessel limits the number of different integrands to a manageable set that can be tabulated in a normalized form and used repeatedly for all structures examined. The forth advantage is that radial boundary conditions between dielectric regions are not required in order to generate characteristic equations for determining steady state frequencies. The polar angle,$\phi$ , dependence is accommodated using a single index Legendre polynomials. The index of the Legendre polynomial is unrestricted in positive values and a series expansion on these functions can be used to represent any polar angle dependent function [35, 36] in spherical space. In the matrix building process the integral of the product of three Legendre polynomials is present. Using a single unrestricted index limits the number of integrand to a manageable set that can be tabulated for use in the analysis of all dielectric structures referenced to a spherical coordinate system. The single index polynomial differs from using two indices and tying one of the indices to the azimuthal order. Using the single index Legendre polynomial greatly simplifies the matrix building process. For the azimuthal coordinate $\phi$, defined in the (x, y) plane, the 2π periodicity is available and thus its contribution can be represented through complex exponential in $\phi$ weighted by integer multiples. Due to the functional form of the basis functions used, we identify the technique as the BLF for Bessel-Legendre-Fourier.

The inverse of the relative dielectric profile and the field components series expansions, related to the basis functions are:

The matrix is chosen to have the following form. The lower subscript of the elements in the 3 by 3 matrix represents the component of the right hand side vector of (1). It also serves to indicate the orthogonal field component that was applied prior to computation domain integration. The letters $\left(R,\text{Θ},\text{Φ}\right)$ represents the vector component on the left hand side of (1). The field related expansion coefficients are listed in the column matrix. In building the matrix it is important to ensure that matrix element positions align with the column vector expansion coefficients. In addition each matrix element is produced through expression involving the summation over the expansion coefficients of the dielectric series. See appendix A for element generating expressions. The eigenvalues are the set, ${\left(\omega /c\right)}^2$, of values returned using a numerical eigenvalue solver such as the eig() function in MATLAB^©. The eigenvectors associated to each eigenvalue are also conveniently obtained using the eig() function. For any dielectric profile which fits the spherical symmetric representation and is examined using the process presented here produces four different categories of modal field characteristics. Localized states normally have high field strengths confined to the central coordinate region and negligible to zero field values at the computational radial edge. These states usually show a strong decaying field towards the edge. Edge states on the other hand show strong field components confined mostly at the edge of the computational domain. These are highly dependent on the computation edge properties and are not normally states examined using this technique. Changing the location of the computational edge will affect the edge state eigen-frequencies and field profiles, while the localized states remain unchanged. Super-states are states with significant field values throughout the entire structure and strongly dependant on the computational domain edge. These states are normally not the required states and can be neglected. Interface states are states that show strong field confinement at the interface between dielectric regions. When the fields for these states decay gradually to zero at the computation radial edge, they may represent acceptable mode profiles. Such states can be returned when searching for whispering gallery modes, bottle, SNAP and ring resonator states for example. In general the mode solver returns one state for each basis function used to express each of the field component in the series in (3). Using 100 different basis functions in the series will return 300 eigenvalues and eigenvectors. The properties of the states can be filtered by examining the dominant field component in the eigenvector, the value of the corresponding wavelength compared to the structure size parameters and the dominant basis function expansion coefficient. The real or complex nature of the eigen-frequency and the relative magnitude of the real and imaginary parts can serve to filter returned states as well. A true measure of the modes properties is determined by reconstructing the field using (3) and plotting the field in different coordinate planes.

2.1 Azimuthal mode order tuning

The orthogonality integration related to the in-plane angle $\phi$ plays an important part in reducing the matrix order and serves to tune the computation process to search for a particular in-plane mode type. The integral is written in (5) and places a restriction on the azimuthal indices combinations for non-zero values.

3. High-dielectric spheres

A single high dielectric constant sphere placed in air is a simple structure to test the accuracy of the Bessel-Legendre-Fourier mode solving technique as published and well established results are available for comparison [24, 25]. The dielectric sphere is selected to have a radius $a$ and relative dielectric constant${\epsilon }_r=9$. For computational purposes, the air region is extended to a radius $R=2a$ with the dielectric sphere placed at the centre. The geometry is shown in Fig. 1(a). Since the BLF technique is configured to solve for steady states, the air region extending to 2a is sufficient to isolate the sphere localized states from the computational edge influence. States of the dielectric sphere have effectively a zero field component on the computation edge when well localized. A similar reasoning is also employed when searching for localized states in single defect containing photonic crystals using the plane wave expansion technique [29]. For the geometry examined here, the relative dielectric profiles has only radial variations, abrupt change at $r=a$, and as such the series of the inverse relative dielectric profile can be represented using a single index, $p$, setting $q$ and $n$ to zero. The series expansion coefficients obtained from (2) using an orthogonal integration process are plotted in Fig. 1(b) versus basis function index$p$. The expansion coefficients obtained contain only real parts as the mediums involved are lossless and $\phi$ coordinate independent. The rapid zero approaching coefficient with increase basis index assures that the dielectric decomposition is well converged. The process of reconstructing the inverse dielectric using Eq. (2) and comparing to the original structure was used to confirm the accuracy of the series representation for dielectric of Fig. 1. This reconstruction test was also used to ensure proper series representation for all dielectric profiles presented later.

 

Fig. 1 (a) – Geometry of the solid-uniform, ${\epsilon }_r=9$, dielectric sphere, radius a, in the spherical computation domain of radius R = 2a. (b) – Real part of the dielectric series expansion coefficients plotted versus spherical Bessel basis function index. Imaginary parts all zero.

Download Full Size | PDF

To populate the matrix (4), the field components of the series in (3) were expressed using 70 spherical Bessel terms (p = 1 to 70) and 71 Legendre terms (n = 0 to 70). Legendre expansion terms are required for the field even though the dielectric is uniform in order to resolve the $\theta$ related modal profile properties. The q index was set to zero which tunes the solver for localized states that are monopole like in the (x, y) plane. The number of basis functions ensures sufficiently converged localized states for the dielectric sphere in air. Convergence test were performed by increasing the number of terms in the field component series and re-computing the fields and state frequencies. To ensure computational accuracy of results presented latter, similar convergence tests were performed. When monopole states are considered in a dielectric medium showing no $\phi$ dependence, the field components decouple into two groups; TE with components $H_r, H_{\theta },E_{\phi }$; and TM with components $E_r, E_{\theta },H_{\phi }$. The decoupling of fields can be confirmed by introducing $q,q^{*}$ and $q_{\text{Ω}}$ equal to zero in the matrix element generating expressions (9) to (17). The decoupling of fields for monopoles in the dielectric sphere environment is similar to the decoupling of the fields observed when the out of plane direction has no dielectric variation as observed in the plane wave analysis of photonic crystals [37]. Population of the matrix elements, for the E field equation, using the generating expressions in Appendix A, results in 4 of the 9 element blocks being all zero, expression (6) below, and enables the $E_{\phi }$ matrix block to be solved independently for its eigenvalues and eigenvectors. A reduction in the matrix order by a factor of 3, in addition to azimuthal mode order tuning, greatly speeds up the computation process when dealing with large matrices.

The eigenvalues returned from (6) and converted to normalized resonance wavenumber, $ka$, for several lower order monopoles are listed in Table 1 for the $E_{\phi }$ (TE) field component. Computed values using BLF are tabulated along with corresponding values extracted from [24, 25] and indicate an excellent agreement between the two different techniques. The bold values correspond to localized states that are plotted in Figs. 2 and 3. Mode nomenclature refers to the number of field maximums observed along the coordinate unit vector directions. Field profiles also match those plotted in [24, 25].

Table 1. Scaled resonance wavenumbers determined using BLF and from [24, 25] for the $E_{\phi }$ , TE polarization. Bold entries have field profiles plotted in Figs. 2 and 3.

View Table | View all tables in this article

 

Fig. 2 The (x, y) plane normalized intensity profile $E_{\phi }$ field component (Left), and the (y, z) plane normalized intensity profile $E_{\phi }$ field component (Right), for the R₁Θ₁Ψ₀ localized state listed in Table 1, computed using BLF,$ka=1.0604$. The field profile and corresponding eigenvalue match those of [24, 25].

Download Full Size | PDF

 

Fig. 3 The (x, y) plane normalized intensity profile $E_{\phi }$ field component (Left), and the (y, z) plane normalized intensity profile $E_{\phi }$ field component (Right), for the R₁Θ₅Ψ₀ localized state listed in Table 1, computed using BLF, $ka=2.7040$. The field profile and corresponding eigenvalue match those of [24, 25].

Download Full Size | PDF

The upper left corner 2 by 2 matrix element blocks of the matrix (6) can be solved for the states containing $\left(E_r, E_{\theta }\right)$ field components, TM. The range of scaled wavenumbers corresponding to the lower set of localized states is listed in Table 2 along with those extracted from [24, 25]. The BLF computed values are computed to an accuracy of ± 0.0001 while those of the reference are cited to only two decimal points. The discrepancy between techniques is less than 5%, however, the cited paper used analytical expressions and many assumptions have been presented in their solutions. The modal profiles for the bold entries are plotted in Figs. 4 to 6 for the various field components and are in excellent agreement with the features of the profiles plotted in the [24, 25].

Table 2. Scaled resonance wavenumbers determined using BLF and from [24, 25] for the $E_r, E_{\theta }$ , TM polarization. Bold entries have field profiles plotted in Figs. 4 to 6 (TM)

View Table | View all tables in this article

 

Fig. 4 The (x, y) plane normalized intensity profile $E_{\theta }$ field component (Left), and the (y, z) plane normalized intensity profile $E_{\theta }$ field component (Right), for R₂Θ₁Ψ₀ localized state listed in Table 2, computed using BLF, $ka=2.4954$.

Download Full Size | PDF

 

Fig. 5 The (x, y) plane normalized intensity profile $E_{\theta }$ field component (Left), and the (y, z) plane normalized intensity profile $E_{\theta }$ field component (Right), for R₂Θ₅Ψ₀ localized state listed in Table 2, computed using BLF, $ka=4.2850$.

Download Full Size | PDF

 

Fig. 6 The (y, z) plane normalized intensity profile $E_r$ field component for R₂Θ₁Ψ₀ localized state listed in Table 2, computed using BLF, $ka=2.4954$. (Left). The (y, z) plane normalized intensity profile $E_r$ field component for R₂Θ₅Ψ₀ localized state listed in Table 2, computed using BLF, $ka=4.2850$ (Right). The (x, y) plane have zero $E_r$ field values.

Download Full Size | PDF

3.1 Whispering gallery modes

The high dielectric sphere in air is known to support high azimuthal order states referred to as whispering-gallery modes. These states generally present a field profile highly confined at the interface between the high dielectric sphere and surrounding air region and a large number of field amplitude cycles within the azimuthal plane sphere’s circumference [13]. The BLF technique is suitable for calculating the WGM for large $ka$ values as in [13], in this computation case WGM of low $ka$ were considered. WGM were computed by populating the matrix using $\left(q,q^{*}\right)=(20,20)$ for azimuthal order 20 and $\left(q,q^{*}\right)=(40,40)$ for azimuthal mode order 40. The field components were represented using 60 for the spherical Bessel and 61 for the Legendre polynomial index maximums. Since the dielectric expansion has non-zero expansion coefficients only for $q_{\text{Ω}}=0$, there is no mixing of the azimuthal mode orders for this dielectric structure and each WGM mode order can be separately solved using a reduced order matrix. The entire populated matrix for each mode must be solved as all 9 generating expressions yield non-zero matrix elements in (4). In intensity profile a WGM of order Ψ₂₀ and $ka=8.3280$, shown in Fig. 7(a), displays 40 field extremes (20 field amplitude cycles) when observed in the (x, y) plane and strong confinement to the (x, y) plane when viewed in the (x, z) plane. In Fig. 7(b) intensity profile a WGM of order Ψ₄₀ and $ka=15.4701$, is shown, similar properties of the mode profile confinement are presented and displays 80 field extremes (40 field amplitude cycles) when observed in the (x, y) plane. The field is highly confined to the high dielectric side of the sphere-air interface with decaying field present on the air side.

 

Fig. 7 Intensity profile for the $E_{\theta }$ field component of the FBL computed whispering-gallery-modes of the solid sphere in air. (a) Top pair – Azimuthal mode order 20 with field R₁Θ₁Ψ₂₀ and $ka=8.3280$. (b) Bottom pair – Azimuthal mode order 40 with field R₁Θ₁Ψ₄₀ and $ka=15.4701$. Left – (x, y) plane. Right – (y, z) plane.

Download Full Size | PDF

3.2 Spherical cavity

A spherical cavity structure can be configured from the dielectric sphere in air by setting the computational domain edge at $R=a$ instead of $2a$. This forces the surface of the dielectric sphere to act as a perfect boundary with zero field values for all components at $r=a$. This is an artificial boundary condition as it is imposed by the spherical Bessel contributions to the basis functions all being zero at the computational edge. However the boundary condition does hold for the$E_{\phi }$, $E_{\theta }$, $H_{\phi }$ and $H_{\phi }$ components at a perfect conductor or perfect magnetic boundary placed on the sphere. The spherical cavity structure presents a single dielectric value within the computational domain and can be represented using only the spherical Bessel part of the basis functions. The constant dielectric profile still requires a series expansion using (2) as the expansion does not contain an averaging or constant offset term. The inverse dielectric was resolved using 100 spherical Bessel basis terms ($p_{\text{Ω}}$ from 1 to 100) with $n_{\text{Ω}}$ and $q_{\text{Ω}}$ set to 0.

Since the boundary conditions are satisfied for the $E_{\phi }$ field component, the matrix was populated using expressions related to the electric field equation. For monopole states in the azimuthal plane, the $E_{\phi }$ matrix element block is decoupled from the other two field component blocks and can be solved independently. The scaled resonance wavenumber, computed using the first 50 spherical Bessel terms and first 51 Legendre polynomials to populate the matrix are presented in Table 3 along with corresponding values extracted using other techniques [38, 39]. Reference values are quoted accurate to 4 decimal places since they can be determined directly from the wave equation without the introduction of any approximations. The FFB quoted results are also presented for well converged numerical values and show excellent agreement with published values. The difference between the two sets of values is less than 0.03%. To further demonstrate the applicability of the FFB technique in determining the states of the cavity, the field profiles for the Table 3 entries in bold are plotted in Figs. 8 and 9 in the two orthogonal planes. It is easy to verify that these field profiles correspond to the field profiles of the [38, 39].

Table 3. Comparison of the scaled resonance wavenumber obtained using BLF and from [38, 39]. Bold-field profiles plotted in Figs. 8 and 9 computed using the BLF technique with eigenvector dominated by the $E_{\phi }$ field component.

View Table | View all tables in this article

 

Fig. 8 The normalized intensity profile $E_{\phi }$ field component for R₁Θ₁Ψ₀ with, $ka=4.4935$. Left – (x, y) plane. Right – (y, z) plane.

Download Full Size | PDF

 

Fig. 9 The normalized intensity profile for $E_{\phi }$ field component for R₂Θ₃Ψ₀ for $ka=10.4181$. Left – (x, y) plane. Right – (y, z) plane.

Download Full Size | PDF

4. Spherical “Bragg” resonators

When solid spheres are coated by alternating layers of thin dielectric materials, several of the resonance states are caused by “total / partial internal reflection” at the shell layer region. The ability to control the resonance properties in the three dimensional domain makes the spherical Bragg shell resonator a structure of considerable interest and suitable for numerous applications. Adding one layer of high refractive index dielectric on a dielectric microsphere results in enhancement of the sensitivity for a WGMs based sensor [40]. Such geometry is also used for higher efficiency back reflector in dye sensitized solar cells [41]. Optical reflect array nano-antenna made of dielectric core and a shell made of a plasmonic material have achieved narrow beamwidth [17]. For a sphere coated by Bragg shells light can be confined strongly, for a range of wavelengths. This is important for controlling spontaneous emission for atoms placed inside the Bragg shells [9]. In Fig. 10, a microsphere placed in the center and coated with Bragg spherical shells of 24 layers is shown. The structure is placed in air and details of the dielectric values and design of the structure are taken from [27]. The radius of the microsphere is 1 µm with relative dielectric constant of 2.25. The Bragg layers are formed from alternating regions with 0.301 µm (${\epsilon }_r=2.10$) and 0.213 µm (${\epsilon }_r=4.20$). The computation domain is extended 5 µm beyond the external layer resulting in $R=12.47$ µm. In [27], these structures were studied using scatting matrix method [26]. In this section the result will be calculated using BLF approach and compared to the cited paper. The calculations demonstrate the validity and accuracy of the BLF technique for these structures. Results are quoted in the frequency domain to match the domain of the cited results.

 

Fig. 10 A dielectric sphere coated by spherical shells dielectric. The shells structure acts as quarter wave Bragg reflector, the relative dielectric constant has real values only. First 100 radial dielectric expansion coefficients for the spherical shell structure.

Download Full Size | PDF

Due to the small width of each dielectric region, the inverse dielectric series were extended to include the first 400 spherical Bessel basis terms ($p_{\text{Ω}}$ from 1 to 400). The rotational symmetry present permits setting $n_{\text{Ω}}$ and $q_{\text{Ω}}$ to 0. In Fig. 10, the first 100 expansion coefficients are plotted indicating sufficient convergence of the radial contribution to the series expansion.

In this computation, the field components were represented using 60 the spherical Bessels and 61 Legendre polynomials. The BLF computed state for $H_{\phi }$ field component at a frequency $f=188$ THz consisting of R₁Θ₁Ψ₀, as shown in Fig. 11. The modal profile and computed frequency match the cited paper and only a 2% difference in frequency values is observed. The localized mode is highly confined around the center, with zero value at the edge of the computation domain. Results of different order and field component were are also obtained and examined.

 

Fig. 11 The normalized intensity profile $H_{\phi }$ field component for R₁Θ₁Ψ₀, computed using BLF. The field profile and corresponding eigenvalue match those of [27]. Left – (x, y) plane. Right – (y, z) plane.

Download Full Size | PDF

In Figs. 12 and 13 mode profiles for $E_{\theta }$ field component for R₁Θ₁Ψ₃₀ at $f=134.12$ THz ($ka=2.8089$) and R₁θ₃Ψ₃₀, at $f=142.09$ THz ($ka=2.9760$) are shown, respectively. The field components were represented using 60 spherical Bessels and 61 Legendre polynomials. The fields are confined at the outer layers of the spherical shells.

 

Fig. 12 The normalized intensity profile $E_{\theta }$ field component for R₁Θ₁Ψ₃₀, computed using BLF,$f=134.12$ THz ($ka=2.8089$). Left – (x, y) plane. Right – (y, z) plane.

Download Full Size | PDF

 

Fig. 13 The normalized intensity profile $E_{\theta }$ field component for R₁Θ₃Ψ₃₀, computed using BLF, $f=142.09$ THz ($ka=2.9760$). Left – (x, y) plane. Right – (y, z) plane.

Download Full Size | PDF

In some instance metal layers are present within the dielectric layered structure [42–44]. The solver presented here is suitable for non-metallic medium. The inclusion of metal properties in the numerical solver is currently being developed for plasmonic geometries of a spherical symmetry.

4.1 Aperiodic reflector resonators

The geometry of the aperiodic reflector is similar to the Bragg layer structure but designed to have a low dielectric fill factor. The high dielectric regions close to the center are thin and gradually increase in width as their radius increases. The outermost layer is designed to be ¼ wavelength at some design value. The parameters for the structure in Fig. 14 are extracted from [45] and consisted of alternating spherical layers of air and nested alumina. The nested alumina used in this structure has real and imaginary parts to the dielectric value. The real part is ${\epsilon }_{r, real}=9.8$, and the imaginary part is ${\epsilon }_{r, imaginary}=0.000098$. The imaginary part is used to represent the absorption loss of the material with loss tangential of 10⁻⁵. In Fig. 14, the real and imaginary expansion space for the dielectric decomposition is shown. Large number (500 spherical Bessel basis terms) is used to ensure accurate representation of the dielectric structure.

 

Fig. 14 A sphere of air coated by aperiodic spherical shells, the dielectric has real and imaginary values. Real radial dielectric expansion coefficients for the aperiodic shell structure (Left). Imaginary radial dielectric expansion coefficients for the aperiodic shell structure (Right).

Download Full Size | PDF

In Fig. 15, the mode profile for the $E_{\phi }$ field components for R₁Θ₁Ψ₀ is plotted using 60 and 61 for the spherical Bessel and Legendre polynomial indices. The wavenumber with real and imaginary values is calculated as $k=2.11\times {10}^2-j 2.13\times {10}^{-5} m^{-1}$. The localized mode profile match the results obtained in the cited paper, the field is decaying outside the center sphere due to the aperiodic Bragg shells. The frequency result also matches the referenced result. The resonance frequency has real and imaginary parts, the real part represent confined energy and the imaginary part represents energy loss due to the absorption.

 

Fig. 15 The normalized intensity profile $E_{\phi }$ field component for R₁Θ₁Ψ₀, computed using BLF. The field profile and corresponding eigenvalue match those of [42]. $f_{real}=10.06$ GHz, $f_{imaginary}=1.0\times {10}^3$ Hz. Left – (x, y) plane. Right – (y, z) plane.

Download Full Size | PDF

In [45] the Q factor was determined using analytical calculations. Using the electric filling factor and multiplying it by the loss tangent of the thin dielectric medium the material loss are computed. For the referenced paper the Q factor is given $4.0\times {10}^6$. Using the BLF method the losses of the material and the Q factor are calculated based on the real and imaginary frequency values:

5. Conclusion

A Bessel-Legendre-Fourier basis (BLF) set is used to cast Maxwell's wave equation into an eigenvalue problem from which the stationary states of dielectric structure that possess spherical symmetry can be determined. The technique is highly efficiency in determining the frequencies, wavelengths, field profiles and Q values for lossy structures that have an underlying spherical geometry. The numerical approach is presented here for the first time and the generating expressions are provided. Under certain conditions related to structure geometry and desired mode properties, the solver can be tuned yielding a reduced order matrix to solve. The application of the technique have been restricted to structures only showing radial variations as published resonator properties for these structures are readily available to compare the accuracy of the BLF results. The technique presented here can be readily applied to other geometries fitting in the spherical boundary containing dielectric variations in any combination of the spherical coordinate directions.