An improved method for calculating resonances of multiple dielectric disks arbitrarily positioned in the plane

We present a numerically improved multipole formulation fo r the calculation of resonances of multiple disks located at a rbi r ry positions in a 2-d plane, and suitable for the accurate computation of t he resonances of large numbers of disks and of high-wavenumber eigenstate s. Using a simple reformulation of the field expansions and boundary co nditions, we are able to transform the multipole formalism into a linea r igenvalue problem, for which fast and accurate methods are available. Observing that the motion of the eigenvalues in the complex plane is analyti c with respect to a two parameter family, we present a numerical algorithm t o compute a range of multiple-disk resonances and field distributions using only two diagonalizations. This method can be applied to photoni c molecules, photonic crystals, photonic crystal fibers, and random lase rs. © 2009 Optical Society of America OCIS codes: (140.3945) Microcavity; (230.5750) Resonators; (230.455 5) Coupled resonators; (000.3860) Mathematical methods in physics. References and links 1. B. Little, S. Chu, H. Haus, J. Foresi, and J. Laine, “Micror ing resonator channel dropping filters,” J. Lightwave Technol.15, 998–1005 (1997). 2. T. Carmon, T. Kippenberg, L. Yang, H. Rokhsari, S. Spillan e, and K. Vahala, “Feedback control of ultra-high-Q microcavities: application to micro-Raman lasers and micr opa ametric oscillators,” Opt. Express 13, 3558–3566 (2005). 3. S. Preu, H. G. L. Schwefel, S. Malzer, G. H. Döhler, L. J. Wa ng, M. Hanson, J. D. Zimmerman, and A. C. Gossard, “Coupled whispering gallery mode resonators in th e Terahertz frequency range,” Opt. Express 16, 7336–7343 (2008). 4. D. S. Wiersma and A. Lagendijk, “Light diffusion with gain and random lasers,” Phys. Rev. E 54, 4256–4265 (1996). 5. H. Cao, J. Y. Xu, D. Z. Zhang, S.-H. Chang, S. T. Ho, E. W. Seel ig, X. Liu, and R. P. H. Chang, “Spatial Confinement of Laser Light in Active Random Media,” Phys. Rev . Lett. 84, 5584–5587 (2000). 6. H. E. Türeci, L. Ge, S. Rotter, and A. D. Stone, “Strong int eractions in multimode random lasers,” Science 320, 643–646 (2008). 7. E. B. Becker, G. F. Carey, and J. T. Oden, Finite Elements (Pentice-Hall, Englewood Cliffs, N.J., 1981). 8. J. B. Davies, “Finite element analysis of waveguides and c avities – a review.” IEEE T. Magn. 29, 1578 (1993). 9. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Mead e, Photonic crystals: Molding the Flow of Light (Princeton University Press, Pinceton, 2008). 10. M. Fujii, C. Koos, C. Poulton, J. Leuthold, and W. Freude, “Nonlinear FDTD analysis and experimental verification of four-wave mixing in InGaAsP-InP racetrack microres onators,” IEEE Photon. Technol. Lett. 18, 361–363 (2006). #110647 $15.00 USD Received 29 Apr 2009; revised 30 Jun 2009; accepted 30 Jun 2009; published 17 Jul 2009 (C) 2009 OSA 20 July 2009 / Vol. 17, No. 15 / OPTICS EXPRESS 13178 11. J. Wiersig, “Boundary element method for resonances in d ielectric microcavities,” J. Opt. Soc. Am. A 5, 53–60 (2003).physics/0206018 . 12. S. V. Boriskina, P. Sewell, T. M. Benson, and A. I. Nosich, “Accurate simulation of two-dimensional optical microcavities with uniquely solvable boundary integral eq uations and trigonometric Galerkin discretization,” J. Opt. Soc. Am. A21, 393–402 (2004). 13. A. B. Movchan, N. V. Movchan, and C. G. Poulton, Asymptotics of dilute and densely packed composites (Imperial College Press, London, 2002). 14. H. E. Türeci, H. G. L. Schwefel, P. Jacquod, and A. D. Ston e, “Modes of wave-chaotic dielectric resonators,” Prog. Opt.47, 75–137 (2005). physics/0308016 . 15. W. T. Perrins, D. R. McKenzie, and R. C. McPhedran, “Trans port Properties of Regular Arrays of Cylinders,” P. Roy. Soc. Lond. A Mat. 369, 207–225 (1979). URLhttp://www.jstor.org/stable/2398611 . 16. B. T. Kuhlmey, T. P. White, G. Renversez, D. Maystre, L. C. Botten, C. M. de Sterke, and R. C. McPhedran, “Multipole method for microstructured optical fibers. II. I mplementation and results,” J. Opt. Soc. B 19, 2331– 2340 (2002). URLhttp://josab.osa.org/abstract.cfm?URI=josab-19-10-2 331. 17. A. Spence and C. Poulton, “Photonic band structure calcu lations using nonlinear eigenvalue techniques,” J. Comput. Phys.204, 65–81 (2005). 18. E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. D. Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK Users’ Guide, 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999). 19. H. E. Tureci and H. G. L. Schwefel, “An efficient Fredholm m ethod for the calculation of highly excited states of billiards,” J. Phys. A40, 13,869–13,882 (2007).


