Pseudo-spectral analysis of radially-diagonalized Maxwell ’ s equations in cylindrical co-ordinates

We present a robust and accuracy enhanced method for analyzing the propagation behavior of EM waves in z−periodic structures in (r,φ,z)−cylindrical co-ordinates. A cylindrical disk, characterized by the radiusa and the periodicity lengthLz, defines the fundamental cell in our problem. The permittivity of the dielectric inside this cell is characterized by an arbitrary, single-valued function ε(r,φ,z) of all three spatial co-ordinates. We consider both open and closed boundary problems. Irrespective of the type of the boundary conditions on the surface = a, our method requires the discretization of the fields in the interior of the disk only. Inside the disk volume, we expand the fields in terms of planewaves on discrete cylindrical surfaces r i = i∆, with ∆ being the discretization step length. The fields on the nested surfaces r i = i∆ in the interior of the simulation domain are interrelated by the application of a simple, yet, powerful finite difference scheme. In free space outside the disk, the fields are expanded in terms of the closed-form eigensolutions of the Maxwell’s equations in cylindrical co-ordinates. In order to uniquely determine the involved unknown coefficients, the solutions in the interiorand exterior domains are matched on the disk’s bounding surface. Our formulation utilizes a radially-diagonalized form of Maxwell’s equations, and merely involves four (out of the six) field components. It is demonstrated that our formulation is perfectly suited, but by no means limited, to cylindrical symmetric problems. © 2003 Optical Society of America OCIS codes: (000.4430) Numerical approximation and analysis, (000.3860) Mathematical methods in physics References and links 1. Joannopoulos J. D., Meade R. D., Winn J. N., ‘‘Photonic Crystals: Molding the Flow of Light,’’ Princeton, September 1995. 2. Knight J. C., Birks T. A., Russell P. St. J., Atkin D. M., ‘‘All-silica single-mode optical fiber with photonic crystal cladding,’’ Opt. Lett.21, 1547-1549 (1996). 3. Cucinotta A., Selleri A., Vincetti L., Zoboli M., ‘‘Holey Fiber Analysis Through the Finite-Element Method,’’ IEEE Photon. Technol. Lett. 3, 147-149 (1991). 4. Chan, C.T, Yu Q. L., Ho K. M., ‘‘Order-N spectral method for electromagnetic waves,’’ Phys. Rev. B 51, 1663516642 (1995). (C) 2003 OSA 17 November 2003 / Vol. 11, No. 23 / OPTICS EXPRESS 3048 #2988 $15.00 US Received Septebmer 04, 2003; Revised October 31, 2003 5. Li Z.-Y., Lin L.-L., ‘‘Photonic band structures solved by a plane-wave-based transfer-matrix method,’’ Phys. Rev. E67, 046607 (2003). 6. Merle Elson J., Tran P., ‘‘Dispersion in photonic media and diffraction from gratings: a different modal expansion for the R-matrix propagation technique, ’’ J. Opt. Soc. Am. A 12, 1765-1771 (1995). 7. Li L., ‘‘Multilayer modal method for diffraction gratings of arbitrary profile, depth, and permittivity,’’ J. Opt. Soc. Am. A,10, 2581-2591 (1993). 8. Johnson S. G., Joannopoulos J.D., ‘‘Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,’’ Opt. Express8, 173-190 (2001) http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173 . 9. Bienstman P., Rigorous and Efficient Modelling of Wavelength Scale Photonic Components , doctoral disseration 2001, Ghent University, Belgium. Available online at http://photonics.intec.rug.ac.be/download/phd 104.pdf 10. Li Z.-Y., Ho K.-M., ‘‘Analytic modal solution to light propagation through layer-by-layer metallic photonic crystals,’’ Phys. Rev. B. 67, 165104 (2003). 11. Varis K., Baghai-Wadji A. R., ‘‘Hybrid Planewave/Finite Difference Transfer Method for Solving Photonic Crystals in Finite Thickness Slabs,’’ EDMO, November 15-16, 2001, Vienna, Austria, pp. 161-166. 12. Varis K., Baghai-Wadhi A.R., ‘‘Z-diagonalized Planewave/FD Approach for Analyzing TE Polarized Waves in 2D Photonic Crystals,’’ ACES, March 24-28, 2003, Monterey, CA, USA. 13. Varis K., Baghai-Wadji A. R., ‘‘A Novel 2D Pseudo-Spectral Approach of Photonic Crystal Slabs,’’ ACES Special Issue, (submitted). 14. Varis K., Baghai-Wadji A.R., ‘‘A Novel 3D Pseudo-spectral Analysis of Photonic Crystal Slabs,’’ ACES Special Issue, (submitted). 15. Liu S., Li L. W., Leong M. S., Yeo T. S., ‘‘Theory of Gyroelectric Waveguides,’’ PIER, 29, 231-259 (2000). 16. Freund R. W., ‘‘A Transpose-Free Quasi-Minimal Residual Algorithm for Non-Hermitian Linear Systems,’’ SIAM Journal of Scientific Computing, 14, 470-482 (1993). 17. Nachtigal N., Freund R. W., Reeb J. C., ‘‘QMRPACK user’s guide,’’ ORNL Technical Report ORNL/TM-12807, August 1994, also available online at ht p://www.cs.utk.edu/ ̃santa/homepage/12807.ps.gz 18. Press W. H., Teukolsky S. A., Wetterling W. T., Flannery B. P., Numerical recipes in C: the Art of Scientific Computing, Cambridge University Press, Cambridge, pp. 493-495.


