SYSTEM AND METHOD FOR FAST SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS, BASED ON AN EIGENPERMITTIVITY MODAL APPROACH

Abstract

A method for providing fast and efficient solution of partial differential equations to calculate the permittivity modes of an arbitrarily complex scatterer geometry using a modal approach, comprising the steps of defining the background geometry and the scatterer's geometry; embedding each scatterer in a simpler geometry; calculating the base transverse modes for each embedding geometry; for each scatterer, calculating the longitudinal modes; calculating the overlap matrix; solving and the resulting eigenvalue equation, using the base transverse modes that have been calculated for each embedding geometry and the longitudinal modes that have been calculated for the each scatterer; if there is more than one scatterer, hybridizing the modes of each pair of scatterers, otherwise, solving the resulting eigenvalue equation for the complete structure; projecting the source is on target modes; substituting the result is in the final equation.

Claims

1. A method for providing fast and efficient solution of partial differential equations to calculate the permittivity modes of an arbitrarily complex scatterer geometry using a modal approach, comprising: a) defining the background geometry and the scatterer's geometry; b) embedding each scatterer in a simpler geometry; c) calculating the base transverse modes for each embedding geometry; d) for each scatterer: calculating the longitudinal modes; calculating the overlap matrix; solving and the resulting eigenvalue equation, using the base transverse modes that have been calculated for each embedding geometry and the longitudinal modes that have been calculated for said each scatterer; e) if there is more than one scatterer, hybridizing the modes of each pair of scatterers, otherwise, f) solving the resulting eigenvalue equation for the complete structure; g) projecting the source is on target modes; and h) substituting the result is in the final equation.

2. A method according to claim 1, wherein a limited set of overlap integrals between the modes of an embedding simple shape are computed, for allowing solving a smaller eigenvalue equation for the modes.

3. A method according to claim 1, wherein the modes of an arbitrarily complex scatterer embedded inside a bigger scatterer of a simple shape are expanded, using the completeness of the permittivity modes inside said scatterer.

4. A method according to claim 1, wherein the eigenvalues of the simpler shapes are computed using an algebraic equation.

5. A method according to claim 1, wherein the entries in the eigenvalue equation are overlap integrals of the known modes over the domain of an arbitrary scatterer.

6. A method according to claim 1, wherein the calculation of modes incorporates naturally the mode discontinuity at the scatterer boundary, for minimizing the size of the eigenvalue problem that has to be solved.

7. A method according to claim 2, wherein the overlap integrals are evaluated by replacing the volume integrals by surface integrals and the surface integrals by line integrals.

8. A method according to claim 1, wherein modes of the target geometry are expanded based on the modes of a simpler geometry.

9. A system for providing fast and efficient solution of partial differential equations to calculate the permittivity modes of an arbitrarily complex scatterer geometry using a modal approach, comprising at least one processor, adapted to: a) define the background geometry and the scatterer's geometry; b) embed each scatterer in a simpler geometry; c) calculate the base transverse modes for each embedding geometry; d) for each scatterer: calculate the longitudinal modes; calculate the overlap matrix; solve and the resulting eigenvalue equation, using the base transverse modes that have been calculated for each embedding geometry and the longitudinal modes that have been calculated for said each scatterer; e) hybridize the modes of each pair of scatterers if there is more than one scatterer, otherwise, f) solve the resulting eigenvalue equation for the complete structure; g) project the source is on target modes; and h) substitute the result is in the final equation.

10. A system according to claim 9, in which modes of the target geometry are expanded based on the modes of a simpler geometry.

11. A system according to claim 9, in which the calculation of modes incorporates naturally the mode discontinuity at the scatterer boundary being boundary-adapted modes, for minimizing the size of the eigenvalue problem that has to be solved.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0047] The above and other characteristics and advantages of the invention will be better understood through the following illustrative and non-limitative detailed description of preferred embodiments thereof, with reference to the appended drawings, wherein:

[0048] FIG. 1 shows a schematic illustration of the re-expansion approachthe actual scatterer, defined by (r), is embedded in a larger simple shape (a circle in the current case, defined by 0(r)); and

[0049] FIG. 2 is a flowchart of the fast solution of partial differential equations using a modal approach, according to an embodiment of the invention.

DETAILED DESCRIPTION OF THE INVENTION

[0050] The present invention proposes a method for fast solution of partial differential equations, using a modal approach. The new proposed method provides an efficient way to calculate the permittivity modes of an arbitrarily complex scatterer geometry, which does not require solving a differential equation by discretizing all space. Instead, the proposed method requires just a computation of a limited set of overlap integrals between the modes of an embedding simple shape and a solution of a small eigenvalue problem. In particular, the completeness of the permittivity modes inside the scatterer is exploited to expand the modes ({tilde over (E)}) of an arbitrarily complex scatterer (defined by (r), shown in FIG. 1) embedded inside a bigger scatterer of a simple shape (e.g., a slab, a circle or a sphere, defined by .sub.0(r)), namely,

[00001]$\begin{array}{cc}\tilde{E}\left(r\right)={.Math.}_n{c}_n{E}_n\left(r\right),& \left(1\right)\end{array}$

where E.sub.n are the modes of the simple shape and c.sub.n is a vector of weights. As the eigenvalues of the latter simple shapes E.sub.n are relatively easy to compute using just an algebraic equation (the so-called dispersion relations), the computational burden is mostly to calculate the weights c.sub.n. A simple substitution into the underlying equations shows that the vector c.sub.n is a solution of an eigenvalue equation whose size is determined by the number of known modes E.sub.n that significantly contribute to the modes of the arbitrary scatterer {tilde over (E)}. The entries in that eigenvalue equation are just overlap integrals of the known modes over the domain of the arbitrary scatterer. Hence, discretization of all space is unnecessary. Discretization is necessary just for the matrix element calculation, the number of which is significantly smaller than what would have been the number of elements necessary to discretize all space.

[0051] FIG. 1 shows a schematic illustration of the re-expansion approach. The actual scatterer, defined by (r), is embedded in a larger simple shape (a circle in the current case, defined by .sub.0(r)).

[0052] The modes of the simple structures have an a-priory known functional form, and the eigenvalue itself is found from the roots of the dispersion relations, which are algebraic rather than differential equations. A highly efficient way to find those (generally complex) roots is the contour methods ([29]). However, an additional set of modes is required, specifically, modes that incorporate the permittivity discontinuity, the so-called longitudinal modes, which are characterized by zero permittivity eigenvalues.

[0053] The present invention proposes an efficient way to compute the longitudinal modes, according to which, the calculation of modes incorporates naturally the mode discontinuity at the scatterer boundary. This way, a minimal set of modes is found, typically 2-3 orders of magnitude smaller compared with previous approaches [10]. This minimizes the size of the eigenvalue problem that has to be solved, and effectively eliminates noise problems and the appearance of both transverse and longitudinal spurious modes.

[0054] The proposed method uses techniques that allow evaluating the resulting overlap integrals by replacing the volume integrals by surface integrals and the surface integrals by line integrals, in order to compute the modes more rapidly. This is an improvement over previous implementations [10], which involve only simple perturbations for which the overlap integrals could be calculated analytically. As a result, this method requires a substantially smaller computational effort, compared to the standard commercially-available techniques, while eliminating convergence problems or noise-induced spurious modes. Also, there are no missing modes.

[0055] The proposed method makes modal methods a powerful computational tool, which is characterized by high speed, high accuracy and high reliability.

[0056] The following examples show a rapid calculation of the overlaps for various non-trivial shapes and overall computation times of modes and scattering problems, which is several order of magnitude faster than available commercial software.

Generalized Normal Mode Expansion (GENOME)

[0057] One goal of the present invention is to solve Maxwell's equations with an arbitrary source J(r),

[00002]$\begin{array}{cc}\left(E\right)-k^2\left(r\right)E=i{}_{0}J,& \left(2\right)\end{array}$

where harmonic e.sup.it time variation and non-magnetic media are being assumed in this case. It is done by first obtaining the Lippmann-Schwinger equation and then, expanding the solution in terms of its eigenmodes. It is assumed that the structure defined by its permittivity profile (r) rests in a background of uniform permittivity .sub.b. This permits the manipulation of (2) to yield

[00003]$\begin{array}{cc}\left(E\right)-k^2{}_{b}E=i{}_{0}J+k^2\left(\left(r\right)-{}_b\right)E.& \left(3\right)\end{array}$

[0058] This can be solved using the Green's function for uniform media,

[00004]$\begin{array}{cc}\left({\left(\stackrel{=}{G}\right)}_0\right)-k^2{}_b{\stackrel{=}{G}}_0=\stackrel{=}{I}{}^3\left(r-r^{}\right),& \left(4\right)\end{array}$

which has a simple known analytic form G.sub.0(|rr|) depending on the dimensionality of the problem. Applying (4) to both terms on the RHS of (3) yields its Green's function solution, the Lippmann-Schwinger equation,

[00005]$\begin{array}{cc}E\left(r\right)=E_0\left(r\right)+k^2{\stackrel{=}{G}}_0\left(.Math.r-r^{}.Math.\right)\left(\left(r^{}\right)-{}_b\right)E\left(r^{}\right){\mathrm{dr}}^{}.& \left(5\right)\end{array}$

[0059] The term E.sub.0(r) is the known radiation pattern of external sources in a uniform background

[00006]$\begin{array}{cc}{E}_0\left(r\right)=i{}_0{\stackrel{=}{G}}_0\left(.Math.r-r^{}.Math.\right)J\left(r^{}\right){\mathrm{dr}}^{},& \left(6\right)\end{array}$

[0060] To solve (5), the appropriate normal modes of the system are being defined, obtained by neglecting source terms in (3). Attention is restricted to the simpler case obtained by assuming that the permittivity of the inclusion is uniform, so that the simulation domain is defined by only two permittivities: an eigenpermittivity .sub.m applicable to the inclusion's interior, and the fixed background .sub.b,