Introduction
Coupled dielectric resonators, particularly in the form of microdisks, form an important building block for a large range of nanophotonic devices.Used as add-drop filters, these structures were among the first commercially viable nanophotonic devices [1], and their ability to concentrate the strength of the electromagnetic field while simultaneously retaining a strong degree of tunability, makes them particularly attractive to experimentalists whose aim is to exploit nonlinearities in a controllable manner [2].Recently, coupled microdisks have been employed as "photonic molecules" at Terahertz frequences [3], and ensembles of microdisks have been proposed for use as novel laser sources and for studying problems in fundamental physics such as Anderson localization of photons [4,5,6].
In general, the calculation of the resonant frequencies and consequent field distributions of these structures must be performed numerically.The numerical algorithms available for this calculation range from finite element (FE) and finite-difference time-domain (FDTD) simulations [7,8,9,10], to boundary element methods [11,12].At high wavenumbers, or with many resonators, the numerical solution becomes computationally challenging: given the innate symmetry of the disks, the most elegant way to compute the resonant frequencies and consequent field distributions is to use the multipole method [13], otherwise known as the T-matrix method or the multiple-scattering method.In this family of algorithms the fields in and around the disks are expanded in a series of cylindrical functions; the resonant wavenumber is then computed by matching the electromagnetic boundary conditions on the edge of the disk.This method is simple enough to be implemented in a short time when results are needed quickly, and yet is sufficiently powerful to allow fast and reliable predictions of resonant frequencies.Usually this difficulty is overcome by obtaining a reasonable estimate for the position of the resonance, and then using an efficient two-dimensional minimum-finding algorithm (such as Powell's method) to find the zero of the systems determinant.Despite these advantages, the finding of resonances using this method can be a complicated process, and is the most time-consuming part of the multipole algorithm.This is because the boundary-matching reduces to solving a nonlinear eigenvalue problem for a dense complex-valued non-Hermitian system, for which there is no option but to perform a search in the complex plane for the resonant wavenumber.This procedure of finding the resonances not only slows the method down considerably, but one also runs the risk of missing a resonance if the initial estimate is too far removed from the resonant value, In this paper we present an alternate approach for finding the resonances of coupled circular resonators that overcomes the limitations of the standard multipole approach.This is made possible by the re-formulation of the multipole algorithm in terms of a generalized linear eigenvalue problem, in the manner of the scattering quantization approach [14].This means that the resonances of the system are readily and efficiently calculated by tracking their evolution in the complex plane; thereby dispensing with the costly root-searching procedure.In addition, the algorithm is easier to implement, and requires no decision parameters to determine whether or not a resonance is "real".The method is particularly efficient for the computation of high wavenumber resonances, and can also be used for calculating the sets of resonant states of extremely large numbers of resonators.

Multipole expansion
For resonators with a thickness below the wavelength, the three dimensional problem can be reformulated into a two dimensional problem with an appropriate wavelength dependent effective index.We consider a two-dimensional Euclidean domain with N c circular disks of radius R i and index of refraction n i (See Fig. 1).Between the i th and j th disk the distance is designated d i j and the angle θ i j ; the disks may be placed arbitrarily so long as they do not overlap.In two di- mensions the full electromagnetic problem can be decomposed into two distinct polarizations, each described by a scalar field ψ(r) and corresponding to the z-component of the electric or magnetic field: for TE polarization ψ = E z while for TM polarization ψ = H z .In each case the scalar field obeys the Helmholtz equation where n(r) = n i within the i th disk, as well as the appropriate boundary conditions on the edge of each disk.k is the wavenumber and R is a length scale which we choose in our calculations to be the radius of the largest disk.We then expand the solution in terms of cylindrical basis functions around each of the disks.The field is given within the i th disk by and outside the disk the field is expanded into In these equations the functions H (1,2) m (z) are Hankel functions of the first and second kind, corresponding to outgoing and incoming waves respectively.J m refers to the standard Bessel function of the first kind.The Hankel functions H (1) m in (2) represent the outgoing field contribution from the i th cylinder, and the non-singular part (containing the J m terms) can then be identified with the sum of the outgoing components from all the other disks in the cluster.We note that the expansion (2) has a domain of convergence extending from the boundary of the i th cylinder to the boundary of the nearest cylinder; if a complete field expansion is desired then usually one must make use of the Wijngaard expansion, as noted in [15].Using the expansion (2) about each of the cylinders and identifying the differing field contributions leads to the Rayleigh-type identity: where the coefficients X ji nm govern the coupling between disks i and j for the different multipole orders n and m: Writing the unknown coefficients as vectors, such that a = [a i m ] T etc., this equation can be written in matrix form as d = Xc.The boundary conditions can then be written in the general form where the quantities B n and B ′ n are square matrices with dimensions corresponding to the vectors a, b, c and d.The problem is then completely specified by noting that the singularity at the center of each disk must vanish, so that a = b.
In the standard multipole formulation, one then eliminates the coefficients a, b and c, thereby obtaining an infinite linear system of the form where u is a vector related to the remaining unknown coefficients d.This system can be truncated to a finite multipole order N trunc , with truncation order determined by refractive index, kR, and the relative distances between the cylinders, as analysed in [16].Zeros of the determinant of this system correspond to electromagnetic resonances of the cluster of disks.This classical approach gives in principle a complete solution to the system.In practice however Eq. ( 6) is a nonlinear eigenvalue problem for the eigenvalue kR, and so cannot be readily solved except using a direct search for zeros of the determinant in the complex plane.This two-dimensional search involves many iterations to find a single solution, and although methods for speeding this process have been proposed [17], the task remains onerous.In addition it is easy to overlook resonances in a given range, especially for high frequencies where the solutions are closely spaced.
In contrast to this usual procedure we reformulate the system by first using the identity (3) to eliminate d, and then re-arrange the equations for the boundary conditions to obtain Setting a = b as before, we then extend the block matrix on the right-hand side to obtain The solution to this equation is the same as the solution to the generalized eigenvalue problem of the form where M(kR) and N(kR) are the block matrices in Eq. ( 8).The condition λ (kR) = 1 + 0i, known as the quantization condition [14], corresponds to a resonance of the system.Contrary to the traditional method outlined above, Eq. ( 9) is a standard generalised linear eigenvalue problem, for which there exist extremely stable and accurate routines [18].Additionally, we will see that the reformulation of the eigenvalue problem leads to regular behavior of the eigenvalue λ with respect to changes in the wavenumber kR.This allows us to derive an new algorithm for the calculation of the resonant states.

New algorithm: description, implementation and results
From the eigenvalue problem (9) we obtain a set of eigenvalues and eigenvectors {λ µ , u µ }, which depend parametrically on the complex wavenumber kR.A generic set of eigenvalues λ µ for kR = 20 is shown in Fig. 2(a).Each eigenvalue represents a potential resonance of the system; in order to find the resonance one must know how to change the complex parameter kR such that the eigenvalue reaches the quantization point λ (k q R) = 1 + 0i.For a single isolated disk each eigenvalue corresponds to a family of resonances with the same angular order.For strongly coupled disks this simple classification breaks down, however even in this case each eigenvalue corresponds to a quantized state for a particular value of kR.By changing kR and tracking the motion of the eigenvalues to the quantization point, a complete set of modes can be obtained.
In order to find the resonances efficiently we must know how to change the parameter kR such that one eigenvalue approaches a true resonance of the system (at λ (k q R) = 1 + 0i) most rapidly.We expect a small change in kR to only shift the eigenvalues slightly; in Figure 3(a) we analyze this behavior of the eigenvalues with respect to a change in Re [kR] and Im [kR].Interestingly, we see that a real positive change in kR induces a counter-clockwise rotation of the eigenvalue in the complex plane.Similarly, a positive imaginary change in kR induces a radial decrease.Figure 3(b) and (c) show that the rate of change of the eigenvalues with respect to kR (this quantity we henceforth refer to as the speed of the eigenvalues) is very nearly constant.The dynamics of the eigenvalues may be understood by analysis of the matrix elements in (9) using the Debye expansion for cylindrical functions [14].
We can now define an efficient numerical algorithm based on this regular behaviour: 1. Perform a diagonalization of M and N at k 0 R and retrieve {λ µ 0 , u µ 0 } 2. Perform a second diagonalization M(k 0 R + δ kR), where δ kR is a small complex number with |δ kR| ≪ kR, and retrieve λ µ Following this algorithm multiple eigenvalues can be generated after only two diagonalizations, compared to a direct search which diagonalizes the matrix for each sample point in the complex kR plane.This can lead to a considerable enhancement of the multipole method.An analysis of the efficiency of a similar root search is shown in detail in Ref. [19].Once an eigenvalue reaches the quantization point, a family of additional resonances can be computed by increasing Re [kR] in such a way that the eigenvalue performs a complete rotation in the complex plane.
In the case of an isolated disk, this corresponds to an increase in the radial quantum number of the resonance.We show this for two coupled disks in Fig. 4. We note that the motion of the eigenvalue does not quite exactly track the unit circle; a small correction to Im [kR] is necessary to obtain the new quantized state.This reflects the physical observation that the higher radial order resonances possess a higher radiative loss.
As the truncation order N trunc increases, the size of the system (9) increases proportionally and more eigenvalues appear very close to the quantization point, that is, at λ = 1 + δ i, where δ ≪ 1 is a small positive number.One might assume that a small decrease in Re [kR] would bring this new state to quantization, however, for these new states the speed goes to zero as the frequency decreases, and so the quantization point cannot ever be reached from above.As the real part of the wavenumber increases, these eigenvalues increase in speed, following the unit circle until their first quantization after one complete circuit.The eigenvalues of the system (9) have no physical meaning when the quantization condition is not fulfilled, since the boundary conditions on the disk edges are not quite satisfied.However, it is a useful observation that the eigenvectors u of ( 9) only change slightly as the eigenvalue traverses the complex plane.To quantify this statement, we define the overlap of two eigenvectors by the standard scalar product for the complex vector u.In Figure 2(b) we trace this overlap as Re [kR] is increased (leading to a rotation of the eigenvalue by approximately π/2 radians).We see that the overlap with the initial eigenvector hardly changes at all, despite the fact that the change in kR can no longer be considered small.This fact can be exploited if one wishes to avoid an extra diagonalization to reconstruct the fields: a good approximation to the resonant fields may be obtained by using the eigenvectors from the initial diagonalisation u 0 , together with the quantized eigenvalue λ q in the field expansion.This can save much effort for very large systems, and can also be used to obtain fields for a large number of the states without using additional matrix computations.
To demonstrate the effectiveness of this method, we have used it to compute the set of resonances for a range of systems (see Fig. 5), including a set of 16 randomly positioned disks (Fig. 5(a)), three strongly-coupled disks with large values of kR (Fig. 5(b)), and a random array of 1024 disks (Fig. 5(c)).This method has also been successfully applied to the calculation of resonances in Terahertz photonic molecules; the predicted resonances achieved excellent agreement with experimental measurements [3].
The algorithm presented here works best for large values of |kR|, in which the eigenvalues behave most regularly, provided that the truncation order is chosen appropriately [16].For low frequency solutions, especially with |nkR| < 5, the eigenvalues do not follow such wellbehaved paths and the predictive algorithm, which assumes constant speeds and motion, that is decoupled between the real and imaginary parts of kR, loses efficacy.However the eigenvalues still follow analytic paths in the complex plane and in principle the resonances can still be calculated, albeit with additional computational expense.This type of complicated motion can be observed in Fig. 5(d); the regularity of the eigenvalue trajectories is disrupted by the presence of anti-crossing-like features [14].However each eigenvalue still possesses long stretches over which the behavior is very regular.

Conclusion
An efficient numerical algorithm for finding resonances in coupled dielectric disk has been presented.With an accelerated complex root search mechanism we are able to reduce the number of necessary diagonalizations to only two.This method will be of particular use in systematically computing the resonances of non-symmetric clusters of disks, where the problem of the determination of resonances cannot be reduced using symmetry considerations.Furthermore, this method should be readily expandable to fibre geometries and coupled rings.

Fig. 2 .
Fig. 2. a) Distribution of the eigenvalues for two coupled disks at kR = 20.0,d = 0.1, and n = 1.5 b) Overlap of some representative eigenfunctions (defined by the scalar product), with respect to kR.