Introduction
In recent years, wavelength scale photonic devices have gained considerable attention.A particular class, in which the geometry is periodic in one, two or three dimensions, is called photonic crystal [1].Devices with periodicity on a plane and light propagation in the perpendicular direction are sometimes referred to as photonic crystal fibers [2].Numerous computational approaches have been developed for solving field distributions and dispersion diagrams characterizing these structures.These include methods based on finite elements [3], finite-difference time-domain [4], transfer matrix method [5,6,7], plane wave expansion [8] and modal expansion [9].
The underlying idea in transfer matrix methods is the characterization of a given slab in terms of a matrix, which interrelates the fields on the two bounding surfaces of the slab.To this end the interface conditions on one surface of the slab are projected onto the second surface utilizing problem-specific expansion functions.Once the transfer matrix has been created, various quantities of practical interest can be computed; such as, the scattered-and transmitted field distributions, or the eigenfrequencies of the problem.
In the modal expansion method, one diagonalizes the Maxwell's equation with respect to a distinguished direction and computes the eigenfunctions of the resulting operator.(The notion of diagonalization will be addressed shortly.)The corresponding eigenvalues represent the propagation constants in the direction of the diagonalization.Once the eigenfunctions have been computed, they can be used to expand the fields.The constructed eigenfunctions are, however, valid basis functions only in those regions, which can be characterized by translationally invariant dielectric functions in the diagonalization direction.For each distinctly invariant layer, a separate set of eigenfunctions is needed.This becomes infeasible if the problem can not be divided into a reasonable number of invariant layers.On the other hand, in semi-infinite spaces the modal expansion method is extraordinarily useful: once the eigenfunctions and expansion coefficients are known, fields are uniquely determined everywhere in the invariant semi-space.This makes the implementation of open boundaries particularly efficient.The distinction between the modal-expansion based methods and the transfer-matrix based methods is not very sharp: indeed eigenfunctions can be, and often are, utilized in the construction of the transfer matrix (see e.g.[10]).
In this work we describe a method, which is particularly suited to problems, which are periodic in two space directions and non-periodic in the third one.The present work is formulated in cylindrical co-ordinates; the adaptation of the method to two dimensional [11,12,13] and three dimensional [14] problems in Cartesian co-ordinates, however, has been thoroughly discussed in previous works.
Our approach is related to the transfer matrix method, the modal expansion method as well as the planewave expansion method.We exploit the advantages of each of these methods and adapt them to the particular case of two dimensional periodicity, allowing arbitrary dielectric functions in the unit cell.The problems we are particularly interested in are periodic both in the z− and θ −directions, and non-periodic in the radial r−direction.(θ −periodicity is a natural consequence of the choise of the cylindrical co-ordinate system.)We expand the fields in terms of planewaves in z− and θ −directions for two reasons: i) the implementation of periodic and Bloch periodic boundary conditions is trivially straighforward, and ii) the spatial derivatives can be evaluated very efficiently.In the radial direction, we utilize a combination of two bases: i) inside the cylinderical simulation domain 0 < r < a, where we assume arbitrary but z−periodic dielectric function, we utilize finite differences, while ii) outside this cylinder, we assume constant dielectric and expand the fields in terms of eigenfunctions, which exist in closed-form.This strategy ensures maximal flexibility in the choise of the dielectric function, combined with the ease of implementation of open boundaries.
Our approach could be formulated in terms of the transfer matrix method with the only unknowns being the fields on two opposite boundaries of the cylinder: the outer boundary r = a and the axis r = 0. Instead, we include the field components on nested cylindrical surfaces r i = i∆, with 0 < r i < a, as explicit unknowns.In this manner we obtain sparse matrices, which can be solved, according to our extensive numerical experiments, with greater robustness and accuracy.