[00007]$\begin{array}{cc}\left(E_m\right)-k^2{}_b{E}_m=\frac{1}{s_m}{k}^2{}_b\left(r\right)E_m.& \left(7\right)\end{array}$

where s.sub.m is the mth eigenvalue,

[00008]$\begin{array}{cc}\frac{1}{s_m}\frac{{}_m-{}_b}{{}_b},& \left(8\right)\end{array}$

(r) is a step function that is unity inside the inclusion and zero elsewhere. Alternatively, the integral form of the eigenvalue equation can be obtained from (5) by neglecting E.sub.0,

[00009]$\begin{array}{cc}{s}_m{E}_m\left(r\right)=k^2{}_b{\stackrel{=}{G}}_0\left(.Math.r-r^{}.Math.\right)\left(r^{}\right)E_m\left(r^{}\right){\mathrm{dr}}^{}.& \left(9\right)\end{array}$

[0061] This provides the modal expansion solution of the Lippmann-Schwinger equation (5) for the simpler case of a uniform inclusion permittivity .sub.i, which after some manipulation is

[00010]$\begin{array}{cc}E\left(r\right)=E_0\left(r\right)+\frac{i}{{}_0}{.Math.}_m{E}_m\left(r\right)\frac{{}_i-{}_b}{\left({}_m-{}_i\right)\left({}_m-{}_b\right)}{E}_m^{}\left(r^{}\right).Math.J\left(r^{}\right){\mathrm{dr}}^{}.& \left(10\right)\end{array}$

[0062] Thus, the solution for any source configuration can be obtained almost immediately once the modes of the structure are available.

[0063] Returning to the defining differential equation (7), there are some properties of the modes. It allows three types of modes, of which one shall use two: transverse and longitudinal electric, referred hereafter simply as longitudinal. Transverse modes have a nonzero eigenvalue .sub.m, so the eigenvalue equation has two forms in the two piecewise uniform regions,

[00011]$\begin{array}{cc}\left(E_m\right)-k^2{}_m{E}_m=0r\mathrm{interior}& \left(11\right)\end{array}$$\left(E_m\right)-k^2{}_b{E}_m=0r\mathrm{exterior}.$

[0064] These modes are divergence-free within the two regions, which can be seen by applying the divergence operator to each of (11), obtaining .Math.E.sub.m=0 in both cases. They are not divergence-free along the boundary between the interior and the exterior, but nevertheless are called transverse modes.

[0065] Typically, only transverse modes are necessary for the expansion (10), since the solution E(r) is also everywhere transverse except along the boundary and locations where the source J(r) is nonzero. These non-transverse parts of E(r) are accounted for by the transverse modes and by E.sub.0(r), respectively. However, since the interest is in perturbing the structure (r) and its interfaces, the locations where the solution E(r) is non-transverse are also perturbed, requiring an additional set of longitudinal modes.

[0066] Longitudinal modes arise from the special set of eigenfunctions of (7) that exist when s.sub.m=1 or .sub.m=0. The eigenvalue equation then reduces to

[00012]$\begin{array}{cc}\left(E_m\right)=0,& \left(12\right)\end{array}$

in the interior. In contrast to the transverse modes, these modes are nonzero only inside the inclusion. The field is identically zero in the exterior since the index contrast with the background is infinite.

[0067] It will be shown that longitudinal electric modes are generated from the subset of (12) where E.sub.m=0, and their non-divergence-free property is also proved.

[0068] The implementation of the proposed Generalized Normal Mode Expansion method requires the calculation of 3D integrals, for the computations of the mode overlaps, as shown in eq. ( . . . ). This can be done, for example, by any software hat is based on the Finite Element Method (such as freefem.orgFreeFEM is a partial differential equation solver for non-linear multi-physics systems in 1D, 2D, 3D and 3D border domains, getfem.orgGetFEM is an open source library based on collaborative development. It aims to offer the most flexible framework for solving potentially coupled systems of linear and nonlinear partial differential equations with the finite element method, etc.). The computation of the target modes and eigenvalues may be performed by software that performs eigenvalue search of given matrices. This can be, for example, using the standard LAPACK (a standard software library for numerical linear algebra) package (implemented, e.g., via MATLAB).

[0069] The size of the simulated domain determines the memory and processor requirements. Those are expected to be significantly smaller compared with the implementation of standard numerical methods.

The Re-Expansion Method

[0070] The main computational effort in any modal expansion method is typically devoted to finding the modes, and there is always a need for a fast, efficient, and reliable method. It is may be achieved by expanding modes of the target geometry, defined by (r), using the modes of a simpler geometry (r) as a basis. The aim is to find the modes of an arbitrary geometry by perturbing the modes of a much simpler geometry. Because the basis is complete, the method is not merely perturbative in the traditional sense, as one can treat perturbations of any depth, successfully generating the modes of any target geometry to any desired accuracy.

[0071] It is required to solve the eigenvalue equation

[00013]$\begin{array}{cc}\left({\tilde{E}}_{}\right)-k^2{}_b{\tilde{E}}_{}=\frac{1}{{\tilde{s}}_{}}{k}^2{}_b\stackrel{}{}\left(r\right){\tilde{E}}_{},& \left(13\right)\end{array}$

where {tilde over (E)}.sub. and {tilde over (s)}.sub. are the perturbed eigenmodes and eigenvalues, and {tilde over ()}(r) differs from (r),

[00014]$\begin{array}{cc}\stackrel{}{}\left(r\right)=\left(r\right)+\left(r\right).& \left(14\right)\end{array}$

[0072] It is assumed that the modes of the geometry defined by (r) are known, while (r) can be arbitrary, as long as it is contained within the so-called interior of (r).

[0073] The solution proposed by the present invention focuses on piecewise uniform structures, so (r) has discontinuous steps. This case is perhaps the most difficult to treat, since the eigenmodes also have discontinuities. However, the formulation can be applied unmodified for (r) with smooth spatial variations, or even tensorial perturbations.

Derivation

[0074] The derivation follows the familiar pattern of expanding the target differential equation using a complete basis, initially with unknown coefficients. Then a project is performed to evaluate the coefficients, resulting in a set of integrals over the perturbation that populate a linear system of equations, which is solved for the final solution.

[0075] The perturbed modes {tilde over (E)}.sub. are expanded using the modes E.sub.m as a basis, since they are completely inside the inclusion (as illustrated in Eq. (1) above),

temporarily suppressing the index for brevity. Then (13) becomes

[00015]$\begin{array}{cc}\left(\tilde{E}\right)-k^2{}_b\tilde{E}-k^2\left(\stackrel{}{}-{}_b\right)\stackrel{}{}\left(r\right)\tilde{E}=0& \left(16\right)\end{array}$${.Math.}_n{c}_n\left[\left(E_n\right)-k^2{}_b{E}_n-k^2\left(\stackrel{}{}-{}_b\right)\left[\left(r\right)+\left(r\right)\right]E_n\right]=0.$

[0076] The first two terms satisfy the eigenvalue equation of the basis modes (7),

[00016]$\begin{array}{cc}{.Math.}_n{c}_n\left[k^2\left({}_n-{}_b\right)\left(r\right)E_n-k^2\left(\stackrel{}{}-{}_b\right)\left[\left(r\right)+\left(r\right)\right]E_n\right]=0.& \left(17\right)\end{array}$

[0077] The unknown coefficients are found by projecting onto the basis modes via multiplication by adjoint modes E.sub.m.sup. and integrating over all space,

[00017]$\begin{array}{cc}{.Math.}_n{c}_n\left[\left({}_n-{}_b\right)E_m^{}.Math.\left(r\right)E_n-\left(\stackrel{}{}-{}_b\right)E_m^{}.Math.\left[\left(r\right)+\left(r\right)\right]E_n\right]\mathrm{dr}=0.& \left(18\right)\end{array}$

[0078] This may be simplified using the orthogonality of the basis modes,

[00018]$\begin{array}{cc}{E}_m^{}.Math.\left(r\right)E_n\mathrm{dr}={}_{\mathrm{mn}},& \left(19\right)\end{array}$

which is described later. This yields

[00019]$\begin{array}{cc}\left({}_m-{}_b\right)c_m=\left(\stackrel{}{}-{}_b\right){.Math.}_n\left[{}_{\mathrm{mn}}+E_m^{}\left(i\right).Math.E_n\mathrm{dr}\right]c_n& \left(20\right)\end{array}$$\stackrel{}{s}{c}_m=s_m{.Math.}_n\left[{}_{\mathrm{mn}}+V_{\mathrm{mn}}\right]c_n,$

where V has been defined with matrix elements

[00020]$\begin{array}{cc}{V}_{\mathrm{mn}}=E_m^{}\left(r\right).Math.\left(r\right)E_n\left(r\right)\mathrm{dr}.& \left(21\right)\end{array}$

[0079] This yields a system of linear equations to be solved for the perturbed modes

[00021]$\begin{array}{cc}\mathrm{diag}\left(s_m\right)\left(V+I\right)c=\stackrel{}{s}c,& \left(22\right)\end{array}$

where I is the identity matrix and

[00022]$\begin{array}{cc}\mathrm{diag}\left(s_m\right)=\left[\begin{array}{cccc}{s}_1& & & \\ & s_2& & \\ & & & \\ & & & s_m\end{array}\right],V=\left[\begin{array}{cccc}{V}_{11}& V_{12}& & \\ V_{21}& V_{22}& & \\ & & & \\ & & & V_{\mathrm{mn}}\end{array}\right],c=\left[\begin{array}{c}{c}_1\\ c_2\\ .Math.\\ c_m\end{array}\right].& \left(23\right)\end{array}$

