OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 11, Iss. 23 — Nov. 17, 2003
  • pp: 3048–3062
« Show journal navigation

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

K. Varis and A. R. Baghai-Wadji  »View Author Affiliations


Optics Express, Vol. 11, Issue 23, pp. 3048-3062 (2003)
http://dx.doi.org/10.1364/OE.11.003048


View Full Text Article

Acrobat PDF (493 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

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 radius a and the periodicity length Lz , 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 r=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 ri =iΔ, with Δ being the discretization step length. The fields on the nested surfaces ri =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 interior- and 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

1. Introduction

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.

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 ri =iΔ, with 0<ri <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 Lz 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 ri =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 ri =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 rN =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:

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

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:

(ϕ,z,μ,ε,ω)[hϕhzeϕez]=r[hϕhzeϕez]
(1)

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 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 is that repeated application of to (1) results in higher-order r-derivatives of the transversal fields. This can be shown fairly easily: Assume (2) is given:
ψ=rψ
(2)

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

(2)ψ=rψ
(3)

Substituting ∂Ψ⃗/∂r for ℒΨ⃗ at the RHS of (3) results in

(2)ψ=2r2ψ.
(4)

A proof by induction shows that (5) is valid for any whole number n:

(n)ψ=nrnψ
(5)

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
{nψ(r,ϕ,z)rn}|r=r0=(n)(ϕ,z)ψ(r0,ϕ,z).
(6)

Consequently, for the Taylor expansion of Ψ⃗(r 0+h,ϕ,z) at r=r 0 we can write

ψ(r0+h,ϕ,z)=n=0hnn!{(n)(ϕ,z)ψ(r0,ϕ,z)}.
(7)

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.

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

×E=jωμH
(8)
×H=jωεH
(9)

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

×f=1rurruϕuzrϕzfrrfϕfz,
(10)

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.

rfϕ=f˜ϕ
(11)

From the symmetry properties of Maxwell’s equations we know that (9) can be obtained from (8) by the replacements EH, 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

1r[ϕezze˜ϕrrez+rzerre˜ϕϕer]=jωμ[hr1rh˜ϕhz].
(12)

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 ez . This fact implies that these latter field components are the “essential” field components in our equations, and, therefore, er has to be eliminated. (iii) Due to the symmetry of our problem we conclude that our diagonalized form will ultimately also involve the variables ϕ and hz , and that the component hr has to be eliminated. The first equation in (12) serves to expressing the undesirable field component hr in terms of the “essential” field components ez and ϕ . Likewise using the second curl equation er can be expressed in terms of the “essential” field components hz and ϕ . More explicitly, we obtain:

er=1jωεrϕhz+1jωεrzh˜ϕ
(13)

Substituting this equation into the second and third equations in (12) and rearranging the terms, we obtain:

z1jωεrϕhz+z1jωεrzh˜ϕjωμrh˜ϕ=rez
(14)
jωμrhzϕ1jωεrϕhz+ϕ1jωεrzh˜ϕ=re˜ϕ
(15)

Constructing the counterparts of (14) and (15) and combining the results, we obtain the desired diagonalized form:

[00A11A1200A21A22B11B1200B21B2200][eze˜ϕhzh˜ϕ]=r[eze˜ϕhzh˜ϕ]
(16)

The explicit forms for the differential operators Aij and Bij in (16) are summarized below:

A11=z1jωεrϕ
(17)
A12=z1jωεrzjωμr
(18)
A21=ϕ1jωεrϕ+jωμr
(19)
A22=ϕ1jωεrz
(20)
B11=ϕ1jωεrz
(21)
B12=1jωμr2zϕ+jωεr
(22)
B21=1jωμr2ϕ2jωεr
(23)
B22=1jωμr2ϕz
(24)

Once the “transversal” field components ϕ , ez , ϕ , and hz , have been determined, the remaining “normal” field components er and hr can be solved from

[erhr]=[01jωμrϕ01jωμrz1jωεrϕ01jωεrz0]T[eze˜ϕhzh˜ϕ].
(25)

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 fm,n (r)

f(r,ϕ,z)=m,nfm,n(r)ejkmzejnϕ.
(26)

Here, km =2πm/Lz +Kz with Kz 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 hz and hϕ on surfaces r i=iΔ, and discretize ez 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:

Δ[A11iA12iA21iA22i][h˜ϕihzi]+[e˜ϕi12ezi12]=[e˜ϕi+12ezi+12]
(27)
Δ[B11i+12B12i+12B21i+12B22i+12][e˜ϕi+12ezi+12]+[h˜ϕihzi]=[h˜ϕi+1hzi+1]
(28)

Here, the operators Apqi and Bpqi+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, ϕ, 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

h˜ϕ00.
(29)

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~ϕ00; 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˜ϕ120 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

lE·dl=jωμSH·dS
(30)

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

E·dl=e˜ϕdϕ.
(31)

On the RHS of (30), we choose to integrate on a plane perpendicular to the z-axis. The integrand becomes

H·dS=hzdS.
(32)

In the finite difference approximation, we take hz 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:

jωμSH·dS=jωμπΔ24hz0.
(33)

Substituting the field expansion given in (26) into (30), and using (31) and (33) result in

02πm,ne˜ϕ,m,n12ejkmzejnϕ=jωμπΔ24m,nhz,m,n12ejkmzejnϕ.
(34)

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

m,ne˜ϕ,m,n122πδ[n]ejkmz.
(35)

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

m,ne˜ϕ,m,n12δ[n]δ[n̂]4π2ejkmz=jωμπΔ24m,nh0,m,n122πδ[nn̂]ejkmz.
(36)

Finally, the equations for individual coefficients can be obtained by multiplying both sides of the corresponding equations by the exponential functions exp(-jkz) for M distinct planewaves, and, consecutively, integrating over the periodicity interval with the length Lz . This procedure results in the desired set of equations:

hz,m̂,n̂0δ[n̂]8jωμΔ2eϕ,m̂,012=0.
(37)

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

15. S. Liu, L. W. Li, M. S. Leong, and T. S. Yeo, “Theory of Gyroelectric Waveguides,” PIER , 29, 231–259 (2000). [CrossRef]

])