1 3 .
Determine the approximate angular and radial speed of eigenvalues λ 4. Assuming constancy of speeds, extrapolate an approximate quantization wavevector k µ q R

Fig. 3 .
Fig. 3. a) Several representative eigenvalues traced in the complex plane for changes of Re [kR] (circular motion) and Im [kR] (radial motion).b) Angular speed (change of φ with respect to kR) of eigenvalues, with respect to the angular position φ in the complex plane.c) Log of the radial speed of the eigenvalues with respect to Im [kR].In all the calculations we used (kR = [20, . . ., 22], d = 0.1, n = 1.5, N trunc = 33).

Fig. 4 .
Fig. 4. (left) false color representation of the a wavefunction with m = 28 at kR = 21.60andone iteration later at kR = 24.74− 0.04i (bottom).On the right the trace of the eigenvalue, in red, through one rotation in the complex plane.Also shown are the eigenvalues of the final diagonalization in black.

Fig. 5 .
Fig. 5. a) False color representation of the absolute value squared wavefunction of 16 disks on a randomized rectangular grid, for the state with kR = 4.502758 − 0.07872i, n = 1.5, N trunc = 10.b) Wavefunction of a high frequency whispering gallery mode in three coupled disks at kR = 40.52586− 0.0005924i, n = 1.5, d = 0.1, R = 1, N trunc = 60.c) Wavefunction in a randomized grid of 1024 dielectric disks of kR ∼ 1, n = 1.5, N trunc = 2. d) Traces of the eigenvalues of Eq. (9) in the 16 disk randomized on a rectangular grid.Six representative eigenvalues are traced with the highest overlap criteria.The 'radial' motion of the eigenvalues is generated by tracing the eigenvalues through a change in Im [kR], kR = [3.69− 0.08i, . . ., 3.69 − 0.78i], the 'rotational' motion by a change in the real component kR = [3.69− 0.08i, . . ., 5.30 − 0.08i], n = 1.5, N trunc = 8.Contrary to Fig. 3 the motion is not any longer completely regular, however the general behavior is good enough for a general prediction.