[0080] Solution of (22) yields the perturbed modes {tilde over (E)}.sub. and their eigenvalues {tilde over (s)}.sub. for the geometry defined by {tilde over ()}(r), expanded in terms of the basis modes E.sub.m. The numerical implementation of (22) first requires preparation of all the basis modes, both transverse and longitudinal, considered further in the description. Next, the matrix elements of (21) must be computed via integrals considered later. Finally, (22) is ready to be solved, using any numerical linear algebra package for small dense systems. For numerical convenience, the system of equations (22) should be symmeterized, a process that is specified later on. After solving (22), a few minor steps complete the process, while enabling the solution of the modal expansion solution (10). The perturbed adjoint modes adjoint modes are necessary, {tilde over (E)}.sub.m.sup. which are considered further in the description. The modes {tilde over (E)}.sub. also require normalization, which is performed later.

Adjoint Form

[0081] The modal expansion solution (10) requires not only the direct modes E.sub.m, but also requires the adjoint modes E.sub.m.sup. for the purposes of projection. All adjoint modes satisfy the same governing eigenvalue equation and have the same eigenvalues as the direct modes. This is a consequence of the symmetry of the free space Green's function (9). Furthermore, the modes are self-adjoint, that is, the adjoint modes E.sub.m.sup. are identical to E.sub.m and the orthogonality relation is

[00023]$\begin{array}{cc}{E}_m^{}.Math.\left(r\right)E_{n}dr={}_{mn}.& \left(36\right)\end{array}$

[0082] This is true unless the structure exhibits a form of symmetry that permits degeneracy. In this case, the set of degenerate modes are sometimes related by a non-standard orthogonality relation. Typically though, the adjoint mode is obtained by reversing the functional dependence along the direction of symmetry, for example e.sup.im.fwdarw.e.sup.im, so the basis is no longer self-adjoint. However, even in this case, a self-adjoint basis can be obtained by using the set {sin(m), cos(m)} for the angular variation.

[0083] According to one embodiment, the analytical modes of the cylinder are employed as a basis, incorporating the complex cylindrical harmonics, {e.sup.im}. Despite the non-self-adjoint nature of this basis, it is preferred due to convenience when evaluating the matrix elements of V, (21). Use of this non-self-adjoint basis means that the expansion of the perturbed adjoint modes may differ from the direct modes,

[00024]$\begin{array}{cc}{\tilde{E}}^{}\left(r\right)={.Math.}_m{c}_m^{}{E}_m^{}\left(r\right),& \left(37\right)\end{array}$

where c.sub.m.sup. labels the coefficients. Explicit computation is used to prove this fact. The derivation for these coefficients is almost identical as before, beginning instead with the adjoint form

[00025]$\begin{array}{cc}\left({\tilde{E}}^{}\right)-k^2{}_b{\tilde{E}}^{}-k^2\left(\stackrel{~}{}-{}_b\right)\stackrel{}{}\left(r\right){\tilde{E}}^{}=0.& \left(38\right)\end{array}$

[0084] Inserting the expansion (37) and eventually projecting onto the direct mode E.sub.n gives

[00026]$\begin{array}{cc}{.Math.}_m{c}_m^{}\left[\left({}_m-{}_b\right)E_m^{}.Math.\left(r\right)E_n-\left(\stackrel{~}{}-{}_b\right)E_m^{}.Math.\left[\left(r\right)+\left(r\right)\right]E_n\right]dr=0.& \left(39\right)\end{array}$

[0085] Due to slight differences in the subscripts, a system of linear equations featuring the transpose of matrix elements V is obtained, defined in (21),

[00027]$\begin{array}{cc}\mathrm{diag}\left(s_m\right)\left(V^{}+I\right)c^{}=\stackrel{}{s}{c}^{}.& \left(40\right)\end{array}$

[0086] It should be noted that in the notation c.sup. is a column vector. Thus, the coefficients of expansion for the adjoint fields, c.sub.n.sup., are identical to c.sub.n if V is symmetric or complex symmetric. The familiar case of complex conjugate adjoints are obtained if V is Hermitean.

[0087] The symmetry of V does not depend on the perturbation geometry (r), since this is simply a scalar function. Instead, it only depends the relationship between the basis modes E.sub.m and their adjoints E.sub.m.sup., which in turn stems from the symmetry properties of the operator G.sub.0(|rr|) of the governing eigenvalue problem (9). It can be seen that because G.sub.0(|rr|) is symmetric, V can also be put into symmetric form, thus providing a simple relationship between c.sub.n and c.sub.n.sup..

Symmetry

[0088] For numerical convenience, the eigenvalue equation (22) may be cast in a more symmetric form by dividing and multiplying by diag({square root over (s.sub.m)}),

[00028]$\begin{array}{cc}\mathrm{diag}\left(\sqrt{s_m}\right)\left(V+I\right)\mathrm{diag}\left(\sqrt{s_m}\right)\left[\mathrm{diag}\left(\frac{1}{\sqrt{s_m}}\right)c\right]=\stackrel{}{s}\left[\mathrm{diag}\left(\frac{1}{\sqrt{s_m}}\right)c\right].& \left(41\right)\end{array}$

[0089] However, an impediment to full symmeterization remains since the matrix V may itself not possess any symmetry properties, such as being equal to its transpose V.sup. or its conjugate transpose V*. In this case, the lack of symmetry in V is due to the unusual orthogonality relationship obeyed by the cylindrical harmonics {exp(im)}, which can be remedied by switching to the trigonometric basis {cos(m),sin(m)}. Thus, full symmeterization can be attained since V.sub.mn bears sufficient similarity to V.sub.nm that they can be symmeterized using a similarity transform S,

[00029]$\begin{array}{cc}SV{S}^{-1}=W=S^{-1}{V}^{}S,& \left(42\right)\end{array}$

where W.sup.=W. The second equality of (42) is a consequence of the first equality and the definitions of S and W. In more concrete terms, a unitary basis transformations of the type

[00030]$\begin{array}{cc}\sqrt{2}\left[\begin{array}{c}\cos \left(m\right)\\ \sin \left(m\right)\end{array}\right]=\frac{\sqrt{2}}{2}\left[\begin{array}{cc}1& 1\\ -i& i\end{array}\right]\left[\begin{array}{c}{e}^{-im}\\ e^{im}\end{array}\right]& \left(43\right)\end{array}$

is required, to be applied to the appropriate pairs of V. The identity J.sub.m(r)=(1).sup.mJ.sub.m(r) can be employed if Bessel functions of negative order require conversion. In this way, the desired similarity transform S can be constructed.

[0090] The eigenvalue equation is now obtained

[00031]$\begin{array}{cc}\mathrm{diag}\left(\sqrt{s_m}\right)\left(SV{S}^{-1}+I\right)\mathrm{diag}\left(\sqrt{s_m}\right)\left[\mathrm{diag}\left(\frac{1}{\sqrt{s_m}}\right)\mathrm{Sc}\right]=\stackrel{}{s}\left[\mathrm{diag}\left(\frac{1}{\sqrt{s_m}}\right)\mathrm{Sc}\right],& \left(44\right)\end{array}$

which after simplification yields the desired symmeterized form

[00032]$\begin{array}{cc}\mathrm{diag}\left(\sqrt{s_m}\right)\left(W+I\right)\mathrm{diag}\left(\sqrt{s_m}\right)b=\stackrel{}{s}b,& \left(45\right)\end{array}$

featuring the rescaled eigenvectors

[00033]$\begin{array}{cc}b=\mathrm{diag}\left(\frac{1}{\sqrt{s_m}}\right)\mathrm{Sc}.& \left(46\right)\end{array}$

[0091] Equivalently, the adjoint eigenvalue equation (40) may be symmeterized, using the rescaled eigenvector

[00034]$\begin{array}{cc}b=\mathrm{diag}\left(\frac{1}{\sqrt{s_m}}\right)S^{-1}{c}^{},& \left(47\right)\end{array}$

again yielding (45). The coefficients {c.sub.m} and adjoint coefficients {c.sub.m.sup.} are now obtained from the same eigenvalue equation.

[0092] This has several advantages: it eliminates any numerical noise arising from independently solving two similar eigensystems, it eases the numerical burden, and enables a more streamlined numerical implementation.

Normalization and Orthogonality

[0093] All modal expansion methods require that the modes be normalized, as this enables projection. Thus, the following condition should be evaluated

[00035]$\begin{array}{cc}{N}_{}^2{\tilde{E}}_{}^{}\left(r\right).Math.\stackrel{}{}\left(r\right){\tilde{E}}_{}\left(r\right)dr=1,& \left(48\right)\end{array}$

where N.sub. is the normalization factor. Thus, the modes produced by numeric solution of the linear system of equations (22) are not necessarily correctly normalized, since most linear algebra packages typically normalize according to linear algebra definitions. However, it was shown that a convenient evaluation of (48) is possible if the basis modes are already normalized, by expanding

[00036]$\begin{array}{cc}\begin{array}{c}{\tilde{E}}_{}^{}\left(r\right)\stackrel{}{}\left(r\right){\tilde{E}}_{}\left(r\right)dr={.Math.}_m{.Math.}_n{c}_m^{}{c}_n{E}_m^{}\left(r\right).Math.\left[\left(r\right)+\left(r\right)\right]E_n\left(r\right)dr\\ ={.Math.}_m{.Math.}_n{c}_m^{}{}_{mn}{c}_n+{.Math.}_m{.Math.}_n{c}_m^{}{V}_{mn}{c}_n\\ ={\left(c^{}\right)}^{}\left(V+I\right)c\\ ={{\stackrel{}{s}}_{}\left(c^{}\right)}^{}\mathrm{diag}\left(s_m^{-1}\right)c.\end{array}& \left(49\right)\end{array}$

[0094] To obtain the result, the orthonormality of the basis modes (19) in the second equality of (49) has been exploited, and the eigenvalue equation (22) in the fourth equality.

