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

2. 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 harmonic 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:

3. 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:

withℒ=ℒ(∂/∂ϕ,∂/∂z). (For notational simplicity material parameters are not shown explicitly in (2).) Applying ℒ to both sides of (2), writing ℒ ⁽²⁾ for ℒℒ at the LHS, and interchanging the order of ℒ and ∂/∂r at the RHS, we obtain:

Substituting ∂Ψ⃗/∂r for ℒΨ⃗ 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 ₀), we can formally build the derivatives of Ψ⃗ with respect to r to any order desired. Thereby, we merely need to know the function Ψ⃗(r ₀,ϕ,z).

Consequently, for the Taylor expansion of Ψ⃗(r ₀+h,ϕ,z) at r=r ₀ 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 time-harmonic oscillations (exp(-jωt)):

One possible way for expressing the curl operator in cylindrical co-ordinates is

with u _r, r u _ϕ, 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 ẽ_ϕ 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 ẽ_ϕ . 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 $\mathcal{A}$_ij and $\mathcal{B}$_ij in (16) are summarized below:

Once the “transversal” field components ẽ_ϕ , 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.)

4. Field expansions

On (z,ϕ)-cylindrical surfaces we expand the fields in terms of planewaves with general (r-dependent) 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.

5. 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 ${\mathbf{A}}_{\mathit{\text{pq}}}^i$ and ${\mathbf{B}}_{\mathit{\text{pq}}}^{i+1/2}$ 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 ẽ _ϕ 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.

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

6.1. 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 ${\stackrel{\mathbf{~}}{\mathbf{e}}}_{\mathrm{\varphi }}^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 ${\tilde{\mathbf{e}}}_{\varphi }^{\frac{1}{2}}\equiv 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 ₀. 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

Interchanging the order of summation and integration at the LHS, and integrating each term individually we obtain:

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(-jn̂ϕ) for 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:

6.2. 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 ${\lambda }_m^2$=ω ² εµ-$k_m^2$ 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.

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

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

8. 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].

8.1. 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)³) operations. In the proposed way, the matrices involved are all diagonal, and the resulting dominating factor O(MNln(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 ∂$\widehat{\epsilon }$ /∂ϕ=∂$\widehat{\epsilon }$ /∂z≡0. On each cylindrical surface, the value of $\widehat{\epsilon }$ 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.

 

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

Download Full Size | PDF

8.2. 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 ₂ norm ∥f∥₂ 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.

9. Numerical results

9.1. 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.5 and n ₂=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.

9.2. 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.5 and the heights of both disks are h ₁=h ₂=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].

 

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.

Download Full Size | PDF

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.

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

 

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.

Download Full Size | PDF

 

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

Download Full Size | PDF

10. 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 fields on 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.