Open Access
Add to BibSonomyAbstract
In this paper, we developed an efficient method for searching the resonant eigenfrequency of dielectric optical microcavities by the boundary element method. By transforming the boundary integral equation to a general eigenvalue problem for arbitrary, symmetric, and multi-domain shaped optical microcavities, we analyzed the regular motion of the eigenvalues against the frequency. The new strategy can predict multiple resonances, increase the speed of convergence, and avoid non-physical spurious solutions. These advantages greatly reduce the computation time in the search process of the resonances. Moreover, this method is not only valuable for dielectric microcavities, but is also suitable for other photonic systems with dissipations, whose resonant eigenfrequencies are complex numbers.
© 2011 Optical Society of America
1. Introduction
Dielectric optical microcavities have recently attracted significant interests, with wide applications in microlasers, sensors, and quantum optics [1]. Furthermore they provide a test bed for quantum and wave chaos [2]. To study the electromagnetic properties of dielectric microcavities, one of the kernel tasks is to solve Maxwell equations, or more specifically to find the quasi-normal resonant modes. For an irregular shaped boundary of such a cavity, including multi-domain cavities, numerical methods, such as the Boundary Element Method (BEM) are the only way to calculate the resonances. The BEM utilizes the Green’s function approach, which reduces the two-dimensional Helmholtz equation to a one-dimensional boundary integral on the boundary. The BEM has been extensively applied in photonics for the theoretical analysis of high quality (Q) factor resonances of deformed microcavities [3–7], the solution of the propagation modes in waveguides [8,9], and the Goos-Hänchen shift of the Gaussian beam reflection at a non-regular boundary [10].
Dielectric optical microcavities are known to be an open system due to their energy dissipation at infinity. Thus, their resonant eigen-frequencies are complex numbers (kR), where the finite quality factor of the resonance is defined as Q = −Re(kR)/2Im(kR), which directly corresponds to the photon lifetime τ = Q/kc in the cavity, where c is the speed of light in vacuum. This means that searching for the resonances take places in the two-dimensional complex plane. Newton’s method is usually used in order to search for the minimum of the determinant of discretized boundary integration matrix [3]. In this method, the initial value (kR0) is very important as kR only converges to the nearest root. In addition, there appear non-physical, so called spurious solutions, which can not be distinguished from high-Q modes with the Newton’s method [3].
In this paper, we present a more efficient approach of finding the resonant eigenfrequencies of an arbitrarily shaped, multi-domain dielectric resonator with the BEM. Firstly, we rewrite the matrix equation into a generalized eigenvalues (GE) problem, similar to Refs. [11–13]. Secondly, we analyze the behavior of the eigenvalues in the complex plane as a function of a small change in the complex wavenumber kR. Thirdly, we use this analytic behavior of the eigenvalues to predict multiple resonant positions near the starting point, and thus reduce the convergence time significantly. Fourthly, spurious solutions can be easily avoided during the search process. This improved BEM root search strategy can be further extended for various applications, such as resonances of photonic molecules, propagation modes in photonic crystal fibers [8, 9], and constant flux modes for microlaser [14, 15].
2. Brief Introduction to the Boundary Element Method
Maxwell equations for two-dimensional dielectric microcavities can be reduced to the Helmholtz equation
where n is the refractive index of cavity, k = 2π/wavelength is the wavenumber, and Ψ is the wavefunction, which can be either the electric or magnetic component of the electromagnetic field. Here, we use the dimensionless wave number kR for convenience, where R is the size of the photonic structure. The corresponding Green’s function is where , and is the zeroth-order Hankel function of the first kind. By employing the Green’s function, the two dimensional Helmholtz equation can be reduced to an integral on the boundary. For multiple regions Ωj filled with a homogeneous dielectric, and bounded by a smooth boundary Γj, the boundary integral equations (BIE) can be written as [3] where B(s′,s) = −2G(s, s′,kR), C(s′,s) = 2∂vG(s, s′,kR) −δ(s − s′), and ϕ(s) = ∂vψ(s) with ∂v being the normal derivative. By dividing the boundary into small elements, the entire set of BIEs can be written in the matrix form [3] where Bj and Cj represent the integral operators in the region Ωj. If the number of elements is l, then M is a 2l × 2l matrix, ϕ and ψ are l-component vectors.The problem of finding the resonant eigen-frequencies reduces thus to find the complex numbers kR to satisfy
In Fig. 1(a), we plot the distribution of the determinant of M(kR) in the two dimensional complex plane. There are several local minima in this figure made visible through the equipotential lines, corresponding to possible eigen-frequency of modes. The locations of the minima are not regularly distributed in the complex plane. Since they are not easily predictable, the root search in the two dimensional space is complex and usually time consuming. A standard numerical procedure for extracting the zeros of this determinant is the Newton’s method [3]. Given a starting point kR0 and step size Δ, it compares det[M(kR0)], det[M(kR0 + Δ)], and det[M(kR0 + Δi)], then approximately estimates the values of kR to satisfy det[M(kR)] = 0. The step is repeated until |kR − kR0| ≤ δ, where δ is the tolerance precision of kR. The drawbacks of this method are obvious: (1) The root strongly depends on the starting point kR0. When there are two modes with frequencies close to each other, it is difficult to distinguish and find both of them. (2) Only one root can be found with one given kR0. (3) It lacks the capability to exclude spurious solutions. For example, as shown with the red circle in Fig. 1(a), several minima near the real axis Im(kR) = 0 are spurious solutions to the interior Dirichlet problem. They are very hard to be distinguished from high-Q modes which also have Im(kR) ≈ 0.
Fig. 1 (a) The distribution of the determinant det[M(kR)] for a circular cavity (n = 1.45). Red and white circles indicate two minima, corresponding to a spurious and physical resonance, respectively. (b) The motion of the generalized eigenvalues λ for the same cavity, with different real parts kR = 14.5+{0.0, 0.02,...,0.50} (blue circles), and different imagine parts kR = 14.5 – {0.0i, 0.02i,...,0.70i} (red circles).
The BEM method is also commonly used to calculate the resonances in billiards, where the eigenfrequencies are real and correspond to the Dirichlet boundary condition. The singular value decomposition (SVD) [16] and eigenvalue decomposition (EVD) [12] can be applied to search resonances, which are very efficient for finding multiple resonances in a single iteration. However, the SVD is not easily adaptable for the dielectric microcavities, as the eigenfrequencies are complex. An adaptation of the EVD for complex eigenfrequencies is presented in this publication.
3. Root Search Strategy
In contrast to the usual procedure based on Newton’s method, we can rearrange Eq. (5) as done in reference [13]
where M′(kR) = M(kR) – N(kR) and u = {ϕ,ψ}T is a 2l-component vector of boundary values. Solving for the eigenvalues λ, the solutions to Eq. (2) require λ(kR) = 1 + 0i, known as the quantization condition [11].Equation (6) is a standard generalized eigenvalue (GE) problem, for which there exist extremely stable and accurate routines. In the Newton’s method, the det[M(kR)] is calculated as a single valued function of kR, and only one resonance can be predicted. However, solving Eq. (6) we can find 2l complex eigenvalues for each kR, and extract more information about the BIE. This extra information can be used to find the resonances more efficiently. By carefully choosing the form of the GE equation, we will see a regular behavior of the eigenvalues λ with respect to a change of kR. Figure 1(b) shows such a motion of λ with respect to a change in kR in the complex plane for the general non-symmetric GE method, which is described in detail in the following. A change in Re[kR] induces a rotation of the eigenvalues, while a change of Im[kR] induces a radial motion. The eigenvalues can be traced by an overlap integral with the corresponding eigenvectors, similar to the description in the improved multipole method and scattering quantization approach [11, 13], but the motion in our case is more complex as the circular traces are not centered around the origin (0, 0). If we can find an analytic function for the motion of λ with respect to kR, then the resonance with λ = 1 becomes predictable. This would allow us to derive an algorithm for the calculation of the resonant states. In addition, we can also distinguish the physical and non-physical solutions through their dynamics in λ.
In the following, we will derive the form of the matrix N from Eq. (6) for dielectric microcavities, including the general, non-symmetric, and the symmetric shape respectively.
3.1. General, Non-Symmetric Condition
For a single circular cavity, the field at the boundary can be represented by Fourier expansion as ϕ = ∑mϕmexp{−imφ}, ψ = ∑mψmexp{−imφ}. By using Graf’s addition theorems [17], we can expand the Hankel function as
and employ the Wronskian identity of integer order Bessel functions [18] to rewrite the boundary integrals in a matrix form Using the asymptotic formulation of the Bessel and the Hankel function, we finally arrive at [11, 19] with , Φ = m(α – tanhα), , and . In comparison with the boundary integral equation of Dirichlet boundary condition in billiards [12], we have a more complex matrix. From the approximation e−2Φ ≈ 0 for m < nkR, Eq. (10) can be simplified as In the matrix, only terms tanβ, tanhα and e−2θ vary with kR. Therefore, for any matrix N, the eigenvalues (λ) of the generalized linear eigenvalues equation Eq. (6) should be a function of kR, which can be expressed by these terms. Among these terms, the change of e−2θ with respect to the change of kR is much faster than others. So, we can expect that the dynamics of λ is predominantly determined by the e−2θ. If we can express the λ by a simple linear function of e−2θ, the motion of λ can be very regular and predictable. After several tries, we found that when the eigenvalue can be written as where I is the unit matrix. From Eq. (13), the trace of λ against kR is a circle centered at ξ, with radius |ρ|, and the speed of motion with respect to kR is κ, which is consistent with a regular motion of λ in Fig. 1(b).3.2. Symmetric Condition
When the cavity boundary shape has a symmetry, for example, it is x-axis symmetric, we can reduce the BIEs to only one half of the boundary. The corresponding size of the matrix will be reduced to only 1/4 of the non-symmetric condition. As a result, the computation complexity is reduced. Therefore, the boundary element matrix M(kR) can be reduced to
where M1 and M2 are l × l matrices, and ũ is a l-component vector of boundary values. sx = 1 (−1) is the symmetric condition corresponding to even (odd) symmetry.By exploring the symmetry, the field at the boundary of a single circular cavity can be represented as ϕ = ∑mϕmcos(mφ), ψ = ∑mψmcos(mφ) for sx = 1, and ϕ = ∑mϕmsin(mφ), ψ = ∑mψmsin(mφ) for sx = −1. Similar to the general case, we can rewrite the boundary integrals in a matrix form
where χ = (sx − 1)π/4. For integer m, the Bessel and Hankel functions satisfy J−m(z) = (−1)mJm(z) and . Then, the matrix for symmetric condition can be simplified, we finally obtain Eq. (9). Similar to the general condition, the regular motion of eigenvalues deduces toThus, for both the symmetric and the non-symmetric condition, we can obtain values for ξ, ρ and κ by repeatedly solving the eigenvalues three times at a starting point kR0,
where λ0 = λ(kR0) and λ± = λ(kR0 ± Δ). Therefore, we can predict the motion of eigenvalues as4. Application
In this section, we apply the GE method to differently shaped dielectric microcavities. Following the above described root search strategy, we choose a starting point kR0, and then solve for the parameters {ξ,ρ,κ}. For a 2l × 2l matrix, there are 2l eigenvalues that can be obtained. Thus, we can get 2l sets of parameters {ξj,ρj,κj}, with j = 1,...,2l. We then can predict 2l possible resonant frequencies as
As a first example, we considered a circular silica microdisk, with refractive index n = 1.45. Figure 2(a) shows the non-symmetric GE method for a number of eigenvalues λ (red circles) with kR = {99.5 – 0.4i,...,100.5 – 0.4i}, and shows the motion with respect to kR. Through solving Eq. (17)–(19), at the given kR0 = 99.5 – 0.4i, we can find the parameters {ξ,ρ,κ} for all λ, and then predict the motion of them by Eq. (13), shown in the blue lines. The predictions are very consistent to the actual motion of eigenvalues.Fig. 2 The motion of eigenvalues λ of a two-dimensional circular microdisk cavity with n = 1.45 for the non-symmetric condition (a) and symmetric condition (c). The circles are the eigenvalues with kR = 99.5–0.4i,...,100.5–0.4i, and the blue lines correspond to the prediction retrieved at kR0 = 99.5 – 0.4i. (b) and (d) show the predicted roots (blue crosses) at starting point kR0 = 100.0 – 0.4i (black dot) and compare them to the analytical results (red circles), for the non-symmetric and symmetric conditions, respectively.
At the starting point kR0 = 100.0 – 0.4i, we can predict the resonances by Eq. (21), shown with blue crosses in Fig. 2(b). Comparing to the exact analytical solution to circular cavities [20], the GE method predicts all resonances around kR0 (black solid dot). However, for some eigenvalues far away from the start point kR0, the predicted roots do not agree well with the exact roots. The reason for this, as analyzed above in Eq. (13), is that the regular motion function is just an approximation, which only works in a small area around kR0. When |kR – kR0| is large, the approximation can lead to large errors. We can define a working distance L = |kRp – kR0|, where the area with radius L around the starting point is within the predictable region. From Fig. 2(b), we can estimate L ≈ 1.0, resulting in a valid prediction of about 100 resonances.
When exploring the symmetric properties of circular microdisks, we can apply the symmetric GE method. In Fig. 2(c) we plot the motion of eigenvalues and their predictions. Comparing with non-symmetric GE method, the results of symmetric GE method are almost the same, the working distance is L ≈ 1.0 in Fig. 2(d).
The GE method shows its power especially when applied to non-circular cavities or multiple coupled cavities, where no analytical solution exists. Different to the circular cavity where the cylindrical harmonics {ϕmeimφ} and {ψmeimφ} are orthogonal, the field in a deformed cavity can be expanded in a superposition of multiple components of different m. According to Eq. (13), we can approximately express λ(kR) ≈ ξ + ∑m ρmekR×κm. For small deformations, the dynamics of lambda is predominantly determined by a single or only a few components. And for different m close to each other, the corresponding κm are nearly the same. Therefore, the approximate function of λ, Eq. (13) is still adequate to describe the regular eigenvalue motion for deformed cavities. When increasing the deformation, more components of m get involved and give rise to much more complicated dynamics and the motion of λ is less regular, thus the work distance L is reduced. Our numerical results indicated that the regular rotation and radial motion of λ is maintained even for squares and full chaotic cavities. As an example, we use GE method for two coupled stadium cavities with aspect ratio 2:1 [5]. In Fig. 3(a), we predict the resonances by the non-symmetric GE method, which agrees well with the exact roots found via the BEM with tolerance precision δ = 10−6, albeit with a reduced working distance of about L < 0.1, resulting in a prediction of 20 resonances. One of the resonances in the coupled cavities is plotted in Fig. 3(b).
Fig. 3 (a) The predicted roots of two coupled stadium cavities (blue crosses) at start point kR0 = 20.025 – 0.07i (black dot), compared to the exact resonant frequencies (calculated via BEM) (red circles). (b) False color representation of the electromagnetic field intensity of a resonance in the two coupled stadium cavities. The refractive index of the cavity is n = 2.0 and the aspect ratio is 2 : 1. The centers of the two cavity are at (3.0, 2.0) and (4.0, 4.1).
5. Discussion
Because the working distance L varies for different cavity boundary shapes, we must know the value of L before applying the GE method to search roots. A simple way to estimate the L is to solve the GE at two separate start points kRa and kRb, and obtain two arrays of predicted roots {kRa1, kRa2,...,kRa2l} and {kRb1, kRb2,...,kRb2l}. Since the predictable region are circles of radius L with center at starting points, we can estimate the unknown L from the overlap of the two arrays. Plot them in a graph and find the correspondence between the predicted roots of the two different start points, similar to Fig. 2(b). The distance between the two starting points should be smaller than L, for example, we can take |kRb – kRa| = 0.05. Through this method, we can estimate the L quickly by solve GE only twice without knowing the exact roots.
In addition, it is worth to note that the GE method only improves the efficiency of the root searching process, not the precision of BEM. As discussed in Ref. [4], the error of the BEM is inevitable due to the finite discrete of the boundary, the precision of BEM can only be improved by increasing the number of elements. However, the L does not increase by increasing of the precision of BEM. The limited L of the GE method does not come from the computation errors, which are intrinsic to the irregular motion of the eigenvalues.
Although the prediction area shrinks for a non-regular cavity, it is still surprising that GE method can predict multiple resonances by one single iteration, thus reducing the time for the roots searching. Taking the two coupled stadium cavities as an example, if we want to find all modes in the regime 19.5 < Re(kR) < 20.5 and −1 < Im(kR) < 0 with BEM, we have to solve all determinants in order to produce a figure like Fig. 1(a). In order to find all resonances, we set the step size as Δ = 0.01 to get sufficient resolution, which means that we need to solve Eq. (5) 10,000 times. In contrast, in the GE method, we can set the step size to Δ = 0.2 to find all resonances, thus we only need to solve the GE problem for 3×25 = 75 times. In the numerical calculation, the computation time for solving a determinant is similar to that of solving the GE problem. Therefore, the GE method is more than two orders of magnitude faster than the previous method which searches the minima of det(M). In addition, when two modes have eigenfrequencies near to each other, with a distance smaller than Δ, it is difficult find both of them by the Newton’s method.
Due to the predictable motion of λ, we can also expect a much faster convergence speed of the root search. In Fig. 4, we compare the Newton’s method and GE method for various cavities. All results show a much better convergence of the GE method. Noting that the starting point of the GE mode is placed 0.05 away from the exact resonance, whereas the starting point for the Newton method is 0.005 away from the exact resonance. The Newton method is limited to the convergence to the root which is closest to the starting point, whereas the GE mode can sustain a starting point further away within the working distance L.
Fig. 4 Comparison of the convergence between the Newton (solid line) and GE (dashed line) method, for a single disk cavity, a single stadium cavities and two coupled stadium cavities. “G” and “S” stand for general, non-symmetric and symmetric condition, respectively.
In addition to the predictability and fast convergence, the GE method can also avoid the non-physical, spurious roots. Because the spurious roots moves much faster than physics roots, we can set a threshold of κ to exclude the spurious roots. In Fig. 1(a), there are some spurious minima in the determinant distribution near the real axis. However, these spurious solutions are avoided in the GE method, as we can see from Fig. 2(b). There are some eigenvalues with very small R0, corresponding to resonances very far away from the starting point, a threshold is also needed to exclude them.
6. Conclusion
In summary, we have presented an efficient numerical algorithm for the boundary element method in order to find resonances in multiple arbitrarily shaped dielectric microcavities. We analyzed the motion of the eigenvalues of the generalized eigenvalue problem in both symmetrically and non-symmetrically shaped microcavities, and demonstrated the capability to predict the positions of multiple resonances. In addition, the new strategy can avoid the non-physical spurious solutions and increase the convergence speed. This method can also be valuable for studying the resonances of lossy photonic structures, such as microcavities, photonic crystal fibers, surface plasmon resonances, and can further be extended to acoustical systems.
We thank Hakan Türeci for fruitful discussions. The work was supported by the National Fundamental Research Program of China under Grant No. 2011CB921200, the Knowledge Innovation Project of Chinese Academy of Sciences, National Natural Science Foundation of China under Grant No. 11004184, and Fundamental Research Funds for the Central Universities. C.-L. Zou acknowledges the support by the scholarship award for excellent doctoral students.
References and links
1. K. Vahala, “Optical Microavities,” Nature 424, 839–845 (2004). [CrossRef]
2. J. U. Nöckel and A. D. Stone, “Ray and wave chaos in asymmetric resonant optical cavities,” Nature 385, 45–47 (1997). [CrossRef]
3. J. Wiersig, “Boundary element method for resonances in dielectric microcavities,” J. Opt. A: Pure Appl. Opt. 5, 53–60 (2003). [CrossRef]
4. C.-L. Zou, Y. Yang, Y.-F. Xiao, C.-H. Dong, Z.-F. Han, and G.-C. Guo, “Accurately calculating high quality factor of whispering-gallery modes with boundary element method,” J. Opt. Soc. Am. B 26, 2050–2053 (2009). [CrossRef]
5. S.-Y. Lee, M. S. Kurdoglyan, S. Rim, and C.-M. Kim, “Resonance patterns in a stadium-shaped microcavity,” Phys. Rev. A 70, 023809 (2004). [CrossRef]
6. J. Wiersig and M. Hentschel, “Combining directional light output and ultralow loss in deformed microdisks,” Phys. Rev. Lett. 100, 033901 (2008). [CrossRef] [PubMed]
7. C.-L. Zou, F.-W. Sun, C.-H. Dong, X.-W. Wu, J.-M. Cui, Y. Yang, G.-C. Guo, and Z.-F. Han, “Mechanism of unidirectional emission of ultrahigh Q whispering gallery mode in microcavities,” http://arxiv.org/abs/0908.3531.
8. H. Cheng, W. Crutchfield, M. Doery, and L. Greengard, “Fast, accurate integral equation methods for the analysis of photonic crystal fibers I: theory,” Opt. Express 12, 3791–3805 (2004). [CrossRef] [PubMed]
9. E. Pone, A. Hassani, S. Lacroix, A. Kabashin, and M. Skorobogatiy, “Boundary integral method for the challenging problems in bandgap guiding, plasmonics and sensing,” Opt. Express 15, 10231–10246 (2007). [CrossRef] [PubMed]
10. L.-M. Zhou, C.-L. Zou, Z.-F. Han, G.-C. Guo, and F.-W. Sun, “Negative Goos-Hänchen shift on a concave dielectric interface,” Opt. Lett. 36, 624–626 (2011). [CrossRef] [PubMed]
11. H. E. Türeci, H. G. L. Schwefel, P. Jacquod, and A. D. Stone, “Modes of wave-chaotic dielectric resonators,” Prog. Opt. 47, 75–137 (2005). [CrossRef]
12. H. E. Türeci and H. G. L. Schwefel, “An efficient Fredholm method for the calculation of highly excited states of billiards,” J. Phys. A: Math. Theor. 40, 13869–13882 (2007). [CrossRef]
13. H. G. L. Schwefel and C. G. Poulton, “An improved method for calculating resonances of multiple dielectric disks arbitrarily positioned in the plane,” Opt. Express 17, 13178–13186 (2009). [CrossRef] [PubMed]
14. H. E. Türeci, A. D. Stone, and B. Collier, “Self-consistent multimode lasing theory for complex or random lasing media,” Phys. Rev. A 74, 043822 (2006). [CrossRef]
15. H. E. Türeci, L. Ge, S. Rotter, and A. D. Stone, “Strong interactions in multimode random lasers,” Science 320, 643–646 (2008). [CrossRef] [PubMed]
16. A. Bäcker, “Numerical aspects of eigenvalue and eigenfunction computations for chaotic quantum systems,” in Lecture Notes in Physics Vol. 618, S. Graffi and M. Degli Esposti, eds. (Spinger, 2003), pp. 91–144. [CrossRef]
17. J. H. Graf, “Über die Addition und Subtraction der Argumente bei Bessel’schen Functionen nebst einer Anwendung,” Math. Ann. 43, 136–144 (1893). [CrossRef]
18. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions (Applied Mathematics Series 55), (National Bureau of Standards1966), Chapter 9, pp. 360, Eq. (9.1.16).
19. H. E. Türeci, “Wave Chaos in Dielectric Resonators: Asymptotic and Numerical Approaches,” Ph.D. thesis (Yale University, New Haven, Connecticut, 2003).
20. J.-W. Ryu, S. Rin, Y.-J. Park, C.-M. Kim, and S. -Y. Lee, “Resonances in a circular dielectric cavity,” Phys. Lett. A 372, 3531–3536 (2008). [CrossRef]