[0095] The normalization integral is thus reduced to a sum,

[00037]$\begin{array}{cc}\frac{1}{N_{}^2}={\stackrel{}{s}}_{}{.Math.}_m\frac{c_m^{}{c}_m}{s_m}.& \left(50\right)\end{array}$

[0096] The evaluation is even simpler in terms of the symmeterized coefficients b.sub.m defined in (46) and (47),

[00038]$\begin{array}{cc}\frac{1}{N_{}^2}={\stackrel{}{s}}_{}{.Math.}_m{b}_m^2.& \left(51\right)\end{array}$

[0097] This is because .sub.m b.sub.m.sup.2=1 is precisely the dot product between the normalized left and right eigenvectors of the symmeterized matrix operator in (45), which is naturally enforced by some linear algebra packages.

[0098] An extension of (49) can be used to test the orthonomality of modes, by numerically confirming the equality

[00039]$\begin{array}{cc}{\left(D^{}\right)}^{}\mathrm{diag}\left(s_m^{-1}\right)\mathrm{Ddiag}\left({\stackrel{}{s}}_{}\right)=I,& \left(52\right)\end{array}$

where D is a matrix of normalized coefficients composed of columns N.sub.c.sub., D.sup. is composed of columns of N.sub.c.sub..sup., and diag({tilde over (s)}.sub.) is a diagonal matrix of all eigenvalues {{tilde over (s)}.sub.}. The identity (52) is numerically satisfied to machine precision regardless of truncation, since the eigenvectors of a complex symmetric matrix (45) obey a complex orthogonality condition from which (52) can be deduced. One interesting consequence is that the set of generated modes {{tilde over (E)}.sub.} is automatically orthogonal. The exception is modes belonging to the same eigenvalue, such as symmetry degenerate modes, which need to be orthogonalized manually.

Transverse Basis Modes

[0099] Perturbation methods uses a basis of known modes of a simple geometry to expand the modes of a complex target geometry. In broad terms, the transverse basis modes reproduce the smooth near-field features and all of the far-field features of the target modes. Typically, modes of circles and spheres are employed, as these modes can be obtained analytically by separation of variables, ultimately reducing to a transcendental equation which can be solved for the eigenvalues. A brief overview of these well known results is provided, pertaining primarily to 2D modes of a circle. Very similar procedures apply in 1D and 3D, and only the form of the mathematical functions differ.

[0100] Applying separation of variables to the Lippmann-Schwinger equation (7) yields the function form of the eigenmodes. For structures that are invariant in the third dimension, two distinct polarizations exist: transverse electric and transverse magnetic. The field profiles for the two cases differ. For the TM modes, it is

[00040]$\begin{array}{cc}{E}_z=H_m\left({}_{b}r\right)e^{\mathrm{im}},r\mathrm{exterior}& \left(53\right)\end{array}$$E_z=J_m\left({}_{m}r\right)e^{\mathrm{im}},r\mathrm{interior}.$

[0101] For the TE modes, it is simpler to specify the magnetic fields

[00041]$\begin{array}{cc}{H}_z=H_m\left({}_{b}r\right)e^{\mathrm{im}},r\mathrm{exterior}& \left(54\right)\end{array}$$H_z=J_m\left({}_{m}r\right)e^{\mathrm{im}},r\mathrm{interior},$

and use the Maxwell curl equation

[00042]$\begin{array}{cc}E=\frac{i}{k}H& \left(55\right)\end{array}$

to obtain the electric fields. The magnetic field representation is preferred when for example, performing integrals over the perturbation, which is described later in the description.

[0102] The propagation constants in (53) and (54) are related to the background permittivity and the as yet unknown eigenpermittivity,

[00043]$\begin{array}{cc}{}_m^2=k^2{}_m,{}_b^2=k^2{}_b,& \left(56\right)\end{array}$

[0103] The eigenpermittivity is found via a field stitching exercise at the boundary between the interior and exterior, defined at radius B. This leads to a form of the well-known dispersion relation of a step-index fiber,

[00044]$\begin{array}{cc}\left(\frac{1}{{}_{m}B}\frac{J_m^{}\left({}_{m}B\right)}{J_m\left({}_{m}B\right)}-\frac{1}{{}_{b}B}\frac{H_m^{}\left({}_{b}B\right)}{H_m\left({}_{b}B\right)}\right)\left(\frac{{}_m}{{}_{m}B}\frac{J_m^{}\left({}_{m}B\right)}{J_m\left({}_{m}B\right)}-\frac{{}_b}{{}_{b}B}\frac{H_m^{}\left({}_{b}B\right)}{H_m\left({}_{b}B\right)}\right)=0.& \left(57\right)\end{array}$

Orthogonality

[0104] Orthogonality is an important property that permits projection. The transverse modes satisfy the orthogonality relation (19), revealing the simple form of the adjoint modes E.sub.m.sup. in the process. However, there are two subtleties associated with symmetry and longitudinal modes that shall be also discussed.

[0105] The proof of orthogonality is similar to other such proofs, proceeding from a construction that involves the operator which generates the eigenmodes,

[00045]$\begin{array}{cc}{E}_n^{}\left(r\right).Math.\left(r\right){\stackrel{=}{G}}_0\left(.Math.r-r^{}.Math.\right)\left(r^{}\right)E_m\left(r^{}\right){\mathrm{drdr}}^{}=E_n^{}\left(r\right).Math.\left(r\right){\stackrel{=}{G}}_0\left(.Math.r-r^{}.Math.\right)\left(r^{}\right)E_m\left(r^{}\right){\mathrm{dr}}^{}\mathrm{dr}=s_m{E}_n^{}\left(r\right).Math.\left(r\right)E_m\left(r\right)\mathrm{dr}.& \left(58\right)\end{array}$

E.sub.n.sup.(r) is a column vector. Considering the alternative way of evaluating the construction,

[00046]$\begin{array}{cc}{E}_n^{}\left(r\right).Math.\left(r\right){\stackrel{=}{G}}_0\left(.Math.r-r^{}.Math.\right)\left(r^{}\right)E_m\left(r^{}\right){\mathrm{drdr}}^{}=\left[{\stackrel{=}{G}}_0^{}\left(.Math.r-r^{}.Math.\right)\left(r\right)E_n^{}\left(r\right)\mathrm{dr}\right].Math.E_m\left(r^{}\right)\left(r^{}\right){\mathrm{dr}}^{}=s_n{E}_n^{}\left(r^{}\right).Math.\left(r^{}\right)E_m\left(r^{}\right){\mathrm{dr}}^{},& \left(59\right)\end{array}$

which follows because the dot product implies transposition, a.Math.ba.sup.b. [G.sub.0.sup.(|rr|)(r)E.sub.n(r)dr] may be simplified with one assumption,

[00047]$\begin{array}{cc}{E}_m^{}=E_m\left(r\right),& \left(60\right)\end{array}$

that adjoint is simply the original field. Then, [G.sub.0.sup.(|rr|)(r)E.sub.n(r)dr] satisfies the defining eigenvalue equation due to the symmetry of the Green's tensor both under transposition and the interchange of r and r,

[00048]$\begin{array}{cc}{\stackrel{=}{G}}_0^{}\left(.Math.r-r^{}.Math.\right)={\stackrel{=}{G}}_0\left(.Math.r-r^{}.Math.\right).& \left(61\right)\end{array}$

[0106] Combining these results, the following is obtained

[00049]$\begin{array}{cc}\left(s_n-s_m\right)E_n^{}\left(r\right).Math.\left(r\right)E_m\left(r\right)\mathrm{dr}=0,& \left(62\right)\end{array}$

so eigenvectors belonging to different eigenvalues are orthogonal. The form of the adjoint necessary to yield orthogonality has been also identified. These properties are guaranteed by the complex symmetric nature of the Green's tensor. However, (62) is silent on the orthogonality of degenerate modes that have the same eigenvalue. Indeed, such modes will not automatically be orthogonal, and it often needs to be enforced by some other means.

[0107] Two cases are relevant to this issue: degeneracy due to symmetry, and the set of longitudinal modes.

[0108] For modes that are degenerate due to symmetry, one way to obtain an orthogonal set is by blind application of say the Gram-Schmidt algorithm, using the integral in (62) as the metric. The resulting eigenmodes will be orthogonal, but typically will not have any recognisable symmetry properties upon visual inspection. This is not problematic for numerical purposes, but can hinder physical interpretation. Alternatively, one can consider the set of symmetry operations that leave the structure invariant. The eigenmodes E.sub.m of Green's tensor in (10) will simultaneously be eigenmodes of these symmetry operations. Just as the properties of the operator (10) impose an orthogonality relation on its eigenmodes, the symmetry operations impose an orthogonality relation on its eigenmodes. While the operator (10) is unable to establish orthogonality among a symmetry degenerate set, the symmetry operations can. In practice, the formal properties of the symmetry operations do not need to be analyzed.

[0109] Instead, orthogonality can be established, after a degenerate set of eigenmodes has been found, by choosing the appropriate linear combination, such that they are also eigenmodes of the symmetry operations. For more complex symmetries, the tools of group theory can be used.

[0110] Several types of symmetry such as continuous rotational symmetry, continuous translational symmetry, and discrete translational symmetry (also known as periodicity) may be considered. The eigenmodes of these symmetry operations are known to have the spatial variation e.sup.im, e.sup.ikz, and e.sup.ik.sup.B.sup.R, with the latter applicable to any lattice point R. However, a problem arises when these spatial variations are inserted into the integral of (62), and integration is performed along the direction of symmetry. The result is zero, for example, e.sup.im e.sup.im d=0, arriving at the paradoxical conclusion that each mode is orthogonal to itself. This can be remedied by using the {cos m, sin m} basis, instead. Since one wish to use the e.sup.im eigenmodes, their properties will be described.

