Efficient and spurious-free integral-equation-based optical waveguide mode solver

Modal analysis of waveguides and resonators by integralequation formulations can be hindered by the existence of sp urious solutions. In this paper, spurious solutions are shown to be eliminated by introduction of a Rayleigh-quotient based matrix singular ity measure. Once the spurious solutions are eliminated, the true modes may be determined efficiently and reliably, even in the presence of degeneracy , by an adaptive search algorithm. Analysis examples that demonstrate the e fficacy of the method include an elliptical dielectric waveguide, two une qual touching dielectric cylinders, a plasmonic waveguide, and a realist ic micro-structured optical fiber. A freely downloadable version of an optical wa veguide mode solver based on this article is available. © 2007 Optical Society of America OCIS codes:(000.4430) Numerical approximation and analysis; (060.0060 ) Fiber optics and optical communications; (230.7370) Waveguides; (230.5750) Resonators; (240.6680) Surface Plasmons; (230.3990) Micro-optical devices; References and links 1. K. Hayata, M. Koshiba, M. Eguchi, and M. Suzuki, “Vectoria l Finite-Element Method Without Any Spurious Solutions for Dielectric Waveguiding Problems Using Transv er e Magnetic-Field Component,” IEEE Trans. Microwave Theory Tech. 34, 1120–1124 (1986). 2. B. Rahman, F. Fernandez, and J. Davies, “Review of finite ele m nt methods for microwave and optical waveguides,” Proc. IEEE79, 1442–1448 (1991). 3. J. T. Chen, T. W. Lin, K. H. Chen, and S. W. Chyuan, “True and s purious eigensolutions for the problems with the mixed-type boundary conditions using BEMs,” Finite Elem. Anal. Des.40, 1521–1549 (2004). 4. C. Vassallo, “1993-1995 optical mode solvers,” Opt. Quant m Electron.29, 95–114 (1997). 5. Y. Leviatan, A. Boag, and A. Boag, “Generalized formulatio ns for electromagnetic scattering from perfectly conducting and homogeneous material bodies-theory and numeri cal solution,” IEEE Trans. on Antennas Propag. 36, 1722–1734 (1988). 6. A. Hochman and Y. Leviatan, “Analysis of strictly bound mode s in photonic crystal fibers by use of a sourcemodel technique,” J. Opt. Soc. Am. A 21, 1073–1081 (2004). 7. W. Schroeder and I. Wolff, “The origin of spurious modes in n umerical solutions of electromagnetic field eigenvalue problems,” IEEE Trans. Microwave Theory Tech. 42, 644–653 (1994). 8. W. Schroeder and I. Wolff, “Full-wave analysis of the influ ence of conductor shape and structure details on losses in coplanar waveguide,” IEEE MTT S International Microwave Symposium Digest 3, 1273–1273 (1995). 9. O. Bucci, G. D’Elia, and M. Santojanni, “Non-redundant imp lementation of method of auxiliary sources for smooth 2D geometries,” Electron. Lett. 41, 22 (2005). 10. B. E. Spielman and R. F. Harrington, “Waveguides of Arbitr ary Cross Section by Solution of a Nonlinear Integral Eigenvalue Equation,” IEEE Trans. Microwave Theory Tech. 20, 578–585 (1972). 11. V. Labay and J. Bornemann, “Matrix singular value decompos ition for pole-free solutions of homogeneous matrix equations as applied to numerical modeling methods,” IEEE Micr owave Guid. Wave Lett. 2, 49–51 (1992). #84481 $15.00 USD Received 25 Jun 2007; revised 26 Aug 2007; accepted 26 Aug 2007; published 18 Oct 2007 (C) 2007 OSA 29 October 2007 / Vol. 15, No. 22 / OPTICS EXPRESS 14431