Φm,n1=[e˜ϕezh˜ϕhz]=[jnHn(λmr)λ2jkmHn(λmr)rωεkmrHn(λmr)0],
(38)
Φm,n2=[e˜ϕezh˜ϕhz]=[jrωμkmrHn(λmr)0nHn(λmr)λ2kmHn(λmr)].
(39)

Here, Hn (·) denotes Hankel function of order n and the first kind, and λm2=ω 2 εµ-km2 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):

Ψ(r,ϕ,z)=m,n[am,nΦm,n1(r)+bm,nΦm,n2(r)]ejkmzejnϕ.
(40)

Here, the 4×1 vector Ψ(r,ϕ,z) comprises the transversal field components in real space, with am,n and bm,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 am,n and bm,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

M(ω,Kz)f=0.
(41)

The parameters ω and Kz 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

12. K. Varis and A.R. Baghai-Wadhi, “Z-diagonalized Planewave/FD Approach for Analyzing TE Polarized Waves in 2D Photonic Crystals,” ACES, March 24–28, 2003, Monterey, CA, USA.

, 13

13. K. Varis and A. R. Baghai-Wadji, “A Novel 2D Pseudo-Spectral Approach of Photonic Crystal Slabs,” ACES Special Issue, (submitted).

, 14

14. K. Varis and A.R. Baghai-Wadji, “A Novel 3D Pseudo-spectral Analysis of Photonic Crystal Slabs,” ACES Special Issue, (submitted).

].

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

11. K. Varis and A. R. Baghai-Wadji, “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

12. K. Varis and A.R. Baghai-Wadhi, “Z-diagonalized Planewave/FD Approach for Analyzing TE Polarized Waves in 2D Photonic Crystals,” ACES, March 24–28, 2003, Monterey, CA, USA.

, 13

13. K. Varis and A. R. Baghai-Wadji, “A Novel 2D Pseudo-Spectral Approach of Photonic Crystal Slabs,” ACES Special Issue, (submitted).