[0111] All symmetry operators are orthogonal operators, satisfying the property .sub.|.sub.=|. This leads to the more familiar case of the adjoint mode being the complex conjugate, .sup.=*, or alternatively .sup.(r)=(r). This conflicts with the definition (60), so a hybrid definition of the adjoint becomes necessary for these three symmetries. For example, if the geometry possesses both continuous rotational symmetry in and continuous translational symmetry in z, the adjoint may be defined as

[00050]$\begin{array}{cc}{E}_m^{}\left(r,,z\right)=E_m\left(r,-,-z\right).& \left(63\right)\end{array}$

Longitudinal Basis Modes

[0112] Longitudinal basis modes are responsible for reproducing the near-field variations of the target modes at locations where the permittivity profile is changing. In particular, they reconstruct any field discontinuities due to step-index changes. Longitudinal basis modes are necessary because transverse basis modes satisfy Gauss' law, .Math.E=0, in regions of uniform permittivity. However, by definition, perturbations cause non-uniformities in the permittivity profile. Only the displacement field remains everywhere transverse, .Math.D=0, and E acquires a longitudinal component, .Math.0, wherever e varies. Thus, it is insufficient to rely solely on transverse basis modes E.sub.m in constructing the perturbed modes {tilde over (E)}.sub.. Since longitudinal modes are seldom used, their properties are discussed, along with a way for obtaining an appropriate minimal set of longitudinal modes.

Definitions and Properties

[0113] The focus is on the first set of modes that arises from (12), whereby E.sub.m=0, producing the longitudinal electric modes. Some of their properties are discussed, with the goal of demonstrating that they are non-divergence-free for at least one point in the inclusion interior. Firstly, rotationless vector fields can be expressed as a potential

[00051]$\begin{array}{cc}{E}_m={}_m.& \left(64\right)\end{array}$

[0114] Since .sub.m=0, there is infinite index contrast between the interior and exterior, and the field is identically zero outside. In order to satisfy the boundary conditions at the interface, the parallel component of the interior electric field must be zero, but the perpendicular component can be arbitrary. The first condition implies that

[00052]$\begin{array}{cc}{}_{.Math.}=E_{.Math.}=0,& \left(65\right)\end{array}$

along the boundary, so .sub.m is constant here. This is directly analogous to the external surface of a perfect conductor being an equipotential surface. Due to the freedom associated with the definition (64), any constant may be added to .sub.m without affecting E.sub.m. Thus, is eas shown that the modes satisfy Dirichlet boundary conditions

[00053]$\begin{array}{cc}{}_m{.Math.}_B=0& \left(66\right)\end{array}$

where B defines the interface between the interior and exterior.

[0115] The boundary conditions (66) are sufficient to demonstrate that these modes are longitudinal. Note that .sup.2.sub.m cannot be everywhere zero, otherwise E.sub.m would be a constant. Since .sup.2 .sub.m=.Math.E.sub.m cannot be everywhere zero, it may be concluded that E.sub.m must have a longitudinal component. Unlike the transverse modes, great freedom exists in the construction of longitudinal modes. Indeed their functional forms can be almost arbitrary, subject only to the two restrictions (64) and (66). Since the longitudinal modes all share the same zero eigenvalue, the modes E.sub.m are not in general orthogonal, as they are not subject to (62). There is also no defined differential operator for .sub.m which would furnish an orthogonality relation. This is a minor inconvenience, as an orthogonalization procedure such as Gram-Schmidt can be applied if necessary, using (19) as the metric.

[0116] An appropriate useful basis of longitudinal modes must be constructed for use in the expansion (15). Firstly, a general orthonormal basis can be constructed via a simple prescription for separable geometries, by demanding that .sub.m be orthogonal along each separable coordinate. This is considered later in the description, and leads to a basis based on generalized Fourier series. However, many basis modes may be required.

[0117] In particular, if sharp interfaces are present, the Fourier-Bessel basis is susceptible to Gibbs phenomenon (the peculiar manner in which the Fourier series of a piecewise continuously differentiable periodic function behaves at a jump discontinuity), resulting in numerical noise in the context of the numerical eigenvalue solution of (22). A more numerically robust minimal set can be obtained by considering the particular perturbation geometry, detailed further in the description. This results in an efficient and reliable basis, exhibiting exponential convergence.

Fourier-Bessel Basis

[0118] For circular and spherical geometries, an orthonormal basis of longitudinal modes can be constructed using the Fourier-Bessel basis, which is capable of representing any longitudinal field. As an example, a recipe for a 2D circular domain of radius B is provided, deriving the basis from first principles by demanding orthogonality. Polar coordinates, (r, ) are employed. In the angular direction, the continuous rotational symmetry furnishes the well known orthogonal set of cylindrical harmonics, e.sup.im. In the radial direction, orthogonality is expressed as

[00054]$\begin{array}{cc}{}_0^B{}_m{}_{n}rdr={}_{mn}& \left(67\right)\end{array}$

subject to the known boundary conditions that {.sub.m} all be zero at r=B and be bounded at r=0. This defines a set of orthogonal functions or a generalized Fourier series, where r is the weight function. The Sturm-Liouville eigenvalue problem associated with (67) is the Bessel differential equation. The overall solution may now be constructed,

[00055]$\begin{array}{cc}{}_{mn}=L_{mn}{J}_m\left(\frac{u_{m,n}}{B}r\right)e^{\mathrm{im}},& \left(68\right)\end{array}$

where two indexes n, m are now employed, the normalization constant is L.sub.mn, and u.sub.m,n is the n-th root of Bessel function J.sub.m(z). In the simple case of a circular geometry, the adjoint modes are given by the complex conjugate,

[00056]$\begin{array}{cc}{}_{mn}^{}={}_{\mathrm{mn}}^{*},& \left(69\right)\end{array}$

since only the angular part of (68) is complex. After normalization, the desired orthonormal set is obtained.

[0119] It was observed that the construction (68) obeys the Helmholtz equation with Dirichlet boundary conditions,

[00057]$\begin{array}{cc}{}^2{}_{mn}+{}_{mn}^2{}_{mn}=0,& \left(70\right)\end{array}$

where .sub.mn=u.sub.m,n/B. Furthermore, the modes defined by (68) can now be shown to be everywhere non-transverse by explicit computation,

[00058]$\begin{array}{cc}.Math.E_{mn}={}^2{}_{mn}=-{}_{mn}^2{}_{mn}0.& \left(71\right)\end{array}$

[0120] Finally, one may normalize the modes according to (19), which can be simplified to a scalar expression via

[00059]$\begin{array}{cc}\begin{array}{c}{E}_{mn}^{}.Math.E_{mn}dr={}_{mn}^{}.Math.{}_{mn}dr\\ =.Math.\left({}_{mn}^{}{}_{mn}\right)dr-{}_{mn}^{}{}^2{}_{mn}dr\\ ={}_{mn}^{}{}_{mn}.Math.\mathrm{dS}+{}_{mn}^2{}_{mn}^{}{}_{mn}dr\\ ={}_{mn}^2{}_{mn}^{}{}_{mn}\mathrm{dr}.\end{array}& \left(72\right)\end{array}$

[0121] The surface integral vanishes due to the boundary condition. The radial part of the final expression can be evaluated analytically using the defining orthogonality relation,

[00060]$\begin{array}{cc}{}_0^1{J}_l\left(u_{l,m}r\right)J_l\left(u_{l,n}r\right)rdr=\frac{{}_{mn}}{2}{J}_{l+1}^2\left(u_{l,m}\right),& \left(73\right)\end{array}$

to yield

[00061]$\begin{array}{cc}{L}_{mn}={\left[\sqrt{}{}_{mn}{\mathrm{BJ}}_{m+1}\left({}_{mn}B\right)\right]}^{-1}.& \left(74\right)\end{array}$

[0122] Unfortunately, a large number of Fourier-Bessel modes (68) is necessary for piecewise uniform structures, an unavoidable fact when a set of smooth functions is used to expand discontinuities. This causes two problems for numerical eigenvalue routines during the solution of (22). Firstly, many spurious modes {tilde over (E)}.sub. are liable to be produced, with eigenvalue {tilde over ()}.sub.=0. These are longitudinal modes of the perturbed structure, so as discussed their modal fields are arbitrary and have little physical relevance for the modal expansion of (10). These can be easily identified and discarded. Unfortunately, the numerical solution of (22) often also results in some spurious modes that do not have {tilde over ()}.sub.=0, probably due to numerical noise. Though perhaps easy to identify by eye via visual inspection of the modal fields, distinguishing such modes through automated methods can be difficult. Inclusion of such modes in the expansion of (10) can have disastrous consequences. Thus, using the Fourier-Bessel longitudinal modes is limited.

[0123] Secondly, even the useful modes {tilde over (E)}.sub. are negatively impacted by numerical noise and sensitivity during the solution of the eigensystem (22). This manifests as highly speckled modal fields {tilde over (E)}.sub., which again is detrimental to the expansion (10). The numerical issues stem from any inaccuracy in evaluating matrix elements (21), but even if these elements have been evaluated analytically, there is also noise due to the numerical evaluation of the eigenvectors of (22). Particularly troublesome is any inaccuracy when calculating the coefficients of Fourier-Bessel modes with high spatial frequencies, as these introduce considerable numerical noise into modal fields {tilde over (E)}.sub.. For these reasons, a more appropriate basis of longitudinal modes is essential for piecewise uniform structures.

Boundary Adapted Longitudinal Mode Basis

[0124] An efficient and reliable basis of longitudinal modes can be derived for piecewise uniform scatterers. As discussed in the introduction, modes are longitudinal only in regions where the permittivity profile is changing, which is restricted to a step change at the boundary of piecewise uniform scatterers. To demonstrate, an explicit expression for the divergence of a mode {tilde over (E)}.sub. was obtained, via