Introduction
Spurious solutions have been known to exist in various techniques for modal analysis of optical waveguides and resonators. Although more well-known in the context of the Finite-Element Method (FEM) [1,2], spurious solutions also appear when the method of moments (MoM) is applied to a surface formulation using local-domain basis functions, like in the Boundary-Element Method (BEM) [3], and also using entire-domain basis functions [4]. As shown in this paper,  [5,6]. Some spurious solutions are the result of the discretization of a continuous operator [7]. Others, however, can be traced to the continuous operator itself, which has, owing to its compactness, a continuum of eigenfunctions of arbitrarily small eigenvalue [8]. In this paper we focus on methods based on integral-equation formulations, such as the MoM and SMT, and show how spurious solutions may be eliminated by a modification of the matrix singularity measure used to determine the eigenvalues.
Methods based on surface and off-surface formulations are among the most successful electromagnetic modelling techniques. Compared with differential methods, the representation of the fields usually requires fewer unknowns, and curved boundaries and unbounded domains are easily accommodated. These advantages are largely a consequence of the basis functions used to approximate the fields. The basis functions are the Maxwellian fields of sources radiating in a homogeneous domain, located either on the media boundaries (as in the MoM), or slightly offset from them (as in the SMT). This choice exploits the tacit assumption that the material is piecewise homogeneous to render a very compact approximation of the fields. Indeed, in the context of scattering problems, it has been shown that the number of unknowns in the SMT is practically optimal and it approaches the number of degrees of freedom of the scattered fields [9].
Despite the advantages of integral-equation-based methods they are not as popular, in the context of mode determination, as the FEM and other differential methods. The main reason is probably that in contrast to the latter, mode determination in the MoM and the SMT leads to a nonlinear eigenvalue problem. The relevant impedance matrix depends nonlinearly on the frequency (and propagation constant, for waveguides), because the basis functions must obey Maxwell's equations for every potential eigen-frequency (and propagation constant). It turns out then, that the advantages of using Maxwellian fields as basis functions are offset by the cost of solving the resultant nonlinear eigenvalue problem.
Clearly, if the nonlinear eigenvalue problem could be solved efficiently and reliably, integralequation methods would become more useful for modal analysis. A number of works have addressed this issue. The simplest approaches consist of sampling the determinant of the impedance matrix in ever increasing sampling rates to try and verify that no modes are missed [10]. The problem with this and similar approaches which search for minima of the smallest singular value [11] is that modes may be degenerate or almost degenerate, and hence, there is no a priori bound on the required sampling rate. To overcome this problem, more sophisticated methods that rely on the so-called orthogonality in the local limit [12] have been proposed [12,13]. The idea is to take advantage of the (near) orthogonality of solution vectors that correspond to (nearly) degenerate modes by continuously correlating sets of candidate solution vectors at different frequencies or propagation constants. In [12], the minima of the smallest singular value are found by keeping track of singular values that correspond to all the candidate solution vectors on a given interval of frequency or propagation constant. Similarly, in [13], an eigenvalue-tracking method is used to track the smallest eigenvalues of the impedance matrix in the complex eigenvalue plane as the frequency is varied. The eigenmodes are determined by detecting real-line crossings of the complex eigenvalues. While these methods are more robust, their implementation is somewhat complicated and the overhead compared with a simple search cannot be entirely overlooked.
The search algorithm proposed in this paper makes use of the same quasi-orthogonality principle, though it does not require continuous tracking of many eigenvalues. Instead, the algorithm first tries to find the modes by simple (yet adaptive) sampling of the reciprocal of the matrix singularity measure introduced in this paper. This function is devoid of the irregularities caused by spurious solutions and typically resembles a sawtooth waveform whose minima correspond to the eigenmodes. The adaptive sampling algorithm exploits this to sample the function sparsely. Only after a reasonably good image is obtained, the finer details are resolved by use of the orthogonality in the local limit.
Although the algorithm assumes that the modes have real frequency or propagation constant, moderately lossy or leaky modes can be found very efficiently by estimating the corresponding imaginary part and starting a search in the complex plane from this estimate. A dielectric waveguide solver based on this paper that includes many examples and a graphical-user interface is freely available for download [14].
The remainder of this paper is organized as follows: In Section 2, we discuss the source and character of the spurious solutions in integral-equation formulations. The difficulties caused by these solutions are described in Section 3. In Section 4, we present the proposed modification of the matrix singularity measure and show how it eliminates the spurious solutions. As explained in Section 5, once free from spurious solutions, the matrix singularity measure can be sampled adaptively to determine the modes. Numerical examples are given in Section 6, and the last section is a summary.

