Unified theory of whispering gallery multilayer microspheres with single dipole or active layer sources

Jonathan M. M. Hall; Tess Reynolds; Matthew R. Henderson; Nicolas Riesen; Tanya M. Monro; Shahraam Afshar

doi:10.1364/OE.25.006192

1. Introduction

Due to their guidance of optical whispering gallery modes (WGMs), dielectric microspheres have attracted a lot of interest in different fields of research, such as remote atmospheric sensing [1, 2], biosensing technologies [3–5], photonic band-gap devices [6, 7], fluorescence spectroscopy [8–10], nonlinear optics [11–16], superscattering [17], and metamaterial perfect absorbers [18,19]. The extensive literature on these applications, together with the numerous theoretical works on light scattering from microspheres [1,2,9,20–26] based on Mie scattering [27], clearly demonstrate a continually growing interest in this field of research. Technological advances have led to new possibilities for using microspheres with several dielectric layers, coatings of active materials, or doped with fluorescence nanoparticles [28,29]. The excitation of WGM resonances of microspheres can also occur in numerous ways, such as via phase-matched optical waveguides [6,30] or fiber tapers [31–34], prism coupling [35], fluorescence emission of incorporated nanoparticles [36–38], or active material coatings [39]. As a result, there is now a necessity to have more advanced models, as well as efficient numerical simulation tools, for describing these resonators and their excitation schemes accurately.

Simulation techniques are required for the interpretation of measured whispering gallery mode (WGM) spectra obtained from optical resonators. They may be used to identify the polarizations of the modes in a given WGM spectrum [40], and provide insight into how they can be used for sensing applications [25,41–43]. Furthermore, the ability to calculate the underlying geometric parameters of a given resonator based solely on its spectrum [44] by scanning over a wide parameter space makes the development of a fast, reliable and general model for resonators of high importance.

The derivation of a model that contains multiple layers of concentric spheres has, until now, not been treated comprehensively. The development of a WGM spectrum model by Chew [20,45] considered spherical dielectric particles with embedded dipole sources, with the motivation being the modeling of Raman and fluorescence scattering. The Chew model was then extended to include a uniform distribution of dipole sources placed on the surface of a sphere [9]. A multilayer variant to the Chew model exists, but contains no derivation of the power spectrum [45]. Meanwhile, a generalization of Mie scattering theory, developed for spherical concentric ‘onion’ resonators, is constrained to external ray excitation [21–23], and does not provide emitted power spectra with which to compare to WGM resonances.

Analytic models typically make use of the transfer-matrix approach [2], which is faster and more convenient to construct mathematically than the multilayer Chew model. However, this approach suffers from numerical instabilities for certain parameter values [46], and needs to be treated carefully. More recent work on single-layer coated microspheres clearly separates the spectrum into TE and TM modes, and calculates the resonance positions [24,26]. Note that, in the case of [2], the labeling of the TE and TM modes has been interchanged. The key result of the present work is the derivation of a model for multilayer microsphere resonators, which may include dipole sources in any layer, or an active layer of sources. The model is able to generate the wavelength positions of the TE and TM modes, and calculate the emitted power spectrum, from a formalism that is general for any excitation strategy, unlike previous works [2].

The format of the article is as follows. In Section 2, the mathematics of the problem is summarized, and the conventions used in the literature are harmonized. In solving the boundary value problem of a multilayer microsphere using the transfer-matrix approach, the resonance positions are calculated. The formulae for the normalized power as a function of wavelength are then derived for two main cases. The first case is mode excitation via a single dipole placed on the surface with any desired orientation, analogous to embedded nanoparticles. The second case is a uniform distribution of dipoles of random orientation, extending the work by Chew [9], which is analogous to a fluorescent dye coating as commonly used in biosensing applications [47, 48]. The distribution of dipoles can occur in any layer of the multilayer resonator. In Section 3, the formulae are shown to apply to specific cases, highlighting the differences between multilayer, single-layer and uncoated microsphere spectra. We showcase the results for active layer coated microspheres, and demonstrate the novelty of the unified multilayer model. A discussion of the implementation of the transfer-matrix algorithm is also included, demonstrating improvements made to the mean execution time. Tests showing that the model converges to the specific cases as the layer thickness becomes vanishingly small are given in the Appendix A.

2. Theory

In this section, we state the basic equations required to find the wavelength positions of the resonances, and the radiation power spectrum of a multilayer microsphere, with either dipole sources or active layers for the excitation method.

2.1. Geometry

We consider a microsphere with an arbitrary number of concentric layers, N. The refractive index distribution and thicknesses are illustrated in Fig. 1. It is assumed that each layer includes a dipole emitter located at position r′_j, where j is the layer number and the prime symbol is specifically used for the position of the sources. We use the radial coordinate r when representing fields at an arbitrary point in space, and r_j to specify unambiguously the radius of the boundary between the j^th and j + 1^th region. The outermost region is N + 1, extended to infinity, and the innermost region is 1, hence the radius of the inner region is r₁ and r_N is the boundary of region N.

Fig. 1 The geometry of a spherical resonator with N layers. (a) A single layer contains a uniform distribution of dipoles, to represent an active layer. (b) One or more individual dipoles can be placed in a given layer, to represent one or more embedded nanoparticles.

Download Full Size | PDF

2.2. Conventions

Considering the spherical symmetry of the problem, we use Vector Spherical Harmonics (VSH). Different conventions exist for the definition of VSH in the literature, e.g. atomic physics and electrodynamics (see [27] and [49]). We modify the definition given by Barrera [50],

Y_{l m} = Y_{l m} (θ, ϕ) \hat{r}, Ψ_{l m} = (\frac{1}{i \sqrt{l (l + 1)}}) r \nabla Y_{l m} (θ, ϕ), Φ_{l m} (θ, ϕ) = (\frac{1}{i \sqrt{l (l + 1)}}) r \times \nabla Y_{l m} (θ, ϕ),

where Y_lm(θ, ϕ) are standard Scalar Spherical Harmonics [51], and l and m are the azimuthal quantum numbers. The vector functions Y_lm(θ, ϕ), Ψ_lm(θ, ϕ), and Φ_lm(θ, ϕ) are orthonormal, and form a complete set (See Eqs. (68) and (69)), i.e. any vector field in spherical coordinates can be expanded based on these functions,

E (r, θ, ϕ) = \sum_{l = 0}^{\infty} \sum_{m = - l}^{m = 1} [E_{l m}^{r} (r) Y_{l m} (θ, ϕ) + E_{l m}^{(1)} (r) ψ_{l m} (θ, ϕ) + E_{l m}^{(2)} (r) Φ_{l m} (θ, ϕ)],

where the coefficients

E_{l m}^{r}

,

E_{l m}^{(1)}

and

E_{l m}^{(2)}

can be found by using orthogonality relations (see Appendix A). Note that Y_lm is in the radial direction and Ψ_lm and Φ_lm are in the transverse plane perpendicular to

\hat{r}

, and hence

E_{l m}^{(1)}

and

E_{l m}^{(2)}

represent the coefficients of the transverse field.

2.3. Transfer matrix method

We solve Maxwell’s Equations (in Gaussian units),

E = \frac{i c}{ω ε} (\nabla \times H), H = - \frac{i c}{ω μ} (\nabla \times E),

taking into account the boundary conditions of the problem. In general, the total electric and magnetic field in each region j can be written as the sum of the fields due to the dipole in the layer, denoted E_jd, H_jd, and those associated with the reflection and transmission from other layers, denoted E_j, H_j. Hence

E_{j}^{total} = E_{j} + E_{j d}

and

H_{j}^{total} = H_{j} + H_{j d}

. The fields (E_j, H_j) can be written as [45]:

E_{j} = \sum_{l, m} (\frac{i c}{n_{j}^{2} ω}) A_{j} \nabla \times [j_{l} (k_{j} r) Φ_{l m} (θ, ϕ)] + (\frac{i c}{n_{j}^{2} ω}) B_{j} \nabla \times [h_{l}^{(1)} (k_{j} r) Φ_{l m} (θ, ϕ)] + C_{j} j_{l} (k_{j} r) Φ_{l m} (θ, ϕ) + D_{j} h_{l}^{(1)} (k_{j} r) Φ_{l m} (θ, ϕ),

H_{j} = \sum_{l, m} - (\frac{i c}{μ_{j} ω}) c_{j} \nabla \times [j_{l} (k_{j} r) Φ_{l m} (θ, ϕ)] - (\frac{i c}{μ_{j} ω}) D_{j} \nabla \times [h_{l}^{(1)} (k_{j} r) Φ_{l m} (θ, ϕ)] + (\frac{1}{μ_{j}}) A_{j} j_{l} (k_{j} r) Φ_{l m} (θ, ϕ) + (\frac{1}{μ_{j}}) B_{j} h_{l}^{(1)} (k_{j} r) Φ_{l m} (θ, ϕ) .

Here, A_j, B_j, C_j and D_j are coefficients that are determined through the boundary conditions. Note that A_j and B_j describe the transverse component of H_j and thus the TM modes, whereas C_j and D_j describe the TE modes, as defined by Jackson [51]. In this notation, the TM modes are defined as the modes for which H_j has no component in the radial direction, whereas the TE modes are those for which E_j has no such radial component. Also, the choice of Bessel functions j_l(k_jr) and

h_{l}^{(1)} (k_{j} r)

ensures that appropriate functions of r can be constructed for any layer, including the innermost and outermost regions. Using the properties of the orthonormal functions Y_lm, Ψ_lm, and Φ_lm (Eqs. (68) through (70)), we rewrite Eqs. (4) and (5) as

E_{j} = \sum_{l, m} [- (\frac{i c}{n_{j}^{2} ω}) \frac{\sqrt{l + (l + 1)}}{i} [A_{j} \frac{1}{r} j_{l} (k_{j} r) + B_{j} \frac{1}{r} h_{l}^{(1)} (k_{j} r)] Y_{l m} (θ, ϕ) - (\frac{i c}{n_{j}^{2} ω}) {A_{j} \frac{1}{r} \frac{d}{d r} [r j_{l} (k_{j} r)] + B_{j} \frac{1}{r} \frac{d}{d r} [r h_{l}^{(1)} (k_{j} r)]} ψ_{l m} (θ, ϕ) + {C_{j} j_{l} (k_{j} r) + D_{j} h_{l}^{(1)} (k_{j} r)} Φ_{l m} (θ, ϕ)],

H_{j} = \sum_{l, m} [(\frac{i c}{μ_{j} ω}) \frac{\sqrt{l + (l + 1)}}{i} [C_{j} \frac{1}{r} j_{l} (k_{j} r) + D_{j} \frac{1}{r} h_{l}^{(1)} (k_{j} r)] Y_{l m} (θ, ϕ) + (\frac{i c}{μ_{j} ω}) {C_{j} \frac{1}{r} \frac{d}{d r} [r j_{l} (k_{j} r)] + D_{j} \frac{1}{r} \frac{d}{d r} [r h_{l}^{(1)} (k_{j} r)]} Ψ_{l m} (θ, ϕ) + (\frac{1}{μ_{j}}) {A_{j} j_{l} (k_{j} r) + B_{j} h_{l}^{(1)} (k_{j} r)} Φ_{l m} (θ, ϕ)] .