[00062]$\begin{array}{cc}.Math.{\tilde{D}}_{}=.Math.\left(\left(r\right){\tilde{E}}_{}\right)=\left(r\right).Math.{\tilde{E}}_{}+{\tilde{E}}_{}.Math.\left(r\right)=0.& \left(75\right)\end{array}$

[0125] Now, assuming (r) is nonzero only where the permittivity profile is discontinuous by , it was found that .Math.{tilde over (E)}.sub. is proportional to the field {tilde over (E)}.sub. and ,

[00063]$\begin{array}{cc}.Math.{\tilde{E}}_{}=-\frac{1}{\left(\right)}{\tilde{E}}_{}.Math.\left(r\right)=-\frac{}{\left(\right)}{\tilde{E}}_{}.Math.n_p\left(r-r_p\right),& \left(76\right)\end{array}$

where r.sub.p defines the perturbation boundary and n.sub.p is the normal vector.

[0126] (76) was introduced for computational purposes, but it is opaque to interpret. Its meaning becomes clear if instead the integral form, {tilde over (D)}.sub..Math.dS=0 is considered. This naturally leads to the condition that {tilde over (E)}.sub. is discontinuous by the ratio of permittivities across a step interface,

[00064]$\begin{array}{cc}{n}_p.Math.\left({}_b{\tilde{E}}_{}^{+}-{}_{}{\tilde{E}}_{\stackrel{}{}}\right)=0,& \left(77\right)\end{array}$

where {tilde over (E)}.sub..sup. and {tilde over (E)}.sub..sup.+ are the fields on the interior and background side of the boundary, respectively. Thus, a localized divergence (76) merely corresponds to a discontinuity in {tilde over (E)}.sub. across the boundary. Furthermore, it is also apparent that the very purpose of longitudinal modes is to reproduce this discontinuity, since the transverse modes are incapable of doing so. Thus, the aim is to introduce a set of intrinsincally discontinuous modes for this purpose. These all necessarily have a longitudinal component similar to (76).

[0127] In principle, only a single longitudinal mode defined by (76) subject to the Dirichlet boundary condition (66) is sufficient. But, neither {tilde over (E)}.sub. nor AE are known prior to the solution of the eigensystem (22). Thus, instead is defined a suitably complete basis of longitudinal modes {E.sub.l} for expanding {tilde over (E)}.sub. so that the appropriate linear combination satisfies (76). Indeed, {E.sub.l} should be capable of expanding any member of the set {{tilde over (E)}.sub.}. Since (76) is only non-zero along the boundary, {E.sub.l} only needs to capture the field variation along the boundary, corresponding to an expansion of the factor multiplying (rr).

[0128] As a simple example, for a 2D inclusion with an interface given by (a.sub.p(),) in polar coordinates, one may use the 1D set

[00065]$\begin{array}{cc}.Math.E_l={}^2{}_l=\frac{1}{2r}\left(r-a_p\left(\right)\right)e^{\mathrm{il}},& \left(78\right)\end{array}$

where l is some integer angular order. This corresponds to Dirac delta function in the radial direction, with a factor related to the Jacobian, while the field variation along the boundary is expanded via a Fourier series in the angular direction. This simple definition is often more than adequate for simple convex structures where a.sub.p() is differentiable to all orders. Rapid, possibly exponential, convergence in l is possible, since Fourier series are known to converge exponentially for smooth functions on periodic domains.

[0129] For more complex structures, a simpler definition for the set {.sub.l} may be used,

[00066]$\begin{array}{cc}.Math.E_l={}^2{}_l=w_l\left(r\right){}_{}\left(r\right),& \left(79\right)\end{array}$

where is defined to be the perturbation boundary, which is a closed surface in 3D and a closed contour in 2D. The special Dirac-delta function (r) was employed, which has the property of having uniform weight over all of . The advantage of uniform weight, and thus (79) over (78), becomes apparent for contours with segments parallel to the direction. The (r) converts integrals over all space to an integral over ,

[00067]$\begin{array}{cc}f\left(r\right){}_{}\left(r\right)\mathrm{dr}={}_{}f\left(r\right)\mathrm{dS},& \left(80\right)\end{array}$

reducing the integration order by one.

[0130] This is the most useful property of (r) for the present invention's purposes, and in practice makes (r) quite convenient to handle, as will be showed. Finally, the set {w.sub.l(r)} should be chosen to obtain a complete and rapidly converging basis along , and indeed only needs to defined along . As an example, a contour with discontinuities can be partitioned into piecewise segments, expanding each segment using Chebyshev polynomials, which are known for their convergence properties.

[0131] For illustration purposes, a concrete evaluation of the somewhat abstract definition (79) is given, exploiting (80) to simplify (79). For example, the fundamental solution of the Laplace operator is defined as the Green's function

[00068]$\begin{array}{cc}{}^2{G}_0\left(r,r^{}\right)=\left(r-r^{}\right),& \left(81\right)\end{array}$

which enables .sub.l(r) to be found via an integral over ,

[00069]$\begin{array}{cc}\begin{array}{c}{}_l\left(r\right)=G_0\left(r,r^{}\right)w_l\left(r^{}\right){}_{}\left(r^{}\right){\mathrm{dr}}^{}\\ ={}_{}{G}_0\left(r,r^{}\right)w_l\left(r^{}\right){\mathrm{dS}}^{}\end{array}.& \left(82\right)\end{array}$

[0132] The evaluation of (82) is simplest if is parameterized. For a 2D domain, the 1D contour can be specified as a vector function of the parameter t

[00070]$\begin{array}{cc}\left(t\right)\left(x\left(t\right),y\left(t\right)\right),& \left(83\right)\end{array}$

as well as the weight, w.sub.l(t). Finally, .sub.l(r) can be obtained from

[00071]$\begin{array}{cc}{}_l\left(r\right)=G_0\left(r,\left(t\right)\right)w_l\left(t\right).Math.\stackrel{.}{}\left(t\right).Math.\mathrm{dt},& \left(84\right)\end{array}$

integrating over the appropriate range of t and using numerical quadrature if necessary. Note that |{dot over ()}(t)| represents the infinitesimal arc length. The numerical details of this procedure are presented later in the description. If is composed of several piecewise smooth segments, then (t) can also be decomposed into n segments {.sub.1(t), .sub.2(t), . . . , .sub.n(t)}, each specified over a different range of t, and (84) can be integrated segment by segment. Also, a decomposition of is numerically expedient if is well approximated by cubic splines.

[0133] Equation (79) defines a Poisson equation which, along with Dirichlet boundary conditions, yields a unique solution for each order 1. Since each mode E.sub.l is inherently discontinuous, it is capable of replacing many hundreds or thousands of Fourier-Bessel basis modes (68), which suffer from slow convergence due to Gibbs phenomenon. An even more important consideration than efficiency is the numerical robustness of {.sub.l} for the solution of (22). Since there are much fewer degrees of freedom than in the Fourier-Bessel basis, this highly-constrained basis has limited freedom for generating numerical noise, greatly ameliorating numerical sensitivity issues. Also critical is that this basis entirely eliminates the problem of spurious modes, since by construction the modes have zero longitudinal component everywhere except the perturbation boundary.

[0134] One of the few disadvantages of boundary adapted longitudinal modes is that the resulting basis {.sub.l} is not in general orthogonal with respect to the inner product (19), thus an orthogonalization procedure is usually first necessary. Lwdin's method was employed, which is briefly summarized later.

[0135] The boundary adapted longitudinal modes are valid for all material parameters and wavelengths, and can even be adapted to inclusions that have both sharp and smooth spatial variations of permittivity. In comparison, the boundary element method requires the frequency and material parameters to be specified, and is applicable only to piecewise uniform structures. Thus, the boundary adapted longitudinal modes are more flexible and tends not to lead to nonlinear eigenvalue problems.

[0136] Numerical implementation are discussed later on in the description, beginning with normalization and perturbation integrals. Then, two methods of generating the boundary adapted longitudinal modes defined by (79) are presented. The first method is conceptually simple, and involves obtaining .sub.l using the Green's function for the Laplace operator and numerically integrating, described further in the description. This is a fast, accurate, and reliable method, but it involves some numerical intracies due to the divergence of the Green's function. It is suitable for implementation in production quality code. The second method involves expanding .sub.l using the Fourier-Bessel basis, described further in the description. Using the Fourier-Bessel basis to expand the boundary adapted longitudinal modes is much more numerically reliable than direct use of the Fourier-Bessel basis, and this fact is indirect proof of the strength of the boundary adapted basis. The expansion is numerically simple, but still suffers from low accuracy due to Gibbs phenomenon and slow speed due to the need to evaluate many integrals. Thus, it should only be used for testing and comparison purposes. An entirely numerical solution of (79) is possible, for example, via the finite element method. Indeed, this would constitute a simple application of the finite element method, and avoids the specialized techniques necessary for the Green's function approach, for example.

[0137] Normalization and perturbation integrals of boundary adapted longitudinal modes Another advantage of boundary adapted longitudinal modes is the simple numerical evaluation of normalization (19) and perturbation integrals (21). For normalization, consider the integral

[00072]$\begin{array}{cc}\begin{array}{c}{N}_{\mathrm{kl}}{E}_k^{}.Math.E_l\mathrm{dr}={}_k^{}.Math.{}_l\mathrm{dr}\\ =.Math.\left({}_k^{}{}_l\right)\mathrm{dr}-{}_k^{}{}^2{}_l\mathrm{dr}\\ ={}_{}{}_k^{}{}_l.Math.\mathrm{dS}-{}_k^{}{w}_l{}_{}\left(r\right)\mathrm{dr}\\ =-{}_{}{}_k^{}{w}_l\mathrm{dS}\end{array}.& \left(85\right)\end{array}$