Outline of the method
Our simulation domain consists of a cylindrical disk characterized by the radius a and the periodicity length L z in the z−direction.Inside the disk we assume an arbitrarily inhomogeneous dielectric.To perform the analysis the interior domain is subdivided into discrete cylindrical surfaces r i = i∆ on which the fields are expanded in terms of planewaves.Here, ∆ is the discretization length.The expansion coefficients for the fields on the interlaced surfaces r i = i∆ are interrelated by the application of a simple, yet powerful, implementation of the finite difference scheme.The boundary conditions on the outermost cylindrical surface r N = a may be open or closed.In case of open boundary, we assume free space outside the disk, and expand the fields in terms of the eigenfunctions of the Maxwell's equations in cylindrical co-ordinates, which we obtain in closed-form.The phase-periodic condition in the z−direction is implicitly accounted for by choosing Bloch waves propagating in this direction.
The method is particularly suitable for analyzing cylindrical symmetric problems.However, since the dielectric function can vary arbitrarily in all three spatial co-ordinates, any z−periodic geometry can be handled straightforwardly.Needless to say that for cylindrical symmetric problems, and in particular for z−independent structures, the number of basis functions required reduces substantially.For problems which are homogeneous in the z−direction, only one har-monic function suffices to describe the propagating mode.For cylindrical symmetric problems, the eigenmodes are Bessel functions and the radial dependence can be described in terms of two appropriately chosen planewaves.
Other noteworthy features of our technique are: • Realistic boundary conditions can be implemented easily and efficiently.
• Iterative solvers are utilized which require modest computational resources.
• The wavenumber and the frequency are both input parameters allowing targeted and efficient computations.
• The method is applicable to eigenmode-as well as field excitation problems.
• In case of open boundary problems, the availability of the fields in closed-form in free space allows us to easily express the fields anywhere in space once the solution in the interior domain is known.