This explicitly shows that the electric and magnetic fields (E_j, H_j in the layer j are written in terms of the orthonormal function set Y_lm(θ, ϕ), Ψ_lm, (θ, ϕ), and Φ_lm (θ, ϕ), in a form consistent with Eq. (2). The boundary conditions at the interfaces of the layers imply that transverse components of the fields are continuous, while there is a discontinuity in the normal components of the fields. We focus on transverse components of the fields, indicated by superscript (1) and (2) in Eq. (2), in layer j. We write them in a matrix form using Eqs. (6) and (7)

E H_{j} = M_{j} (r) A_{j},

where the block diagonal matrix

M_{j} (r) = (\begin{matrix} M_{j}^{T M} & 0 \\ 0 & M_{j}^{T E} \end{matrix})

and the vectors EH_j and A_j are defined as

M_{j} (r) = \frac{1}{k_{j} r} (\begin{matrix} - (\frac{i}{n_{j}}) ψ_{l}^{'} (k_{j} r) & - (\frac{i}{n_{j}}) χ_{l}^{'} (k_{j} r) & 0 & 0 \\ (\frac{1}{μ_{j}}) ψ l (k_{j} r) & (\frac{1}{μ_{j}}) χ_{l} (k_{j} r) & 0 & 0 \\ 0 & 0 & ψ_{l} (k_{j} r) & χ_{l} (k_{j} r) \\ 0 & 0 & (\frac{i n_{j}}{μ_{j}}) ψ_{l}^{'} (k_{j} r) & (\frac{i n_{j}}{μ_{j}}) χ_{l}^{'} (k_{j} r) \end{matrix}),

E H_{j} (r) = (\begin{matrix} E_{j}^{(1)} (r) \\ H_{j}^{(2)} (r) \\ E_{j}^{(2)} (r) \\ H_{j}^{(1)} (r) \end{matrix}), A_{j} = (\begin{matrix} A_{j} \\ B_{j} \\ C_{j} \\ D_{j} \end{matrix}) .

In arriving at Eq. (9) we have used c/n_jω = 1/k_j, Riccati-Bessel and Riccati-Hankel functions, ψ_l(k_jr) = k_jrj_l(k_jr) and $χ_{l} (k_{j} r) = k_{j} r h_{l}^{(1)} (k_{j} r),$ and their derivatives with respect to their arguments, $ψ_{l}^{'} (k_{j} r)$ and $χ_{l}^{'} (k_{j} r)$ . Note that the determinant of 2 × 2 blocks, $M_{j}^{T M}$ and $M_{j}^{T E}$ , are given by

\det (M_{j}^{T M} (r)) = \frac{i}{μ_{j} n_{j} k_{j}^{2} r^{2}} W_{k}_{_{j} r} [ψ_{l} (k_{j} r), χ_{l} (k_{j} r)] = - \frac{1}{μ_{j} n_{j} k_{j}^{2} r^{2}},

\det (M_{j}^{T E} (r)) = \frac{1}{k_{j}^{2} r^{2}} (\frac{i n_{j}}{μ_{j}}) W_{k} {_{_{j}}}_{r}_{_{j}} [ψ_{l} (k_{j} r), χ_{l} (k_{j} r)] = - \frac{n_{j}}{μ_{j} k_{j}^{2} r^{2}},

where the Wronskian is defined as

W_{x} [f (a x), g (a x)] \equiv f (a x) g' (a x) - f' (a x) g (a x), for the derivative with respect to x .

Using Eq. (19) in [2], it can be shown that

W_{k_{j}}_{r} [ψ_{l} (k_{j} r), χ_{l} (k_{j} r)] = i

. We write the fields E_jd(r), H_jd(r) due to a dipole located in the layer j and at position r′_j in a similar way as Eq. (8), and hence we modify the field equations in [20] to write

E H_{j d} (r) = θ (r_{j}^{'} - r) M_{j} (r) a_{j}_{L} + θ (r - r_{j}^{'}) M_{j} (r) a_{j H} .

Here, θ(r) is a step function [θ(r) = 0 for r < 0 and θ(r) = 1 for r ≥ 0] that ensures the correct Bessel function j_l(k_jr) or

h_{l}^{(1)} (k_{j} r)

is selected to evaluate the fields at point r either when

r < r_{j}^{'}

(shown by subscript L in a_jL) or

r > r_{j}^{'}

(shown by subscript H in a_jH), respectively. The vectors EH_jd, a_jL and a_jH are given by

E H_{j d} = (\begin{matrix} E_{j d}^{(1)} \\ H_{j d}^{(2)} \\ E_{j d}^{(2)} \\ H_{j d}^{(1)} \end{matrix}), a_{j L} = (\begin{matrix} a_{j E L} \\ 0 \\ a_{j M L} \\ 0 \end{matrix}), a_{j H} = (\begin{matrix} 0 \\ a_{j E H} \\ 0 \\ a_{j M H} \end{matrix}),

where the coefficients a_jEL, a_jML, a_jEH and a_jML take the form [20]:

a_{j E L} (r_{j}^{'}) = 4 π k_{j}^{2} \sqrt{\frac{μ_{j}}{ϵ_{j}}} P . \nabla_{j}^{'} \times [h_{l}^{(1)} (k_{j} r_{j}^{'}) Φ_{l m}^{*} (θ_{j}^{'}, ϕ_{j}^{'})], a_{j M L} (r_{j}^{'}) = 4 π i k_{j}^{3} \frac{1}{ϵ_{j}} h_{l}^{(1)} (k_{j} r_{j}^{'}) P . Φ_{l m}^{*} (θ_{j}^{'}, ϕ_{j}^{'}),

a_{j E H} (r_{j}^{'}) = 4 π k_{j}^{2} \sqrt{\frac{μ_{j}}{ϵ_{j}}} P . \nabla_{j}^{'} \times [j_{l} (k_{j} r_{j}^{'}) Φ_{l m}^{*} (θ_{j}^{'}, ϕ_{j}^{'})], a_{j M H} (r_{j}^{'}) = 4 π i k_{j}^{3} \frac{1}{ϵ_{j}} j_{l} (k_{j} r_{j}^{'}) P . Φ_{l m}^{*} (θ_{j}^{'}, ϕ_{j}^{'}) .

The dipole coefficients a_jEL and a_jEH contribute to the TM modes, whereas a_jML and a_jMH contribute to the TE modes [27,51], which has caused confusion in the literature [2]. Here, we use the subscripts E and M of the dipole coefficients to indicate their origin in the electric and magnetic multipole expansions, respectively. The

\nabla_{j}^{'}

symbol indicates derivatives with respect to the position

r_{j}^{'}

, and P is the dipole moment vector. Establishing field components in the layer j, we can construct the total field as the sum of fields in Eqs. (8) and (14), i.e.,

E H_{j}^{T} (r) = M_{j} (r) A_{j} + θ (r_{j}^{'} - r) M_{j} (r) a_{j L} + θ (r - r_{j}^{'}) M_{j} (r) a_{j H},

where the superscript T indicates total fields in layer j. The continuity of the transverse components of the electric and magnetic fields at the interface of regions j and j + 1 leads to

E H_{j}^{T} (r_{j}, θ, ϕ) = E H_{j + 1}^{T} (r_{j}, θ, ϕ),

M_{j} (r_{j}) [A_{j} + a_{j H}] = M_{j + 1} (r_{j}) [A_{j + 1} + a_{j}_{+ 1 L}],

A_{j + 1} = M_{j + 1}^{- 1} (r_{j}) M_{j} (r_{j}) A_{j} + M_{j + 1}^{- 1} (r_{j}) M_{j} (r_{j}) a_{j H} - a_{j + 1 L} .

Equation (20) can be used recursively to connect the coefficients in the outermost region, N + 1, to the innermost region, 1, leading to

A_{N + 1} = T (N + 1, 1) A_{1} + D,

where matrix T(N + 1, j) is defined by

\begin{array}{l} T (N + 1, j) = M_{N + 1}^{- 1} (r_{N}) M_{N} (r_{N}) M_{N}^{- 1} (r_{N - 1}) M_{N - 1} (r_{N - 1}) M_{N - 1}^{- 1} (r_{N - 2}) M_{N - 2} (r_{N - 2}) \dots \\ M_{j + 2}^{- 1} (r_{j + 1}) M_{j + 1} (r_{j + 1}) M_{j + 1}^{- 1} (r_{j}) M_{j} (r_{j}), for 1 \leq j < N + 1, \end{array}

T (N + 1, j) = I_{4 \times 4}, for j = N + 1 .

Note that T(N + 1, j) is always composed of repetitive blocks of matrices in the form of

M_{j + 1}^{- 1} (r_{j}) M_{j} (r_{j})

. Explicitly writing this matrix, it can be found that

M_{j + 1}^{- 1} (r_{j}) M_{j} (r_{j}) = \frac{i}{n_{j} + 1 μ_{j}} \frac{k_{j + 1}}{k_{j}} G (j + 1, j) = \frac{i}{n_{j} + 1 μ_{j}} \frac{k_{j + 1}}{k_{j}} {(\begin{matrix} G_{11} & G_{12} & 0 & 0 \\ G_{21} & G_{22} & 0 & 0 \\ 0 & 0 & G_{33} & G_{34} \\ 0 & 0 & G_{43} & G_{44} \end{matrix})}_{(j + 1, j)} .

The sub-matrices

G^{T M} = (\begin{matrix} G_{11} & G_{12} \\ G_{21} & G_{22} \end{matrix})

and

G^{T E} = (\begin{matrix} G_{33} & G_{34} \\ G_{43} & G_{44} \end{matrix}) .

take the form

(\begin{matrix} G_{L} ψ_{l}^{'} (k_{j} r_{j}) χ_{l} (k_{j}_{+ 1} r_{j}) - G_{R} ψ_{l} (k_{j} r_{j}) χ_{l}^{'} (k_{j + 1} r_{j}) & G_{L} χ_{l}^{'} (k_{j} r_{j}) χ_{l} (k_{j}_{+ 1} r_{j}) - G_{R} χ_{l} (k_{j} r_{j}) χ_{l}^{'} (k_{j}_{+ 1} r_{j}) \\ G_{L} ψ_{l}^{'} (k_{j}_{+ 1} r_{j}) ψ_{l} (k_{j} r_{j}) - G_{R} ψ_{l} (k_{j}_{+ 1} r_{j}) ψ_{l}^{'} (k_{j} r_{j}) & G_{L} ψ_{l}^{'} (k_{j}_{+ 1} r_{j}) χ_{l} (k_{j} r_{j}) - G_{R} ψ_{l} (k_{j}_{+ 1} r_{j}) χ_{l}^{'} (k_{j} r_{j}) \end{matrix})

for

G_{L} = μ_{j} n_{j + 1}^{2} / n_{j}

and G_R=µ_j₊₁n_j₊₁ in the case of G^TM, and G_L= µ_j₊₁n_j and G_L= µ_jn_j₊₁ in the case of G^TE.

The constant vector $D$ is given by

D = \sum_{j = 1}^{N + 1} T (N + 1, j) (1 - δ_{j, N + 1}) a_{j H} - T (N + 1, j) (1 - δ_{j, 1}) a_{j L} = (D_{1}, D_{2}, D_{3}, D_{4}) .

The sum in Eq. (27) is effectively over the regions that contain dipoles, since terms associated with regions with no dipoles are zero. Equation (22) can be inverted to obtain A₁ in terms of A_N₊₁ as follows,

A_{1} = T^{- 1} (N + 1, 1) A_{N + 1} + T^{- 1} (N + 1, 1) D .

Note that the constant vector

D

, contains information about the structure through the scattering matrix T, and information about the dipole sources, through a_jH and a_jL. However, matrices T(N + 1, 1) and T⁻¹(N + 1, 1) are independent of any source, and depend only on the parameters of the structure. Thus they represent the scattering matrix of the entire structure. To avoid confusion between the elements of matrices T(N + 1, 1) and T(N + 1, j), from now on, we use T ≡ T(N + 1, 1) and T^j ≡ T(N + 1, j) for j ≠ 1. The matrices of the form T(N + 1, j) are block diagonal matrices, since they have been built based on block diagonal matrices M. We define matrix

S \equiv (\begin{matrix} S^{T M} & 0 \\ 0 & S^{T E} \end{matrix}) = T^{- 1} (N + 1, 1) = T^{- 1}

and hence

S = {(\begin{matrix} T_{11} & T_{12} & 0 & 0 \\ T_{21} & T_{22} & 0 & 0 \\ 0 & 0 & T_{33} & T_{34} \\ 0 & 0 & T_{43} & T_{44} \end{matrix})}^{- 1} = (\begin{matrix} S_{11} & S_{12} & 0 & 0 \\ S_{21} & S_{22} & 0 & 0 \\ 0 & 0 & S_{33} & S_{34} \\ 0 & 0 & S_{43} & S_{44} \end{matrix}) .

In the innermost region, the coefficients of Hankel functions in Eq. (4) and (5) are zero, and hence the vector A₁ has the form A₁ = (A₁, 0, C₁, 0). Using this, Eq. (28) can be solved to find A_N₊₁ in terms of A₁ and C₁ as

A_{N + 1} = D_{1} + \frac{S_{22}}{(- S_{21} S_{12} + S_{11} S_{22})} A_{1},

B_{N + 1} = D_{2} - \frac{S_{21}}{(- S_{21} S_{12} + S_{11} S_{22})} A_{1},

C_{N + 1} = D_{3} - \frac{S_{44}}{(S_{43} S_{34} - S_{33} S_{44})} C_{1},

D_{N + 1} = D_{4} + \frac{S_{43}}{(S_{43} S_{34} - S_{33} S_{44})} C_{1} .

The coefficients A_N₊₁ = (A_N₊₁, B_N₊₁, C_N₊₁, D_N₊₁) determine the fields in the region N + 1, i.e., the outermost region, and hence are the scattering coefficients of the whole system of the microsphere and its sources. These coefficient are the same as a_n and b_n coefficients of Mie scattering, and represent different magnetic and electric dipole moments (B_N₊₁ corresponds to a_n, and D_N₊₁ corresponds to b_n). Once the values of the field coefficients in the outermost layer are known from Eqs. (30) through (33), the coefficients in any other layer, A_j, can be calculated recursively from Eq. (21). Then, Eqs. (6) and (7) may be used to calculate the total fields, E_j and H_j.

Equations (30) through (33) are general, and can be used for a range of scenarios:

If there are no sources (i.e. no dipoles in the structure and no incident wave) then $D = 0$ , and from the above equations one can find the ratios B_N+1/A_N+1 and D_N+1/C_N+1, which determine all the far field scattering coefficients. From these ratios one can find the TM and TE resonances of the structure, respectively (see Sections 7.1 and 7.2).
If there are dipole sources in the structure, then, without loss of generality, we can choose A_N+₁ = C_N₊₁ = 0, solve for A₁ and C₁ from Eqs. (30) and (32) and then find B_N₊₁ and D_N₊₁ from Eqs. (31) and (33), respectively. This is due to the fact that the total emitted power must be the same in both the near and far fields.
If there is only an incident field, then $D = 0$ , and A_N₊₁ and C_N₊₁ are known from the incident wave expansion, and in turn A₁ and C₁ are known, which means that B_N₊₁ and D_N₊₁ are also known.

2.4. Structure resonances

We consider the case where there are no sources-neither plane waves nor dipoles. Thus $D = 0$ in Eqs. (31) through (32) and we can find

B_{N + 1} = - \frac{S_{21}}{S_{22}} A_{N + 1},

D_{N + 1} = - \frac{S_{43}}{S_{44}} C_{N + 1} .

In the electric field of Eq. (4)B_N+₁ is the coefficient of the

h_{l}^{(1)} (k r)

term and A_N₊₁ is the coefficient of the j_l(kr) term in the outermost region N + 1. Hence the ratio of B_N₊₁/A_N₊₁ should approach infinity near a resonance with a transverse magnetic component only. Similarly, in the magnetic field of Eq. (5)D_N₊₁ is the coefficient of the

h_{l}^{(1)} (k r)

term and C_N₊₁ is the coefficient of the j_l(kr) term, in the outer layer N + 1. In this case, the ratio of D_N₊₁/C_N₊₁ should approach infinity near a resonance with a transverse electric component only. As a result, both S₂₁/S₂₂ and S₄₃/S₄₄ → ∞, which means the TM and TE resonances of the structure can be found by setting

T_{11} = 0 for TM resonances,

T_{33} = 0 for TE resonances .

These equations are in general multivalued, and for a given azimuthal quantum number l, the solutions to Eqs. (36) and (37) form the fundamental radial modes and their harmonics. Since the higher order harmonics have larger values of the Bessel function argument (kr) at their roots than the fundamental modes, all modes that appear within a wavelength range, including higher order modes, can be obtained. The numerical code associated with this paper (http://www.photonicsimulation.net) uses the optical and geometrical properties of any structure, to find the T matrix, and numerically solves Eqs. (36) and (37) for any given wavelengths range.

2.5. Scattered power in the outer region

We are interested in calculating the total radiated power of the system. Thus, we need to calculate the fields as r → ∞. In the outermost region N + 1 and for r ≫ r′_j, the total transverse parts of the fields in Eq. (18) are given by

E H_{N + 1}^{T} (r) = M_{N + 1} (r) A_{N + 1} + θ (r - r_{N + 1}^{^{'}}) M_{N + 1} (r) a_{N + 1 H} (r_{N + 1}^{^{'}}) .

In the limit of r → ∞, this equation leads to the following forms for the scattered fields,

\begin{array}{l} E_{s c} = - \sum_{l, m} (\frac{i}{n_{N + 1}}) \frac{1}{k_{N + 1} r} χ_{l}^{^{'}} (k_{N + 1} r) Ψ_{l m} (θ, ϕ) [B_{N + 1} + a_{N + 1 E H} (r_{N + 1}^{^{'}})] \\ + \sum_{l, m} Φ_{l m} (θ, ϕ) \frac{1}{k_{N + 1} r} χ_{l} (k_{N + 1} r) [D_{N + 1} + a_{N}_{+ 1 M H} (r_{N + 1}^{^{'}})], \end{array}

\begin{array}{l} H_{s c} = \sum_{l, m} (\frac{1}{μ_{N + 1}}) Φ_{l m} (θ, ϕ) \frac{1}{k_{N + 1} r} χ_{l} (k_{N + 1} r) [B_{N + 1} + a_{N}_{+ 1 E H} (r_{N + 1}^{^{'}})] \\ + \sum_{l, m} (\frac{i n_{N + 1}}{μ_{N + 1}}) \frac{1}{k_{N + 1} r} χ_{l}^{^{'}} (k_{N + 1} r) Ψ_{l m} [D_{N + 1} + a_{N + 1 M H} (r_{N + 1}^{^{'}})] . \end{array}

Note that in arriving at the above equation, we have used Eqs. (2), (9), and (10) and j_l(kr) → 0 as r → ∞. The total scattered power through a sphere of radius r can then be calculated by

\begin{array}{l} P_{total} = r^{2} \int S_{s c} . \hat{r} d Ω = \frac{c}{8 π μ_{N + 1}} r^{2} \int (E_{s c} \times μ_{N + 1} H_{s c}^{*}) . \hat{r} d Ω \\ = \frac{c}{8 π μ_{N + 1}} \sum_{l, m} (\frac{- i}{n_{N + 1} k_{N + 1}^{2}}) χ_{l}^{^{'}} k_{N + 1} r χ_{l}^{*} (k_{N + 1} r) {| B_{N + 1} + a_{N + 1}_{E H} (r_{N + 1}^{^{'}}) |}^{2} \\ + i \frac{n_{N + 1}}{k_{N + 1}^{2}} χ_{l} (k_{N + 1} r) χ_{l}^{*^{'}} (k_{N + 1} r) {| D_{N + 1} + a_{N + 1}_{M H} (r_{N + 1}^{^{'}}) |}^{2}, \end{array}

where we have used the orthonormal properties of Ψ_lm(θ, ϕ) and Φ_lm(θ, ϕ) functions given in Eq. (68) through (70). Note that in the limit of r → ∞,

χ_{l} (z) \to z h_{l}^{(1)} (z) \to i^{- l - 1}

exp(iz) and

χ_{l}^{^{'}} (z) \to i^{- l}

exp(iz) [2], and hence

P_{total} = \frac{c}{8 π} \sqrt{\frac{ϵ_{N + 1}}{μ_{N + 1}}} \frac{1}{k_{N + 1}^{2}} \sum_{l, m} [(\frac{1}{n_{N + 1}^{2}}) {| B_{N + 1} + a_{N + 1 E H} (r_{N + 1}^{^{'}}) |}^{2} + {| D_{N + 1} + a_{N + 1 M H} (r_{N + 1}^{'}) |}^{2}] .

Equation (42) is general, where B_N₊₁, D_N₊₁,

a_{N + 1 E H} (r_{N + 1}^{^{'}})

and

a_{N + 1 M H} (r_{N + 1}^{^{'}})

can be calculated based on Eqs. (31), (33), (15), and (16) respectively. Note that if there is no dipole in the outermost region, then

a_{N + 1 E H} (r_{N + 1}^{^{'}}) = a_{N + 1 M H} (r_{N + 1}^{^{'}}) = 0

.

2.6. A dipole in one layer

In this section, we assume that there exists only one dipole in the layer j, where j = 1,…, N + 1. Then, according to Scenario 2 of Section 2.3,

\begin{array}{l} A_{1} = - \frac{(- S_{21} S_{12} + S_{11} S_{22})}{S_{22}} D_{1}, & C_{1} = \frac{(S_{43} S_{34} - S_{33} S_{44})}{S_{44}} D_{3}, \\ B_{N + 1} = D_{2} + \frac{S_{21}}{S_{22}} D_{1}, & D_{N + 1} = D_{4} + \frac{S_{43}}{S_{44}} D_{3} . \end{array}

Then, using Eqs. (15), (27), and (29), we get

\begin{array}{l} D = T (N + 1, j) (1 - δ_{j, N + 1}) a_{j H} - T (N + 1, j) (1 - δ_{j, 1}) a_{j L} \\ = (\begin{array}{l} T_{12}^{j} (1 - δ_{j, N + 1}) a_{j E H} - T_{11}^{j} (1 - δ_{j, 1}) a_{j E L} \\ T_{22}^{j} (1 - δ_{j, N + 1}) a_{j E H} - T_{21}^{j} (1 - δ_{j, 1}) a_{j E L} \\ T_{34}^{j} (1 - δ_{j, N + 1}) a_{j M H} - T_{33}^{j} (1 - δ_{j, 1}) a_{j M L} \\ T_{44}^{j} (1 - δ_{j, N + 1}) a_{j M H} - T_{43}^{j} (1 - δ_{j, 1}) a_{j M L} \end{array}), \end{array}

based on which we can find

\begin{array}{l} B_{N + 1} + a_{N + 1 E H} (r_{N + 1}^{^{'}}) = \\ (T_{22}^{j} + \frac{S_{21}}{S_{22}} T_{12}^{j}) (1 - δ_{j, N + 1}) a_{j E H} (r_{j}^{^{'}}) - (T_{21}^{j} + \frac{S_{21}}{S_{22}} T_{11}^{j}) (1 - δ_{j, 1}) a_{j E L} (r_{j}^{^{'}}) + δ_{j, N + 1} a_{j E H} (r_{j}^{^{'}}) \\ = α_{l} a_{j E H} (r_{j}^{^{'}}) - β_{l} a_{j E L} (r_{j}^{^{'}}), \end{array}

where

α_{l} = (T_{22}^{j} + \frac{S_{21}}{S_{22}} T_{12}^{j}) (1 - δ_{j, N + 1}) + δ_{j, N + 1} and β_{l} = (T_{21}^{j} + \frac{S_{21}}{S_{22}} T_{11}^{j}) (1 - δ_{j, 1})

and

\begin{array}{l} D_{N + 1} + a_{N + 1 M H} (r_{N + 1}^{^{'}}) = \\ (T_{44}^{j} + \frac{S_{43}}{S_{44}} T_{34}^{j}) (1 - δ_{j, N + 1}) a_{j M H} (r_{j}^{^{'}}) - (T_{43}^{j} + \frac{S_{43}}{S_{44}} T_{33}^{j}) (1 - δ_{j, 1}) a_{j M L} (r_{j}^{^{'}}) + δ_{j, N + 1} a_{j M H} (r_{j}^{^{'}}) \\ = γ_{l} a_{j M H} (r_{j}^{^{'}}) - ζ_{l} a_{j M L} (r_{j}^{^{'}}), \end{array}

where

γ_{l} = (T_{44}^{j} + \frac{S_{43}}{S_{44}} T_{34}^{j}) (1 - δ_{j, N + 1}) + δ_{j, N + 1} and ζ_{l} = (T_{43}^{j} + \frac{S_{43}}{S_{44}} T_{33}^{j}) (1 - δ_{j, 1})

Note that α_l and β_l correspond to the contributions from the TM modes, whereas γ_l and ζ_l correspond to the contributions from the TE modes. Equations (45) and (47) are general, and can be applied to a dipole in any layer j, including the innermost layer 1 or the outermost layer N + 1. Considering Eqs. (16) and (17), we can rewrite $B_{N + 1} + a_{N + 1 E H} (r_{N + 1}^{^{'}})$ and $D_{N + 1} + a_{N + 1 M H} (r_{N + 1}^{^{'}})$ , which appear in the total scattered power in Eq. (42), as

B_{N + 1} + a_{N + 1 E H} (r_{N + 1}^{^{'}}) = 4 π k_{j}^{2} \sqrt{\frac{μ_{j}}{ϵ_{j}}} P . \nabla_{j}^{^{'}} \times {[α_{l} j_{l} (k_{j} r_{j}^{^{'}}) - β_{l} h_{l}^{(1)} (k_{j} r_{j}^{^{'}})] Φ_{l m}^{*} (θ_{j}^{^{'}}, ϕ_{j}^{^{'}})},

D_{N + 1} + a_{N + 1 M H} (r_{N + 1}^{^{'}}) = 4 π i k_{j}^{3} \frac{1}{ϵ_{j}} [γ_{l} j_{l} (k_{j} r_{j}^{^{'}}) - ζ_{l} h_{l}^{(1)} (k_{j} r_{j}^{^{'}})] P . Φ_{l m}^{*} (θ_{j}^{^{'}}, ϕ_{j}^{^{'}}) .

Using the properties of orthonormal functions Y_lm, Ψ_lm, and Φ_lm(θ, ϕ), and Eq (70), we note the forms of $B_{N + 1} + a_{N + 1 E H} (r_{N + 1}^{^{'}}) = P . [f_{l} (r_{j}^{^{'}}) Y_{l m} + g_{l} (r_{j}^{^{'}}) Ψ_{l m}]$ and $D_{N + 1} + a_{N + 1 M H} (r_{N + 1}^{^{'}}) = P . k_{l} (r_{j}^{^{'}}) Φ_{l m}^{*}$ , where $f_{l} (r_{j}^{^{'}})$ , $g_{l} (r_{j}^{^{'}})$ , and $k_{l} (r_{j}^{^{'}})$ are functions of $r_{j}^{'}$ , as in Eqs. (49) and (50). Hence we can use the properties of the dyadic products of Y_lm, Ψ_lm, and Φ_lm, Eq. (71), to perform the summation over m and simplify $\sum_{l, m} {| B_{N + 1} + a_{N + 1 E H} (r_{N + 1}^{'}) |}^{2}$ and $\sum_{l, m} {| D_{N + 1} + a_{N + 1 M H} (r_{N + 1}^{'}) |}^{2}$ , which appear in Eq. (42),

\begin{array}{l} \sum_{l, m} | B_{N + 1} + a_{N + 1 E H} (r_{N + 1}^{'}) |^{2} = \\ 16 π^{2} k_{j}^{6} (\frac{μ_{j}}{ϵ_{j}}) \sum_{l} {\frac{2 l + 1}{4 π} l (l + 1) \frac{|] α_{l} j_{l} (k_{j} r_{j}^{'}) - β_{l} h_{l}^{(1)} (k_{j} r_{j}^{'}) [|^{2}}{k_{j}^{2} r^{'}^{2}} | P_{r} |^{2} \\ + \frac{2 l + 1}{8 π} \frac{| {α_{l} \frac{d}{d r_{j}^{'}}] r_{j}^{'} j_{l} (k_{j} r_{j}^{'}) [- β_{l} \frac{d}{d r_{j}^{'}}] r_{j}^{'} h_{l}^{(1)} (k_{j} r_{j}^{'}) [} |^{2}}{k_{j}^{2} r_{j}^{' 2}} (| P_{θ} |^{2} + | P_{ϕ} |^{2})}, \end{array}

\begin{array}{l} \sum_{l, m} {| D_{N + 1} + a_{N + 1 M H} (r_{N + 1}^{'}) |}^{2} = \\ 16 π^{2} k_{j}^{6} (\frac{1}{ϵ_{j}^{2}}) \sum_{l} \frac{2 l + 1}{8 π} {| [γ_{l m} j_{l} (k_{j} r_{j}^{'}) - ζ_{l m} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2} ({| P_{θ} |}^{2} + {| P_{ϕ} |}^{2}) . \end{array}

Here, P_r, P_θ, and P_ϕ are the polar components of the polarization vector P. Based on the above equation, we can find the total scattered power from a sphere as the sum of powers due to normal and transverse components of P,

\begin{array}{l} P_{t o t a l} = P_{⊥} + P_{∥} = \frac{c}{2} \sqrt{\frac{ϵ_{N + 1}}{μ_{N + 1}}} \frac{k_{j}^{4} n_{j}^{2}}{n_{N + 1}^{2}} \frac{1}{ϵ_{j}^{2}} \sum_{l} (2 l + 1) {(\frac{n_{j}^{2}}{n_{N + 1}^{2}}) l (l + 1) \frac{{| [α_{l} j_{l} (k_{j} r_{j}^{'}) - β_{l} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2}}{k_{j}^{2} r_{j}^{' 2}} {| P_{r} |}^{2} \\ + [(\frac{n_{j}^{2}}{n_{N + 1}^{2}}) \frac{{| {α_{l} \frac{d}{d r_{j}^{'}} [r_{j}^{'} j_{l} (k_{j} r_{j}^{'})] - β_{l} \frac{d}{d r_{j}^{'}} [r_{j}^{'} h_{l}^{(1)} (k_{j} r_{j}^{'})]} |}^{2}}{k_{j}^{2} r_{j}^{' 2}} + {| [γ_{l} j_{l} (k_{j} r_{j}^{'}) - ζ_{l} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2}] [\frac{{| P_{θ} |}^{2} + {| P_{ϕ} |}^{2}}{2}]} . \end{array}

One can normalize the powers P_⊥ and P_∥ to powers radiated by a dipole in a bulk material with (n_j, ϵ_j, µ_j), i.e. $P_{⊥}^{0} = c k_{j}^{4} {| P_{r} |}^{2} / (3 ϵ_{j} n_{j})$ and $P_{∥}^{0} = c k_{j}^{4} ({| P_{θ} |}^{2} + {| P_{ϕ} |}^{2}) / (3 ϵ_{j} n_{j})$ to obtain

\frac{P_{⊥}}{P_{⊥}^{0}} = \frac{1}{2} \sqrt{\frac{ϵ_{N + 1}}{μ_{N + 1}}} \frac{n_{j}^{2}}{n_{N + 1}^{2}} \frac{3 n_{j}}{ϵ_{j}} \sum_{l} (\frac{n_{j}^{2}}{n_{N + 1}^{2}}) (2 l + 1) l (l + 1) \frac{{| [α_{l} j_{l} (k_{j} r_{j}^{'}) - β_{l} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2}}{k_{j}^{2} r_{j}^{' 2}}

\begin{array}{l} \frac{P_{∥}}{P_{∥}^{0}} = \frac{1}{4} \sqrt{\frac{ϵ_{N + 1}}{μ_{N + 1}}} \frac{n_{j}^{2}}{n_{N + 1}^{2}} \frac{3 n_{j}}{ϵ_{j}} \times \\ \sum_{l} (2 l + 1) {[(\frac{n_{j}^{2}}{n_{N + 1}^{2}}) \frac{{| {α_{l} \frac{d}{d r_{j}^{'}} [r_{j}^{'} j_{l} (k_{j} r_{j}^{'})] - β_{l} \frac{d}{d r_{j}^{'}} [r_{j}^{'} h_{l}^{(1)} (k_{j} r_{j}^{'})]} |}^{2}}{k_{j}^{2} r_{j}^{' 2}} + {| [γ_{l} j_{l} (k_{j} r_{j}^{'}) - ζ_{l} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2}]} \end{array}

2.7. One active layer

In this section, we consider a multilayer structure where one of the layers consists of active material. In this context, we add randomly oriented and uniformly distributed dipoles, with density $ρ (r_{j}^{'}) = 1$ , into that layer. As a result, one must integrate Eqs. (54) and (55) with respect to $r_{j}^{'}$ , which is located within the layer j. Since dipoles are randomly oriented in the layer, we can write $〈 P_{t o t a l} / P^{0} 〉 = \frac{1}{3} 〈 P_{⊥} / P_{⊥}^{0} 〉 + \frac{2}{3} 〈 P_{∥} / P_{∥}^{0} 〉$ , where

〈 \frac{P_{⊥}}{P_{⊥}^{0}} 〉 = \frac{1}{2} \sqrt{\frac{ϵ_{N + 1}}{μ_{N + 1}}} \frac{n_{j}^{2}}{n_{N + 1}^{2}} \frac{3 n_{j}}{ϵ_{j}} \sum_{l} (\frac{n_{j}^{2}}{n_{N + 1}^{2}}) (2 l + 1) l (l + 1) \int_{j} \frac{{| [α_{l} j_{l} (k_{j} r_{j}^{'}) - β_{l} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2}}{k_{j}^{2} r_{j}^{' 2}} d^{3} r_{j}^{'} / \int_{j} d^{3} r_{j}^{'}

= \frac{1}{2} \sqrt{\frac{ϵ_{N + 1}}{μ_{N + 1}}} \frac{n_{j}^{2}}{n_{N + 1}^{2}} \frac{3 n_{j}}{k_{j}^{2} ϵ_{j} V_{j shell}} 4 π \sum_{l} (\frac{n_{j}^{2}}{n_{N + 1}^{2}}) l (l + 1) I_{l}^{(1)},

I_{l}^{(1)} = (2 l + 1) \int_{j} {| [α_{l} j_{l} (k_{j} r_{j}^{'}) - β_{l} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2} d r_{j}^{'},

where the volume of the shell is

V_{j shell} = 4 π (r_{j}^{' 2} - r_{j - 1}^{' 2})

. If j = 1, then by convention, r₀ is set to zero, as in that case the volume is simply the sphere bounded by the innermost radius. Similarly,

\begin{array}{l} 〈 \frac{P_{∥}}{P_{∥}^{0}} 〉 = \frac{1}{4} \sqrt{\frac{ϵ_{N + 1}}{μ_{N + 1}}} \frac{n_{j}^{2}}{n_{N + 1}^{2}} \frac{3 n_{j}}{k_{j}^{2} ϵ_{j}} (\frac{n_{j}^{2}}{n_{N + 1}^{2}}) 4 π \sum_{l} (2 l + 1) {[\int_{j} (\frac{n_{j}^{2}}{n_{N + 1}^{2}}) | {α_{l} \frac{d}{d r_{j}^{'}} [r_{j}^{'} j_{l} (k_{j} r_{j}^{'})] \\ {- β_{l} \frac{d}{d r_{j}^{'}} [r_{j}^{'} h_{l}^{(1)} (k_{j} r_{j}^{'})]} |}^{2} + k_{j}^{2} r_{j}^{' 2} {| [γ_{l} j_{l} (k_{j} r_{j}^{'}) - ζ_{l} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2} d r_{j}^{'} / \int_{j} d^{3} r_{j}^{'}]}, \end{array}

= \frac{1}{4} \sqrt{\frac{ϵ_{N + 1}}{μ_{N + 1}}} \frac{n_{j}^{2}}{n_{N + 1}^{2}} \frac{3 n_{j}}{K_{j}^{2} ϵ_{j} V_{j shell}} 4 π \sum_{l} (\frac{n_{j}^{2}}{n_{N + 1}^{2}}) I_{l}^{(2)} + I_{l}^{(3)},

I_{l}^{(2)} = (2 l + 1) \int_{j shell} {| {α_{l} \frac{d}{d r_{j}^{'}} [r_{j}^{'} j_{l} (k_{j} r_{j}^{'})] - β_{l} \frac{d}{d r_{j}^{'}} [r_{j}^{'} h_{l}^{(1)} (k_{j} r_{j}^{'})]} |}^{2} d r_{j}^{'},

I_{l}^{(3)} = (2 l + 1) \int_{j shell} k_{j}^{2} r_{j}^{' 2} {| [γ_{l} j_{l} (k_{j} r_{j}^{'}) - ζ_{l} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2} d r_{j}^{'} .

The total power is then

\begin{array}{l} 〈 \frac{P_{t o t a l}}{P^{0}} 〉 = \frac{1}{3} 〈 \frac{P_{⊥}}{P_{⊥}^{0}} 〉 + \frac{2}{3} 〈 \frac{P_{∥}}{P_{∥}^{0}} 〉 \\ = \frac{1}{2} \sqrt{\frac{ϵ_{N + 1}}{μ_{N + 1}}} \frac{n_{j}^{2}}{n_{N + 1}^{2}} \frac{n_{j}}{k_{j}^{2} ϵ_{j} V_{j shell}} 4 π \sum_{l} [(\frac{n_{j}^{2}}{n_{N + 1}^{2}}) l (l + 1) I_{l}^{(1)} + (\frac{n_{j}^{2}}{n_{N + 1}^{2}}) I_{l}^{(2)} + I_{l}^{(3)}] . \end{array}

Defining a functional form

Ψ [p_{l} q_{l}] (z) = \int p_{l} (z) q_{l} (z) d z = constant + \frac{1}{2} {[z - \frac{l (l + 1)}{z}] p_{l} q_{l} - \frac{1}{2} (p_{l} q_{l}^{'} + p_{l}^{'} q_{l}) + z p_{l}^{'} q_{l}^{'}},

where p_l(z) and q_l(z) can be any of ψ(x) = x j_l(x) or

χ (x) = x h_{l}^{(1)} (x)

, one can then find

\begin{array}{l} [l (l + 1) I_{l}^{(1)} + I_{l}^{(2)}] = \\ \frac{1}{k_{j}} {{| α_{l} |}^{2} (l + 1) (Ψ [{| ψ_{l - 1} |}^{2}] (k_{j} r_{j}) - Ψ [{| ψ_{l - 1} |}^{2} (k_{j} r_{j - 1})) + {| α_{l} |}^{2} l (Ψ [{| ψ_{l + 1} |}^{2}] (k_{j} r_{j}) - Ψ [{| ψ_{l + 1} |}^{2}] (k_{j} r_{j - 1})) \\ + {| β_{l} |}^{2} (l + 1) (Ψ [{| χ_{l - 1} |}^{2}] (k_{j} r_{j}) - Ψ [{| χ_{l - 1} |}^{2}] (k_{j} r_{j - 1})) + {| β_{l} |}^{2} l (Ψ [{| χ_{l + 1} |}^{2}] (k_{j} r_{j}) - Ψ [{| χ_{l + 1} |}^{2}] (k_{j} r_{j - 1})) \\ - (l + 1) (α_{l} β_{l}^{*} (Ψ [ψ_{l - 1} χ_{l - 1}^{*}] (k_{j} r_{j}) - Ψ [ψ_{l - 1} χ_{l - 1}^{*}) (k_{j} r_{j - 1})] \\ + α_{l}^{*} β_{l} (Ψ [ψ_{l - 1}^{*} χ_{l - 1}] (k_{j} r_{j}) - Ψ [ψ_{l - 1}^{*} χ_{l - 1}] (k_{j} r_{j - 1}))) \\ - l (α_{l} β_{l}^{*} (Ψ [ψ_{l + 1} χ_{l + 1}^{*}] (k_{j} r_{j}) - Ψ [ψ_{l + 1} χ_{l + 1}^{*}] (k_{j} r_{j - 1})) \\ + α_{l}^{*} β_{l} (Ψ [ψ_{l + 1}^{*} χ_{l + 1}] (k_{j} r_{j}) - Ψ [ψ_{l + 1}^{*} χ_{l + 1}] (k_{j} r_{j - 1})))}, \end{array}

where r_j and r_j₋₁ are the radii of the upper and lower interfaces of the region j, respectively. Similarly, we can calculate

\begin{array}{l} I_{l}^{(3)} = (2 l + 1) \int_{j shell} k_{j}^{2} r_{j}^{' 2} {| [γ_{l} j_{l} (k_{j} r_{j}^{'}) - ζ_{l} h_{l}^{(1)} (k_{j} r_{j}^{'})] |}^{2} d r_{j}^{'} = \\ \frac{1}{k_{j}} (2 l + 1) {{| γ_{l} |}^{2} (Ψ [{| ψ_{l} |}^{2}] (k_{j} r_{j}) - Ψ [{| ψ_{l} |}^{2}] (k_{j} r_{j - 1})) + {| ζ_{l} |}^{2} (Ψ [{| χ_{l} |}^{2}] (k_{j} r_{j}) - Ψ [{| χ_{l} |}^{2}] (k_{j} r_{j - 1})) \\ - (γ_{l} ζ_{l}^{*} (Ψ [ψ_{l} χ_{l}^{*}] (k_{j} r_{j}) - Ψ [ψ_{l} χ_{l}^{*}] (k_{j} r_{j - 1})) + γ_{l}^{*} ζ_{l} (Ψ [ψ_{l}^{*} χ_{l}] (k_{j} r_{j}) - Ψ [ψ_{l}^{*} χ_{l}] (k_{j} r_{j - 1})))} . \end{array}

With

[l (l + 1) I_{l}^{(1)} + I_{l}^{(2)}]

and

I_{l}^{(3)}

now known, the total averaged power can be calculated from Eq. (63).

3. Demonstration and discussion

We now consider several scenarios of interest that can be uniquely treated using this model. First, we examine the behavior of the WGM spectrum of a silica microsphere (with dispersion included [52]) coated with a single high refractive index layer (n₂ = 1.7), surrounded by water (n₃ = 1.33), as the thickness of the coating d is changed from 5 nm to 15 nm. This scenario corresponds to N = 2, where N is the number of layers, including the central silica sphere, and the high index layer. Figure 2 shows the results for an electric dipole placed just outside the surface, with an orientation perpendicular or parallel to the surface of the sphere. A range of wavelengths 0.59–0.61 µm is simulated, and the outer diameter is kept fixed at 25 µm. It is found that there is a systematic shift in the prominent WGM peaks towards higher wavelengths, as the thickness of the layer is increased. The free spectral range, however, remains largely unchanged over this range of wavelength values. In the limit d → 0, the results match the simple case of the Chew model, as anticipated [20].

Fig. 2 Spectra obtained from a simulated silica microsphere, incorporating dispersion, coated with a high refractive index layer (n₂ = 1.7) with a diameter of 25 µm, surrounded by water. A single electric dipole is oriented (a) in the radial direction and (b) in the tangential direction.

Download Full Size | PDF

The thin vertical lines marking the position of the resonances, with the corresponding mode numbers and labels, are obtained from the structure resonance positions in Eqs. (36) and (37) for a two-layer microsphere, by setting N = 2 in our general formalism. This can be done in a straightforward way in the code associated with this paper (see pg. 10). The thick vertical lines indicate the resonance positions obtained from the Arnold formalism [24,26], which agree exactly. Note that in the case of the parallel excitation in Fig. 2(b), there is a small contribution from the TM mode as well, as expected from Eq. (55).

A silica microsphere of the same size, 25 µm, is then modeled as being coated with a polymer (PMMA) with dispersion incorporated through the Sellmeier equations for both materials [52,53], as shown in Fig. 3. Across the wavelength range of 0.59 to 0.61 µm, the refractive index of silica varies from 1.4584 to 1.4577, wheres the index of PMMA varies from 1.4913 to 1.4902, respectively. PMMA is simulated in this example because it can act as an active layer [47] and is a straightforward way of testing the functionality of the model. At the diameter considered, higher order modes are not strongly coupled to, and so both the TE and TM modes remain distinct, as can be seen in Fig. 3(a). The refractive index contrast between silica and PMMA is relatively small, and so the dependence of the mode positions on the thickness of the coating is more mild. A close-up view of the plot is shown in Fig. 3(b). The vertical lines indicate the full-width at half-maximum (FWHM) positions for a variety of peaks, from which the Q-factors can be extracted. The Q-factor is estimated directly from the width of a spectral peak δλ, centred around a wavelength λ₀, through Q = λ₀/δλ. Q-factors corresponding to a selection of TM and TE modes are also shown above each peak. Note that for the TM modes (the left three peaks) there is a decrease in the Q-factor as d increases, but for the TE modes there is an apparent increase in the Q-factor.

Fig. 3 Spectra for a silica microsphere coated with a polymer layer (PMMA), both of which include dispersion. The polymer layer functions as an active layer. (a) Both TE and TM modes are excited. (b) A zoomed in plot showing the FWHM of several peaks and their Q-factors.

Download Full Size | PDF

The sensitivity of the WGM peaks can be examined by varying the refractive index of the surrounding medium, n₃. In Fig. 4(a), d = 10 nm, the mean TE peak shift leads to a sensitivity of S ≡ dλ/dn₃ = 60.0 nm/R.I.U. For a surrounding index of 1.36, the TE and TM modes are broad and overlapping, and an extraction of the Q-factor using the FWHM approach will necessarily include contributions from both mode polarizations. This is a consequence of the fact that the Q-factors obtained from the spectrum include both intrinsic, non-radiative contributions from the geometry of the structure, and contributions from the relative coupling of the active layer to the modes. This makes it more difficult to separate nearby modes and define independent Q-factors. We also calculate a figure of merit (FOM), defined as Q.S [43], to assess the sensing performance of the microspheres. For a given sensitivity S, the value of Q chosen is the mean value of the peak as it shifts to its new position, as a function of n₃. The inset of Fig. 4(a) shows a decrease in the FOM as a function of λ. In Fig. 4(b), d = 50 nm, both the sensitivity and the FOM are larger, with S = 63.3 nm/R.I.U.

Fig. 4 The sensitivity of silica microspheres coated with PMMA as a function of the surrounding refractive index, for two example layer thicknesses. (a) d = 10 nm, (b) d = 50 nm. Inset: the figure of merit (FOM), Q.S, in units of 10⁵ nm/R.I.U., as a function of λ.

Download Full Size | PDF

The summation $\sum_{l = 1}^{\infty}$ , in the expressions for P/P⁰(λ), is calculated to an upper bound, l_max, determined by a prescribed tolerance τ, so that

\frac{| {(P / P^{0})}_{l = l_{\max + 1}} - {(P / P^{0})}_{l = l_{\max}} |}{{(P / P^{0})}_{l = l_{\max + 1}}} < τ .

This prescription is sufficient so long as the behavior of P/P⁰ is convergent, which is usually the case, except for unstable parameter regions, described below. At each value of l > 1, the spherical Bessel and Hankel functions are calculated using the recursion relations, and function calls are minimized to improve the efficiency of the calculation.

Examples of the scaling behavior of the execution times (T) for the functions P/P⁰(λ), with respect to wavelength, for numbers of layers N = 1, 2 and 3 are shown in Fig. 5 for a fixed outer diameter of 25 µm. The results are fairly insensitive to the layer thickness, allowing Fig. 5 to be a fairly accurate measure of the execution time for a given number of layers and prescribed tolerance, τ. It was found that the implementation of the recursion relations resulted in an improvement of approximately one order of magnitude in the execution time compared to function-call methods.

Fig. 5 The execution time T of the formulae P/P⁰(λ) as a function of wavelength. The results are shown for a numbers of layers N = 1, 2 and 3. The results for a single-dipole excitation oriented (a) parallel, and (b) perpendicular to the surface of the sphere are similar in magnitude. (c) The results of a single uniform distribution of dipoles within the center of the sphere begin to plateau as the wavelength becomes small. The tolerance selected is τ = 1 × 10¹².

Download Full Size | PDF

4. Conclusion

A method for the modeling of whispering gallery modes in optical resonators, including various excitation scenarios that closely mirror experimental techniques, is an important step toward facilitating the design and analysis of novel resonator architectures. We developed an algorithm based on the solutions of the boundary value problem for multilayer spherical resonators, with improved execution times compared to standard functional methods. The model is able to handle an arbitrary number of concentric, spherical dielectric layers, and extract the resonance positions from the characteristic equation. Formulae for the power spectrum are derived for the case of a single dipole or an active layer source, and the behavior of the mode positions are examined for several different scenarios. Estimates of the computation time for the normalized emitted power at a single wavelength, and improvements to stability issues inherent in the transfer-matrix approach, have been discussed. The spectrum simulator reported represents an important step towards a general, fast and efficient method for extracting the underlying parameters and properties of a given resonator purely from its spectrum, simulating spectra over a wide parameter space, and predicting features of novel resonator designs for photonic band-gap devices and biosensing applications.

5. Appendix A: Properties of vector spherical harmonics (VSH)

The VSH are orthogonal in the usual three-dimensional sense,

Y_{l m} \cdot Ψ_{l m} = Φ_{l m} \cdot Ψ_{l m} = Φ_{l m} \cdot Y_{l m} = 0,

and also orthonormal in the Hilbert space

\begin{array}{l} \int Y_{l m} \cdot Y_{l^{'} m^{'}}^{*} d Ω = \int Ψ_{l m} \cdot Ψ_{l^{'} m^{'}}^{*} d Ω = \int Φ_{l m} \cdot Φ_{l^{'} m^{'}}^{*} d Ω = δ_{l l^{'}} δ_{m m^{'}}, \\ \int Y_{l m} \cdot Ψ_{l^{'} m^{'}}^{*} d Ω = \int Ψ_{l m} \cdot Φ_{l^{'} m^{'}}^{*} d Ω = \int Y_{l m} \cdot Φ_{l^{'} m^{'}}^{*} d Ω = 0, \end{array}

In addition, it can be shown that

\begin{array}{l} \nabla \times (f (r) Y_{l m}) = \frac{1}{r} f (r) Φ_{l m}; \nabla \times (f (r) Ψ_{l m}) = (\frac{d f}{d r} + \frac{1}{r} f (r)) Φ_{l m}, \\ \nabla \times (f (r) Φ_{l m}) = - \frac{\sqrt{l (l + 1)}}{i r} f Y_{l m} - (\frac{d f}{d r} + \frac{1}{r} f) Ψ_{l m} . \end{array}

The following summation rules apply to dyadic products of the VSH:

\begin{array}{l} \sum_{m = - l}^{m = l} Φ_{l m} Φ_{l m}^{*} = \sum_{m = - l}^{m = l} Ψ_{l m} Ψ_{l m}^{*} = \frac{2 l + 1}{8 π} (e_{θ} e_{θ} + e_{ϕ} e_{ϕ}); \sum_{m = - l}^{m = l} Y_{l m} Y_{l m}^{*} = \frac{2 l + 1}{4 π} (e_{r} e_{r}); \\ \sum_{m = - l}^{m = l} Φ_{l m} Y_{l m}^{*} = \sum_{m = - l}^{m = l} Φ_{l m} Ψ_{l m}^{*} = \sum_{m = - l}^{m = l} Y_{l m} Ψ_{l m}^{*} = 0 . \end{array}

6. Appendix B: Relation between conventions

Depending on the application, different definitions have been used for the VSH. The following functions have been used in the literature

X_{l m} (θ, ϕ) = (\frac{1}{i}) (\frac{1}{\sqrt{l (l + 1)}}) r \times \nabla Y_{l m} (θ, ϕ) [23, 51]

X_{l l m} (θ, ϕ) = (\frac{1}{i}) (\frac{1}{\sqrt{l (l + 1)}}) r \times \nabla Y_{l m} (θ, ϕ) Eq . (5.9.14) of [54]

\begin{array}{l} Y_{L}^{(m)} = (\frac{1}{i}) (\frac{1}{\sqrt{l (l + 1)}}) r \times \nabla Y_{l m} (θ, ϕ); Y_{L}^{(e)} = - (\frac{1}{i}) (\frac{1}{\sqrt{l (l + 1)}}) r \nabla Y_{l m} (θ, ϕ); \\ (\frac{- 1}{i}) Y_{l m} (θ, ϕ) \hat{r} = Y_{L}^{(o)} [2] \end{array}

Comparing these equations with the orthonormal functions used here, Eqs. (1), we find,

Y_{l l m} = X_{l m} = Y_{L}^{(m)} = Φ_{l m}, Y_{L}^{(e)} = - Ψ_{l m}, (\frac{1}{i}) Y_{L}^{(o)} = Y_{l m} .

7. Appendix C: Examples

7.1. Example 1: A microsphere

For a microsphere, we set N = 1, and then using Eq. (23), the T matrix is

T (2, 1) = M_{2}^{- 1} (r_{1}) M_{1} (r_{1}),

and hence we can use Eq. (25) and (26) with j = 1 to find

T_{11} = \frac{i}{n_{2} μ_{1}} G_{11} = \frac{i}{n_{2} μ_{1}} \frac{k_{2}}{k_{1}} \frac{n_{2}}{n_{1}} [n_{2} μ_{1} ψ_{l}^{^{'}} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - n_{1} μ_{2} ψ_{l} (k_{1} r_{1}) χ_{l}^{^{'}} (k_{2} r_{1})],

T_{33} = \frac{i}{n_{2} μ_{1}} \frac{k_{2}}{k_{1}} G_{33} = \frac{i}{n_{2} μ_{1}} \frac{k_{2}}{k_{1}} [n_{1} μ_{2} ψ_{l}^{^{'}} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - n_{2} μ_{1} ψ_{l} (k_{1} r_{1}) χ_{l}^{^{'}} (k_{2} r_{1})] .

By assuming µ₁ = µ₂, the TM and TE resonance conditions, T₁₁ = 0 and T₃₃ = 0, lead to

n_{2} \frac{ψ_{l}^{^{'}} (k_{1} r_{1})}{ψ_{l} (k_{1} r_{1})} = n_{1} \frac{χ_{l}^{^{'}} (k_{2} r_{1})}{χ_{l} (k_{2} r_{1})} TM resonance condition,

n_{1} \frac{ψ_{l}^{^{'}} (k_{1} r_{1})}{ψ_{l} (k_{1} r_{1})} = n_{2} \frac{χ_{l}^{^{'}} (k_{2} r_{1})}{χ_{l} (k_{2} r_{1})} TE resonance condition .

Eqs. (79) and (80) exactly match Eqs. (19) and (13) in [25] and Eq. (4.53) in [27], respectively. It should be noted that in Eq. (33) of [2] the TE and TM modes are interchanged.

7.1.1. Special case: A dipole in the outer region

Let us assume a dipole source is located in the outermost region, i.e., j = N + 1 = 2. Then $D$ from Eqs. (27) and (23) becomes

D = - T^{2} (2, 2) a_{2 L} = - I_{4 \times 4} a_{2 L}

Using Eq. (45) and (47), we can find

B_{2} + a_{2 E H} (r_{j}^{^{'}})

and

D_{2} + a_{2 M H} (r_{j}^{^{'}})

in

B_{2} + a_{2 E H} (r_{2}^{'}) = α_{l} α_{2 E H} - β_{l} α_{2 E H L}; where α_{l} = 1 and β_{l} = (T_{21}^{2} + \frac{S_{21}}{S_{22}} T_{11}^{2}) = \frac{S_{21}}{S_{22}}

D_{2} + a_{2 M H} (r_{2}^{'}) = γ_{l} α_{2 M H} - ζ_{l} α_{2 M L}; where γ_{l} = 1 and ζ_{l} = (T_{43} + \frac{S_{43}}{S_{44}} T_{33}) = \frac{S_{43}}{S_{44}}

Having identified the coefficient α_l, β_l, γ_l, ζ_l, we can find the power due to normal and transverse components using Eq. (54) and (55) as

\frac{P_{⊥}}{P_{⊥}^{0}} = \frac{3}{2} \sum_{l} (2 l + 1) l (l + 1) \frac{| [j_{l} (k_{2} r_{2}^{'}) + \frac{S_{21}}{S_{22}} h_{l}^{(1)} (k_{2} r_{2}^{'})] |^{2}}{k_{2}^{2} r_{2}^{' 2}},

\frac{P_{| |}}{P_{| |}^{0}} = \frac{3}{4} \sum_{l} (2 l + 1) {[\frac{{| {\frac{d}{d r_{2}^{'}} [r_{2}^{'} j_{l} (k_{2} r_{2}^{'})] + \frac{S_{21}}{S_{22}} \frac{d}{d r_{2}^{'}} [r_{2}^{'} h_{l}^{(1)} (k_{2} r_{2}^{'})]} |}^{2}}{k_{2}^{2} r_{2}^{' 2}} + {| [j_{l} (k_{2} r_{2}^{'}) + \frac{S_{43}}{S_{44}} h_{l}^{(1)} (k_{2} r_{2}^{'})] |}^{2}]} .

Now since S = T⁻¹, Eq. (29), then

\frac{S_{21}}{S_{22}} = - \frac{T_{21}}{T_{11}}

and

\frac{S_{43}}{S_{44}} = - \frac{T_{43}}{T_{33}}

In addition,

T (2, 1) = M_{2}^{- 1} (r_{1}) M_{1} (r_{1})

and hence we can use Eq. (25) and (26) to find

\frac{S_{21}}{S_{22}} = - \frac{[n_{1} μ_{2} ψ_{l}^{'} (k_{2} r_{1}) ψ_{l} (k_{1} r_{1}) - n_{2} μ_{1} ψ_{l} (k_{2} r_{1}) ψ_{l}^{'} (k_{1} r_{1})]}{[n_{2} μ_{1} ψ_{l}^{'} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - n_{1} μ_{2} ψ_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1})]}

= \frac{[(ε_{2} / n_{2}) ψ_{l} (k_{2} r_{1}) ψ_{l}^{'} (k_{1} r_{1}) - (ε_{1} / n_{1}) ψ_{l}^{'} (k_{2} r_{1}) ψ_{l} (k_{1} r_{1})]}{[(ε_{2} / n_{2}) ψ_{l}^{'} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - (ε_{1} / n_{1}) ψ_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1})]},

\frac{S_{43}}{S_{44}} = - \frac{n_{2} μ_{1} ψ_{l}^{'} (k_{2} r_{1}) ψ_{l} (k_{1} r_{1}) - n_{1} μ_{2} ψ_{l} (k_{2} r_{1}) ψ_{l}^{'} (k_{1} r_{1})}{n_{1} μ_{2} ψ_{l}^{'} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - n_{2} μ_{1} ψ_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1})}

= \frac{(ε_{1} / n_{1}) ψ_{l} (k_{2} r_{1}) ψ_{l}^{'} (k_{1} r_{1}) - (ε_{2} / n_{2}) ψ_{l}^{'} (k_{2} r_{1}) ψ_{l} (k_{1} r_{1})}{(ε_{1} / n_{1}) ψ_{l}^{'} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - (ε_{2} / n_{2}) ψ_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1})} .

Having found the coefficients

\frac{S_{21}}{S_{22}}

and

\frac{S_{43}}{S_{44}}

, we can find

\frac{P_{⊥}}{P_{⊥}^{0}}

and

\frac{P_{| |}}{P_{| |}^{0}}

. Equations (86) and (88) are the same as Eqs. (6) and (7) of Chew [20].

7.1.2. Special case: A dipole in the inner region

Let us assume a dipole in the innermost region, then N = 1, j = 1 and from Eqs. (15), (23) and (27) we find $D = T (2, 1) a_{j_{H}} = (T_{12} a_{1}_{E H}, T_{22} a_{1}_{E H}, T_{34} a_{1}_{M H}, T_{44} a_{1}_{M H})$ . Using Eq. (45) and (47), we can find $B_{2} + a_{2 E H} (r_{1}^{'})$ and $D_{2} + a_{2 M H} (r_{1}^{'})$ ,

B_{N + 1} + a_{N + 1 E H} (r_{2}^{'}) = (T_{22} + \frac{S_{21}}{S_{22}} T_{12}) a_{j E H}

= α_{l} a_{1}_{E H} - β_{l} a_{1 E L}; where α_{l} = (T_{22} + \frac{S_{21}}{S_{22}} T_{12}) and β_{l} = 0

D_{N + 1} + a_{N + 1 M H} (r_{2}^{'}) = (T_{44} + \frac{S_{43}}{S_{44}} T_{34}) a_{1 M H}

= γ_{l} a_{1 M H} - ζ_{l} a_{1 M L}; where γ_{l} = (T_{44} + \frac{S_{43}}{S_{44}} T_{34}) and ζ_{l} = 0

According to Eq. (27)S = T⁻¹(2, 1) hence (S₂₂T₂₂ + S₂₁T₁₂) = (S₄₄T₄₄ + S₄₃T₃₄) = 1 and as a result,

α_{l} = \frac{1}{S_{22}}

and

γ_{l} = \frac{1}{S_{44}}

Now, since S = T⁻¹(2, 1), then S₂₂ = det(T^TM)⁻¹T₁₁ and S₄₄ = det(T^TE)⁻¹T₃₃ and hence we can use Eqs. (25) and (26) to find

S_{22} = {[\det {(M_{2})}^{- 1} \det (M_{2})]}^{- 1} \frac{i}{n_{2} μ_{1}} \frac{k_{2}}{k_{1}} \frac{n_{2}}{n_{1}} [n_{2} μ_{1} ψ_{l}^{'} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - n_{1} μ_{2} ψ_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1})],

= - i \frac{k_{1} n_{1}}{ϵ_{2} ϵ_{1}} \sqrt{\frac{ϵ_{2}}{μ_{2}}} \frac{n_{1}}{n_{2}} [ϵ_{1} r_{1} j_{l} (k_{1} r_{1}) χ' (k_{2} r_{1}) - ϵ_{2} ψ_{l}^{'} (k_{1} r_{1}) r h_{l}^{(1)} (k_{2} r_{1})] .

Similarly, we find

S_{44} = {[\det {(M_{2})}^{- 1} \det (M_{1})]}^{- 1} \frac{i}{n_{2} μ_{1}} \frac{k_{2}}{k_{1}} [n_{1} μ_{2} ψ_{l}^{'} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - n_{2} μ_{1} ψ_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1})],

= i (\frac{n_{1}}{μ_{2}}) [\frac{n_{2}}{ϵ_{2}} ψ_{l}^{'} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - \frac{n_{1}}{ϵ_{1}} ψ_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1})] .

Thus, we can find

\frac{P_{⊥}}{P_{⊥}^{0}}

and

\frac{P_{| |}}{P_{| |}^{0}}

as

\frac{P_{| |}}{P_{| |}^{0}} = \frac{1}{4} \sqrt{\frac{ϵ_{2}}{μ_{2}}} \frac{3 n_{1} ϵ_{1}}{{(k_{1} r_{1})}^{2}} \sum_{l} (2 l + 1) {\frac{{| \frac{d}{d r_{1}^{'}} [r_{j}^{'} j_{l} (k_{1} r_{1}^{'})] |}^{2}}{k_{1}^{2} r_{1}^{' 2} | D_{l} |^{2}} + \frac{μ_{1} μ_{2}}{ϵ_{2} ϵ_{1}} \frac{{| j_{l} (k_{1} r_{1}^{'}) |}^{2}}{| D_{l}^{'} |^{2}}},

\frac{P_{⊥}}{P_{⊥}^{0}} = \frac{1}{2} \sqrt{\frac{ϵ_{2}}{μ_{2}}} \frac{3 n_{1} ϵ_{1}}{{(k_{1} r_{1})}^{2}} \sum_{l} (2 l + 1) l (l + 1) \frac{| j_{l} (k_{1} r_{1}^{'}) |^{2}}{k_{1}^{2} r_{1}^{' 2} | D_{l} |^{2}},

for

D_{l} = [ϵ_{1} j_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1}) - ϵ_{1} ψ_{l}^{'} (k_{1} r_{1}) h_{l}^{(1)} (k_{2} r_{1})], D_{l}^{'} = D_{l} (ϵ \leftrightarrow μ) .

These are exactly the same as Eqs. (1, 3a) and Eq. (2,3) of Chew [9].

7.1.3. Special case: A microsphere with an active inner region

Now assume that the inner region is filled with randomly oriented dipoles. Then, we can use Eqs. (63), (65), and (66) together with $α_{l} = \frac{1}{S_{22}}$ , β_l = 0, ζ = 0, and $γ_{l} = \frac{1}{S_{44}}$ , calculated in the previous section, to find

\frac{1}{(2 l + 1)} [l (l + 1) I_{l}^{(1)} + I_{l}^{(2)}] = \frac{α_{l}^{2}}{(2 l + 1) k_{1}} {(l + 1) Ψ [ψ_{l - 1}^{2}] (k_{1} r_{1}) + l^{2} Ψ [ψ_{l + 1}^{2}] (k_{1} r_{1})} .

Note that in the case of a microsphere, r₀is set to zero, as well as all functionals of the form Ψ[·](k₁r₀) = 0. Similarly, we can also find

I_{l}^{(3)} = \frac{1}{k_{1}} (2 l + 1) γ_{l}^{2} Ψ [ψ_{l}^{2}] (k_{1} r_{1})

Having

\frac{1}{(2 l + 1)} [l (l + 1) I_{l}^{(1)} + I_{l}^{(2)}]

and

I_{l}^{(3)}

, we find

〈 \frac{P_{t o t a l}}{P^{0}} 〉 = \frac{1}{3} 〈 \frac{P_{⊥}}{P_{⊥}^{0}} 〉 + \frac{2}{3} 〈 \frac{P_{| |}}{P_{| |}^{0}} 〉 = \frac{3}{2} \sqrt{\frac{ϵ_{2}}{μ_{2}}} \frac{n_{1}^{2}}{n_{2}^{2}} \frac{n_{1}}{k_{1}^{2} ϵ_{1} r_{1}^{3}} \sum_{l} [l (l + 1) I_{l}^{(1)} + I_{l}^{(2)} + I_{l}^{(3)}] .

7.2. Example 2: A shell

If N = 2, i.e., a microsphere coated by a single layer, then using Eqs. (23), (25), and (26), we have the following T(3, 1) matrix:

T = M_{3}^{- 1} (r_{2}) M_{2} (r_{2}) M_{2}^{- 1} (r_{1}) M_{1} (r_{1}) = \frac{i k_{3}}{n_{2} μ_{2}} \frac{i k_{2}}{n_{1} μ_{1}} G (3, 2) G (2, 1) .

One can find the TM and TE resonances of the structure by setting T₁₁ = 0 and T₃₃ = 0, respectively, from Eqs. (36) and (37)

T_{11} = \frac{i k_{3}}{n_{2} μ_{2}} \frac{i k_{2}}{n_{1} μ_{1}} [G_{11} (3, 2)] G_{11} (2, 1) + G_{12} (3, 2) G_{21} (2, 1)] = 0

T_{33} = \frac{i k_{3}}{n_{2} μ_{2}} \frac{i k_{2}}{n_{1} μ_{1}} [G_{33} (3, 2)] G_{33} (2, 1) + G_{34} (3, 2) G_{43} (2, 1)] = 0

which results in the following resonances

\frac{μ_{2} n_{3} χ_{l} (k_{3} r_{2})}{n_{2} μ_{3} χ_{l}^{'} (k_{3} r_{2})} = \frac{\frac{C}{D} ψ_{l} (k_{2} r_{2}) + χ_{l} (k_{2} r_{2})}{\frac{C}{D} ψ_{l}^{'} (k_{2} r_{2}) + χ_{l}^{'} (k_{2} r_{2})}} (TM)

and

\frac{n_{2} μ_{3} χ_{l} (k_{3} r_{2})}{n_{3} μ_{2} χ_{l}^{'} (k_{3} r_{2})} = \frac{\frac{E}{F} ψ_{l} (k_{2} r_{2}) + χ_{l} (k_{2} r_{2})}{\frac{E}{F} ψ_{l}^{'} (k_{2} r_{2}) + χ_{l}^{'} (k_{2} r_{2})} (TE)

where

\frac{C}{D} = \frac{[n_{2} μ_{1} ψ_{l}^{'} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - n_{1} μ_{2} ψ_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1})]}{[n_{1} μ_{2} ψ_{l}^{'} (k_{2} r_{1}) ψ_{l} (k_{1} r_{1}) - n_{2} μ_{1} ψ_{l} (k_{2} r_{1}) ψ_{l}^{'} (k_{1} r_{1})]}

\frac{E}{F} = \frac{[n_{1} μ_{2} ψ_{l}^{'} (k_{1} r_{1}) χ_{l} (k_{2} r_{1}) - n_{2} μ_{1} ψ_{l} (k_{1} r_{1}) χ_{l}^{'} (k_{2} r_{1})]}{[n_{2} μ_{1} ψ_{l}^{'} (k_{2} r_{1}) ψ_{l} (k_{1} r_{1}) - n_{1} μ_{2} ψ_{l} (k_{2} r_{1}) ψ_{l}^{'} (k_{1} r_{1})]} .

Eqs. (107) and (108) exactly match Eq. (7) in [26] and Eq. (10) in [24], respectively.

A numerical comparison is shown in Fig. 6. The multilayer model for a single layer converges to the microsphere case for a vanishingly small layer coating, or a vanishingly small internal sphere size. It was found that both the limits converge to within numerical precision.

Fig. 6 Demonstration that the multilayer model reproduces the microsphere model results, for D = 6 µm, n₁ = 1.59, and n₂ = 1.33. (a) Spectra for single-dipole excitation in both tangential and radial orientations. (b) Spectra for a uniform distribution of dipoles are also shown.

Download Full Size | PDF

Funding

Australian Research Council (ARC) (FL130100044, CE140100003).

References and links

1. H. Chew, P. J. McNulty, and M. Kerker, “Model for Raman and fluorescent scattering by molecules embedded in small particles,” Phys. Rev. A 13, 396–404 (1976). [CrossRef]

2. A. Moroz, “A recursive transfer-matrix solution for a dipole radiating inside and outside a stratified sphere,” Annals of Physics 315, 352–418 (2005). [CrossRef]

3. F. Vollmer, D. Braun, A. Libchaber, M. Khoshsima, I. Teraoka, and S. Arnold, “Protein detection by optical shift of a resonant microcavity,” Appl. Phys. Lett. 80, 4057–4059 (2002). [CrossRef]

4. S. Arnold, M. Khoshsima, I. Teraoka, S. Holler, and F. Vollmer, “Shift of whispering-gallery modes in microspheres by protein adsorption,” Opt. Lett. 28, 272–274 (2003). [CrossRef] [PubMed]

5. H. Quan and Z. Guo, “Simulation of whispering-gallery-mode resonance shifts for optical miniature biosensors,” J. Quant. Spectrosc. Radiat. Transfer 93, 231–243 (2005). [CrossRef]

6. Z. Guo, H. Quan, and S. Pau, “Near-field gap effects on small microcavity whispering-gallery mode resonators,” J. Phys. D: Appl. Phys. 39, 5133 (2006). [CrossRef]

7. K. Imakita, H. Shibata, M. Fujii, and S. Hayashi, “Numerical analysis on the feasibility of a multi-layered dielectric sphere as a three-dimensional photonic crystal,” Opt. Express 21, 10651–10658 (2013). [CrossRef] [PubMed]

8. R. Chance, A. Prock, and R. Silbey, “Molecular fluorescence and energy transfer near interfaces,” Adv. Chem. Phys. 37, 65 (1978).

9. H. Chew, “Radiation and lifetimes of atoms inside dielectric particles,” Phys. Rev. A 38, 3410–3416 (1988). [CrossRef]

10. K. Kolmakov, V. N. Belov, J. Bierwagen, C. Ringemann, V. Muller, C. Eggeling, and S. W. Hell, “Red-emitting rhodamine dyes for fluorescence microscopy and nanoscopy,” Chem. Eur. J. 16, 158–166 (2010). [CrossRef]

11. M. Li, X. Wu, L. Liu, and L. Xu, “Kerr parametric oscillations and frequency comb generation from dispersion compensated silica micro-bubble resonators,” Opt. Express 21, 16908–16913 (2013). [CrossRef] [PubMed]

12. D. Farnesi, A. Barucci, G. C. Righini, G. N. Conti, and S. Soria, “Generation of hyper-parametric oscillations in silica microbubbles,” Opt. Lett. 40, 4508–4511 (2015). [CrossRef] [PubMed]

13. N. Riesen, S. Afshar V., A. François, and T. M. Monro, “Material candidates for optical frequency comb generation in microspheres,” Opt. Express 23, 14784–14795 (2015). [CrossRef] [PubMed]

14. G. Kozyreff, J. Dominguez-Juarez, and J. Martorell, “Nonlinear optics in spheres: from second harmonic scattering to quasi-phase matched generation in whispering gallery modes,” Laser Photon. Rev. 5, 737–749 (2011). [CrossRef]

15. D. Farnesi, A. Barucci, G. C. Righini, S. Berneschi, S. Soria, and G. Nunzi Conti, “Optical frequency conversion in silica-whispering-gallery-mode microspherical resonators,” Phys. Rev. Lett. 112, 093901 (2014). [CrossRef] [PubMed]

16. W. Liang, A. A. Savchenkov, Z. Xie, J. F. McMillan, J. Burkhart, V. S. Ilchenko, C. W. Wong, A. B. Matsko, and L. Maleki, “Miniature multioctave light source based on a monolithic microcavity,” Optica 2, 40–47 (2015). [CrossRef]

17. Z. Ruan and S. Fan, “Superscattering of light from subwavelength nanostructures,” Phys. Rev. Lett. 105, 013901 (2010). [CrossRef] [PubMed]

18. Y. He, H. Deng, X. Jiao, S. He, J. Gao, and X. Yang, “Infrared perfect absorber based on nanowire metamaterial cavities,” Opt. Lett. 38, 1179–1181 (2013). [CrossRef] [PubMed]

19. C. Qu and E. C. Kinzel, “Polycrystalline metasurface perfect absorbers fabricated using microsphere photolithography,” Opt. Lett. 41, 3399–3402 (2016). [CrossRef] [PubMed]

20. H. Chew, “Transition rates of atoms near spherical surfaces,” J. Chem. Phys. 87, 1355–1360 (1987). [CrossRef]

21. Y. Xu, W. Liang, A. Yariv, J. G. Fleming, and S.-Y. Lin, “High-quality-factor bragg onion resonators with omnidirectional reflector cladding,” Opt. Lett. 28, 2144–2146 (2003). [CrossRef] [PubMed]

22. Y. Xu, W. Liang, A. Yariv, J. G. Fleming, and S.-Y. Lin, “Modal analysis of bragg onion resonators,” Opt. Lett. 29, 424–426 (2004). [CrossRef] [PubMed]

23. W. Liang, Y. Xu, Y. Huang, A. Yariv, J. Fleming, and S.-Y. Lin, “Mie scattering analysis of spherical Bragg “onion” resonators,” Opt. Express 12, 657–669 (2004). [CrossRef] [PubMed]

24. I. Teraoka and S. Arnold, “Enhancing the sensitivity of a whispering-gallery mode microsphere sensor by a high-refractive-index surface layer,” J. Opt. Soc. Am. B 23, 1434–1441 (2006). [CrossRef]

25. I. Teraoka and S. Arnold, “Theory of resonance shifts in TE and TM whispering gallery modes by nonradial perturbations for sensing applications,” J. Opt. Soc. Am. B 23, 1381–1389 (2006). [CrossRef]

26. I. Teraoka and S. Arnold, “Whispering-gallery modes in a microsphere coated with a high-refractive index layer: polarization-dependent sensitivity enhancement of the resonance-shift sensor and TE–TM resonance matching,” J. Opt. Soc. Am. B 24, 653–659 (2007). [CrossRef]

27. C. F. Bohren and D. R. Huffman, Absorption and Scattering by a Sphere (Wiley-VCH Verlag GmbH, 2007), pp. 82–129.

28. J. J. Kaufman, G. Tao, S. Shabahang, E.-H. Banaei, D. S. Deng, X. Liang, S. G. Johnson, Y. Fink, and A. F. Abouraddy, “Structured spheres generated by an in-fibre fluid instability,” Nature 487, 463–467 (2012). [CrossRef] [PubMed]

29. G. Tao, J. J. Kaufman, S. Shabahang, R. RezvaniNaraghi, S. V. Sukhov, J. D. Joannopoulos, Y. Fink, A. Dogariu, and A. F. Abouraddy, “Digital design of multimaterial photonic particles,” Proc. Nat. Acad. Sci. USA 113, 6839–6844 (2016). [CrossRef] [PubMed]

30. V. N. Astratov, J. P. Franchak, and S. P. Ashili, “Optical coupling and transport phenomena in chains of spherical dielectric microresonators with size disorder,” Appl. Phys. Lett. 85, 5508–5510 (2004). [CrossRef]

31. D. K. Armani, T. J. Kippenberg, S. M. Spillane, and K. J. Vahala, “Ultra-high-Q toroid microcavity on a chip,” Nature 421, 925–928 (2003). [CrossRef] [PubMed]

32. D. Xiao-Wei, L. Shao-Hua, F. Su-Chun, X. Ou, and J. Shui-Sheng, “All-fibre micro-ring resonator based on tapered microfibre,” Chin. Phys. B 17, 1029 (2008). [CrossRef]

33. F. Monifi, S. Odemir, J. Friedlein, and L. Yang, “Encapsulation of a fiber taper coupled microtoroid resonator in a polymer matrix,” IEEE Photon. Technol. Lett. 25, 1458–1461 (2013). [CrossRef]

34. N. Riesen, T. Reynolds, A. François, M. R. Henderson, and T. M. Monro, “Q-factor limits for far-field detection of whispering gallery modes in active microspheres,” Opt. Express 23, 28896–28904 (2015). [CrossRef] [PubMed]

35. M. L. Gorodetsky and V. S. Ilchenko, “Optical microsphere resonators: optimal coupling to high-Q whispering-gallery modes,” J. Opt. Soc. Am. B 16, 147–154 (1999). [CrossRef]

36. E. Nuhiji, F. G. Amar, H. Wang, N. Byrne, T.-L. Nguyen, and T. Lin, “Whispering gallery mode emission generated in tunable quantum dot doped glycerol/water and ionic liquid/water microdroplets formed on a superhydrophobic coating,” J. Mater. Chem. 21, 10823–10828 (2011). [CrossRef]

37. J. Geng, R. W. Ziolkowski, R. Jin, and X. Liang, “Numerical study of the near-field and far-field properties of active open cylindrical coated nanoparticle antennas,” IEEE Photon. J. 3, 1093–1110 (2011). [CrossRef]

38. N. V. Kryzhanovskaya, M. V. Maximov, and A. E. Zhukov, “Whispering-gallery mode microcavity quantum-dot lasers,” IEEE J. Quant. Electron. 44, 189–200 (2014). [CrossRef]

39. A. François, T. Reynolds, and T. M. Monro, “A fiber-tip label-free biological sensing platform: A practical approach toward in-vivo sensing,” Sensors 15, 1168–1181 (2015). [CrossRef] [PubMed]

40. T. C. Preston and J. P. Reid, “Determining the size and refractive index of microspheres using the mode assignments from Mie resonances,” J. Opt. Soc. Am. A 32, 2210–2217 (2015). [CrossRef]

41. Y. Yang, J. Ward, and S. N. Chormaic, “Quasi-droplet microbubbles for high resolution sensing applications,” Opt. Express 22, 6881–6898 (2014). [CrossRef] [PubMed]

42. J. M. M. Hall, S. Afshar V., M. R. Henderson, A. François, T. Reynolds, N. Riesen, and T. M. Monro, “Method for predicting whispering gallery mode spectra of spherical microresonators,” Opt. Express 23, 9924–9937 (2015). [CrossRef] [PubMed]

43. T. Reynolds, M. R. Henderson, A. François, N. Riesen, J. M. M. Hall, S. Afshar V., S. J. Nicholls, and T. M. Monro, “Optimization of whispering gallery resonator design for biosensing applications,” Opt. Express 23, 17067–17076 (2015). [CrossRef] [PubMed]

44. J. M. M. Hall, A. François, S. Afshar V., N. Riesen, M. R. Henderson, T. Reynolds, and T. M. Monro, “Determining the geometric parameters of microbubble resonators from their spectra,” J. Opt. Soc. Am. B 34(1), 2699–2706 (2017). [CrossRef]

45. H. Chew, M. Kerker, and P. J. McNulty, “Raman and fluorescent scattering by molecules embedded in concentric spheres,” J. Opt. Soc. Am. 66, 440–444 (1976). [CrossRef]

46. W. Yang, “Improved recursive algorithm for light scattering by a multilayered sphere,” Appl. Opt. 42, 1710–1720 (2003). [CrossRef] [PubMed]

47. A. François, K. J. Rowland, S. Afshar V., M. R. Henderson, and T. M. Monro, “Enhancing the radiation efficiency of dye doped whispering gallery mode microresonators,” Opt. Express 21, 22566–22577 (2013). [CrossRef] [PubMed]

48. A. François, K. J. Rowland, and T. M. Monro, “Highly efficient excitation and detection of whispering gallery modes in a dye-doped microsphere using a microstructured optical fiber,” Appl. Phys. Lett. 99, 141111 (2011). [CrossRef]

49. J. Blatt and V. Weisskopf, Theoretical Nuclear Physics, Dover Books on Physics Series (Dover Publications, 1991).

50. R. G. Barrera, G. A. Estevez, and J. Giraldo, “Vector spherical harmonics and their application to magnetostatics,” Eur. J. Phys. 6, 287 (1985). [CrossRef]

51. J. Jackson, Classical electrodynamics (Wiley, 1975).

52. I. H. Malitson, “Interspecimen comparison of the refractive index of fused silica,” J. Opt. Soc. Am. 55, 1205–1209 (1965). [CrossRef]

53. S. N. Kasarova, N. G. Sultanova, C. D. Ivanov, and I. D. Nikolov, “Analysis of the dispersion of optical plastic materials,” Opt. Mater. 29, 1481–1490 (2007). [CrossRef]

54. A. Edmonds, Angular Momentum in Quantum Mechanics, Investigations in Physics Series (Princeton University, 1996).

Unified theory of whispering gallery multilayer microspheres with single dipole or active layer sources

Abstract

1. Introduction

2. Theory

2.1. Geometry

2.2. Conventions

2.3. Transfer matrix method

2.4. Structure resonances

2.5. Scattered power in the outer region

2.6. A dipole in one layer

2.7. One active layer

3. Demonstration and discussion

4. Conclusion

5. Appendix A: Properties of vector spherical harmonics (VSH)

6. Appendix B: Relation between conventions

7. Appendix C: Examples

7.1. Example 1: A microsphere

7.1.1. Special case: A dipole in the outer region

7.1.2. Special case: A dipole in the inner region

7.1.3. Special case: A microsphere with an active inner region

7.2. Example 2: A shell

Funding

References and links

Cited By

Figures (6)

Equations (110)

Optics Express