[0138] The surface integral after the third equality vanishes since all longitudinal modes are zero along the boundary by definition, while the final integral is due to the definition of (r), (80). The final integral is compatible with standard numerical quadrature techniques. Due to the Dirac-delta function, the order of the integral is always reduced by one. The normalization integral is obtained for k=l, while the general case of kl yields modal overlaps.

[0139] The matrix of overlap integrals N calculated in (85) allows normalization and orthogonalization using Lwden's method, which treats all basis vectors equally and preserves any symmetry, unlike the Gram-Schmidt procedure. Despite its simplicity, it is usually only treated in the computational chemistry field and may not be commonly known, so a brief exposition is provide. It exploits the structure of N, which is Hermitian and positive definite, so long as the basis is linearly independent. The normalized longitudinal modes D.sub.m are obtained by transforming the vectors by the matrix square root of N,

[00073]$\begin{array}{cc}{D}_m={.Math.}_k{E}_k{N}_{\mathrm{km}}^{-\frac{1}{2}}.& \left(86\right)\end{array}$

[0140] An explicit computation demonstrates the desired properties,

[00074]$\begin{array}{cc}\begin{array}{c}{D}_m^{}.Math.D_n\mathrm{dr}={.Math.}_k{.Math.}_l{N}_{\mathrm{km}}^{*-\frac{1}{2}}{N}_{\ln }^{-\frac{1}{2}}{E}_k^{}.Math.E_l\mathrm{dr}\\ ={.Math.}_k{.Math.}_l{N}_{\mathrm{km}}^{*-\frac{1}{2}}{N}_{\mathrm{kl}}{N}_{\ln }^{-\frac{1}{2}}\\ ={}_{\mathrm{mn}}\end{array},& \left(87\right)\end{array}$

using the symmetry of

[00075]$N^{-\frac{1}{2}},$

which is inherited from N.

[0141] (85) avoids treating the discontinuity in E.sub.l, since the result involves .sub.l, which is everywhere continuous. This discontinuity is first evaluated at all points along the perturbation boundary. This can be derived directly from the defining equation (79) by integrating over an infinitesimal Gaussian pillbox that straddles the boundary,

[00076]$\begin{array}{cc}\begin{array}{c}{}^2{}_l\mathrm{dr}=w_l\left(r\right){}_{}\left(r\right)\mathrm{dr}\\ {}_l.Math.\mathrm{dn}=w_l\left(r\right)\mathrm{dS}\\ \left({}_l^{+}-{}_l^{-}\right).Math.\hat{n}\mathrm{dS}=w_l\left(r\right)\mathrm{dS}\\ \left(E_l^{+}-E_l^{-}\right).Math.\hat{n}=w_l\left(r\right)\end{array}& \left(88\right)\end{array}$

where E.sub.l.sup.+E.sub.l.sup. is the difference E.sub.l across the boundary, and {circumflex over (n)} is the unit vector normal to the interface. Therefore, the magnitude of the discontinuity is given by the magnitude of the Dirac-delta function. As for the parallel component of E.sub.l, the intuitive result that its discontinuity is zero can be derived using similar arguments starting from the identity

[00077]$\begin{array}{cc}{}_l=0.& \left(89\right)\end{array}$

the matrix element (21) is now treated,

[00078]$\begin{array}{cc}\begin{array}{c}{V}_{\mathrm{kl}}{E}_k^{}.Math.\left(r\right)E_l\mathrm{dr}={}_k^{}.Math.{}_l\mathrm{dr}\\ =.Math.\left({}_k^{}{}_l\right)\mathrm{dr}-{}_k^{}{}^2{}_l\mathrm{dr}\\ ={}_{}{}_k^{}{E}_l.Math.\mathrm{dS}\end{array},& \left(90\right)\end{array}$

where the surface integral is now along the interior boundary of (r). This time, the integral involving .sup.2.sub.l, vanishes, since .sup.2.sub.l=0 within the region (r), which is an open set. This means there is a need to evaluate the integral along the interior boundary of the perturbation so E.sub.l by E.sub.l.sup. are being replaced. Then by invoking (88), a particularly convenient result follows,

[00079]$\begin{array}{cc}\begin{array}{c}{V}_{\mathrm{kl}}={}_{}{}_k^{}\left(E_l\mathrm{ave}-\frac{E_l}{2}\right).Math.\mathrm{dS}\\ ={}_{}{}_k^{}{E}_l\mathrm{ave}.Math.\mathrm{dS}-\frac{1}{2}{}_{}{}_k^{}{w}_l\mathrm{dS}\end{array},& \left(91\right)\end{array}$

where

[00080]$E_l\mathrm{ave}=\frac{1}{2}\left(E_l^{+}+E_l^{-}\right)$

is the sum of the two discontinuous fields. It will be shown now that an efficient technique for obtaining E.sub.lave exists in the following parts of the description. Conveniently, the second term in (91) is precisely the overlap integral obtained in (85).

Fourier-Bessel Expansion of Boundary Adapted Longitudinal Modes

[0142] The Fourier-Bessel basis is used to expand (79),

[00081]$\begin{array}{cc}{}_l={.Math.}_{n,m}{c}_{\mathrm{mn}}{}_{\mathrm{mn}}.& \left(92\right)\end{array}$

[0143] The coefficients can be evaluated via expansion and projection,

[00082]$\begin{array}{cc}\begin{array}{c}{.Math.}_{n,m}{c}_{\mathrm{mn}}{}^2{}_{\mathrm{mn}}={}^2{}_l\\ -{.Math.}_{n,m}{c}_{\mathrm{mn}}{}_{\mathrm{mn}}^2{}_{\mathrm{mn}}=w_l\left(r\right){}_{}\left(r\right)\\ -{.Math.}_{n,m}{c}_{\mathrm{mn}}{}_{\mathrm{mn}}^2{}_{\mathrm{kl}}^{}{}_{\mathrm{mn}}\mathrm{dr}={}_{\mathrm{kl}}^{}{w}_l\left(r\right){}_{}\left(r\right)\mathrm{dr}\\ c_{\mathrm{mn}}=-{}_{}{w}_l{}_{\mathrm{mn}}^{}\mathrm{dS}\end{array}& \left(93\right)\end{array}$

where the overlap integral between modes was evaluated using (72) assuming the normalization (19). A similar derivation can be performed for the adjoint modes .sub.l.sup., obtaining expansion coefficients c.sub.mn.sup..

[0144] Normalization of .sub.l is simple if the basis .sub.mn is already normalized. Instead of (85), the integral can be obtained from a simple sum,

[00083]$\begin{array}{cc}\begin{array}{c}{}_l^{}.Math.{}_l\mathrm{dr}={.Math.}_{\mathrm{kl}}{.Math.}_{\mathrm{mn}}{c}_{\mathrm{kl}}^{}{c}_{\mathrm{mn}}{}_{\mathrm{kl}}^{}.Math.{}_{\mathrm{mn}}\mathrm{dr}\\ ={.Math.}_{\mathrm{mn}}{c}_{\mathrm{mn}}^{}{c}_{\mathrm{mn}}\end{array}.& \left(94\right)\end{array}$

[0145] If desired, the integrals of the matrix elements (21) involving .sub.l can also be evaluated in the basis of .sub.mn, for example

[00084]$\begin{array}{cc}{}_l^{}.Math.\left(r\right){}_t\mathrm{dr}={.Math.}_{k,l}{.Math.}_{n,m}{c}_{\mathrm{kl}}^{}{c}_{\mathrm{mn}}{}_{\mathrm{kl}}^{}.Math.\left(r\right){}_{\mathrm{mn}}\mathrm{dr}.& \left(95\right)\end{array}$

[0146] Expanding boundary adapted longitudinal modes using the Fourier-Bessel has some disadvantages. Firstly, the resulting boundary adapted longitudinal modes are susceptible to Gibbs phenomenon. This limits the accuracy of the matrix elements (95), and ultimately the accuracy of the perturbed modes {tilde over (E)}.sub.. One possible visible manifestation of this inaccuracy is a halo in the {tilde over (E)}.sub. field at r=B, the basis boundary. This halo is a numerical artifact resulting from incorrect contributions among the various {.sub.l}. Better methods exist for finding the boundary adapted longitudinal modes.

Green's Function Solution of Boundary Adapted Longitudinal Modes

[0147] Further details of evaluating (79) are now supplied, directly using the Green's function for the Laplace operator, (81). In particular, the explicit form of the 2D Green's function when subject to the Dirichlet boundary condition G.sub.D(r=B,r)=0 can be found by the method of images,

[00085]$\begin{array}{cc}{G}_D\left(r,r^{}\right)=\frac{1}{4}\log \left(\frac{{.Math.r-r.Math.}^2{B}^2}{{.Math.r-r.Math.}^{2}r{}^2}\right),& \left(96\right)\end{array}$

where r is the location of the source and r is the location of the image source,

[00086]$\begin{array}{cc}{r}^{}=\frac{B^2}{r{}^2}{r}^{}.& \left(97\right)\end{array}$

[0148] To obtain .sub.l, integration is being performed over all source terms, as in (82).

[0149] The evaluation of E.sub.l(r) is also necessary, for the matrix elements (90), for example.

[0150] The gradient of G.sub.D(r,r) is thus derived. To simplify the notation, the following two vectors are introduced

[00087]$\begin{array}{cc}{R}^{}r-r^{},R^{}r-r^{}.& \left(98\right)\end{array}$

[0151] By chain rule, the following is obtained

[00088]$\begin{array}{cc}\begin{array}{c}{G}_D\left(r,r^{}\right)={}_R{G}_D\left(r,r^{}\right)R^{}+{}_R{G}_D\left(r,r^{}\right)R^{}\\ =\frac{1}{2}\left[\frac{r-r}{{.Math.r-r.Math.}^2}-\frac{r-r}{{.Math.r-r.Math.}^2}\right]\end{array},& \left(99\right)\end{array}$$\mathrm{since}$$\begin{array}{cc}{}_R{G}_D\left(r,r^{}\right)=\frac{1}{R},R^{}=\frac{R}{R}.& \left(100\right)\end{array}$