Spurious solutions of integral equations
Wolff and Schroeder have pointed out two different types of spurious solutions in electromagnetic modelling schemes. The more well-known type was discussed thoroughly in [7]. These solutions are due to the discretization process; they solve the discrete approximation of a continuous operator equation but do not correspond to any of its solutions. In effect, the spurious solutions solve the equation in a weighted average sense, but if the error of the solution is evaluated by use of slightly different weighting functions a much larger error is found. In the context of the SMT, these spurious solutions were also recognized in [15]. They can be avoided at a marginal extra cost by overdetermination, i.e., by using more weighting functions than basis functions (about twice as many).
The other type of spurious solutions was alluded to in [8], where the BEM was used to analyze coplanar waveguides. In contrast to the previous type of spurious solutions, these do not appear at discrete frequencies or propagation constants. Instead they form a continuum of spurious solutions which exist at each and every frequency or propagation constant. These spurious solutions do correspond to eigenfunctions of the continuous operator.
The integral equations that are discretized in the BEM and SMT are usually Fredholm equations of the first kind. For example, the integral equation for the TM modes of a perfectly conducting cylindrical waveguide oriented along the z direction is given by (1) Here, H (2) 0 (·) is the Hankel function of zero order and second kind, J z is the unknown current distribution, and the radial wave number k ρ obeys the separation equation where k 0 is the free-space wave number, equal to ω/c. The assumed and suppressed time and longitudinal variation are exp( jωt) and exp(− jβ z), respectively, and c is the speed of light in vacuum. The radius vector r points to the waveguide boundary contour, C, and the equation must be satisfied at all points r ∈ C. The integration is carried out on the contour C ′ , and the radius vector that points to C ′ is r ′ . The integration variable, dl ′ , is a differential line element. is placed on C ′ which is a slightly dilated version of C. A non-trivial solution to Eq. (1) exists only for specific pairs of ω and β for which the integral operator has a zero eigenvalue.
In the context of integral equations, it is well-known that the compact integral operator, induced by a smooth or weakly-singular symmetric kernel, K(r, r ′ ), has an infinite number of real eigenvalues, {λ j }, which tend to zero [16, p. 38]. This is a direct consequence of the lowpass nature of the kernel function; rapidly varying eigenfunctions have, consequently, small eigenvalues. For example, in a circular geometry such an eigenfunction would be a cylindrical current sheet, of radius ρ 0 , carrying a unit-amplitude, circularly-harmonic current, exp ( jnφ ).
The longitudinal electric field, E z , in the region ρ ≤ ρ 0 due to this current is given by where J n is the nth order Bessel function and ε 0 is the free-space permittivity. When n ≫ k ρ ρ 0 , the spatial variation of the current is fast on a radial-wavelength scale, and E z can be approximated by This shows that the field inside of a unit amplitude current sheet can be made arbitrarily small by using a sufficiently rapidly varying current. As n is increased, the field right on the current sheet decreases and it also decays more rapidly as the observation point moves away from the current sheet. When the operator L is discretized, the finite number of basis functions, N, limits the spatial variation of the solutions. Hence, when N is small, the error of the spurious solution is large, and assuming the true solutions have a smaller error they will be easily discernable. On the other hand, as N is increased the error of spurious solutions decreases and they may hide a true solution. In the SMT this problem is more severe owing to the factor ( ρ ρ 0 ) n which appears in Eq. (5). As the current is displaced from the waveguide boundary, this factor strongly reduces the field there, thus reducing the error of the spurious solution.
Alternative surface formulations, such as Müller's [17], which lead to Fredholm equations of the second kind have also been proposed as a way of avoiding ill-conditioning (which is due to the spurious solutions). Applications of this approach to mode determination of optical waveguides can be found in [18,19].