Diagonalization of the Maxwell's equations
Our aim is to diagonalize the Maxwell's equation with respect to the radial co-ordinate.To be more specific, we wish to write the Maxwell's equations in the form: This representation has several distinguished properties including: • Only the ''transversal'' field components, i.e., those involved in the interface conditions on cylindrical surfaces (r = constant) enter our calculations.
• The operator L involves spatial derivatives ∂ /∂φ and ∂ /∂ z, and material parameters which are defined on r = constant surfaces.
• The derivative of the transversal fields with respect to r appears at the RHS only.
• In view of (1) we realize that once we have determined the transversal field distribution on a certain cylindrical surface (r = constant), we progressively can compute the transversal fields on any other neighboring cylindrical surface.
• A further property of L is that repeated application of L to (1) results in higher-order r−derivatives of the transversal fields.This can be shown fairly easily: Assume (2) is given: with L = L (∂ /∂φ, ∂ /∂ z).(For notational simplicity material parameters are not shown explicitly in (2).)Applying L to both sides of (2), writing L (2) for LL at the LHS, and interchanging the order of L and ∂ /∂ r at the RHS, we obtain: Substituting ∂ ψ/∂ r for L ψ at the RHS of (3) results in A proof by induction shows that ( 5) is valid for any whole number n: Therefore, keeping r constant (r = r 0 ), we can formally build the derivatives of ψ with respect to r to any order desired.Thereby, we merely need to know the function ψ(r 0 , φ , z).
• Taylor Expansion: From the above we have Consequently, for the Taylor expansion of ψ(r 0 + h, φ , z) at r = r 0 we can write We could explore these ideas further, and discuss their implications for designing a host of robust and efficient algorithms.In this contribution, however, we limit ourselves to the basic form given in (2).It will be shown that this formula allows a powerful implementation of the central finite difference method.
The reasons for choosing the diagonalized form is threefold: i) The diagonalized form separates the non-periodic co-ordinate from the periodic ones and allows us to easily use different bases in different directions.ii) The diagonalized form only involves transversal field components, which are continuous across interfaces in radial direction.This makes the implementation of the finite difference algorithm particularly easy.iii) The eigenfunctions of (2) can be solved in closed-form for translationally invariant media.As will be shown further below, the fields expressed in the modal basis can be matched with the fields in the finite difference basis simply by equating the sets of expansion coefficients.
Finally, and merely for the sake of completeness, it should be pointd out that the diagonalization procedure, if carried out successfully, factorizes a given differential operator into commutative constituents.The rich properties of commutative operators, and the corresponding algebras and representations can be systematically utilized to create fast and robust algorithms.In contructing our basis functions used in this work we have tacitly made use of these properties.
To further develop the underlying theory we start with the Maxwell's curl equations for timeharmonic oscillations (exp(− jωt)): One possible way for expressing the curl operator in cylindrical co-ordinates is with u r , ru φ , and u z , respectively, denoting the unit vectors in the r−, φ −, and z−directions.It is advantageous to use the variable substitution defined in (11) for all φ −directional field components.This substitution not only simplifies the manipulations, but also enhances the accuracy of the computations.
From the symmetry properties of Maxwell's equations we know that ( 9) can be obtained from ( 8) by the replacements E ↔ H, and ε ↔−µ.Therefore, it is sufficient to perform our manipulations on (8) only; results associated with ( 9) can be obtained by the aforementioned replacements.
Using (10) for the evaluation of the curl in ( 8) we obtain Remembering that we are seeking a diagonalization with respect to r, a moment's reflection on Eqs. ( 12) reveals the following facts: (i) The first equation in (12) does not involve any r−derivatives.Therefore, it will not appear in our diagonalized form explicitly. (ii) The second and third equations, respectively, involve the r−derivative of e φ and e z .This fact implies that these latter field components are the ''essential'' field components in our equations, and, therefore, e r has to be eliminated.(iii) Due to the symmetry of our problem we conclude that our diagonalized form will ultimately also involve the variables h φ and h z , and that the component h r has to be eliminated.The first equation in (12) serves to expressing the undesirable field component h r in terms of the ''essential'' field components e z and e φ .Likewise using the second curl equation e r can be expressed in terms of the ''essential'' field components h z and h φ .More explicitly, we obtain: Substituting this equation into the second and third equations in (12) and rearranging the terms, we obtain: Constructing the counterparts of ( 14) and ( 15) and combining the results, we obtain the desired diagonalized form: The explicit forms for the differential operators A ij and B ij in ( 16) are summarized below: Once the ''transversal'' field components e φ , e z , h φ , and h z , have been determined, the remaining ''normal'' field components e r and h r can be solved from Recall that these equations are the first equations in the curl equations, which do not contain any r−derivatives.(In (25) the superscript T denotes transposition.)

Field expansions
On (z, φ )-cylindrical surfaces we expand the fields in terms of planewaves with general (rdependent) expansion coefficients f m,n (r) Here, k m = 2πm/L z + K z with K z denoting the Bloch phasing factor in the z−direction.Details concerning the discretization in r−variable will be explained below.

Discretization in the radial direction
In ( 16) the r−derivative of the e−field is given by the h−field components alone, and vice versa.This property suggests the discretization of the electric-and magnetic fields on interlaced cylindrical surfaces.Stated more precisely, we discretize h z and h φ on surfaces r i = i∆, and discretize e z and e φ on surfaces r i+1/2 =(i + 1/2) ∆, where ∆ stands for the discretization step length in the radial direction.In the following we will refer to a specific r−position by providing an appropriate multiplier of ∆.
We establish a finite-difference relationship between the coefficients on adjacent cylindrical surfaces by using Taylor series expansions of the fields, and keeping first-order (polynomial) terms only: Here, the operators A i pq and B i+1/2 pq are discrete analogues corresponding to the differential operators defined in ( 17), ( 18), ( 19), ( 20), ( 21), ( 22), ( 23) and ( 24).The vectors h z , h φ , e z , and e φ comprise expansion coefficients of dimension M × N of the corresponding fields.
It should be pointed out that in terms of the discrete matrix operators in ( 27) and (28) we can compute the r−derivative of the Fourier coefficients on a given cylindrical surface.In this manner we can formulate the r−dependence of the Fourier coefficients in terms of finite differences, exactly the same way as we would operate with real domain fields.

Boundary conditions
For a complete specification of our boundary value problem, we need to formulate the boundary conditions.We need to define two boundary conditions for fields in the proximity of the axis of the simulation disk (r = 0), and two conditions on outer bounding surface (r = a).

Axis of the cylinder: r = 0
In view of ( 27) and (28) we can define the magnetic fields on cylindrical surfaces r = i∆, including r = 0. Based on the variable substitution in (11) we require that This requirement constitutes our first boundary condition.
The formulation of the second boundary condition on a cylindrical surface, which ''tightly'' embraces the r = 0 axis, is slightly more complicated.Naively, it would be tempting to set e 0 φ ≡ 0; however, the transversal electric fields in the discrete system are defined on the layers r =(i + 1/2) ∆, rather than on the center axis of the simulation disk.However, setting e φ ≡ 0 would mean assuming a fictitious impedance, and therefore, would alter the physics of our problem.
To remedy these difficulties, we propose the following: We use Faraday's law in integral form to generate the missing boundary condition.Next we approximately evaluate the integrals in (30) by substituting field expansions in (26) with a priori unknown coefficients.Utilizing the orthogonality of the harmonic basis functions we obtain an adequate number of equations for the determination of the unknowns.On the LHS of (30), we choose the integration path to be the periphery of the circle: [l : r = ∆/2].On this path, the involved scalar product becomes On the RHS of (30), we choose to integrate on a plane perpendicular to the z−axis.The integrand becomes In the finite difference approximation, we take h z to be constant over the surface: 0 < r < ∆/2, 0 < φ < 2π, z = z 0 .Upon this assumption, the integral at the RHS assumes the form: Substituting the field expansion given in (26) into (30), and using (31) and (33) result in Here δ [n] is the Kronecker delta symbol.The orthogonality of the φ −dependent harmonics can be used to create a sufficient number of equations.We multiply both sides of (34) by exp (− j nφ ) all N harmonics in the basis and integrate over the periodicity interval with the length 2π.These steps result in the following set of equations for various n Finally, the equations for individual coefficients can be obtained by multiplying both sides of the corresponding equations by the exponential functions exp (− jk m z) for M distinct planewaves, and, consecutively, integrating over the periodicity interval with the length L z .This procedure results in the desired set of equations: φ , m,0 = 0. (37)