[0152] Correspondingly,

[00089]$\begin{array}{cc}{E}_l\left(r\right)={}_l\left(r\right)={}_{}{G}_D\left(r,r^{}\right)w_l\left(r^{}\right){\mathrm{dS}}^{}.& \left(101\right)\end{array}$

[0153] In principle, (82) and (101) are sufficient to obtain the modes and perform all necessary integrals, but in practice, their evaluation is often problematic if attempted using naive methods. Difficulties arise if r sits on the perturbation boundary, such that r and r coincide, since the Green's function kernel (96) has a logarithmic singularity. Yet, it is vital to obtain .sub.l and E.sub.l along the perturbation boundary, as the most efficient methods for evaluating integrals such as (85) rely on these values. Fortunately, efficient and reliable quadrature methods exist for the evaluation integrable singularities. Meanwhile, (99) has a simple pole singularity, so (101) is not integrable and can only defined in the Cauchy principal value sense. One might then wonder whether or not (101) yields a physical quantity. Actually, the Cauchy principal value represents the average of the two discontinuous values across the boundary E ave(r), required for evaluating (91). This result is a consequence of the Plemelj formula. Furthermore, since the difference is also known from (88), the individual values E.sup.+ and E.sup. can be deduced. Efficient and reliable numerical evaluation of these integrals is treated in later in the description.

[0154] Finally, the winding number of a closed curve around a point r is introduced, which is useful for distinguishing whether r is interior or exterior. While this quantity is not essential for numerical implementation, it is of great aid for sanity checks. It can be expressed as

[00090]$\begin{array}{cc}W=\frac{1}{2}{}_{}d=\frac{1}{2}{}_{}\frac{r-r}{{.Math.r-r.Math.}^2}.Math.d\hat{n},& \left(102\right)\end{array}$

where is the angle between rr and some reference vector, such as {circumflex over (x)}. The obvious similarity to (99) hints at the deep underlying connection. The winding number is applicable to any curve, though simpler, faster criteria are available for simpler curves.

Numerical Quadrature and Singular Integrands

[0155] Numerical quadrature and singular integrands treats the quadrature techniques for singular integrals along 1D intervals. Evaluation is simplest if a parameterization (t) of the interface is known, as given by (83), and furthermore, if this parameterization is regular. This reduces (82) and other relevant integrals to 1D integrals along t in the form of (84), for example. Three quadrature techniques are presented, leading to three types of integrands: smooth functions without singularities, functions with logarithmic singularities, and Cauchy principle value integrals. All are variations of Gaussian quadrature, which are among the most efficient techniques.

[0156] As a preliminary step, all integrals along the curve are converted to an integral in t. This was done for (84), and for a general scalar integrand one has

[00091]$\begin{array}{cc}{}_{}f\left(r\right)\mathrm{ds}=f\left(\left(t\right)\right).Math.\stackrel{.}{}\left(t\right).Math.\mathrm{dt}.& \left(103\right)\end{array}$

[0157] Meanwhile, a surface integral of the form (90) becomes

[00092]$\begin{array}{cc}{}_{}f\left(r\right).Math.d\hat{n}=f\left(\left(t\right)\right).Math.\left(\stackrel{.}{}\left(t\right)\hat{z}\right)\mathrm{dt},& \left(104\right)\end{array}$

where {circumflex over (n)} is outward normal to the curve and {circumflex over (z)} is normal to the plane of the curve. d{circumflex over (n)}={dot over ()}(t){circumflex over (z)} dt is obtained for a positively oriented curve because {dot over ()}(t) is always tangential to a normally parameterized curve (t). These integrals are readily evaluated by Gaussian quadrature if the integrands are nowhere singular, so

[00093]$\begin{array}{cc}{}_{-1}^{1}f\left(t\right)\mathrm{dt}={.Math.}_{n=1}^N{w}_{n}f\left(t_n\right),& \left(105\right)\end{array}$

where t.sub.n is the nth root of the Legendre polynomial P.sub.N(t) and the weights w.sub.n are

[00094]$\begin{array}{cc}{w}_n=\frac{2}{\left(1-t_n^2\right){\left(P{}_N\left(t_n\right)\right)}^2}.& \left(106\right)\end{array}$

[0158] For the Cauchy principal value integral of (101), with a singularity at t=t that lies within the domain of integration, a modified Gaussian quadrature rule exists,

[00095]$\begin{array}{cc}0=--00=--00=--00=--0{}_{-1}^1\frac{f\left(t\right)}{t-t}\mathrm{dt}-\frac{2f\left(t\right)Q_N\left(t\right)}{P_N\left(t\right)}+{.Math.}_{n=1}^N\frac{v_{n}f\left(t_n\right)}{t_n-t},& \left(107\right)\end{array}$

where Q.sub.N(x) is the Legendre function of the second kind, and the weight functions are now

[00096]$\begin{array}{cc}{v}_n=-\frac{2}{\left(N+1\right)P{}_N\left(t_n\right)P_{N+1}\left(t_n\right)}.& \left(108\right)\end{array}$

f(t) may be viewed as the residue of the function g(t)f(t)/(tt) at t=t.

[0159] For the integration of (101), this residue may be found via series expansion. It is sufficient to consider only the singular term (rr)/|rr|.sup.2, since the non-singular term (rr)/|rr|.sup.2 can be efficiently integrated by standard Gaussian quadrature. Using the expansion (t)=(t)+{dot over ()}(t)(tt)+{umlaut over ()}(t)(tt).sup.2/2+ . . . , the following is obtained

[00097]$\begin{array}{cc}\frac{\left(t\right)-\left(t^{}\right)}{{.Math.\left(t\right)-\left(t^{}\right).Math.}^2}\frac{\stackrel{.}{}\left(t^{}\right)}{{.Math.\stackrel{.}{}\left(t^{}\right).Math.}^2}\frac{1}{t-t^{}},& \left(109\right)\end{array}$

where higher orders terms in (tt) have been neglected. This is always possible since {dot over ()}(t) is never zero for a regular parameterization of a curve. Thus, the integration of the singular part of (99) over some segment proceeds by inserting into (103) and (107) to yield

[00098]$\begin{array}{cc}{}_{}\frac{r-r^{}}{{.Math.r-r^{}.Math.}^2}{\mathrm{ds}}^{}2\frac{\stackrel{.}{}\left(t\right)}{.Math.\stackrel{.}{}\left(t\right).Math.}\frac{Q_N\left(t\right)}{P_N\left(t\right)}+{.Math.}_{n=1}^N{v}_n\frac{\left(t\right)-\left(t_n\right)}{{.Math.\left(t\right)-\left(t_n\right).Math.}^2}.Math.\stackrel{}{}\left(t_n\right).Math.,& \left(110\right)\end{array}$

where (t) is parameterization of extending from t=1 to 1 and (t) is the detector point located on .

[0160] A modified Gaussian quadrature rule also exists for integrands with a logarithmic singularity at an endpoint of the integration interval. In our case, the logarithmic singularity is usually within the integration interval, so the interval must first be split into two segments on which the quadrature rule is separately applied. Unlike the above quadrature rules, closed form expressions for the weights and nodes are no longer available, and they must be obtained via a numerical procedure. Fortunately, values have already been computed and available in the literature, so there is no need to perform this numerical procedure.

Evaluation of Perturbation Integrals

[0161] Once all basis modes have been prepared, including the transverse modes and longitudinal modes, the matrix elements of (21),

[00099]$\begin{array}{cc}{E}_s^{}.Math.E_t\left(r\right)\mathrm{dr},& \left(111\right)\end{array}$

are ready to be evaluated. These integrals over the perturbation typically require numerical integration which can be computationally intensive, especially for higher-order modes with high spatial frequencies. To reduce the computational burden, it is possible to use vector identities to reduce the order of integration by one for piecewise uniform perturbations. This technique was already utilized in (90), for example. Such techniques were unnecessary in previous implementations of the perturbation method because only cases where the integrals could be calculated analytically were considered.

[0162] The simplest case to consider is the integral between transverse modes and longitudinal modes. Let E.sub.s be the longitudinal mode and E.sub.t the transverse. Firstly, the electric field E.sub.t of the transverse mode is expressed in terms of its magnetic field using the Maxwell curl equation (55), while the longitudinal mode is E.sub.s=.sub.s. Since (r) is assumed to be a piecewise uniform function, one may decompose the domain of integration into one or more disjoint subdomains that is considered separately, and apply the vector identity

[00100]$\begin{array}{cc}\begin{array}{c}{E}_s^{}.Math.H_{t}dr=.Math.\left(H_t{}_s^{}\right)dr+H_t.Math.{}_s^{}dr\\ =\left(H_t{E}_s^{}\right).Math.\mathrm{dn}\end{array}.& \left(112\right)\end{array}$

[0163] This integral is simple to treat because the second term on the RHS of (112) is identically zero.

[0164] For 2D geometries, the integral can be split into two cases and be further simplified. If the transverse field is of TM polarization, then its in-plane electric field components E.sub. are zero, while the only nonzero electric field components of longitudinal modes are in-plane. Thus, these integrals are always identically zero. If the transverse field is of TE polarization, then H=H.sub.z{circumflex over (z)}, so

[00101]$\begin{array}{cc}\begin{array}{c}\left(H_{z,t}\stackrel{}{z}{E}_s^{}\right).Math.\mathrm{dn}=-\left(H_{z,t}\stackrel{}{z}{E}_s^{}\right).Math.\stackrel{}{z}ds\\ =-H_{z,t}{E}_s^{}.Math.\mathrm{ds}\end{array}.& \left(113\right)\end{array}$