, 14

14. K. Varis and A.R. Baghai-Wadji, “A Novel 3D Pseudo-spectral Analysis of Photonic Crystal Slabs,” ACES Special Issue, (submitted).

].)

M(ω,Kz)f=ρ(ω,Kz)
(42)

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

16. R. W. Freund, “A Transpose-Free Quasi-Minimal Residual Algorithm for Non-Hermitian Linear Systems,” SIAM Journal of Scientific Computing , 14, 470–482 (1993). [CrossRef]

], do not usually converge fast enough without utilizing a suitable pre-conditioner. Our approach is to use a system matrix , 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 , the direct factorization of which is computationally affordable.

Once the routines which perform matrix products with M, and a routine which computes -1 g for a given vector g are available, the preconditioned system can be solved as described in [17

17. N. Nachtigal, R.W. Freund, and J. C. Reeb, “QMRPACK user’s guide,” ORNL Technical Report ORNL/TM-12807, August 1994, also available online at http://www.cs.utk.edu/~santa/homepage/12807.ps.gz

].

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.

8.2. Eigenproblems

The eigenmodes of a system correspond to those ordered values (ω,Kz ), 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 (ω,Kz )-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 ∥f2 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

18. W. H. Press, S. A. Teukolsky, W. T. Wetterling, and B. P. Flannery, Numerical recipes in C: the Art of Scientific Computing, Cambridge University Press, Cambridge, pp. 493–495.

].

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=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 Lz =1. Outside the fiber is free space.

Our results are compared with the planewave method solutions [8

8. S. G. Johnson and J.D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express 8, 173–190 (2001) http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173. [CrossRef] [PubMed]

]. 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 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/Lz and compare the results with those obtained from the eigenmode expansion technique [9

9. P. Bienstman, 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

].

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.

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 kz =0.4673591×2π/Lz , while for the eigenmode expansion method being kz =0.4673595×2π/Lz , 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 rf =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 rc =0.5L (for some length scale L), and the radii of the air holes are ra =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 kz =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.
Fig. 4. The real part of the z-directional eigenmode field pattern in a fiber with circular air holes.

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 Lz , 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 ri =iΔ, with Δ being the discretization step length. The fields on the nested surfaces ri =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.

Acknowledgments

One of the authors (K. V.) thanks the Finnish Graduate School of Modern Optics and Photonics, and the Ausseninstitut (Extension Institute) at Vienna University of Technology for providing research scholarships which made possible a visit to Vienna, Austria for finalizing the current project.

The work was partially done while the second author was visiting the Institute for Mathematical Sciences, National Univesrity of Singapore and Institute of High Performance Computing (IHPC) in 2003. The visit was supported by the Institute and IHPC.

References and links

1.

J. D. Joannopoulos, R. D. Meade, and J. N. Winn, “Photonic Crystals: Molding the Flow of Light,” Princeton, September 1995.

2.

J. C. Knight, T. A. Birks, P. St. J. Russell, and D. M. Atkin, “All-silica single-mode optical fiber with photonic crystal cladding,” Opt. Lett. 21, 1547–1549 (1996). [CrossRef] [PubMed]

3.

A. Cucinotta, A. Selleri, L. Vincetti, and M. Zoboli, “Holey Fiber Analysis Through the Finite-Element Method,” IEEE Photon. Technol. Lett. 3, 147–149 (1991).

4.

C.T Chan, Q. L. Yu, and K. M. Ho, “Order-N spectral method for electromagnetic waves,” Phys. Rev. B 51, 16635–16642 (1995). [CrossRef]

5.

Z.-Y. Li and L.-L. Lin, “Photonic band structures solved by a plane-wave-based transfer-matrix method,” Phys. Rev. E 67, 046607 (2003). [CrossRef]

6.

Elson J. Merle and P. Tran, “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). [CrossRef]

7.

L. Li, “Multilayer modal method for diffraction gratings of arbitrary profile, depth, and permittivity,” J. Opt. Soc. Am. A , 10, 2581–2591 (1993). [CrossRef]