The effect of the continuum of spurious solutions on mode determination
Since the first type of spurious solutions can be handled quite easily by overdetermination, we focus now on the second type and show how it renders mode determination problematic. The difficulties caused by the spurious solutions are manifest even in the simplest analysis cases. This enables us to discuss these difficulties while avoiding undue complications, by determining the TM modes of a circular perfectly conducting waveguide. More challenging examples are given in Section 6. The method of analysis is the SMT, and we will look at the effect of varying the number of filamentary sources and their distance from the boundary. It is assumed that the behavior of conventional MoM can be inferred from the behavior of the SMT in the limit that the number of sources increases and they approach the waveguide boundary.
In the SMT, TM modes are approximated by the fields of an array of N electrical current filaments, uniformly distributed on a circle outside the waveguide, concentric with the waveguide.
The filaments carry a longitudinally varying current, whose dependence on the longitudinal coordinate z, is exp(− jβ z) (this dependence is suppressed). We denote the waveguide radius by R. The radius of the circle on which the sources are located is larger than R by a factor α (α > 1). This geometry is depicted in Fig. 1 N-tuple column vector, denoted by I. They are found by requiring that the longitudinal electric field, E z , be zero at a set of M testing points, uniformly distributed on the waveguide boundary. This leads to a homogeneous matrix equation  [15], which, however, is not directly applicable to non-square matrices and is difficult to compute for even moderately large square matrices. Typically, the singularity measure is sampled on a sufficiently fine grid of k ρ R, yielding a plot like the one shown in Fig. 2a, or Fig. 2 in [20]. The cut-off wave numbers, which are at the first zeros of the Bessel functions, can be clearly seen as peaks of the matrix condition number.
It is interesting to note that the matrix condition number is quite high between the peaks. The high condition number is of course due to the continuum of spurious solutions. In the discrete version, these solutions consist of sources of equal magnitude and alternating sign, that generate fields strongly confined to the vicinity of the sources. Because their E z field is small on the waveguide boundary, they approximately solve Eq. (6), and this is reflected by the high condition number.
The difference between a true mode and a spurious one, is that while both have vanishingly small fields on the boundary, the fields of the spurious solution are also vanishingly small within the boundary. To distinguish between the two, we define a normalized error measure, ∆E, as where [Z] is an impedance matrix that maps the amplitude vector I to the values of E z at a number of sampling points inside the waveguide. Also in Eq. (7), the 2-norm, or root-meansquare value, is denoted by rms(·). The sampling points can be arranged, for example, on a uniform grid of reasonable density. When analyzing dielectric structures, the sampling points can be the testing points already distributed on the media boundaries, since the modal fields are in general not zero there. This is a convenient choice because [Z] can then be obtained directly from [Z] after sign reversal of some of its entries. To evaluate ∆E, the vector that solves Eq. (6) in the least-squares sense is found by a singular value decomposition. The value of ∆E for I found at every k ρ R of Fig. 2a, is shown in Fig. 2b, and it confirms that although the matrix is rather singular between the peaks, there are no modes there since their normalized error is high. We now turn to take a closer look at the singularities shown in Fig 2. Zooming on the first of these singularities, we see in Fig. 3, that it is very sharp and discontinuous. This behavior has undesirable consequences. Because of the discontinuity, there can be no efficient search algorithm to find the singularities since there is no indication that a point near a singularity is in fact close, unless it is closer than the discontinuous edges. The only option is to sample the measure of singularity on a very fine grid. The width of the singularity determines the appropriate sampling resolution, and as shown in Fig. 3, the width depends on the number and location of the sources. Note that a small change in these parameters is enough to shrink the singularity considerably. If the sampling grid is kept constant while the number of sources is increased, the singularity may shrink to the point where it falls between grid points and goes undetected. As it is often essential to verify that the results are correct by increasing the number of sources or changing their positions slightly, this lack of reliability in the detection is a serious drawback.
The behavior seen in Fig. 3 is easily explained by existence of spurious solutions. If for a fixed α, the number of sources is increased, it would be like increasing n in Eq. (5), and if for a fixed N, α is increased, it would be like increasing ρ 0 in Eq. (5). In both cases, the absolute error of the spurious solution would decrease, contracting the singularity.  wide enough, one should therefore use as few sources as possible and they should be placed as close as possible to the boundary. While this facilitates the detection of the singularities, the accuracy of the field approximation is bound to suffer. As we show in the next section, the need to balance these conflicting goals is obviated by the proposed method.
In a conventional surface-formulation, discretized by the MoM, the sources are placed on the boundary, i.e., α = 1. Therefore, the singularities in the MoM are wider than in the SMT for an equal number of unknowns. Nevertheless, the problem will occur in the MoM as the number of basis functions is increased (Usually, a MoM solution will require more unknowns than an SMT solution [9]). Indeed, the continuum of spurious solutions was identified in [8] in the context of the BEM, which can be classified as a surface-formulation.

A modified singularity measure
In [8], a Tikhonov regularization was proposed to deal with the continuum of spurious solutions. The penalty function proposed there was proportional to the norm of the second derivative of the unknown fields at the interfaces between homogeneous regions. The method proposed here has the advantage that it does not involve derivatives of the fields and it also does not require calibration of the penalty function coefficient.
As explained in Section 3, it is the normalized error, ∆E, which takes into account field values inside the waveguide that reliably indicates the existence of a true mode. Hence, instead of evaluating ∆E for the least-squares solution of Eq. (6), it would be better to find the vector that minimizes ∆E for a given k ρ R. The square of normalized error, (∆E) 2 , is a ratio of two positive definite quadratic forms known as a generalized Rayleigh quotient. Direct differentiation shows that the stationary values of a generalized Rayleigh quotient are the generalized eigenvalues of the following generalized eigenvalue problem The generalized eigenvector I min , which corresponds to the minimum generalized eigenvalue, ξ min , yields the minimum ∆E, which is simply ξ min . In essence, we have modified the matrix singularity measure; it can now be defined as the reciprocal of the square root of the small- est generalized eigenvalue of Eq. (8). The main advantage of this scheme is that I min and ξ min change continuously when moving from one singularity to the next, thus allowing reliable detection of the singularities.
The generalized eigenvalue decomposition can be carried out by a number of methods quite efficiently, i.e., without a significant increase in the computation time relative to any of the other common measures of singularity. Since we are only interested in the smallest eigenvalue, Arnoldi methods are most suitable [21]. They require, however, explicit multiplication of the matrices [Z] and [Z] by their Hermitian conjugates, and this operation is prone to round-off errors. In the vast majority of cases, however, this method yields excellent results and is used throughout this paper. A more robust, though slower, alternative would be to use Van Loan To demonstrate the effectiveness of the spurious-free formulation, the plots of Figs. 2 and 3 were recalculated with the same N and α using the proposed matrix singularity measure. The results are shown in Figs. 4 and 5, respectively. The slope from which the singularities protruded in Fig. 2 has disappeared in Fig. 4, and the minima of the curve can be found with far less iterations. This is the subject of the next section. As shown in Fig. 5, the high sensitivity of the width of the singularity to the location and number of sources has been completely eliminated.