Free space boundary condition
For constructing the remaining two boundary conditions on the surface r = a we utilize the free space eigenvectors of the diagonalized form given in (16).The eigenvectors of ( 16) in Fourier domain, which satisfy the Sommerfeld's radiation condition in free space, are (e.g.Ref. [15]) Here, H n (•) denotes Hankel function of order n and the first kind, and λ 2 m = ω 2 εµ − k 2 m is the square of the propagation constant in free space.In order to ensure decaying behavior (integrability), we choose the branch with ℜ(λ m )=0 and ℑ(λ m ) > 0.
The transversal fields outside the cylindrical simulation domain (r > a) can be expanded in terms of the eigenvectors in ( 38) and (39): Here, the 4 × 1 vector Ψ(r, φ , z) comprises the transversal field components in real space, with a m,n and b m,n being expansion coefficients.Note that the basis functions are harmonic in the variables φ and z, and thus orthogonal in (φ , z)−domain.Utilizing this property, we can establish relationships between a m,n and b m,n and the expansion coefficients of the fields in the interior domain.There is, however, an important detail which we should be aware of while implementing these ideas: Inside the simulation disk, the magnetic-and electric fields are defined on interlaced rather than immediate neighboring surfaces.
In order to account for this circumstance, and thus to adequately ''correct'' the electric field coefficients in (40), we shift r by the distance ∆/2 in the radial direction.Note that in our formulation this can be accomplished fairly easily, since the functional form of the eigensolutions in the r−direction is known.

Closed boundary
The implementation of the closed boundary is trivial.In case of electrically-or magnetically conducting boundary, respectively, one requires the electric-or magnetic fields, to vanish on the outermost cylindrical surface.

Constructing the system matrix
Before we embark on solving the electromagnetic wave propagation in the assumed z−corrugated fiber, it is instructive to collect all the occurring expansion coefficients into one vector, say, f.Next we set up the system matrix M, by evaluating ( 27) and (28) at individual (discrete) cylindrical surfaces in the simulation domain, and matching the solutions.We obtain the following homogeneous system The parameters ω and K z have been written explicitly in order to emphasize their role as input parameters.
It is important to point out that M is very well structured.This fact is instrumental in numerical calculations: The equations resulting from ( 27) and (28) create two unity side diagonals, having opposite signs, and being offset from the main diagonal by ±2M × N elements, with M and N denoting the number of planewaves in the z− and φ −directions, respectively.Between these two diagonals the matrix operators A m,n and B m,n generate three block diagonals of the size M × N. The equations originating from free space eigenvectors can be written in a form which preserves this structure.Unfortunately, however, the equations resulting from the boundary conditions in the proximity of the center axis, do not have this desirable structure.Nonetheless, the matrix elements arising from these equations are sparse, and structured systematically in a specific way.

Solving the system of equations
Solutions of the homogeneous system (41) are eigenmodes of our corrugated fiber.We thoroughly have discussed the efficient solution of similar matrices elsewhere.Therefore, here, we restrict ourselves to a few useful hints and guidelines.For more details we refer the interested reader to Refs.[12,13,14].

Excitation problems
The equation system (41) can be easily modified to include excitation problems, involving electric-or magnetic current distributions in transversal directions.Solving excitation problems does not require any modification of the system matrix: The zero vector at the RHS can be replaced by a vector describing the current distribution.The modified system equation is given in (42).(For details considering problems in the Cartesian co-ordinate system refer to [11,12,13,14].) The numerical solution to this system is best obtained by means of iterative solvers.Most iterative solvers do not perform any direct factorization; they basically require a software routine to furnish matrix-vector products.This fact is of significance since the system matrices for three dimensional problems can become prohibitively large.
Instead of discretizing the operators introduced in ( 16) directly, we propose an alternative procedure which enables us to apply these operators to a trial vector f numerically.While, real space derivatives can be effectively evaluated in Fourier domain, multiplications by functions, e.g.ε can be preferably carried out in real domain.Therefore, we will Fourier and inverse Fourier transform the coefficient vector, in order to move back and forth between the two spaces whenever necessary.Direct discretization would create MN by MN complex-valued matrices and multiplications involving them would require O((MN) 3 ) operations.In the proposed way, the matrices involved are all diagonal, and the resulting dominating factor O(MN ln(MN)) stems from the Fast Fourier Transform.Proceeding along these lines, we can construct the matrix product piece by piece.The actual matrix M is never constructed explicitly.
Iterative solvers, such as Transpose Free Quasi Minimal Residual method [16], do not usually converge fast enough without utilizing a suitable pre-conditioner.Our approach is to use a system matrix M, which corresponds to a simplified problem defined by ∂ ε/∂φ = ∂ ε/∂ z ≡ 0. On each cylindrical surface, the value of ε can be obtained by averaging the ε of the original problem over the corresponding cylindrical surface.The resulting simplified problem, leads to a sparse and well-structured matrix M, the direct factorization of which is computationally affordable.
Once the routines which perform matrix products with M, and a routine which computes M −1 g for a given vector g are available, the preconditioned system can be solved as described in [17].
After the solution is obtained, the real domain fields can be recovered by Fourier transforming appropriate coefficient sets.The expansion coefficients characterizing free space exterior domain together with (38), (39) and (40) can be used to compute the fields in r > a.