8.

S. G. Johnson and J.D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express 8, 173–190 (2001) http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173. [CrossRef] [PubMed]

9.

P. Bienstman, 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.

Z.-Y. Li and K.-M. Ho, “Analytic modal solution to light propagation through layer-by-layer metallic photonic crystals,” Phys. Rev. B. 67, 165104 (2003). [CrossRef]

11.

K. Varis and A. R. Baghai-Wadji, “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.

K. Varis and A.R. Baghai-Wadhi, “Z-diagonalized Planewave/FD Approach for Analyzing TE Polarized Waves in 2D Photonic Crystals,” ACES, March 24–28, 2003, Monterey, CA, USA.

13.

K. Varis and A. R. Baghai-Wadji, “A Novel 2D Pseudo-Spectral Approach of Photonic Crystal Slabs,” ACES Special Issue, (submitted).

14.

K. Varis and A.R. Baghai-Wadji, “A Novel 3D Pseudo-spectral Analysis of Photonic Crystal Slabs,” ACES Special Issue, (submitted).

15.

S. Liu, L. W. Li, M. S. Leong, and T. S. Yeo, “Theory of Gyroelectric Waveguides,” PIER , 29, 231–259 (2000). [CrossRef]

16.

R. W. Freund, “A Transpose-Free Quasi-Minimal Residual Algorithm for Non-Hermitian Linear Systems,” SIAM Journal of Scientific Computing , 14, 470–482 (1993). [CrossRef]

17.

N. Nachtigal, R.W. Freund, and J. C. Reeb, “QMRPACK user’s guide,” ORNL Technical Report ORNL/TM-12807, August 1994, also available online at http://www.cs.utk.edu/~santa/homepage/12807.ps.gz

18.

W. H. Press, S. A. Teukolsky, W. T. Wetterling, and B. P. Flannery, Numerical recipes in C: the Art of Scientific Computing, Cambridge University Press, Cambridge, pp. 493–495.

OCIS Codes
(000.3860) General : Mathematical methods in physics
(000.4430) General : Numerical approximation and analysis

ToC Category:
Research Papers

History
Original Manuscript: September 4, 2003
Revised Manuscript: October 27, 2003
Published: November 17, 2003

Citation
Karri Varis and A. Baghai-Wadji, "Pseudo-spectral Analysis of Radially-Diagonalized Maxwell's Equations in Cylindrical Co-ordinates," Opt. Express 11, 3048-3062 (2003)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-11-23-3048


Sort:  Journal  |  Reset  

References

  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). [CrossRef] [PubMed]
  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, 16635- 16642 (1995). [CrossRef]
  5. Li Z.-Y., Lin L.-L., �??�??Photonic band structures solved by a plane-wave-based transfer-matrix method,�??�?? Phys. Rev. E 67, 046607 (2003). [CrossRef]
  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). [CrossRef]
  7. Li L., �??�??Multilayer modal method for diffraction gratings of arbitrary profile, depth, and permittivity,�??�?? J. Opt. Soc. Am. A, 10, 2581-2591 (1993). [CrossRef]
  8. Johnson S. G., Joannopoulos J.D., �??�??Block-iterative frequency-domain methods for Maxwell�??s equations in a planewave basis,�??�?? Opt. Express 8, 173-190 (2001)<a href=" http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173</a> [CrossRef] [PubMed]
  9. </a>Bienstman P., Rigorous and Efficient Modelling of Wavelength Scale Photonic Components, doctoral dissertation 2001, Ghent University, Belgium. Available online at <a href="http://photonics.intec.rug.ac.be/download/phd 104.pdf">http://photonics.intec.rug.ac.be/download/phd 104.pdf</a>
  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). [CrossRef]
  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). [CrossRef]
  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). [CrossRef]
  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 <a href="http://www.cs.utk.edu/�?santa/homepage/12807.ps.gz">http://www.cs.utk.edu/�?santa/homepage/12807.ps.gz</a>
  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

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1. Fig. 2. Fig. 3.
 
Fig. 4.
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited