Edee et al.
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
667
Mode solver based on Gegenbauer polynomial expansion for cylindrical structures with arbitrary cross sections Kofi Edee,1,2,* Mira Abboud,1,2 Gérard Granet,1,2 Jean Francois Cornet,1,2 and Nikolay A. Gippius1,2 1
Clermont Université, Université Blaise Pascal, Institut Pascal, BP 10448, F-63000 Clermont-Ferrand, France 2 CNRS, UMR 6602, Lasmea, F-63177 Aubière, France *Corresponding author: kofi.edee@univ‑bpclermont.fr Received November 15, 2013; revised January 17, 2014; accepted January 24, 2014; posted January 27, 2014 (Doc. ID 201081); published March 3, 2014
We present a modal method for the computation of eigenmodes of cylindrical structures with arbitrary cross sections. These modes are found as eigenvectors of a matrix eigenvalue equation that is obtained by introducing a new coordinate system that takes into account the profile of the cross section. We show that the use of Hertz potentials is suitable for the derivation of this eigenvalue equation and that the modal method based on Gegenbauer expansion (MMGE) is an efficient tool for the numerical solution of this equation. Results are successfully compared for both perfectly conducting and dielectric structures. A complex coordinate version of the MMGE is introduced to solve the dielectric case. © 2014 Optical Society of America OCIS codes: (000.3860) Mathematical methods in physics; (230.7370) Waveguides. http://dx.doi.org/10.1364/JOSAA.31.000667
1. INTRODUCTION The modal method by Gegenbauer expansion (MMGE) was recently introduced for plane-wave diffraction by a lamellar grating [1,2], and the key point of this method is that it allows transverse boundary conditions to be rigorously enforced, which results in fast convergence of the numerical scheme. In the case of a lamellar grating, the permittivity function that characterizes the grating is described by a piecewise function, but physical principles impose the continuity of some electromagnetic-field components at material boundaries. Describing this using a Fourier basis can be tricky because such basis functions do not efficiently take into account localized phenomena, unlike the MMGE. The MMGE is based on a decomposition of the study area into homogeneous areas. In each region, functions representing each electromagneticfield component can be expanded as a linear combination of Gegenbauer polynomials, and thereby the continuity conditions of certain electromagnetic-field components are explicitly enforced. This method has been successfully extended to nonperiodic structures by the introduction of complex coordinates, which allow the simulation of perfectly matched layers (PMLs) [3]. Until now, the MMGE has only been applied to Maxwell’s equations in Cartesian-coordinate form. We discuss in this paper a generalization of the method to cylindrical structures with arbitrary cross sections. We illustrate the method by calculating the eigenmodes of waveguides with arbitrary cross sections. The structure can be solid, dielectric, metallic, or hollow with a perfectly conducting wall. Analytical solutions for the eigenmodes of an object having axial invariance exist only in specific cases such as a cylindrical rod with a circular cross section. In this simple case, specific solutions are sought as the zeros of a transcendental equation. Without going into details, we think that it is easier to search for the modes of a specific structure by solving a matrix eigenvalue equation rather than searching 1084-7529/14/040667-10$15.00/0
for the zeros of a transcendental equation, especially if these solutions are in the complex plane. Almost all programming software now includes fairly robust and sophisticated algorithms for the fast computation of eigenvalues and eigenvectors of a matrix. In this article, we establish a modal equation for these structures by rigorously taking into account the boundaries of the structure. These equations consist of Maxwell’s equations in their tensorial form and are written in a coordinate system adapted to the contour of the waveguide. This system of coordinates has been used previously in [4]. Plane-wave scattering by a perfectly conducting elliptic cylinder was first considered by Sieger [5], who provided little numerical evaluation. In 1938, Morse and Rubenstein [6] proposed a numerical evaluation of the scattering of a perfectly conducting ellipse with a minor diameter equal to zero, that is, a perfectly conducting thin strip. It was not until 1963 that Yeh [7] found a solution to the diffraction of waves by a dielectric elliptical cylinder or by a dielectric ribbon for normal incidence. Yeh solved the problem in an elliptic cylindrical coordinate system and applied the orthogonality properties of the Mathieu functions. The method was extended to oblique incidence a year later [8]. A numerical evaluation of Mathieu functions can lead to numerical instabilities, and to remove these difficulties, some authors [9–12] have suggested the use of a scale transformation that involves finding a coordinate system in which the solution to the elliptic cylinder or ellipsoidal target scattering is well known. This transformation generally converts elliptic cylinder or ellipsoidal targets into a circle or sphere, and the inverse transform to the original coordinate system allows a solution to the scattering problem to be computed. The method presented in this paper is based on the curvilinear coordinate method (C method) [13–15], which relies on converting a rough surface into a plane through a change of coordinates. In the present method, a coordinate system is © 2014 Optical Society of America
668
J. Opt. Soc. Am. A / Vol. 31, No. 4 / April 2014
defined from the profile of the arbitrary cross section of the waveguide. This system of coordinates also converts the profile of the waveguide cross section into a circle as a scale transformation. However, the coordinate system used here is different from that of the scale transformation method, and the numerical solutions (i.e., the eigenmodes of the waveguide) are evaluated in this new coordinate system. In Section 2, we begin by presenting the covariant form of Maxwell’s equations and show how the introduction of Hertz potentials leads elegantly to a modal equation. We then present the coordinate transform and the metric tensor generated by this transformation. In doing this, we also briefly introduce a useful transformation (a complex coordinate transform) that allows the simulation of PMLs in MMGE [3]. While keeping in mind the basic principles of the MMGE, we show in Section 3 how this numerical method is suitable for solving the problem of a waveguide with an arbitrary cross section. All the matrix relations are presented in this section. In Section 4, we present the numerical results for the case of a circular perfectly hollow conducting waveguide to show the influence of the number of Gegenbauer polynomials on the convergence of the numerical scheme. We also study the impact of the PMLs through an analysis of the eigenmodes of a
(
Edee et al.
p ξijk ∂j H k iωε ggij E j ; p ξijk ∂j E k −iωμ ggij H j
i; j; k 1; 2; 3;
where ω is the angular frequency and gij represent contravariant components of the metric tensor obtained by inverting the matrix of covariant components gij , with g being the determinant of gij . We assume in the following that the ε and μ of the materials in which these equations are solved are both constant. If the elements of these matrices do not depend on the variable x3 and g33 is constant, then it can be shown that: 1. The E 3 and H 3 components satisfy the following propagation equation: p p ∂i ggij ∂j F k2 gF 0;
3
where F j Ej and Gj iZH j or F j H j and Gj iE j ∕Z, and Z is the wave impedance. At this point, we introduce the Hertz potentials Πe (electric potential) and Πh (magnetic potential) for reasons that will become clear later:
2. STATEMENT OF PROBLEM AND COVARIANT FORM OF MAXWELL’S EQUATIONS Let us consider, in the usual cylindrical coordinates system ρ; θ; z, an object with an arbitrary cross section, the border of which C is described by a function ρ Pθ. This object is assumed to be invariant along the z axis. Because of this invariance, it can support eigenmodes that can be computed as eigenfunctions of a matrix equation. Our goal is to establish this eigenvalue matrix equation, and to do so, we need to write rigorous continuity equations of some of the electromagneticfield components along the border C. Therefore, we propose a resolution of Maxwell’s equations in a coordinate system that is adapted to the object profile described by the function P. We start from the covariant form of Maxwell’s equations, since in this form these equations do not depend on any coordinate system. Consider the three-dimensional vector space R3 subtended by the basis vectors e1 ; e2 ; e3 . At any point xx1 ; x2 ; x3 in a medium whose electromagnetic properties are a permittivity of ε and a permeability of μ, Maxwell’s equations can be written as
(2)
p where F refers to E3 or H 3 , and k ω με is the wavenumber. 2. All other components of the field can be written in terms of E 3 , H 3 , and their derivatives ∂3 E 3 and ∂3 H 3 :
p ∂23 k2 g33 F 1 ∂1 ∂3 k2 g13 F 3 − k gg21 ∂1 g22 ∂2 g23 ∂3 G3 ; p ∂23 k2 g33 F 2 ∂2 ∂3 k2 g23 F 3 k gg11 ∂1 g12 ∂2 g13 ∂3 G3
circular dielectric waveguide. This preliminary work also allows us to compare our results with those in the literature. The case of an elliptical structure completes the presentation of the numerical results.
(1)
E3 ∂23 k2 g33 Πe H 3 ∂23 k2 g33 Πh :
4
In terms of the potentials, Eq. (2) can be written as follows: p p ∂i ggij ∂j k2 gΠe;h 0;
(5)
and Eq. (3) as p 8 E 1 ∂1 ∂3 k2 g13 Πe − iZk gg21 ∂1 g22 ∂2 g23 ∂3 Πh > > > < E ∂ ∂ k2 g Π iZkpgg11 ∂ g12 ∂ g13 ∂ Π 2 2 3 23 e 1 2 3 h 2 g Π i k p 21 ∂ g22 ∂ g23 ∂ Π > g g H ∂ ∂ k > 1 1 3 13 h 1 2 3 e Z > : p H 2 ∂2 ∂3 k2 g23 Πh − i Zk gg11 ∂1 g12 ∂2 g13 ∂3 Πe : 6
3. COORDINATE TRANSFORM AND EQUATION OF PROPAGATION We introduce here the coordinate transform that is adapted to the profile C. The transform is from the classical cylindrical coordinate system x1 ; x2 ; x3 ρ; θ; z as follows:
Edee et al.
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
669
the field including the radiation conditions at infinity in the case of an isolated object.
4. FRAMEWORK OF THE MMGE: APPLICATION TO CYLINDRICAL COORDINATES
Fig. 1. Geometry of the problem. Example of the coordinate transform of the profile described by Pθ into concentric circles.
( x1′ u ρ0 ρ Px2 x2′ x2 θ x3′ x3 z:
7
Recalling that P denotes the function describing the profile C, the coordinate surfaces x1′ constant are cylinders of a normalized radius constant∕ρ0 (Fig. 1), and the profile C is converted into a cylinder of a unit normalized radius in the coordinate system x1′ ; x2′ ; x3′ . This coordinate transform produces the following metric tensor: 2 gi′ j′
f 2 x2′
6 2′ 1′ _ 2′ 6 4 f x f x x 0
f_ x2′ f x2′ x1′ 2′ f x x1′ 2 f_ x2′ x1′ 2 0
0
3
7 07 5; (8) 1
where f P∕ρ0 and f_ df ∕dθ. We can now express the conp ′ ′ travariant elements ggi j that appear in Eq. (5) as 2 6 p i′ j′ gg 6 4
_ 2′ x1′ 1 a_ 2 x2′ −ax
3
0
_ 2′ −ax
1 x1′
0
0
0
f x2′ f x2′ x1′
7 7; (9) 5
_ with a_ f_ ∕f P∕P and P_ dP∕dθ. Using this expression, Eq. (5) takes the following form: _ 2′ x1′ ∂1′ ∂22′ k2 f f x1′ x1′ Πe;h _ 1′ ∂1′ x1′ ∂1′ − ∂2′ a_ a∂ 1 a_ ax −f f x1′ x1′ ∂23′ Πe;h ;
(10)
which may be rewritten as (11)
From Eqs. (6) and (11), the electromagnetic-field components that are tangential to the coordinate surfaces x1′ constant are written as 8 > E −ΔT k2 Πe > > 3′ > > < H ′ −ΔT k2 Πh 3 ; > E 2′ ∂2′ Π′e ikZ∂~ 1′ a_ 2 ∂~ 1′ − a∂ _ 2′ Πh > > > > : H ∂ Π′ − i k ∂~ a_ 2 ∂~ − a∂ _ Π ; 2′
h
Z
1′
1′
2′
12
e
where ∂~ 1′ x1 ∂1′ and Π′e;h ∂3′ Πe;h . In the following section, we propose a numerical solution to these equations based on the MMGE. In particular, we establish and solve an eigenvalue equation that takes into account the boundary conditions of ′
~ uu χ − iηu − u1 u1 ;
u ∈ u1 ; u2 :
(13)
The parameter χ controls the scale factor in the interval u1 ; u2 , while the parameter η controls the attenuation of the field in this interval. A reader who desires further explanation on these PMLs should refer to [3]. In addition, all the calculations related to the complex variable are presented in Appendix A. Taking into account this change of variables can be done intelligently through the definition of an appropriate inner product. B. Second Step: Polynomial Expansion and Undersized Matrix Relations The field functions depend on the three variables u; θ; z, and the z dependence was described in Section 2. Before discussing the second stage of the MMGE, which concerns only the radial variable u, we briefly discuss the θ dependence of these _ and a_ [symbolifunctions. The profile functions, namely P, P, cally denoted by Θθ], are periodic with respect to θ: Θθ 2πm Θθ, where m is an integer. This allows the following Fourier series expansion: M X
jΘi
ΔT Πe;h ∂23′ Πe;h 0:
2′
A. First Step: Partitioning of the Study Area The first step of the MMGE involves partitioning the computational region Ω into homogeneous subregions Ωj . Without any loss of generality, we illustrate the method using a canonical model that consists of two homogeneous layers Ωj , i ∈ f1; 2g, of which the electromagnetic properties are εj and μj (Fig. 1). The two media are bounded by concentric cylinders of radii u1 and u2 (u2 > u1 ), with borders denoted by ω1 and ω2 . For the problem under consideration, an electric- or magnetic-wall condition is always imposed along the boundary of the computational domain, namely ω2 . If we need to include radiation conditions at infinity in the case of an isolated object filling region Ω1 , then PMLs may be introduced in region Ω2 through a complex coordinate transform:
Θm jem i;
(14)
m−M
where em θ e−imθ . Hence, the Hertz potentials and the electromagnetic-field components have the same periodicity with respect to θ. The second step of the MMGE consists in approximating the restriction of the eigenfunctions, namely jΠj θ; zi in each homogeneous region Ωj , by a linear combination of Gegenbauer polynomials C Λn : jΠj zi e−iγz
N M X X
Πjnm jBΛnm i;
(15)
n1 m−M
where jBΛnm i jC Λn i ⊗ jem i. Here, the symbol ⊗ denotes the ′ Kronecker product, and Πj represents both Πje;h and Πje;h . Equations (14) and (15) are included in the partial differential equation of Eq. (10), and the resulting equation is projected
670
J. Opt. Soc. Am. A / Vol. 31, No. 4 / April 2014
Edee et al.
onto Qj × 2M 1 basis functions BΛnm , where Qj is the number of polynomial functions and 2M 1 is the number of Fourier basis functions. At this stage, one must keep in mind that the solution to the problem at hand needs to be that of a system of equations consisting of not only the partial differential equation of Eq. (10) but also the equations resulting from the boundary conditions of the electromagnetic field along the azimuthal (θ) direction. In preparation for the inclusion of the boundary conditions equations, Qj is chosen such that Qj < N, where N denotes the number of polynomials C Λn . For a given domain Ωj , the number Qj depends on the number of boundary conditions equations and also on the number of borders. In the two-media example (Ω1 and Ω2 ) chosen here, Ω1 has only one border ω1 , and Eq. (10) is projected onto N − 1 × 2M 1 basis functions, whereas Ω2 is limited by two borders ω1 and ω2 , and the partial differential equation solved in this medium is projected onto N − 2 × 2M 1 functions. The right-hand side of Eq. (10) takes the form R G2 ⊗ P2 ;
(16)
regardless of the medium, and the left-hand side is given by _ 2 − D ⊗ _a2 Dθ Dθ a_ 2 L DG−1 0 D ⊗ Iθ a G0 ⊗ Dθ 2 k2 εG2 ⊗ P2 ;
(17)
where Gp hC Λm ; xp C Λn i, p 0; 2; D hC Λm ; xdC Λn ∕dxi; Dθ diag−im; a_ and P are Toeplitz matrices formed by the Fourier coefficients of the function a_ and P, respectively; and Iθ is the identity matrix with respect to θ. Both Eqs. (16) and (17) show that all the matrices on the left-hand side of the Kronecker product are related to the radial variable u, while those on the right-hand side are related to θ. The matrices on the left-hand side should strictly be assigned with the superscript j, but for the sake of simplicity, this superscript is omitted. For a medium Ωj , Eq. (10) can be rewritten as a block matrix of dimension Qj × 2M 1 × N × 2M 1 that involves both the potentials Πje and Πjh and their derivatives Π′je and Π′jh : 2
Lj 0 0 0
32
Πje
3
2
0 0 Rj 0
32
Πje
3
6 76 7 6 76 7 6 0 0 0 Rj 7 6 Πj 7 6 0 L j 0 0 7 6 Πj 7 6 76 h 7 6 76 h 7 76 7 γ 6 j 6 76 7; 6 I 0 0 0 76 Π′je 7 6 0 0 Ij 0 76 Π′je 7 4 54 5 4 54 5 0 0 0 Ij 0 Ij 0 0 Π′jh Π′jh
γ
Φ21∶N;1∶2M1 L2left " 1 #" 1 # Φ1∶N;1∶2M1 Lright 0 0
L2right
Φ21∶N;1∶2M1
Ljleft
8 2 ′ > > > ∂3′ E3′ −ΔT k g33 Πe > > < ∂ ′ H ′ −Δ k2 g Π′ m T 33 3 3 : > ∂3′ E2′ −∂2′ ΔT Πe ikZ∂~ 1′ a_ 2 ∂~ 1′ − a∂ _ 2′ Π′m > > > > : ∂ ′ H ′ −∂ ′ ΔT Πm − i k ∂~ ′ a_ 2 ∂~ ′ − a∂ _ ′ Π′e 3
2
2
Z
1
1
20
2
This is possible because the operator ∂3′ is a constant. These equations then allow us to express Φ1N;1∶2M1 ; Φ2N−1;N;1∶2M1 t in terms of Φ11∶N−1;1∶2M1 ; Φ21∶N−2;1∶2M1 t as "
#
Φ1N;1∶2M1
" T
Φ2N−1;N;1∶2M1
Φ11∶N−1;1∶2M1
# :
Φ21∶N−2;1∶2M1
(21)
Furthermore, the vectors Φ11∶N;1∶2M1 ; Φ21∶N;1∶2M1 t can be expressed in terms of Φ11∶N−1;1∶2M1 ; Φ21∶N−2;1∶2M1 t through the following matrix relation:
(18)
Φ11∶N;1∶2M1
#
Φ21∶N;1∶2M1
" C
Φ11∶N−1;1∶2M1 Φ21∶N−2;1∶2M1
# :
(22)
Finally, we obtain an expression for computing the eigenvalues γ and associated eigenvectors Φ11∶N−1;1∶2M1 ; Φ21∶N−2;1∶2M1 t of a matrix of dimension dim 4N − 12M 1 4N − 22M 1 × 4N − 12M 1 4N − 22M 1: " Lleft C
Φ11∶N−1;1∶2M1 Φ21∶N−2;1∶2M1
"
# γLright C
Φ11∶N−1;1∶2M1 Φ21∶N−2;1∶2M1
# : (23)
; The construction of the matrices C and T has been presented in detail in [2] for the case of a lamellar grating.
where Φj Πje ; Πjh ; Π′je ; Π′jh t ;
C. Third Step: Boundary Conditions in the Radial Direction and Eigenvalues Equation We describe in this subsection the boundary conditions that complement the undersized system established in Section 4.B. These equations express the continuity of the tangential components of the electromagnetic field across the coordinate surfaces u constant (u1 and u2 in the current example). The four tangential components are expressed in term of the Hertz potentials in Eq. (12). The informed reader will note that writing the continuity of the four tangential components will not satisfactorily express the coefficients Π1pq;Nm and Π2pq;N−1;Nm in terms of the other coefficients for both Πje;h functions and their derivatives ∂3′ Πje;h −iγΠe;h Π′je;h . Thus, we add the continuity of the following functions to these equations:
"
where Ij Gj0 ⊗ Iθ . For Ω1 and Ω2 , Q1 N − 1 and Q2 N − 2, respectively, and we obtain the following undersized matrix relation: " 1 #" 1 # Φ1∶N;1∶2M1 Lleft 0 0
are vectors whose elements are the components of the potentials in the tensorial basis jBΛnm i [refer to Eq. (15)]. In the next subsection, the expansion coefficients Π1pq;Nm and Π2pq;N−1;Nm are expressed in terms of lower-order coefficients using the boundary conditions.
(19)
is the matrix on the left-hand side of Eq. (18), while and Ljright is the matrix on the right-hand side. Both Πje;h and Π′je;h
5. NUMERICAL RESULTS In this section, we study the performance of the numerical method in terms of convergence for various values of
Edee et al.
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
N (number of polynomials) and 2M 1 (number of Fourier harmonics). Before we do this, we briefly outline the main steps in the numerical modeling. For an arbitrary object of cross section C described in the cylindrical coordinate system ρ; θ; z: • We define a coordinate system u; θ; z as in Eq. (7) that converts C into a circle. • The functions of the problem then depend on two independent variables u and θ. The u behavior of these functions is approximated by a linear combination of Gegenbauer polynomials, while the periodicity of the problem with respect to the variable θ needs an expansion in a Fourier basis. • The boundary conditions in the radial direction u are then taken into account through the MMGE. One can expect that the convergence of the polynomial series will depend on the radius of the cylinder, while the choice of the number of Fourier harmonics will be highly dependent on the differentiability of the profile C. Thus, we first investigate the convergence with respect to N by considering a perfectly conducting cylinder with well-known eigenmodes. A. Influence of the Number of Polynomials on the Convergence: Case of a Hollow Perfectly Conducting Circular Waveguide The propagation of waves in a hollow perfectly conducting circular waveguide provides a good illustration of the influence of the radial variable u on the convergence of the numerical solutions. Conventionally, solutions to this problem are computed by considering a scalar potential Πnm solution of the form imθ iγ nm z Πnm u; θ; z J m knm e ; ρ ue
(24)
where J m is the mth-order Bessel function. For a given value of m, J m and its derivative J ′m have a denumerably infinite number of zeros that are ordered with respect to the subscript n. Furthermore, knm and γ satisfy the separation equation: ρ 2 2 2 knm ρ γ nm k :
(25)
Some tangential components of the electromagnetic field, namely, E3 in the transverse magnetic (TM) case and E 2 in the transverse electric (TE) case, must vanish at the conducting walls u ρ0 . Hence, we must have J nm knm ρ ρ0 0 and J ′nm knm ρ ρ0 0 in the TM and TE cases, respectively, from which knm may be determined from the zeros of the ρ Bessel functions or their derivatives. The propagation constant γ nm is then computed using Eq. (25). In the numerical approach presented in this paper, we first calculate γ N nm as an eigenvalue of Eq. (23) for a given value of N. Then, using Eq. (25), we compute the propagation constant in the radial direction kNnm , and so the zeros uN ρ nm of the Bessel functions and the zeros u′N nm of their derivatives are given by the formula q ′N 2 N 2 uN nm or unm ρ0 k − γ nm :
(26)
In the particular case of a hollow perfectly conducting circular waveguide, the dimension of the eigenvalue equation is
671
Table 1. First Five Zeros uref n0 of the Zeroth-order Bessel Function and Its Numerical Derivative u′ref n0 Computed from the Eigenvalues γ and for N 30a uref n0
n
uAbra n0
u′ref n0
u′Abra n0
1 2.4048255577 2.4048255577 3.8317059702 3.8317059702 2 5.5200781103 5.5200781103 7.0155866698 7.0155866698 3 8.6537279129 8.6537279129 10.1734681352 10.1734681351 4 11.7915344390 11.7915344391 13.3236919363 13.3236919363 5 14.9309177085 14.9309177086 16.4706300509 16.4706300509 a These values are used as reference values to study the convergence of the ′Abra are comparative values from [16] numerical results. Note that uAbra n0 and un0 (pp. 409–411). Values were calculated with ρ0 λ, ε 3.52 and ∂θ 0.
Table 2. First Five Zeros uref n2 of the Second-order Bessel Function and Its Numerical Derivative u′ref n2 Computed from the Eigenvalues γ and for N 30a n
uref n2
uAbra n2
u′ref n2
u′Abra n2
1 2 3 4 5
5.1356223018 8.4172441404 11.6198411721 14.7959517823 17.9598194950
5.13562 8.41724 11.61984 14.79595 17.95982
3.0542369282 6.7061331942 9.9694678230 13.1703708560 16.3475223183
3.05424 6.70613 9.96947 13.17037 16.34752
a These values are used as reference values to study the convergence of and u′Abra are comparative values the numerical results. Note that uAbra n2 n2 from [16] (pp. 409–411). Values were calculated with u1 λ, ε 3.52 and ∂θ −2i.
N − 1 × N − 1, and in this case, it can be shown that the TM mode, completely characterized by the magnetic potential Πh , and a TE mode described by Πe are not coupled by the boundary conditions. It is therefore not necessary to use the derivatives ∂3 Πe and ∂3 Πh of these functions, and so the second-order eigenvalue equation, Eq. (10) with an unknown γ 2 . In addition, a_ 0, which means that there is no coupling of the θ variable. Equation (10) can be solved for one and only one fixed value of the integer m ∂2 −im. To check the convergence of the method, we began by setting the values calculated using the MMGE for a sufficiently high N of 30 and listed in Tables 1 and 2 as reference values. In the following, these values are denoted by the superscript ref. Tables 1 and 2 list the first five zeros uref nm of the mth Bessel function and its derivative u′ref nm as computed from the eigenvalues γ N nm (N 30) and for m 0 and m 2, respectively. We note that there is excellent agreement between the values calculated with the MMGE and those reported in the ′Abra literature, uAbra [16]. To check the relative error nm and unm N ′N in unm and unm , we compare the number of common digits between these values and the reference values uref nm and u′ref nm ; the results are shown in Figs. 2(a), 2(b), 3(a), and 3(b). It can be seen that our results converge rather rapidly with respect to N to the reference values for both TE and TM polarizations. For given geometrical and physical parameters, namely a radius ρ0 and a relative permittivity ε1 (here ρ0 λ and ε1 3.52 ), we note that the convergence of the results strongly depends on the order n of the zero of the Bessel function and its derivative. This can be easily explained by noting that the argument of the Bessel functions depends inter alia on the propagation constant kNnm along the radial direction ρ [see Eq. (24)]. A higher value of n results in a higher number
672
J. Opt. Soc. Am. A / Vol. 31, No. 4 / April 2014
Edee et al.
0
n=1 n=2 n=3 n=4 n=5
−6
0.6
−8
0.4
−10
0.2
−12
0 −0.2
−14 0
5
10
15
20
25
0
30
0
uref=2.404825557691049 10
0
0.2
0.4
0.6
0.8
1
u
N−1
1.5
0
n=1 n=2 n=3 n=4 n=5
−2 ′ ref Log10|1−u′n0/un0 |
10
N=3 N=4 N=5
0.8
n0
−4
J (urefu/ρ ) 1
10
n0
Log |1−u /uref|
−2
1.2
−4
ref 50
J (u u/ρ )
ref
u50 =14.930917708487028
0
0
N=11 N=13 N=16
1
−6 0.5
−8 −10
0
−12 −14 0
−0.5
5
10
15
20
25
30
0
0.2
0.4
0.6
′N Fig. 2. Relative error in uN nm and unm , that is, the number of digits in ′ref common with the reference values uref nm and unm in Table 1 (∂θ 0).
0
n=1 n=2 n=3 n=4 n=5
−2
1
−4
Fig. 4. Comparison of the zeroth-order Bessel functions and the eigenfunctions computed from the eigenvalue problem for different numbers of polynomials.
of oscillations in the Bessel functions. Therefore, their description requires polynomials of a high degree. Figures 4(a) and 4(b) compare the zeroth-order Bessel function
ref
Log10|1−un2/un2 |
0.8
u
N−1
−6 −8
J0
−10 −12
ref unm u ρ0
27
with the eigenfunction computed from the eigenvalue problem
−14 −16 0
5
10
15
20
25
30
ΠN nm u
N−1
N X
Λ Πnm p C p u;
(28)
p1
0
n=1 n=2 n=3 n=4 n=5
Log |1−u′ /u′ ref| 10 n2 n2
−2 −4 −6 −8 −10 −12 −14 −16 0
5
10
15
20
25
30
N−1 ′N Fig. 3. Relative error in uN nm and unm , that is, the number of digits in ′ref common with the reference values uref nm and unm in Table 2 (∂θ 2i).
where Πnm p are the components of the eigenvector associated N 2 N 2 2 with the eigenvalue γ N nm such that γ nm unm k , and ref N unm is the limit of unm as N tends to infinity. In Fig. 4, ref uref nm 2.4048, and there is little variation in J 0 unm ∕ρ0 u, and so a small number of polynomials allows the series in Eq. (28) to converge quickly to Eq. (27). In Fig. 4(b), uref nm 14.9309, and the rapid oscillations of this function requires that polynomials of a higher order than the previous case are introduced into Eq. (28). B. Influence of the PMLs Function: Eigenfunctions of a Circular Dielectric Waveguide To highlight the impact of certain numerical parameters such as the scaling factor χ on the convergence of the results, we consider the case of a dielectric circular waveguide
Edee et al.
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
Table 3. Influence of the PMLs Functiona N max
N
χ1
χ2
χ3
20 36 52 68 84 100 116 132 148 164 180 196 212 228
4 6 8 10 12 14 16 18 20 22 24 26 28 30
1.383663803 1.237927700 1.206438154 1.201374359 1.200626785 1.200520140 1.200505270 1.200503244 1.200502973 1.200502938 1.200502933 1.200502933 1.200502933 1.200502936
1.000000000 1.325224960 1.235877091 1.209556137 1.202746613 1.201051758 1.200635895 1.200534541 1.200510070 1.200504232 1.200502857 1.200502537 1.200502463 1.200502446
1.000000000 1.387117537 1.280626926 1.227661681 1.209194780 1.203235020 1.201357454 1.200769413 1.200585410 1.200528037 1.200510272 1.200504817 1.200503157 1.200502655
a The effective index of the TE01 mode (γ ref 01 ∕k0 1.200502) of a circular dielectric waveguide with a radius of λ is computed with the MMGE with respect to the size of the eigenvalue matrix for three different values of the scaling factor χ. Values were calculated with u1 0.5λ, u2 4λ, ε1 2.5, and ε2 1.
673
Table 5. Convergence of the Effective Index γ 1 ∕k0 1.24235 of an Elliptical Dielectric Waveguide (Radii of ra 0.5λ and rb 0.45λ) with Respect to Both N and M max 2M 1a M max 3 5 7 9 11 13
N 11
N 13
N 15
N 17
1.242027987 1.242659231 1.242659231 1.242683061 1.242683061 1.242683175
1.241744579 1.242372635 1.242372635 1.242397679 1.242397679 1.242397772
1.241704827 1.242332430 1.242332430 1.242357753 1.242357753 1.242357849
1.241699366 1.242326903 1.242326903 1.242352297 1.242352297 1.242352394
Values were calculated with u1 0.5λ, u2 4r 1 , ε1 2.5, and ε2 1.
a
explained by the fact that using a scaling factor χ corresponds to enlarging the PMLs region by a factor of exactly 2χ. Consequently, the number of polynomials needed to describe the Hertz potentials increases.
(u1 0.5λ, ε1 2.5, and ε2 1). This case provides semianalytical solutions that can be useful for a comparison. Similar to the case of the circular perfectly conducting waveguide, solutions can be subscripted with m, n, where m is related to the θ dependence (∂θ im), and n is obtained from the zeros of a transcendental equation. Since a_ 0, it is possible to show that there is no coupling along the θ direction, and hence Eq. (10) can be solved for a fixed value of the integer m. The effective index of the TE01 mode computed with the semi-analytical method is equal to γ ref 01 ∕k0 1.200502 (k0 is the wavenumber of the vacuum), and Table 3 lists the results computed with the MMGE with respect to the size of the eigenvalue matrix for three different values of χ. In the current study, the computational domain is a circle of radius u2 4λ, and the scaling factor χ allows this domain to be enlarged from u2 to 2χu2 . This parameter works as a real PMLs [3]. In regard to these results, we first remark that the convergence with respect to N defines an N max value: the number of polynomial needed to describe the hole structure for any value of χ. The confinement of the field in the core of the waveguide suggests that it is quite possible to obtain satisfactory results without necessarily requiring large values of χ. An accuracy of 10−6 is obtained for χ 1 and N 20, and the same accuracy is obtained for χ 2 and χ 3 for large values of N or N max: N 24 and N 30, respectively. This can be
Table 4. Convergence of the Effective Index γ1 ∕k0 1.14623 of an Elliptical Dielectric Waveguide (Radii of ra 0.5λ and rb 0.45λ) with Respect to Both N and M max 2M 1a M max 3 5 7 9 11 13
N 11
N 13
N 15
N 17
1.169860902 1.147973728 1.147973728 1.148409721 1.148409721 1.148426112
1.169590596 1.146103877 1.146103877 1.146645363 1.146645363 1.146647215
1.169551905 1.145732438 1.145732438 1.146295794 1.146295794 1.146297765
1.169546524 1.145664811 1.145664811 1.146232876 1.146232877 1.146234887
Values were calculated with u1 0.5λ, u2 4r 1 , ε1 2.5, and ε2 1.
a
Fig. 5. Spatial distribution of the Ez [Fig. (5a)] and Hz [Fig. (5b)] components for γ 1 ∕k0 1.14623 in the case of an elliptical dielectric waveguide with radii of r a 0.5λ and r b 0.45λ.
674
J. Opt. Soc. Am. A / Vol. 31, No. 4 / April 2014
Edee et al.
C. General Case of a Differentiable Profile Function: Eigenfunctions of a Dielectric Waveguide with an Elliptical Cross Section In this subsection, we investigate the impact of the number of Fourier basis functions (M) on the convergence of the numerical results. This is performed within the framework of a dielectric waveguide with a cross section defined by an ellipse: b f θ p ; 1 − e2 cos2 θ
(29)
where the eccentricity e is linked to the axial radii r a and r b by the following relation:
e2 1 −
r 2b : r 2a
(30)
We assume that a_ does not vanish by considering the following numerical parameters: r a 0.5λ and r b 0.45λ. The electric and magnetic parameters remain unchanged. Since a_ does not vanish, the eigenfunctions are coupled along the θ direction through the boundary conditions at u u1 written with the components E 2′ and H 2′ [Eq. (12)] but also with the functions ∂3′ E2′ and ∂3′ H 2′ [Eq. (20)]. The computations of these components need the use of a finite number (2M 1) of _ Tables 4 and 5 show the convergence Fourier coefficients of a. with respect to both N and 2M 1 for two effective indices: γ 1 ∕k0 1.14623 and γ 2 ∕k0 1.24235, respectively. These results compared well with those obtained using the COMSOL software based on the finite element method. With this software we obtain γ 1 ∕k0 1.14622, γ 2 ∕k0 1.24235, and the computation area (circle of radius 50λ) is meshed with 17,632 elements. The results presented in Tables 4 and 5 show remarkable stability in terms of the convergence with respect to M and N. However, the rate of convergence of the eigenvalues differs according to the eigenmode that is observed. For γ 1 ∕k0 , the fifth digit (1.14623) stabilizes as soon as N reaches 17 and M 5, while values of N 15 and M 4 are necessary to reach the same accuracy for γ 2 ∕k0 . The field maps in Figs. 5(a), 5(b), 6(a), and 6(b) represent the spatial distribution of the E3 and H 3 components and help understand this finding.
6. CONCLUSION
Fig. 6. Spatial distribution of the Ez [Fig. (6a)] and Hz [Fig. (6b)] for γ 2 ∕k0 1.242352 in the case of an elliptical dielectric waveguide with radii of r a 0.5λ and r b 0.45λ.
We have presented a mode solver dedicated to waveguides with arbitrary cross sections. We used a change of coordinates to map the profile of the waveguide onto a circle and then introduced a Hertz potential to formulate the eigenvalue problem. We succeed in deriving a matrix eigenvalue equation by implementing the MMGE. The MMGE has once again shown its capability to take into account the continuity conditions of the electromagnetic field and allows us to build an operator incorporating these conditions. The method was successfully applied to calculate the eigenvectors and eigenvalues of cylindrical waveguides with arbitrary cross sections. Both the case of a hollow perfectly conducting waveguide and a dielectric waveguide were treated to demonstrate how the various numerical parameters are involved in the convergence of the results. It appears from these numerical studies that in the case of structures with circular cross sections, the number of polynomials describing the field components in the radial direction and the number of Fourier harmonics describing the azimuthal behavior of the field operate separately in the convergence of the electromagnetic field calculation. More specifically, once the derivative of the waveguide profile vanishes, a study of the convergence of the modes with respect to the number of Fourier harmonics is incongruous because in this case there is no coupling along the θ direction. However, when the derivative of the profile does not vanish, the convergence of the modes of the structure under consideration depends strongly on the differentiability of the profile for a given number of polynomials. In the present implementation, the numerical behavior of the MMGE, according to the differentiability of the profile of the waveguide cross section, is not very different from that of the traditional C method. The differentiability of the profile function determines the convergence of the Fourier series of the derivative function of the profile and thus the numerical scheme according to the azimuthal
Edee et al.
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
variable. At this stage, the real challenge is to extend the method to more general nondifferentiable profiles by introducing a parametric formulation for the azimuthal variable.
Consider the variables ζ belonging to −1; 1 and u ∈ u1 ; u2 defined as follows: uζ
APPENDIX A: SOME USEFUL OPERATORS IN THE INTERVAL −1;1 AND A CHANGE OF VARIABLES Gegenbauer polynomials are always defined in the interval −1; 1. Generally, the computation of all inner products presented in this paper needs a change of variables, and this directly affects the inner products. Before presenting expressions for the change of variables, we present the recursive relationships that allow these inner products to be calculated in the interval ζ ∈ −1; 1. For all ζ ∈ −1; 1, we get ζ
~ uu χ − iηu − u1 u1 ;
(A1)
~ uζ Aζ B;
(A6)
u2 − u1 ; A χ − iη 2
(A7)
and B
However, nC Λn ζ 2n Λ − 1ζC Λn−1 ζ − n 2Λ − 2C Λn−2 ζ;
u2 u1 2u1 2u1 − : u2 − u1 u2 − u1 χ − iηu2 − u1
u~ n 1 n − 1 2Λ Λ jC Λ i jC n−1 i: 2n Λ n1 2n Λ
d d ζ B ; du~ dζ
2 Λ Λ Λ Λ ~ Λn iuu hC Λm juC ~ 1 ;uu ~ 2 A hC m jζC n i BhC m jC n i > > > : Λ 2 Λ 3 Λ 2 Λ Λ Λ 2 Λ Λ hC m ju~ C n iuu ~ 1 ;uu ~ 2 A hC m jζ C n i 2BhC m jζC n i B hC m jC n i:
hC Λm jζC Λn i
A10
ACKNOWLEDGMENTS n 1 n − 1 2Λ Λ Λ hC Λ jC Λ i hC m jC n−1 i; 2n Λ m n1 2n Λ (A2)
This work was funded by grants from the French program investissement d’avenir managed by the National Research Agency (ANR), the European Commission (Auvergne FEDER funds), and the Region Auvergne in the framework of the LabEx IMobS3 (ANR-10-LABX-16-01).
n 1 n − 1 2Λ Λ hC Λ jζC Λn1 i hC m jζC Λn−1 i: 2n Λ m 2n Λ (A3)
REFERENCES
we get hC Λm jζ 2 C Λn i
(A9)
and
hD E E D Ei 8D Λ Λ dC Λn Λ ~ dC n Λ Λ dC n > j u A C jζ j B C C > m m m du~ uu dζ dζ > ~ 1 ;uu ~ 2
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
667
Mode solver based on Gegenbauer polynomial expansion for cylindrical structures with arbitrary cross sections Kofi Edee,1,2,* Mira Abboud,1,2 Gérard Granet,1,2 Jean Francois Cornet,1,2 and Nikolay A. Gippius1,2 1
Clermont Université, Université Blaise Pascal, Institut Pascal, BP 10448, F-63000 Clermont-Ferrand, France 2 CNRS, UMR 6602, Lasmea, F-63177 Aubière, France *Corresponding author: kofi.edee@univ‑bpclermont.fr Received November 15, 2013; revised January 17, 2014; accepted January 24, 2014; posted January 27, 2014 (Doc. ID 201081); published March 3, 2014
We present a modal method for the computation of eigenmodes of cylindrical structures with arbitrary cross sections. These modes are found as eigenvectors of a matrix eigenvalue equation that is obtained by introducing a new coordinate system that takes into account the profile of the cross section. We show that the use of Hertz potentials is suitable for the derivation of this eigenvalue equation and that the modal method based on Gegenbauer expansion (MMGE) is an efficient tool for the numerical solution of this equation. Results are successfully compared for both perfectly conducting and dielectric structures. A complex coordinate version of the MMGE is introduced to solve the dielectric case. © 2014 Optical Society of America OCIS codes: (000.3860) Mathematical methods in physics; (230.7370) Waveguides. http://dx.doi.org/10.1364/JOSAA.31.000667
1. INTRODUCTION The modal method by Gegenbauer expansion (MMGE) was recently introduced for plane-wave diffraction by a lamellar grating [1,2], and the key point of this method is that it allows transverse boundary conditions to be rigorously enforced, which results in fast convergence of the numerical scheme. In the case of a lamellar grating, the permittivity function that characterizes the grating is described by a piecewise function, but physical principles impose the continuity of some electromagnetic-field components at material boundaries. Describing this using a Fourier basis can be tricky because such basis functions do not efficiently take into account localized phenomena, unlike the MMGE. The MMGE is based on a decomposition of the study area into homogeneous areas. In each region, functions representing each electromagneticfield component can be expanded as a linear combination of Gegenbauer polynomials, and thereby the continuity conditions of certain electromagnetic-field components are explicitly enforced. This method has been successfully extended to nonperiodic structures by the introduction of complex coordinates, which allow the simulation of perfectly matched layers (PMLs) [3]. Until now, the MMGE has only been applied to Maxwell’s equations in Cartesian-coordinate form. We discuss in this paper a generalization of the method to cylindrical structures with arbitrary cross sections. We illustrate the method by calculating the eigenmodes of waveguides with arbitrary cross sections. The structure can be solid, dielectric, metallic, or hollow with a perfectly conducting wall. Analytical solutions for the eigenmodes of an object having axial invariance exist only in specific cases such as a cylindrical rod with a circular cross section. In this simple case, specific solutions are sought as the zeros of a transcendental equation. Without going into details, we think that it is easier to search for the modes of a specific structure by solving a matrix eigenvalue equation rather than searching 1084-7529/14/040667-10$15.00/0
for the zeros of a transcendental equation, especially if these solutions are in the complex plane. Almost all programming software now includes fairly robust and sophisticated algorithms for the fast computation of eigenvalues and eigenvectors of a matrix. In this article, we establish a modal equation for these structures by rigorously taking into account the boundaries of the structure. These equations consist of Maxwell’s equations in their tensorial form and are written in a coordinate system adapted to the contour of the waveguide. This system of coordinates has been used previously in [4]. Plane-wave scattering by a perfectly conducting elliptic cylinder was first considered by Sieger [5], who provided little numerical evaluation. In 1938, Morse and Rubenstein [6] proposed a numerical evaluation of the scattering of a perfectly conducting ellipse with a minor diameter equal to zero, that is, a perfectly conducting thin strip. It was not until 1963 that Yeh [7] found a solution to the diffraction of waves by a dielectric elliptical cylinder or by a dielectric ribbon for normal incidence. Yeh solved the problem in an elliptic cylindrical coordinate system and applied the orthogonality properties of the Mathieu functions. The method was extended to oblique incidence a year later [8]. A numerical evaluation of Mathieu functions can lead to numerical instabilities, and to remove these difficulties, some authors [9–12] have suggested the use of a scale transformation that involves finding a coordinate system in which the solution to the elliptic cylinder or ellipsoidal target scattering is well known. This transformation generally converts elliptic cylinder or ellipsoidal targets into a circle or sphere, and the inverse transform to the original coordinate system allows a solution to the scattering problem to be computed. The method presented in this paper is based on the curvilinear coordinate method (C method) [13–15], which relies on converting a rough surface into a plane through a change of coordinates. In the present method, a coordinate system is © 2014 Optical Society of America
668
J. Opt. Soc. Am. A / Vol. 31, No. 4 / April 2014
defined from the profile of the arbitrary cross section of the waveguide. This system of coordinates also converts the profile of the waveguide cross section into a circle as a scale transformation. However, the coordinate system used here is different from that of the scale transformation method, and the numerical solutions (i.e., the eigenmodes of the waveguide) are evaluated in this new coordinate system. In Section 2, we begin by presenting the covariant form of Maxwell’s equations and show how the introduction of Hertz potentials leads elegantly to a modal equation. We then present the coordinate transform and the metric tensor generated by this transformation. In doing this, we also briefly introduce a useful transformation (a complex coordinate transform) that allows the simulation of PMLs in MMGE [3]. While keeping in mind the basic principles of the MMGE, we show in Section 3 how this numerical method is suitable for solving the problem of a waveguide with an arbitrary cross section. All the matrix relations are presented in this section. In Section 4, we present the numerical results for the case of a circular perfectly hollow conducting waveguide to show the influence of the number of Gegenbauer polynomials on the convergence of the numerical scheme. We also study the impact of the PMLs through an analysis of the eigenmodes of a
(
Edee et al.
p ξijk ∂j H k iωε ggij E j ; p ξijk ∂j E k −iωμ ggij H j
i; j; k 1; 2; 3;
where ω is the angular frequency and gij represent contravariant components of the metric tensor obtained by inverting the matrix of covariant components gij , with g being the determinant of gij . We assume in the following that the ε and μ of the materials in which these equations are solved are both constant. If the elements of these matrices do not depend on the variable x3 and g33 is constant, then it can be shown that: 1. The E 3 and H 3 components satisfy the following propagation equation: p p ∂i ggij ∂j F k2 gF 0;
3
where F j Ej and Gj iZH j or F j H j and Gj iE j ∕Z, and Z is the wave impedance. At this point, we introduce the Hertz potentials Πe (electric potential) and Πh (magnetic potential) for reasons that will become clear later:
2. STATEMENT OF PROBLEM AND COVARIANT FORM OF MAXWELL’S EQUATIONS Let us consider, in the usual cylindrical coordinates system ρ; θ; z, an object with an arbitrary cross section, the border of which C is described by a function ρ Pθ. This object is assumed to be invariant along the z axis. Because of this invariance, it can support eigenmodes that can be computed as eigenfunctions of a matrix equation. Our goal is to establish this eigenvalue matrix equation, and to do so, we need to write rigorous continuity equations of some of the electromagneticfield components along the border C. Therefore, we propose a resolution of Maxwell’s equations in a coordinate system that is adapted to the object profile described by the function P. We start from the covariant form of Maxwell’s equations, since in this form these equations do not depend on any coordinate system. Consider the three-dimensional vector space R3 subtended by the basis vectors e1 ; e2 ; e3 . At any point xx1 ; x2 ; x3 in a medium whose electromagnetic properties are a permittivity of ε and a permeability of μ, Maxwell’s equations can be written as
(2)
p where F refers to E3 or H 3 , and k ω με is the wavenumber. 2. All other components of the field can be written in terms of E 3 , H 3 , and their derivatives ∂3 E 3 and ∂3 H 3 :
p ∂23 k2 g33 F 1 ∂1 ∂3 k2 g13 F 3 − k gg21 ∂1 g22 ∂2 g23 ∂3 G3 ; p ∂23 k2 g33 F 2 ∂2 ∂3 k2 g23 F 3 k gg11 ∂1 g12 ∂2 g13 ∂3 G3
circular dielectric waveguide. This preliminary work also allows us to compare our results with those in the literature. The case of an elliptical structure completes the presentation of the numerical results.
(1)
E3 ∂23 k2 g33 Πe H 3 ∂23 k2 g33 Πh :
4
In terms of the potentials, Eq. (2) can be written as follows: p p ∂i ggij ∂j k2 gΠe;h 0;
(5)
and Eq. (3) as p 8 E 1 ∂1 ∂3 k2 g13 Πe − iZk gg21 ∂1 g22 ∂2 g23 ∂3 Πh > > > < E ∂ ∂ k2 g Π iZkpgg11 ∂ g12 ∂ g13 ∂ Π 2 2 3 23 e 1 2 3 h 2 g Π i k p 21 ∂ g22 ∂ g23 ∂ Π > g g H ∂ ∂ k > 1 1 3 13 h 1 2 3 e Z > : p H 2 ∂2 ∂3 k2 g23 Πh − i Zk gg11 ∂1 g12 ∂2 g13 ∂3 Πe : 6
3. COORDINATE TRANSFORM AND EQUATION OF PROPAGATION We introduce here the coordinate transform that is adapted to the profile C. The transform is from the classical cylindrical coordinate system x1 ; x2 ; x3 ρ; θ; z as follows:
Edee et al.
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
669
the field including the radiation conditions at infinity in the case of an isolated object.
4. FRAMEWORK OF THE MMGE: APPLICATION TO CYLINDRICAL COORDINATES
Fig. 1. Geometry of the problem. Example of the coordinate transform of the profile described by Pθ into concentric circles.
( x1′ u ρ0 ρ Px2 x2′ x2 θ x3′ x3 z:
7
Recalling that P denotes the function describing the profile C, the coordinate surfaces x1′ constant are cylinders of a normalized radius constant∕ρ0 (Fig. 1), and the profile C is converted into a cylinder of a unit normalized radius in the coordinate system x1′ ; x2′ ; x3′ . This coordinate transform produces the following metric tensor: 2 gi′ j′
f 2 x2′
6 2′ 1′ _ 2′ 6 4 f x f x x 0
f_ x2′ f x2′ x1′ 2′ f x x1′ 2 f_ x2′ x1′ 2 0
0
3
7 07 5; (8) 1
where f P∕ρ0 and f_ df ∕dθ. We can now express the conp ′ ′ travariant elements ggi j that appear in Eq. (5) as 2 6 p i′ j′ gg 6 4
_ 2′ x1′ 1 a_ 2 x2′ −ax
3
0
_ 2′ −ax
1 x1′
0
0
0
f x2′ f x2′ x1′
7 7; (9) 5
_ with a_ f_ ∕f P∕P and P_ dP∕dθ. Using this expression, Eq. (5) takes the following form: _ 2′ x1′ ∂1′ ∂22′ k2 f f x1′ x1′ Πe;h _ 1′ ∂1′ x1′ ∂1′ − ∂2′ a_ a∂ 1 a_ ax −f f x1′ x1′ ∂23′ Πe;h ;
(10)
which may be rewritten as (11)
From Eqs. (6) and (11), the electromagnetic-field components that are tangential to the coordinate surfaces x1′ constant are written as 8 > E −ΔT k2 Πe > > 3′ > > < H ′ −ΔT k2 Πh 3 ; > E 2′ ∂2′ Π′e ikZ∂~ 1′ a_ 2 ∂~ 1′ − a∂ _ 2′ Πh > > > > : H ∂ Π′ − i k ∂~ a_ 2 ∂~ − a∂ _ Π ; 2′
h
Z
1′
1′
2′
12
e
where ∂~ 1′ x1 ∂1′ and Π′e;h ∂3′ Πe;h . In the following section, we propose a numerical solution to these equations based on the MMGE. In particular, we establish and solve an eigenvalue equation that takes into account the boundary conditions of ′
~ uu χ − iηu − u1 u1 ;
u ∈ u1 ; u2 :
(13)
The parameter χ controls the scale factor in the interval u1 ; u2 , while the parameter η controls the attenuation of the field in this interval. A reader who desires further explanation on these PMLs should refer to [3]. In addition, all the calculations related to the complex variable are presented in Appendix A. Taking into account this change of variables can be done intelligently through the definition of an appropriate inner product. B. Second Step: Polynomial Expansion and Undersized Matrix Relations The field functions depend on the three variables u; θ; z, and the z dependence was described in Section 2. Before discussing the second stage of the MMGE, which concerns only the radial variable u, we briefly discuss the θ dependence of these _ and a_ [symbolifunctions. The profile functions, namely P, P, cally denoted by Θθ], are periodic with respect to θ: Θθ 2πm Θθ, where m is an integer. This allows the following Fourier series expansion: M X
jΘi
ΔT Πe;h ∂23′ Πe;h 0:
2′
A. First Step: Partitioning of the Study Area The first step of the MMGE involves partitioning the computational region Ω into homogeneous subregions Ωj . Without any loss of generality, we illustrate the method using a canonical model that consists of two homogeneous layers Ωj , i ∈ f1; 2g, of which the electromagnetic properties are εj and μj (Fig. 1). The two media are bounded by concentric cylinders of radii u1 and u2 (u2 > u1 ), with borders denoted by ω1 and ω2 . For the problem under consideration, an electric- or magnetic-wall condition is always imposed along the boundary of the computational domain, namely ω2 . If we need to include radiation conditions at infinity in the case of an isolated object filling region Ω1 , then PMLs may be introduced in region Ω2 through a complex coordinate transform:
Θm jem i;
(14)
m−M
where em θ e−imθ . Hence, the Hertz potentials and the electromagnetic-field components have the same periodicity with respect to θ. The second step of the MMGE consists in approximating the restriction of the eigenfunctions, namely jΠj θ; zi in each homogeneous region Ωj , by a linear combination of Gegenbauer polynomials C Λn : jΠj zi e−iγz
N M X X
Πjnm jBΛnm i;
(15)
n1 m−M
where jBΛnm i jC Λn i ⊗ jem i. Here, the symbol ⊗ denotes the ′ Kronecker product, and Πj represents both Πje;h and Πje;h . Equations (14) and (15) are included in the partial differential equation of Eq. (10), and the resulting equation is projected
670
J. Opt. Soc. Am. A / Vol. 31, No. 4 / April 2014
Edee et al.
onto Qj × 2M 1 basis functions BΛnm , where Qj is the number of polynomial functions and 2M 1 is the number of Fourier basis functions. At this stage, one must keep in mind that the solution to the problem at hand needs to be that of a system of equations consisting of not only the partial differential equation of Eq. (10) but also the equations resulting from the boundary conditions of the electromagnetic field along the azimuthal (θ) direction. In preparation for the inclusion of the boundary conditions equations, Qj is chosen such that Qj < N, where N denotes the number of polynomials C Λn . For a given domain Ωj , the number Qj depends on the number of boundary conditions equations and also on the number of borders. In the two-media example (Ω1 and Ω2 ) chosen here, Ω1 has only one border ω1 , and Eq. (10) is projected onto N − 1 × 2M 1 basis functions, whereas Ω2 is limited by two borders ω1 and ω2 , and the partial differential equation solved in this medium is projected onto N − 2 × 2M 1 functions. The right-hand side of Eq. (10) takes the form R G2 ⊗ P2 ;
(16)
regardless of the medium, and the left-hand side is given by _ 2 − D ⊗ _a2 Dθ Dθ a_ 2 L DG−1 0 D ⊗ Iθ a G0 ⊗ Dθ 2 k2 εG2 ⊗ P2 ;
(17)
where Gp hC Λm ; xp C Λn i, p 0; 2; D hC Λm ; xdC Λn ∕dxi; Dθ diag−im; a_ and P are Toeplitz matrices formed by the Fourier coefficients of the function a_ and P, respectively; and Iθ is the identity matrix with respect to θ. Both Eqs. (16) and (17) show that all the matrices on the left-hand side of the Kronecker product are related to the radial variable u, while those on the right-hand side are related to θ. The matrices on the left-hand side should strictly be assigned with the superscript j, but for the sake of simplicity, this superscript is omitted. For a medium Ωj , Eq. (10) can be rewritten as a block matrix of dimension Qj × 2M 1 × N × 2M 1 that involves both the potentials Πje and Πjh and their derivatives Π′je and Π′jh : 2
Lj 0 0 0
32
Πje
3
2
0 0 Rj 0
32
Πje
3
6 76 7 6 76 7 6 0 0 0 Rj 7 6 Πj 7 6 0 L j 0 0 7 6 Πj 7 6 76 h 7 6 76 h 7 76 7 γ 6 j 6 76 7; 6 I 0 0 0 76 Π′je 7 6 0 0 Ij 0 76 Π′je 7 4 54 5 4 54 5 0 0 0 Ij 0 Ij 0 0 Π′jh Π′jh
γ
Φ21∶N;1∶2M1 L2left " 1 #" 1 # Φ1∶N;1∶2M1 Lright 0 0
L2right
Φ21∶N;1∶2M1
Ljleft
8 2 ′ > > > ∂3′ E3′ −ΔT k g33 Πe > > < ∂ ′ H ′ −Δ k2 g Π′ m T 33 3 3 : > ∂3′ E2′ −∂2′ ΔT Πe ikZ∂~ 1′ a_ 2 ∂~ 1′ − a∂ _ 2′ Π′m > > > > : ∂ ′ H ′ −∂ ′ ΔT Πm − i k ∂~ ′ a_ 2 ∂~ ′ − a∂ _ ′ Π′e 3
2
2
Z
1
1
20
2
This is possible because the operator ∂3′ is a constant. These equations then allow us to express Φ1N;1∶2M1 ; Φ2N−1;N;1∶2M1 t in terms of Φ11∶N−1;1∶2M1 ; Φ21∶N−2;1∶2M1 t as "
#
Φ1N;1∶2M1
" T
Φ2N−1;N;1∶2M1
Φ11∶N−1;1∶2M1
# :
Φ21∶N−2;1∶2M1
(21)
Furthermore, the vectors Φ11∶N;1∶2M1 ; Φ21∶N;1∶2M1 t can be expressed in terms of Φ11∶N−1;1∶2M1 ; Φ21∶N−2;1∶2M1 t through the following matrix relation:
(18)
Φ11∶N;1∶2M1
#
Φ21∶N;1∶2M1
" C
Φ11∶N−1;1∶2M1 Φ21∶N−2;1∶2M1
# :
(22)
Finally, we obtain an expression for computing the eigenvalues γ and associated eigenvectors Φ11∶N−1;1∶2M1 ; Φ21∶N−2;1∶2M1 t of a matrix of dimension dim 4N − 12M 1 4N − 22M 1 × 4N − 12M 1 4N − 22M 1: " Lleft C
Φ11∶N−1;1∶2M1 Φ21∶N−2;1∶2M1
"
# γLright C
Φ11∶N−1;1∶2M1 Φ21∶N−2;1∶2M1
# : (23)
; The construction of the matrices C and T has been presented in detail in [2] for the case of a lamellar grating.
where Φj Πje ; Πjh ; Π′je ; Π′jh t ;
C. Third Step: Boundary Conditions in the Radial Direction and Eigenvalues Equation We describe in this subsection the boundary conditions that complement the undersized system established in Section 4.B. These equations express the continuity of the tangential components of the electromagnetic field across the coordinate surfaces u constant (u1 and u2 in the current example). The four tangential components are expressed in term of the Hertz potentials in Eq. (12). The informed reader will note that writing the continuity of the four tangential components will not satisfactorily express the coefficients Π1pq;Nm and Π2pq;N−1;Nm in terms of the other coefficients for both Πje;h functions and their derivatives ∂3′ Πje;h −iγΠe;h Π′je;h . Thus, we add the continuity of the following functions to these equations:
"
where Ij Gj0 ⊗ Iθ . For Ω1 and Ω2 , Q1 N − 1 and Q2 N − 2, respectively, and we obtain the following undersized matrix relation: " 1 #" 1 # Φ1∶N;1∶2M1 Lleft 0 0
are vectors whose elements are the components of the potentials in the tensorial basis jBΛnm i [refer to Eq. (15)]. In the next subsection, the expansion coefficients Π1pq;Nm and Π2pq;N−1;Nm are expressed in terms of lower-order coefficients using the boundary conditions.
(19)
is the matrix on the left-hand side of Eq. (18), while and Ljright is the matrix on the right-hand side. Both Πje;h and Π′je;h
5. NUMERICAL RESULTS In this section, we study the performance of the numerical method in terms of convergence for various values of
Edee et al.
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
N (number of polynomials) and 2M 1 (number of Fourier harmonics). Before we do this, we briefly outline the main steps in the numerical modeling. For an arbitrary object of cross section C described in the cylindrical coordinate system ρ; θ; z: • We define a coordinate system u; θ; z as in Eq. (7) that converts C into a circle. • The functions of the problem then depend on two independent variables u and θ. The u behavior of these functions is approximated by a linear combination of Gegenbauer polynomials, while the periodicity of the problem with respect to the variable θ needs an expansion in a Fourier basis. • The boundary conditions in the radial direction u are then taken into account through the MMGE. One can expect that the convergence of the polynomial series will depend on the radius of the cylinder, while the choice of the number of Fourier harmonics will be highly dependent on the differentiability of the profile C. Thus, we first investigate the convergence with respect to N by considering a perfectly conducting cylinder with well-known eigenmodes. A. Influence of the Number of Polynomials on the Convergence: Case of a Hollow Perfectly Conducting Circular Waveguide The propagation of waves in a hollow perfectly conducting circular waveguide provides a good illustration of the influence of the radial variable u on the convergence of the numerical solutions. Conventionally, solutions to this problem are computed by considering a scalar potential Πnm solution of the form imθ iγ nm z Πnm u; θ; z J m knm e ; ρ ue
(24)
where J m is the mth-order Bessel function. For a given value of m, J m and its derivative J ′m have a denumerably infinite number of zeros that are ordered with respect to the subscript n. Furthermore, knm and γ satisfy the separation equation: ρ 2 2 2 knm ρ γ nm k :
(25)
Some tangential components of the electromagnetic field, namely, E3 in the transverse magnetic (TM) case and E 2 in the transverse electric (TE) case, must vanish at the conducting walls u ρ0 . Hence, we must have J nm knm ρ ρ0 0 and J ′nm knm ρ ρ0 0 in the TM and TE cases, respectively, from which knm may be determined from the zeros of the ρ Bessel functions or their derivatives. The propagation constant γ nm is then computed using Eq. (25). In the numerical approach presented in this paper, we first calculate γ N nm as an eigenvalue of Eq. (23) for a given value of N. Then, using Eq. (25), we compute the propagation constant in the radial direction kNnm , and so the zeros uN ρ nm of the Bessel functions and the zeros u′N nm of their derivatives are given by the formula q ′N 2 N 2 uN nm or unm ρ0 k − γ nm :
(26)
In the particular case of a hollow perfectly conducting circular waveguide, the dimension of the eigenvalue equation is
671
Table 1. First Five Zeros uref n0 of the Zeroth-order Bessel Function and Its Numerical Derivative u′ref n0 Computed from the Eigenvalues γ and for N 30a uref n0
n
uAbra n0
u′ref n0
u′Abra n0
1 2.4048255577 2.4048255577 3.8317059702 3.8317059702 2 5.5200781103 5.5200781103 7.0155866698 7.0155866698 3 8.6537279129 8.6537279129 10.1734681352 10.1734681351 4 11.7915344390 11.7915344391 13.3236919363 13.3236919363 5 14.9309177085 14.9309177086 16.4706300509 16.4706300509 a These values are used as reference values to study the convergence of the ′Abra are comparative values from [16] numerical results. Note that uAbra n0 and un0 (pp. 409–411). Values were calculated with ρ0 λ, ε 3.52 and ∂θ 0.
Table 2. First Five Zeros uref n2 of the Second-order Bessel Function and Its Numerical Derivative u′ref n2 Computed from the Eigenvalues γ and for N 30a n
uref n2
uAbra n2
u′ref n2
u′Abra n2
1 2 3 4 5
5.1356223018 8.4172441404 11.6198411721 14.7959517823 17.9598194950
5.13562 8.41724 11.61984 14.79595 17.95982
3.0542369282 6.7061331942 9.9694678230 13.1703708560 16.3475223183
3.05424 6.70613 9.96947 13.17037 16.34752
a These values are used as reference values to study the convergence of and u′Abra are comparative values the numerical results. Note that uAbra n2 n2 from [16] (pp. 409–411). Values were calculated with u1 λ, ε 3.52 and ∂θ −2i.
N − 1 × N − 1, and in this case, it can be shown that the TM mode, completely characterized by the magnetic potential Πh , and a TE mode described by Πe are not coupled by the boundary conditions. It is therefore not necessary to use the derivatives ∂3 Πe and ∂3 Πh of these functions, and so the second-order eigenvalue equation, Eq. (10) with an unknown γ 2 . In addition, a_ 0, which means that there is no coupling of the θ variable. Equation (10) can be solved for one and only one fixed value of the integer m ∂2 −im. To check the convergence of the method, we began by setting the values calculated using the MMGE for a sufficiently high N of 30 and listed in Tables 1 and 2 as reference values. In the following, these values are denoted by the superscript ref. Tables 1 and 2 list the first five zeros uref nm of the mth Bessel function and its derivative u′ref nm as computed from the eigenvalues γ N nm (N 30) and for m 0 and m 2, respectively. We note that there is excellent agreement between the values calculated with the MMGE and those reported in the ′Abra literature, uAbra [16]. To check the relative error nm and unm N ′N in unm and unm , we compare the number of common digits between these values and the reference values uref nm and u′ref nm ; the results are shown in Figs. 2(a), 2(b), 3(a), and 3(b). It can be seen that our results converge rather rapidly with respect to N to the reference values for both TE and TM polarizations. For given geometrical and physical parameters, namely a radius ρ0 and a relative permittivity ε1 (here ρ0 λ and ε1 3.52 ), we note that the convergence of the results strongly depends on the order n of the zero of the Bessel function and its derivative. This can be easily explained by noting that the argument of the Bessel functions depends inter alia on the propagation constant kNnm along the radial direction ρ [see Eq. (24)]. A higher value of n results in a higher number
672
J. Opt. Soc. Am. A / Vol. 31, No. 4 / April 2014
Edee et al.
0
n=1 n=2 n=3 n=4 n=5
−6
0.6
−8
0.4
−10
0.2
−12
0 −0.2
−14 0
5
10
15
20
25
0
30
0
uref=2.404825557691049 10
0
0.2
0.4
0.6
0.8
1
u
N−1
1.5
0
n=1 n=2 n=3 n=4 n=5
−2 ′ ref Log10|1−u′n0/un0 |
10
N=3 N=4 N=5
0.8
n0
−4
J (urefu/ρ ) 1
10
n0
Log |1−u /uref|
−2
1.2
−4
ref 50
J (u u/ρ )
ref
u50 =14.930917708487028
0
0
N=11 N=13 N=16
1
−6 0.5
−8 −10
0
−12 −14 0
−0.5
5
10
15
20
25
30
0
0.2
0.4
0.6
′N Fig. 2. Relative error in uN nm and unm , that is, the number of digits in ′ref common with the reference values uref nm and unm in Table 1 (∂θ 0).
0
n=1 n=2 n=3 n=4 n=5
−2
1
−4
Fig. 4. Comparison of the zeroth-order Bessel functions and the eigenfunctions computed from the eigenvalue problem for different numbers of polynomials.
of oscillations in the Bessel functions. Therefore, their description requires polynomials of a high degree. Figures 4(a) and 4(b) compare the zeroth-order Bessel function
ref
Log10|1−un2/un2 |
0.8
u
N−1
−6 −8
J0
−10 −12
ref unm u ρ0
27
with the eigenfunction computed from the eigenvalue problem
−14 −16 0
5
10
15
20
25
30
ΠN nm u
N−1
N X
Λ Πnm p C p u;
(28)
p1
0
n=1 n=2 n=3 n=4 n=5
Log |1−u′ /u′ ref| 10 n2 n2
−2 −4 −6 −8 −10 −12 −14 −16 0
5
10
15
20
25
30
N−1 ′N Fig. 3. Relative error in uN nm and unm , that is, the number of digits in ′ref common with the reference values uref nm and unm in Table 2 (∂θ 2i).
where Πnm p are the components of the eigenvector associated N 2 N 2 2 with the eigenvalue γ N nm such that γ nm unm k , and ref N unm is the limit of unm as N tends to infinity. In Fig. 4, ref uref nm 2.4048, and there is little variation in J 0 unm ∕ρ0 u, and so a small number of polynomials allows the series in Eq. (28) to converge quickly to Eq. (27). In Fig. 4(b), uref nm 14.9309, and the rapid oscillations of this function requires that polynomials of a higher order than the previous case are introduced into Eq. (28). B. Influence of the PMLs Function: Eigenfunctions of a Circular Dielectric Waveguide To highlight the impact of certain numerical parameters such as the scaling factor χ on the convergence of the results, we consider the case of a dielectric circular waveguide
Edee et al.
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
Table 3. Influence of the PMLs Functiona N max
N
χ1
χ2
χ3
20 36 52 68 84 100 116 132 148 164 180 196 212 228
4 6 8 10 12 14 16 18 20 22 24 26 28 30
1.383663803 1.237927700 1.206438154 1.201374359 1.200626785 1.200520140 1.200505270 1.200503244 1.200502973 1.200502938 1.200502933 1.200502933 1.200502933 1.200502936
1.000000000 1.325224960 1.235877091 1.209556137 1.202746613 1.201051758 1.200635895 1.200534541 1.200510070 1.200504232 1.200502857 1.200502537 1.200502463 1.200502446
1.000000000 1.387117537 1.280626926 1.227661681 1.209194780 1.203235020 1.201357454 1.200769413 1.200585410 1.200528037 1.200510272 1.200504817 1.200503157 1.200502655
a The effective index of the TE01 mode (γ ref 01 ∕k0 1.200502) of a circular dielectric waveguide with a radius of λ is computed with the MMGE with respect to the size of the eigenvalue matrix for three different values of the scaling factor χ. Values were calculated with u1 0.5λ, u2 4λ, ε1 2.5, and ε2 1.
673
Table 5. Convergence of the Effective Index γ 1 ∕k0 1.24235 of an Elliptical Dielectric Waveguide (Radii of ra 0.5λ and rb 0.45λ) with Respect to Both N and M max 2M 1a M max 3 5 7 9 11 13
N 11
N 13
N 15
N 17
1.242027987 1.242659231 1.242659231 1.242683061 1.242683061 1.242683175
1.241744579 1.242372635 1.242372635 1.242397679 1.242397679 1.242397772
1.241704827 1.242332430 1.242332430 1.242357753 1.242357753 1.242357849
1.241699366 1.242326903 1.242326903 1.242352297 1.242352297 1.242352394
Values were calculated with u1 0.5λ, u2 4r 1 , ε1 2.5, and ε2 1.
a
explained by the fact that using a scaling factor χ corresponds to enlarging the PMLs region by a factor of exactly 2χ. Consequently, the number of polynomials needed to describe the Hertz potentials increases.
(u1 0.5λ, ε1 2.5, and ε2 1). This case provides semianalytical solutions that can be useful for a comparison. Similar to the case of the circular perfectly conducting waveguide, solutions can be subscripted with m, n, where m is related to the θ dependence (∂θ im), and n is obtained from the zeros of a transcendental equation. Since a_ 0, it is possible to show that there is no coupling along the θ direction, and hence Eq. (10) can be solved for a fixed value of the integer m. The effective index of the TE01 mode computed with the semi-analytical method is equal to γ ref 01 ∕k0 1.200502 (k0 is the wavenumber of the vacuum), and Table 3 lists the results computed with the MMGE with respect to the size of the eigenvalue matrix for three different values of χ. In the current study, the computational domain is a circle of radius u2 4λ, and the scaling factor χ allows this domain to be enlarged from u2 to 2χu2 . This parameter works as a real PMLs [3]. In regard to these results, we first remark that the convergence with respect to N defines an N max value: the number of polynomial needed to describe the hole structure for any value of χ. The confinement of the field in the core of the waveguide suggests that it is quite possible to obtain satisfactory results without necessarily requiring large values of χ. An accuracy of 10−6 is obtained for χ 1 and N 20, and the same accuracy is obtained for χ 2 and χ 3 for large values of N or N max: N 24 and N 30, respectively. This can be
Table 4. Convergence of the Effective Index γ1 ∕k0 1.14623 of an Elliptical Dielectric Waveguide (Radii of ra 0.5λ and rb 0.45λ) with Respect to Both N and M max 2M 1a M max 3 5 7 9 11 13
N 11
N 13
N 15
N 17
1.169860902 1.147973728 1.147973728 1.148409721 1.148409721 1.148426112
1.169590596 1.146103877 1.146103877 1.146645363 1.146645363 1.146647215
1.169551905 1.145732438 1.145732438 1.146295794 1.146295794 1.146297765
1.169546524 1.145664811 1.145664811 1.146232876 1.146232877 1.146234887
Values were calculated with u1 0.5λ, u2 4r 1 , ε1 2.5, and ε2 1.
a
Fig. 5. Spatial distribution of the Ez [Fig. (5a)] and Hz [Fig. (5b)] components for γ 1 ∕k0 1.14623 in the case of an elliptical dielectric waveguide with radii of r a 0.5λ and r b 0.45λ.
674
J. Opt. Soc. Am. A / Vol. 31, No. 4 / April 2014
Edee et al.
C. General Case of a Differentiable Profile Function: Eigenfunctions of a Dielectric Waveguide with an Elliptical Cross Section In this subsection, we investigate the impact of the number of Fourier basis functions (M) on the convergence of the numerical results. This is performed within the framework of a dielectric waveguide with a cross section defined by an ellipse: b f θ p ; 1 − e2 cos2 θ
(29)
where the eccentricity e is linked to the axial radii r a and r b by the following relation:
e2 1 −
r 2b : r 2a
(30)
We assume that a_ does not vanish by considering the following numerical parameters: r a 0.5λ and r b 0.45λ. The electric and magnetic parameters remain unchanged. Since a_ does not vanish, the eigenfunctions are coupled along the θ direction through the boundary conditions at u u1 written with the components E 2′ and H 2′ [Eq. (12)] but also with the functions ∂3′ E2′ and ∂3′ H 2′ [Eq. (20)]. The computations of these components need the use of a finite number (2M 1) of _ Tables 4 and 5 show the convergence Fourier coefficients of a. with respect to both N and 2M 1 for two effective indices: γ 1 ∕k0 1.14623 and γ 2 ∕k0 1.24235, respectively. These results compared well with those obtained using the COMSOL software based on the finite element method. With this software we obtain γ 1 ∕k0 1.14622, γ 2 ∕k0 1.24235, and the computation area (circle of radius 50λ) is meshed with 17,632 elements. The results presented in Tables 4 and 5 show remarkable stability in terms of the convergence with respect to M and N. However, the rate of convergence of the eigenvalues differs according to the eigenmode that is observed. For γ 1 ∕k0 , the fifth digit (1.14623) stabilizes as soon as N reaches 17 and M 5, while values of N 15 and M 4 are necessary to reach the same accuracy for γ 2 ∕k0 . The field maps in Figs. 5(a), 5(b), 6(a), and 6(b) represent the spatial distribution of the E3 and H 3 components and help understand this finding.
6. CONCLUSION
Fig. 6. Spatial distribution of the Ez [Fig. (6a)] and Hz [Fig. (6b)] for γ 2 ∕k0 1.242352 in the case of an elliptical dielectric waveguide with radii of r a 0.5λ and r b 0.45λ.
We have presented a mode solver dedicated to waveguides with arbitrary cross sections. We used a change of coordinates to map the profile of the waveguide onto a circle and then introduced a Hertz potential to formulate the eigenvalue problem. We succeed in deriving a matrix eigenvalue equation by implementing the MMGE. The MMGE has once again shown its capability to take into account the continuity conditions of the electromagnetic field and allows us to build an operator incorporating these conditions. The method was successfully applied to calculate the eigenvectors and eigenvalues of cylindrical waveguides with arbitrary cross sections. Both the case of a hollow perfectly conducting waveguide and a dielectric waveguide were treated to demonstrate how the various numerical parameters are involved in the convergence of the results. It appears from these numerical studies that in the case of structures with circular cross sections, the number of polynomials describing the field components in the radial direction and the number of Fourier harmonics describing the azimuthal behavior of the field operate separately in the convergence of the electromagnetic field calculation. More specifically, once the derivative of the waveguide profile vanishes, a study of the convergence of the modes with respect to the number of Fourier harmonics is incongruous because in this case there is no coupling along the θ direction. However, when the derivative of the profile does not vanish, the convergence of the modes of the structure under consideration depends strongly on the differentiability of the profile for a given number of polynomials. In the present implementation, the numerical behavior of the MMGE, according to the differentiability of the profile of the waveguide cross section, is not very different from that of the traditional C method. The differentiability of the profile function determines the convergence of the Fourier series of the derivative function of the profile and thus the numerical scheme according to the azimuthal
Edee et al.
Vol. 31, No. 4 / April 2014 / J. Opt. Soc. Am. A
variable. At this stage, the real challenge is to extend the method to more general nondifferentiable profiles by introducing a parametric formulation for the azimuthal variable.
Consider the variables ζ belonging to −1; 1 and u ∈ u1 ; u2 defined as follows: uζ
APPENDIX A: SOME USEFUL OPERATORS IN THE INTERVAL −1;1 AND A CHANGE OF VARIABLES Gegenbauer polynomials are always defined in the interval −1; 1. Generally, the computation of all inner products presented in this paper needs a change of variables, and this directly affects the inner products. Before presenting expressions for the change of variables, we present the recursive relationships that allow these inner products to be calculated in the interval ζ ∈ −1; 1. For all ζ ∈ −1; 1, we get ζ
~ uu χ − iηu − u1 u1 ;
(A1)
~ uζ Aζ B;
(A6)
u2 − u1 ; A χ − iη 2
(A7)
and B
However, nC Λn ζ 2n Λ − 1ζC Λn−1 ζ − n 2Λ − 2C Λn−2 ζ;
u2 u1 2u1 2u1 − : u2 − u1 u2 − u1 χ − iηu2 − u1
u~ n 1 n − 1 2Λ Λ jC Λ i jC n−1 i: 2n Λ n1 2n Λ
d d ζ B ; du~ dζ
2 Λ Λ Λ Λ ~ Λn iuu hC Λm juC ~ 1 ;uu ~ 2 A hC m jζC n i BhC m jC n i > > > : Λ 2 Λ 3 Λ 2 Λ Λ Λ 2 Λ Λ hC m ju~ C n iuu ~ 1 ;uu ~ 2 A hC m jζ C n i 2BhC m jζC n i B hC m jC n i:
hC Λm jζC Λn i
A10
ACKNOWLEDGMENTS n 1 n − 1 2Λ Λ Λ hC Λ jC Λ i hC m jC n−1 i; 2n Λ m n1 2n Λ (A2)
This work was funded by grants from the French program investissement d’avenir managed by the National Research Agency (ANR), the European Commission (Auvergne FEDER funds), and the Region Auvergne in the framework of the LabEx IMobS3 (ANR-10-LABX-16-01).
n 1 n − 1 2Λ Λ hC Λ jζC Λn1 i hC m jζC Λn−1 i: 2n Λ m 2n Λ (A3)
REFERENCES
we get hC Λm jζ 2 C Λn i
(A9)
and
hD E E D Ei 8D Λ Λ dC Λn Λ ~ dC n Λ Λ dC n > j u A C jζ j B C C > m m m du~ uu dζ dζ > ~ 1 ;uu ~ 2
