1. 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 non-linearities 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-dimensionalminimum-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, or if large numbers of resonances are closely spaced.

 

Fig. 1. Schematic of disk geometry.

Download Full Size | PDF

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 themultipole algorithmin 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.

2. 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 θ_ij; the disks may be placed arbitrarily so long as they do not overlap. In two dimensions 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 ⁽¹⁾ _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 theWijngaard 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 differentmultipole orders n and m:

Writing the unknown coefficients as vectors, such that a=[aⁱ_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 involvesmany iterations to find a single solution, and althoughmethods 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

whereM(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.

 

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.

Download Full Size | PDF

3. 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_qR)=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 wemust knowhow to change the parameter kR such that one eigenvalue approaches a true resonance of the system (at λ(k_qR)=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₀R and retrieve {λ^µ₀,u^µ₀}

2. Perform a second diagonalization M(k ₀ R+δkR), where δkR is a small complex number with |δkR|≪kR, and retrieve λm ₁

3. Determine the approximate angular and radial speed of eigenvalues λ

4. Assuming constancy of speeds, extrapolate an approximate quantization wavevector k^µ_q R

Following this algorithmmultiple 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.

 

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).

Download Full Size | PDF

 

Fig. 4. (left) false color representation of the a wavefunction with m=28 at kR=21.60 (top) and one 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.

Download Full Size | PDF

The eigenvalues of the system (9) have no physicalmeaning 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 ₀, 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 well-behaved 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.

 

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.

Download Full Size | PDF

4. 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.