Spurious-free mode determination
The first phase of the mode determination scheme consists of adaptive sampling of the singularity measure. The objective of this phase is to sample the singularity measure, ∆E, on a nonuniform sampling grid which reveals just enough detail for the next phase to detect all the minima in a prescribed interval of propagation constant or frequency. As the next phase utilizes the orthogonality in the local limit, closely spaced modes need not be resolved, but otherwise all the minima should be detected. In the description that follows it will be assumed that a waveguide is analyzed, and therefore the singularity measure depends both on the frequency, ω, and on the propagation constant, β . In the analysis of a resonator it would depend, of course, only on the frequency.

Adaptive sampling algorithm
The adaptive sampling algorithm is based on the observation that ∆E(β ) for a given frequency ω, resembles a sawtooth waveform. This can be readily observed in Fig. 6a, where ∆E(β ) is plotted on a linear scale (as opposed to the logarithmic scale used in the previous figures). As shown in Fig. 6b, if β is kept fixed and the frequency is allowed to vary, a similar result is obtained for ∆E(ω). The sawtooth-like shape is obtained when the proposed singularity measure is used. Of course, if Fig. 2 were plotted on a linear scale it would still be discontinuous and completely unsuitable for an efficient search algorithm. The algorithm, shown in Fig. 7, exploits the shape of ∆E by assuming that monotonic parts of the function are approximately straight lines. For a given interval of β , the algorithm attempts to ascertain as quickly as possible whether the function ∆E(β ) is monotonic on the interval. As explained in the next paragraph, this is where the shape of ∆E(β ) is utilized. If ∆E(β ) is indeed monotonic, the algorithm ends with empty output. If, on the other hand, a minimum or a maximum is detected, its location is determined to higher accuracy by a standard golden section search [23]. The interval is then subdivided to intervals to the left and to the right of the extremum found, and the algorithm is called recursively with these subintervals as input. The output then is a concatenation of this extremum and those found in the two subintervals. At the end of the recursive process the output is a list of all the extrema encountered. Ascertaining whether a function is monotonic on a given interval is of course a difficult problem in the general case. However, by exploiting the sawtooth form of the function, a simple and efficient algorithm may be proposed (see Fig. 8). The algorithm begins with three samples: at the interval endpoints and at its midpoint. It then adds samples as long as the sample series is monotonic and the number of samples, n, is smaller than a prescribed number, n max , which can be determined empirically as will be shown below. The samples are added as follows. Assuming the function has been sampled n times, a straight line is fit through every three consecutive samples by linear regression. For each fit, the correlation coefficient, r, is evaluated by where β i are the sample points, ∆E f (β i ) is the value of the fit at these points, and the angle brackets denote the mean value. A correlation coefficient close to 1 indicates that the three points are approximately collinear. The algorithm finds the three consecutive points, β m−1 , β m , and β m+1 , which correspond to the smallest value of r, i.e., the ones for which the straight line fit is of worst quality. This deviation from a straight line possibly indicates an extremum, and therefore, the following samples are taken to be (β m−1 + β m )/2 and (β m + β m+1 )/2. If the sample series remains monotonic when the maximum number of samples is reached, it is concluded that ∆E(β ) is monotonic. If an extremum is encountered during the sampling process, ∆E(β ) is obviously not monotonic.
As illustrated in Fig. 9, the sampling points tend towards regions of high curvature where an extremum is likely to be found. However, this sampling scheme will have difficulty detecting both minima of a nearly degenerate mode pair. When the first minimum of the pair is found, while ∆E is monotonic and n n max find the three less collinear points, β m−1 , β m , and β m+1 the algorithm will be called with a new subinterval which includes the undetected mode near one of its endpoints. If the undetected mode is sufficiently close to the endpoint, the function will appear to be monotonic on the subinterval and the mode will be missed. This case can also be seen in Fig. 9, where a mode has been detected at the left endpoint in the previous call to the search algorithm. The mode just to the right of the left endpoint may go undetected, even after further subdivision of the intervals because the number of samples required for its detection could be more than n max . The probability of this happening may be reduced (at a relatively small expense) by choosing the first two samples very close to the endpoints. However, since minima may be arbitrarily close, a mode may still be missed. Such close modes are resolved by the next phase of the algorithm. The first phase is summarized in the animation shown in Fig. 10.

Resolving degenerate and nearly-degenerate modes
It is well-known that eigenmodes of waveguides and resonators obey certain orthogonality relations. In general, however, it is not true that the solution vectors of Eq. (6) that correspond to different modes obey an orthogonality relation. In fact, the number of modes is infinite whereas the number of linearly independent solution vectors is bounded by their length, N. However, the generalized eigenvectors of Eq. (8) corresponding to different eigenvalues (which are the con- A high resolution plot of ∆E is shown in (a) for the same step-index fiber of Fig. 6, but with λ 0 = 2R. In the first step of the sampling scheme (b), the endpoints and midpoint are sampled. Since the sample series is monotonic, two new samples are added at (β 0 + β 1 )/2 and (β 1 + β 2 )/2, as shown in (c). Since the sample series is still monotonic, the three less collinear consecutive points are found (β 3 , β 4 , and β 5 ) and two points are added in between them (d). This last refinement reveals a minimum and a maximum. tinuity condition errors) are orthogonal. So the solution vectors that correspond to degenerate modes are orthogonal. For a generalized eigenvalue problem, Similarly, if the modes are nearly degenerate, their solution vectors will be nearly orthogonal [12]. This quasi-orthogonality can be used to resolve degenerate and nearly degenerate modes. Suppose the previous phase of the algorithm has found one of the minima of two closely spaced minima. The solution vector at this minima is a by-product of the evaluation of the matrix singularity measure ∆E; let this vector be denoted by I 1 . Now, instead of [Z] † [Z] in Eq. (8), substitute the matrix [A] given by   6). The eigenvalues are also the same except for ξ 1 , the eigenvalue of I 1 , which is now ξ 1 + ξ 0 instead. By setting ξ 0 to a high enough value (say, ξ 0 = 1), the minimum which corresponds to ξ 1 can be removed from the ∆E(β ) function which is defined as the smallest generalized eigenvalue at a given β . An important feature of this process is that it does not create any spurious minima. When the modes are only nearly-degenerate, the eigenvectors and eigenvalues of all the modes will be in general altered by the proposed substitution. Nevertheless, as can be seen in Fig. 11, this alteration is negligible when the modes are close enough. After removal of second minimum Fig. 11. The resolution of two close modes in the step-index fiber of Fig. 9. Assuming the sampling process missed the first minimum but found the second minimum, the neighborhood of the second minimum can be searched again, this time with the matrix [A] of Eq. (12). The first minimum can then be easily detected.
The input of the second phase of the algorithm is the list of extrema found in the previous phase. The second phase has to check each minimum found for a nearby mode that was missed. To this end, the sampling algorithm could be run again for each interval defined by two adjacent maxima (or a maximum and an endpoint), with the minimum between them removed. While this is probably the most robust possibility, a faster alternative is to remove the minimum and then check the value of ∆E(β ) at the point of the removed minimum. If the value is small enough, it indicates that there is in fact another minimum nearby and that another search is warranted.
To determine whether another search is warranted the following criterion is used. The value of ∆E(β ) with the minimum removed is checked against the values of ∆E(β ) at the maxima on both sides of the minimum, with the minimum present. If it is smaller than the smallest of these values, another search between the two maxima is conducted. After the new minimum is found the criterion is tested again, this time with both minima removed. This process continues until it is concluded that there are no more minima in between the two maxima, and the search moves on to the next minimum found in the previous phase. The above criterion will be always correct when the modes are degenerate. As the modes move further apart, it may become incorrect and a mode could be missed. The probability of this happening is controlled by the resolution of the previous phase, which is determined by the parameter n max . By varying this parameter in numerous study-cases we have found, empirically, that n max = 7 leads to no modes being missed. However, the time complexity scales slower than linearly with n max , so, for added robustness, n max can be increased while incurring only a small increase in computation time.

Finding moderately lossy modes in the complex plane
It has been assumed until now that the modes are lossless and hence the minima are on the real β axis. For moderately lossy (or leaky) modes, the minima will be found in the complex plane, near the real axis. Usually, only the modes with the smallest losses are of interest, and these are very close to the real axis. The determination of these small losses can be quite challenging as the imaginary part of β , β i , can be several orders of magnitude smaller than the real part, β r . When this is the case, the effect of the losses is to blunt the minima of ∆E(β r ) when it is evaluated on the real line. It is then possible to estimate β i from the shape of the blunted curve. Furthermore, this estimate can be used as the starting point for standard quasi-Newton methods (we used the MATLAB implementation of Broyden's algorithm [25]) which can converge very rapidly to the complex minimum given a good first estimate.
In the previous phases of the algorithm, and especially in the golden section searches, the function ∆E(β r ) is sampled many times in the neighborhood of a minimum. These samples can be used to estimate β i . In analogy to the lossless case, it is assumed that ∆E(β ) depends linearly on the distance from the minimum point, β 0 , i.e., When evaluated on the real β axis, [∆E(β r )] 2 is given by, Where β 0 = β 0r + jβ 0i . This is a second degree polynomial in β r with coefficients p i , i.e., The polynomial which best fits the non-uniformly sampled data is easily obtained by the solving the corresponding Vandermonde matrix equation [26, p.119-124]. By comparing the coefficients in Eq. (13) and Eq. (15), the following estimate for β 0 can be obtained:  An example of the use of this method is shown in Fig. 12. As shown in the figure, the estimate is very close to the real minimum, and a search starting from the estimate should converge quickly. The constant ∆E contours plotted in the bottom part of the figure are equally spaced in ∆E. Their nearly equal spacing in the complex β plane, as can be observed in the figure, indicates that Eq. (13) is a good approximation. To find degenerate modes in the complex plane, a minimum found can be eliminated, and the quasi-Newton method can be started again from the same estimate. It should then converge to a nearby minimum, if present.

Numerical results
To demonstrate the various aspects of the proposed algorithm, a few optical waveguide analysis examples are given in this section. Some of the results shown have been obtained previously by other methods; these are used to validate the code. New data, however, is also presented.

Two touching cylinders
The geometry analyzed consists of two touching dielectric cylinders of unequal radii R 1 and R 2 , surrounded by an air cladding, as shown in the inset of Fig. 13. This geometry has been analyzed previously in [27] by a Rayleigh formulation, and our results are validated against the six digit results for the effective index given there. We also give the dispersion curves (shown in Fig. 13) of the first modes of the two touching cylinders, which were not shown in [27]. The x coordinate of the figure is the normalized frequency, V , given by, where ε rc is the relative permittivity of the cladding medium, and ε rw is the relative permittivity of the two dielectric waveguides. In Fig. 13, the birefringence of the fundamental mode pair, given by, is also shown. Here, n x and n y are the effective indices of the x-and y-polarized modes. Similarly to an elliptical waveguide, the birefringence of the fundamental mode pair tends to zero at the high and low normalized frequency limits [28]. The convergence of the effective index with the number of sources used is shown in Table 1. To find hybrid modes, both electric and magnetic current filaments are used [15], and an equal N β /k 0 ∆E Ref.
[27] 1.463292 Table 1. Convergence of the effective index. The normalized frequency, V , is 3.5. All the rest of the parameters are given in the caption of Fig. 13.
number of each type are placed inside and outside of each homogeneous region. In Table 1, and in all the tables and figures that follow, we use N to denote the number of sources of each type and on each side of a material boundary, so the total number of sources is actually 4N. We find that our solver converges to the same value given in [27]. two cylinders, and penetrates quite moderately into the other cylinder. This is a consequence of the unequal radii, and it would result in only moderate coupling between the two cylinders, even though they are touching.

An almost-circular dielectric waveguide
The next example demonstrates the reliability of the method in finding numerous modes. An elliptical waveguide with an aspect ratio of 1.05 was chosen as a challenging example, as it is expected to have pairs of very closely spaced modes. The modes of the waveguide were found for three normalized frequencies V = 1, 2, 3. In this case, the normalized frequency is again given by Eq. (17), with R 1 the semiminor axis of the ellipse. The number of modes found at each normalized frequency as a function of the number of sources is shown in Fig. 15(a). As To verify the that no modes are missed owing to possible deficiencies of the sampling algorithm (as opposed to insufficient number of sources), symmetry considerations [6, 29] were used, as follows. In an elliptical waveguide, modes can be classified into four different classes according to their symmetry. Since nearly-degenerate modes belong to different symmetry classes, it is easier to find the total number of modes by searching for each class separately. The total number of modes found in this way matched those found without exploiting symmetry. To verify the results further, the propagation constants found were compared with those of a circular dielectric waveguide, for which the characteristic equation (the solution of which yields the propagation constants) is well-known [30, p. 296]. This test also showed that the total number of modes found was correct, when enough sources were used. The interested reader is referred to the Appendix for the propagation constants of all the modes found.
The error in continuity conditions was averaged over all the modes found at a given frequency and number of sources. The average error, ∆E, is shown in Fig. 15(b). An exponential decrease of the average error with the number of basis functions, typical of integral equation methods, can be readily observed. Higher order modes require in general more sources owing to their faster spatial variation. This explains why the average error increases with frequency, for a given number of sources, or conversely, why higher frequencies require more sources to attain a prescribed average error.
In Fig. 15(c), the computation time per mode is shown for the three normalized frequencies studied. The three plots fall practically one on top of the other, even though, as the frequency increases, the average density of the modes increases, and their minimum separation decreases (for V = 1 it is 1.9 × 10 −3 , for V = 2, it is 1.8 × 10 −5 , and for V = 3, it is 4 × 10 −8 ). This is an appealing feature of the sampling algorithm. It implies that, for a given number of sources, the computation time would grow only linearly with the number of modes found, regardless of their average density and minimum spacing.
As usual in methods based on integral equations, the memory requirements are quite modest. Practically all the memory allocated is used up by the impedance matrix. For the calculations of Fig. 15, the number of testing points used was greater than the number of sources by a factor of 1.5, and each complex entry of the matrix required 16 bytes of storage. Therefore, the memory required, in kilobytes, is 0.375N 2 , which is 600 KBytes for the largest matrix used in the calculations of Fig. 15.

Circular and elliptical plasmonic nano-wire
The following example demonstrates complex eigenvalue determination. Plasmonic nano-wires have attracted considerable interest recently, owing to their unique wave-guiding (and light scattering) properties [31-33]. These lossy waveguides are made of a noble metal, such as silver, which is usually modelled by assuming a plasma-like frequency dependent permittivity. In some works (see [34] for example), a simple Drude model is assumed, whereas in others the permittivity is interpolated from measured values. We follow the latter practice and use the measured data given by Johnson and Christy [35].
The characteristic equation for a lossy circular nano-wire is the same as that of a lossless one, with a complex permittivity substituted for the real permittivity. The search for the effective indices, however, must be carried out in the complex plane, and the square root branch in Eq. (2) should be chosen according to the desired behavior at infinity. If, as in the following examples, proper modes are sought, the radial wave number outside the waveguide should have negative real and imaginary parts. For a circular silver nano-wire, the convergence of the effective index to the value obtained from numerical solution of the characteristic equation is shown in Table 2 also appears to converge to the exact value. It worth noting that when the number of sources is small, the estimate can err by a few orders of magnitude. However, the results improve greatly after the complex plane is searched, starting from this estimate. The reason for this is that when too few sources are used, the minimum of ∆E on the real line is blunted more by the numerical error than by the distance of the true minimum from the real line. Thus, using enough sources is important for ensuring rapid convergence to the minimum closest to the real line, especially if many minima are present. The dispersion with semimajor axis, a, of the TM 0 mode of an elliptical nano-wire made of silver is shown in Fig. 16 As the cross-sectional area of the nano-wire is decreased, the power flow density in the z direction concentrates in the metal and this is accompanied by more rapid dissipation and slower phase velocity.

A micro-structured optical fiber
Micro-structured optical fibers usually have quite complex cross-sections. To easily distribute the sources and testing points, it is most convenient to have an analytical parametric representation of the material boundaries. It is then easy to distribute the sources on contracted and dilated versions of the boundaries. The required parametric representation may be obtained, for example, by fitting spline curves to the material boundaries obtained a from thresholded SEM image. This approach was used in [38], where further details about it can be found. The fiber analyzed in this example was proposed for gas sensing applications in [39]. The geometry of the fiber, superimposed on the calculated Re(S z ) of the fundamental mode, is shown in Fig. 17. A parameter of interest in the design of the gas sensor is the fraction of power  observed, for small core diameters, the fraction of power in air can reach 50%.
The dispersion curves and birefringence of the fundamental mode pair are shown in Fig. 19. The two modes are roughly polarized in the x and y directions, and the birefringence is again given by Eq. (18). It is interesting to note that, in contrast to the two-touching cylinders, the birefringence does not vanish in the low frequency limit. This occurs also in micro-structured fibers with elliptical veins [28]. As explained in [28], in contrast to a waveguide surrounded by air, the cladding of the micro-structured fiber is itself birefringent, and therefore, even though the field extends into the cladding in the low frequency limit, the birefringence does not vanish.
A common practice (and one we adopted) when analyzing micro-structured fibers, is to assume that the dielectric extends to infinity. This assumption implies that the modes are necessarily leaky, and their confinement losses can be determined by searching the complex β plane [40]. However, since the fiber has a solid core, the fundamental core mode is guided by total internal reflection. It is thus evanescent in air, and consequently, in a real fiber, any transverse radiation would be reflected back at the outer interface between the fiber and the surrounding air [6]. Confinement losses were therefore not calculated. On the other hand, the air jacket was not modelled either as its effect on the real part of the effective index is anticipated to be negligible away from cut-off.

Summary
We have expanded on the nature of spurious solutions in integral equation formulations, and described a simple method for their elimination. An adaptive sampling algorithm that takes advantage of the resulting form of the matrix singularity measure has been presented. Methods for resolving degenerate and nearly degenerate modes and for determining complex-valued eigenvalues have also been described.