Eigenproblems
The eigenmodes of a system correspond to those ordered values (ω, K z ), which make the system matrix M singular.There are a variety of methods in our disposal for determining the singular points of a matrix.The common feature of most of these ''traditional'' methods is that they rely on the direct factorization of the system matrix.Here, we propose an alternative way, which is based on the field solution to a randomly, yet, suitably chosen excitation.To illuminate the details consider a (ω, K z )−point for which the matrix M in (42) approaches a singularity.(For the sake of simplicity we limit ourselves to non-degenerate cases.)Then, one of the eigenvalues of M approaches zero.Assume that y is the corresponding eigenvector.Under these conditions, f, the solution of (42) approaches y, unless the algebraic scalar product y T • ρ is vanishingly small.At the same time, the L 2 norm ||f|| 2 approaches infinity.This theorem can be easily verified if both the matrix product Mf and the excitation ρ are thought to be expanded in terms of the eigenvectors of M, as is done in [18].
Guided by this idea, we can search for the eigenmodes by varying ω for a given K z (or vice versa), solve fields in response to a random excitation, and use the solution norm as a measure for the singularity.

Bragg fiber
Consider a z−periodic Bragg fiber with the following properties: one unit cell is constructed by stacking two cylindrical elements on top of each other.In addition, a sector of 90 degrees has been removed from both elements as sketched in Fig. 1.The refractive indices of the elements are n 1 = 1.5 and n 2 = 2.The radius of the fiber is r = 0.5, the height of both pieces constituting the unit cell is d = 0.5 and the height of the unit cell is L z = 1.Outside the fiber is free space.
Our results are compared with the planewave method solutions [8].Computations are carried out by using 16 grid points for each of the three spatial variables, in our method, and, using a 128 × 128 × 16 grid for a supercell of size 8 × 8 × 1 in the planewave method.The computed Brillouin dispersion diagrams are shown in Fig. 2. The match is encouraging.The minor discrepancies are most likely due to the coarse grid.

Bragg fiber with cylindrical symmetry
Our second test case is a Bragg fiber which consists of two periodically stacked disks with relative dielectric constants ε r1 = 2.25 and ε r2 = 4.The radius of the fiber is r = 0. Curves with the marker ''o'' are computed with our method while curves with the marker ''x'' are obtained using the planewave method.Inset shows a zoom of the bands near the Brillouin zone edge.heights of both disks are h 1 = h 2 = 0.5.The medium outside the fiber is free space.Due to the cylindrical symmetry of the problem, we only use 4 planewaves in the φ −direction.This allows us to expand the modes in terms of the Bessel functions up to the first order.We solve one mode only with ω = 0.4 × 2πc/L z and compare the results with those obtained from the eigenmode expansion technique [9].
In the reference method the structure of interest is enclosed inside a metallic cylinder.In the immediate proximity of the metallic surface a Perfectly Matched Layer (PML) has been assumed.Ideally, a PML should emulate an open boundary, but in practice there will be reflections and we have to assume free space between the fiber and the PML.In our numerical experiment we surrounded the fiber with 4.9 units of free space and 0.1 units of PML material.The convergence of both methods as a function of discretization is shown in Fig. 3.
Both methods seem to converge to the same value, the most accurate solution with our method being k z = 0.4673591 × 2π/L z , while for the eigenmode expansion method being k z = 0.4673595 × 2π/L z , with the relative difference compared with the latter being 8.068 × 10 −7 .It should be pointed out that due to the intrinsic differences between the two methods, the discretization sizes are not directly comparable.

