Chapter 2 The Finite Element Method
2.3 The Hybrid Edge/Nodal Element
2.3 The Hybrid Edge/Nodal Element
Various types of finite element methods have been developed for the he existence of spurious modes.
Sp
ean solve all
electromagnetic problems, because spurious modes would arise in the olution of the vector wave equation if the wrong differential form is used to pproximate the electric-field vector. Early thinking about spurious modes attributed this problem to a deficiency in imposing the solenoidal nature of the field in the approximation process. A series of papers, beginning with Konrad [A. Konrad, 1976] and followed by [M. Hara, 1983] expounded this idea. Many researchers have been influenced by the notion that spurious modes are caused by the non-solenoidal nature of finite element approximation procedures. Yet, the early thinking is wrong. The true reason of spurious modes is the incorrect approximation of the null space of the curl operator [M. Hano, 1984]. It has been shown that the hybrid edge/nodal ector elements with triangular shape imposing the continuity of the tangential field but leave the normal component discontinuous are very useful for eliminating the spurious solutions.
In the method employed in this thesis, the word “hybrid,” means the order of interpolation functions is mixed. As for the term “edge/nodal,”, the former indicates a set of vector interpolation functions locate at the edges of the elements and are responsible for the transverse field interpolation, and similarly, the latter indicates another set of vector interpolation functions locate at the nodes of the elements and are responsible for the longitudinal field interpolation. Recently, curvilinear hybrid edge/nodal elements are introduced in the simulation of photonic crystal fibers [M. Koshiba et al.].
full-vectorial analysis of guided-wave problems. An important issue for full-vectorial finite element analysis is t
urious modes are numerical solutions of the vector wave equation without physical m ing. The scalar finite elements are not sufficient to
s a
v
The virtue of this kind element lies in the fact that they can match the curved boundary better than rectangular ones with more accuracy and fewer unknowns.
In our study, the CT/LN (constant tangential, linear normal and linear nodal) rectilinear element [M. Koshiba et al., 2000] is used for the
imulation work. Fig.1 shows the CT/LN element which is composed of an s
edge element with three tangential variables, φt1 to φt3, based on constant gential and linear normal vector interpolation functions, and a linear no
tan
dal element with three axial variables, φz1 to φz3 . The tangential component of a specific CT/LN interpolation function is constant along one edge of the triangle element and is zero along the other two edges, while the normal component is a linear function along the three edges.
For elements with 2D triangular shapes, the Cartesian coordinate, x and y, in each element can be approximated with the linear local coordinate functions Li( i =1,2,3), as shown in Fig.2.
3 each element. Note that the relation between the local coordinates is defined as L1 +L2 +L3=1. For 2D problem, L1 , L2 are usually selected as independent variables. The transformation for differentiation is given by
[ ]
Where [J] is the Jacobian matrix. The transformation relation for integration of a function f(x,y) is given by
1
where |J| is the determinat of the Jacobian matrix and is called Jacobian.
The following numerical integration can be applied directly to above integration
Where subscript i denotes the quantity associated with the sampling point i
Edge
Table 1. Interpolation functions
{ }
φt ix{U}+ {V} NodalTable 2. Values of weighting coefficients and local coordinates Wi L1i L2i L2i 2 0.13239415
a 3 a2 a 470142063 a3 =0. 3 0.13239415
a 3 a 3 a2 a4 =0.79742669 4 0.13239415
a4 a 5 a 101286515 a5 =0. 5 0.12593918
6 0.12593918 0.12593918
a 5 a4 a 5
7 a 5 a 5 a4
Fig. 1 The basis functions of the CT/LN element.
Local coordinate (L1,L L )
1
6
2 5 3
4
Point 2, 3
1 (1,0,0) 2 (0,1,0) 3 (0,0,1)
4 (0.5, 0.5, 0)
5 (0, 0.5, 0.5)
6 (0.5, 0, 0.5)
Fig. 2 The local coordinate of a element.
Chapter 3 The Perfectly Matched Layers
3.1 Introduction to PML
The perfectly matched layer (PML) is an artificial medium which serves as an absorbing boundary condition (ABC). This absorbing boundary condition holds great promise for truncating the mesh efficiently in the numerical computation for EM wave problems. It can absorb radiation wave almost without reflection at the absorber interface for arbitrary angle, wavelength and polarization of incidence. In addition,
linear lossy or anisotropic media eixeira, 1998]. Because of the it i ed as “perfe matched layer,”.
s originally proposed by Berenger for FDTD imulation [J. P. Berenger, 1994]. He used the so called “split field,”
ch, for example, Hz is decomposed into Hzx and Hzy. This leads to a modified version of Maxwell equations, where the introduction of the split fields provides extra degrees of freedom that can be used to achieve a perfect reflectionless match at the absorber interface. This is fairly a revolution, since it was quickly shown that the PML outperforms other previous known boundary conditions.
In the subsequent years, there are many d retations about the physical aning of PML . Chew et al., indicated that the Berenger’s PML
an be derived from a more general way base on the concept of complex
coordinate s can be
regarded as a regular isotropic medium complex thickness. Sacks et al., [Z. Sacks, 1995] revealed that th also can be considered as an anisotropic medium. In fact, Chew’s a d Sacks’ statements are equivalent mathematically [F. L. Teixeira, 1998].
Here, we adopt anisotropic PML n our simulation because in finite element analysis, it is convenient to specify the material parameters (i.e.
it is also valid for any [F. L. T
ctly reflectionless, s term
The idea of PML wa s
formalism, in whi
ifferent interp me
c
tretching [W. C. Chew et al., 1997], in which the PML but with a
e PML n i
permittivity & permeability) of the PM
3.2 The Concept of PML
for a plane wave
coordinate transformation to z L.
The concept of PML is shown as follow
0nz
e
-jkOne can apply the following
z d ) z ( s z~
z
0
z ′ ′
=
→
∫
where sz (z′) is so called complex stretched variable [W. C. Chew et al., 1997].
Considering the following structure
z
We can let sz(z′)=1 in non-PML region, and sz(z′)= 1- j*c in PML region. In consequence, z~ = z (real) in non-PML region, and z~ = z - j*c*z (complex) in PML region.
z
This implies that the propagation wave would be absorbed in the PML region. Base on this idea, after a series of transformations, the permittivity and permeability tensors in the PML can be expressed as [F. L. Teixe 1998]:
ira,
] [
)
(detS 1 S S
PML = − ⋅ε⋅
ε
µPML =(detS)−1[S⋅µ⋅S]
Note that the intrinsic impedance of the PML region is equal to that of the
non-PML region, namely, impedance match. So no reflection occurs at the absorber interface.
1) The specification of are shown in Fig. 3, where
nuation of the field in PML regions can be controlled by choosing th
ρ is stance from the PML interface, and t is the thickness of PML.
The atte
e value of αi appropriately.
Fig. 3 The specification of s in the PML.
PML Region PML
arameter
1 2 3 4 5 6 7 8
P
s x s1 s2 1 1 s1 s2 s1 s2
s y 1 1
1 1 1 1 1 1 1 1
s 3 s4 s 3 s 3 s4 s4
sz
t4 t2
t3
4
2
8 3
5 6
1
7
t1
3.3 Solution Issues
In order to reliably predict the features of HFs, the full vector formulation is applied. Because of the losses resulted from the air holes and the finite transverse extent of the confining structure, the effective index is a complex value and its imaginary part is related to the losses. This is termed as a leaky mode. In order to evaluate propagation losses of leaky modes, an anisotropic perfectly matched layer (PML) is employed as an absorption boundary condition at the computational window edges.
How to solve the large sparse generalized eigen value problem to get the leaky mode solution is an important issue in finite element procedure.
Recently, Koshiba et al., employed the imaginary-distance beam propagation method (ID-BPM) based on the FEM to deal with the leaky-mode problems [K. Saitoh et al., 2002]. On the other hand, Selleri et al., used the Arnoldi method based complex modal solver t tain sev eigen- des directly [S. Selle t al., 2001]. The ID-B starts propagation with an initial approximate field and the field evolves into the exact eigen-mode after propagating a long distance. The drawback of this method is that it’s time consuming to achieve high accuracy and only one mode can be obtained each time. On the contrary, solving the eigen-value quation directly is much simpler than the ID-BPM method. For this reason, we use a modified Matlab built-in eigen-solver based on Arnoldi algorithm
eat this problem
propag n los f le ode is defined as
Fig. 4 shows the cross section of the HF surrounded by the PML regions
with thickness se plane, z is the
ropagation direction, and and are the half computational window
=
d . Here x and y are the axes of the transverρ
Wx Wy p
size respectively. The PML parameter S is complex for the leaky mode as:
analysis, which is given
α j
-=1 S
The parameter s controls the attenuation of the field in the PML region through the choice of the appropriate value of α with the parabolic profile
α =
2
max ⎟⎟
⎠
⎞
⎜⎜
⎝
⎛ dp
α ρ ;
whereρis the distance in the PML region from its inner interface.
Consider a PCF with two ring air holes with hexagonal (or triangular ) lattice arrangement, as shown in Fig. 4, the hole pitch Λ =2.3 mµ , silica index = 1.45, air filling ratio =0.5
Λ
d , dρ =2.3µm, 7wx =wy = µm, and operating wavelength λ =1.5µm. Because of the six-fold symmetry of this PCF, for the fundamental mode, one-fourth of the fiber cross section combined with proper boundary conditions is taken into computational region. Fig.5 to Fig.8 are the E field distributions for x component, y component, z component and transverse component respectively. Fig.9 and Fig.10 are the effective indices and propagation losses compared with the published results respectively [K. Saitoh et al., 2002].
Fig. 4 The computational window of the PCF.
d
ρw
xy
x
w
yd
Λ
0 2 4 6 8 10 0
1 2 3 4 5 6 7 8
Ex
Wavelength :1.5 um a/d=0.5
a
d
um um
Fig. 5 The Ex field.
0 2 4 6 8 10
0 1 2 3 4 5 6 7 8
Ey
Wavelength :1.5 um a/d=0.5
um um
Fig. 6 TheEy field.
0 2 4 6 8 10 0
1 2 3 4 5 6 7 8
Ez
Wavelength :1.5 um a/d=0.5
um
Fig. 7 TheEz field.
0 2 4 6 8 10
0 1 2 3 4 5 6 7 8
Et
Wavelength :1.5 um a/d=0.5
um um
Fig. 8 The transverse field.
(a)
(b)
0 0.5 1 1.5 2
1.4 1.41 1.42 1.43 1.44 1.45
Two-ring air holes
Wavelength um
Effective index
0.5 0.6
d/a=0.5
d/a=0.6 d : Hole diameter
a : Lattice constant
Fig. 9 The effective index of two-ring air holes.
(a) Published result in [K. Saitoh et al., 2002].
(b) Our simulation result.
(a)
(b)
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
10-8 10-6 10-4 10-2 100 102 104
P ro p a g a ti o n L o s s fo r Tw o -ri n g a i r h o l e
Wa ve l e n g th u m
Loss
d B /m
d /a = 0 .5
Fig. 10
. Saitoh et al., 2002].
(b) Our simulation result.
The propagation loss of two-ring air holes.
(a) Published result in [K
Chapter 4 Optical Properties of PCFs
4.1 Endless Single Mode
The picture in fig. 11 shows the structure of a typical holey fiber. In [T. A.
Brinks, 1997], the fiber was reported to be single mode over a remarkably wide wavelength range, form 458 to 1550 nm at least, and it is confirmed numerically that HFs are endless single mode for d/Λ<=0.43 [M. Koshiba, 2002]. The cladding effective index, which is a very important design parameter for realizing a single-mode HF, is defined as the effective index of the infinite photonic crystal cla ing if the core is absent. Unlike
nal fibers, the effective index of the HF cladding is very sensitive to the wavelength, and has to be estimated for different frequencies by finding the fundamental space filling mode (FSM) of a cladding unit cell, which is shown in Fig. 12.
Just like conventional fibers, the effective index of the HF guided mode is between the core index and the effective index of the cladding. Fig. 13 (b) shows our simulation results for the fundamental mode of the PCF, and Fig.
13(a) is the published results with the same structure parameters as in (b).
Fig.14 (a) shows the published result of the effective index for the fundamental guided mode and cladding with pitch = 2.3
dd conventio
µm, hole diameter
= 0.6 mµ . Fig.14 (b) shows the simulation result. The lowest solid line is the fective index, the middle is the fundamental mode effective index, and the user is the core index. We can find that this structure support single mode propagation because only one effective index is found in the region between the core index and the cladding effective index. Fig.15 shows the dispers on curve for different air hole diameters with a fixed pitch of 2.3 cladding ef
i µm.
One can see that, as the hole diameter increasing, the dispersion zero poin
shifts to the near part is
approximately preser
t higher frequency range, and the slope of their li
ved.
Here, the group velocity dispersion is defined as ] dλ Re[n
d D=−cλ⋅ 22 ⋅
eff
Fig. 11 The endless single mode fiber.
Unit cell for FSM Unit cell for FSM Unit cell for FSM
Fig. 12 Unit cell for FSM estimation.
(a)
(b)
Fig. 13 Contour plot of the fundamental mode in the PCF (a) Published results [M. Koshiba, 2002]
(b) Our simulation results )
(a
(b)
2 3 4 5 6 7 8 9 10
1.44 1.442 1.444 1.446 1.448 1.45 1.452 1.454 1.456 1.458 1.46
neff of cladding
neff of guided mode
core index
Neff
pitch=2.3 um
pitch/wavelength
Fig. 1 e PCF
( hiba, 2002]
(b) Our simulation result
4 Effective index of the guided mode of th a) Published result [M. Kos
0.5 1 1.5 2 2.5 3 -120
-100 -80 -40
-60 -20 0 20 40
a=0.2um a=0.25um a=0.3um a=0.35um a=0.4um a=0.45um Dg
ps/nm/km
wavelength (um) geometrical dispersion curves
pitch=2.3 um
Fig 15. Our simulation of dispersion curves for a fixed pitch with different air filling ratio a
4.2 Mode Classification
When higher order modes or polariza on properties are considered, the full vector approach is crucial for assessing the true behavior of electromagnetic waves in complex wave g structures such as PCF. It is well known that the PCF is often strongly multimode in the visible and near-infrared regions when the filling factor is large enough. This may lead to a number of intermodally phase-matched nonlinear processes. As a consequence, it is necessary to investigate the modal properties of PCFs including their degeneracy, classification, and so on.
Some studies have discussed the mode properties of PCFs [M. J. Steel et al., 2001], and the classification and degeneracy properties of higher-order modes were discussed further in [R. Guobin et al., 2003]. It is shown that the mode classification of a PCF is similar to that of a step-index fiber, except for modes with the same symmetry as the PCF. When the doublet of the degenerate pairs both have the same symmetry as the PCF, they will be split into two non-degenerate modes.
In the scalar approximation, the characteristics of polarization in the PCF are hidden. With the full vector aroach, the vector modal behavior of the
PCF ca e
symmetry of a l important characteristics of the odes of the waveguide [T. P. White et al., 2002]. Determining the
mmetry type of a particular waveguide enables one to classify the possible modes in terms of mode classes, a ode degeneracies between mode classes. Further, from the azimuthal symmetries of modal electromagnetic fields of a mode class inimum sector of waveguide cross section, together with associated boundary conditions for this sector, which is necessary and su icient for completely determining all the modes of that mode class.
If a uniform wave-guide has n metry and also possesses precisely n reflection planes spaced azimuthally by
ti -guidin
n be predicted. From the theory of group representations, th waveguide controls severa
m sy
nd to predict the m , one can specify a m ff
-fold sym n
π radians, it is said to have symmetry. For the triangular lattice PCF, it has six-fold rotation
ymmetry and Cnv
6
π reflection symmetry, so the point group is C6v. PCFs s
with the C6v symmetry possess the following properties: all the modes can e divid o eight classes according to the minimum sector and boundary co
of these modes. Fig. 24 to
b ed int
nditions. As shown in Fig. 16, class p = 1,2,7,8 are non-degenerate, which exhibit full waveguide symmetry, i.e. C6v symmetry, while class p=3, 4 and p=5, 6 are degenerate pairs, and they exhibit full waveguide symmetry in combination.
Table 3 lists the first 14 modes of PCF (with pitch =2.3 um, air filling ratio =0.8, and wavelength = 0.633 um) with effective index, mode class, degeneracy, computation error, and label published in [R. Guobin et al., 2003]. Fig. 17 to Fig. 23 shows the field distribution
Fig. 28 show the fundamental mode and some higher order modes simulated respectively by using the minimum sector with specific boundary conditions belonging to the mode class listed in Table 3
Fi
ble 3 Mode classes of the triangular-lattice PCF [R. Guobin et al., 2003].
g. 16 Minimum sectors for waveguides withC6vsymmetry, the modes of waveguides are classified into eight classes (p=1, 2, 3, …8). Solid lines indicate PEC boundary condition, and dashed lines indicate PMC boundary condition [R. Guobin et al., 2003].
Ta
(a) (b)
Fig. 17 mode [R. Guobin et al., 2003].
(a) (b) HE11
Fig. 18 (a) TE01 mode (b) HE21 mode [R. Guobin, 2003].
(a) (b)
Fig. 19 (a) HE21 mode (b) TM mode [R. Guobin, 01 2003].
(a) (b)
Fig. 20 (a) HE311 mode (b) EH11 mode [R. Guobin, 2003].
(a) (b)
Fig. 21 (a) EH11 mode (b) HE312 mode [R. Guobin, 2003].
(a) (b)
Fig. 22 (a) HE12 mode (b) HE12 mode [R. Guobin, 2003].
(a) (b)
Fig. 23 (a) EH21 mode (b) EH21 mode [R. Guobin, 2003].
−0.5 0 0.5 1 1.5 0
0.5 1 1.5
Fig. 24 mode.
2
HE11
−0.5 0 0.5 1 1.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Fig. 25 TE01 mode.
−1 −0.5 0 0.5 1 1.5 0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Fig. 26 HE21 mode.
−0.5 0 0.5 1 1.5
0 0.5 1 1.5 2
Fig. 27 HE21 mode.
−0.5 0 0.5 1 1.5 0
0.5 1 1.5 2
Fig. 28 EH21 mode
4.3 Dispersion Flattening
Controlling the chromatic dispersion in optical fibers is a very important problem for communication systems. In both linear and nonlinear regimes, or for any optical systems using ultra short soliton-pulse propagation. In all cases, the almost-flattened fiber-dispersion behavior becomes a crucial issue.
The most appealing feature of photonic crystal fibers is their high flexibility based on the particular geometry of their refractive index distribution. This fact allows us to manipulate the geometrical parameters of the fiber to generate enormous variety of different configurations.
The form of dispersion relation of guided mode for PCFs is very sensitive to the 2D photonic crystal cladding. For this reason, one can expect to control, at least to some extent, the dispersion properties of guided modes by manipulating the geometry of the cladding. It was soon realized that the PCF exhibited dispersion properties very different than thos of conventional fibers. As an som configurations presenting
ero dispersion point bellow that of slica at 1.3
e example, e PCF
µm
z [D. Mogilevtsev, 1998;
. J Bennet, 1999], and some others showing flattened dispersion (one point f zero third-order dispersion) or near zero ultra-flattened dispersion (one oint of zero fourth-order dispersion) profiles [A. Ferrando et al., 2000].
ince the number of different photonic crystal configurations is significant, ne can deduce that it must be possible to elaborate a procedure to tailor the ispersion of the PCF modes in an efficient way. A systematic approach to esign the dispersion properties of the PCF using a systematic procedure has een already suggested in [A. Ferrando et al., 2001]. The analogous design etails of dispersion flattened or dispersion compensation for triangular
ttice PCFs and high-index-core bragg fibers were also proposed in [A.
errando et al., 2000], and the dispersion properties of square lattice PCF ere discussed recently [A. H. Bouk 2004].
Fig. 29 and Fig. 30 show two sets of PCF geometrical dispersion curves:
. Different pitches with a fixed air filling ratio.
. Different air filling ratios with a fixed pitch.
set A, for different pitches, the dispersion zero point would shift and the lope of each curve in the linear region is different. In set B, for different air P
filling factors, the curves are shifted and the slope of their linear part is pproximately preserved.
D ~= Dg - (-Dm)
set to be pa
a
Fig. 31 shows the idea of dispersion flattening design. The line with open circle is the opposite sign of the material dispersion curve, the line with open diamond is the geometrical dispersion curve, and the triangle is the total dispersion, which is defined as
The key factor to achieve the flattened dispersion curve is the control of the slope of the linear part Dg. The sign changed material dispersion -Dm is a smooth and almost linear curve in most of the infrared region. It is clear from Fig. 31 that if the linear part of geometrical dispersion can be
rallel with the material dispersion, therefore the total dispersion will achieve an ideal perfect flattened behavior.
The strategy to obtain such a behavior is then straight forward. It can be started by determining the slope of material dispersion curve at some specific wavelength. In the region where the material dispersion curve is smooth, the slope is approximately the same for a reasonably wide neighborhood around the specific wavelength. Once the slope of the Dm is fixed, we can change the pitch and the air hole diameter with a fixed air filling factor to obtain a Dg curve having a linear region with the same given slope of Dm, as shown in Fig. 29. If both the linear part of material and geometrical dispersion curves overlap in a specific wavelength region, then the dispersion flattened curve can be obtained in the overlapping region.
After that, we can fix the pitch and change the air hole diameter to get a shifted curve, as shown in Fig. 30. Then different widths of flattened region can be obtained.
6
Fig. 29 Dispersion curves with fixed air filling ratio.
1
Fig. 30 Dispersion curves with fixed pitch.