Fields in a fiber with air holes
Our final test case consists of a fiber with six air holes arranged symmetrically around the center.The radius of the fiber is r f = 1, the assumed relative dielectric constant is ε r = 4, and outside the fiber is free space.The centers of the air holes are located on a circle with r c = 0.5L (for some length scale L), and the radii of the air holes are r a = 0.1L.There is no variation in the z−direction.We compute the lowest-order eigenmode for the frequency ω = 0.3 × 2πc/L and find its propagation constant along the z−axis to be k z = 0.33035 × 2π/L.The eigenmode is found by the alternative method outlined above: solve the field pattern associated with a random excitation and seek the maximum of the solution norm.In this way the eigenmode field pattern is obtained immediately by inverse Fourier transform of the resulting set of expansion coefficients.Since the simulation is two dimensional we have used a single planewave in the  3.The converge of the eigenmode propagation constant as a function of discrete basis size.The curve with circular tags is computed with our method, while the curve with cross tags is obtained with the eigenmode expansion method.For our method, the abscissa denotes the discretization step in the radial direction and also the number of planewave components in the z−direction, which are both equal here.For the eigenmode expansion method the horizontal axis denotes the number of eigenmodes used.z−direction only.In the r−direction, we discretized the fields on 64 concentric cylindrical surfaces, while in the φ −direction we have used 64 planewaves.The resulting z−directional fields are shown in Fig. 4.

Conclusion
We presented a robust and performance-enhanced method for numerically analyzing the propagation behavior of electromagnetic waves in structures which are periodic in two directions in space and arbitrarily non-periodic in the remaining direction.In particular, we considered z−periodic structures in (r, φ , z)−cylindrical co-ordinates.A cylindrical disk, characterized by the radius a and the periodicity length L z , defined the fundamental cell in our problem.The permittivity of the dielectric inside this cell was characterized by an arbitrary, single-valued function ε(r, φ , z) of all three spatial co-ordinates.We addressed both open and closed boundary problems.The presented method relied on the diagonalization of Maxwell's equations in the radial (r−) direction.The presented method requires the discretization of the fields in the interior of the disk only.Inside the disk volume, we expanded the fields in terms of planewaves on discrete cylindrical surfaces r i = i∆, with ∆ being the discretization step length the nested surfaces r i = i∆ in the interior of the simulation domain were interrelated by the application of a simple, yet, powerful finite difference scheme.In the free space outside the disk, the fields were expanded in terms of closed-form eigensolutions of the Maxwell's equations in cylindrical co-ordinates.In terms of three examples we presented the general applicability of our method.
Our present research efforts are focused on the development of Wannier function like and other localized analyzing functions for the efficient modeling of defects in photonic crystals and devices.

φφ
,m,n e jk m z e jnφ dφ = jωµ π∆ ,n e jk m z e jnφ .(34)Interchanging the order of summation and integration at the LHS, and integrating each term individually we obtain: ,m,n 2πδ[n]e jk m z .

Fig. 1 .
Fig. 1.The geometry of test case 1. Image shows one unit cell, which is periodically replicated in the z−direction.

Fig. 2 .
Fig. 2. Dispersion diagram for the four guided modes in the geometry shown in Fig. 1.Curves with the marker ''o'' are computed with our method while curves with the marker ''x'' are obtained using the planewave method.Inset shows a zoom of the bands near the Brillouin zone edge.
Fig.3.The converge of the eigenmode propagation constant as a function of discrete basis size.The curve with circular tags is computed with our method, while the curve with cross tags is obtained with the eigenmode expansion method.For our method, the abscissa denotes the discretization step in the radial direction and also the number of planewave components in the z−direction, which are both equal here.For the eigenmode expansion method the horizontal axis denotes the number of eigenmodes used.

Fig. 4 .
Fig. 4. The real part of the z−directional eigenmode field pattern in a fiber with circular air holes.