Photonic crystal fibers ( or holey fibers) consisting of a central defect region surrounded by multiple air holes running parallel to the fiber axis have attracted a lot of research interest in recent years. Due to the array-like arrangement of the air-hole cladding, the holey structure of PCFs can provide more design flexibility than conventional fibers. It can be tailored for wide-band single-mode transmission, for high nonlinearity with small core area, or for zero or flattened dispersion in the optical communication window. Because of the large index difference between the cladding and the core, the scalar approximation for weakly guiding is not applicable and the full vector formalism is needed. Further more, owing to the curved structure of air holes, a numerical algorithm with high accuracy is also needed.
The Finite Element Method (FEM) has become the trend of numerical simulation for studying the mode properties and propagation characteristics of waveguides with arbitrary cross section shapes. In contrast, for the Finite Difference Method (FDM), since the element mesh is rectangular, for structures with curved shapes one has to increase the number of grid points in order to achieve the specified accuracy. If the structure is large or complicated, this will result in much computational efforts. On the other hand, the finite element method permits users to choose the shape of mesh (e.g. triangular or curvilinear) and the order of interpolation functions according to the requirements. When the FEM is adopted, the triangular meshes can be utilized to match the curved boundary better than the rectangular meshes utilized in finite difference. Besides, one can refine the mesh in a specific region rather than in the whole region. The user also can impose high order interpolation functions to reach fast convergence of the solution. All of the above procedures of the finite element method can be introduced to attain higher accuracy with fewer unknowns. These are what FDM cannot easily achieve.
In this thesis, the CT/LN edge element is applied in the simulation. The mathematical formulation will be described first. Then, the dispersion
property, leak be studied
and some sim results.
c Crystal Fibers
radial la
fiber, also named as holey fiber, is similar to th
age loss, and the birefringence of holey fibers will ulation works will be compared with the published
1.2 Introduction to The Photoni
Photonic crystal fibers (PCFs) have in recent years attracted much scientific research and technological development interest. Generally speaking, PCFs may be defined as the optical fibers in which the core and/or the cladding regions consist of micro structured air holes rather than homogeneous materials. The most common type of PCFs, which were first fabricated in 1996 [J. C. Knigh, 1996], consists of a pure silica fiber with an array of air holes extending along the longitudinal axis. Later on, the PCFs fabricated from other host materials [K. M. Kaing, 2002] or with incorporated sections of doped materials have been demonstrated. A considerable amount of modeling and experimental efforts have also been put into the design and fabrication of circularly symmetric PCFs with
yers of alternating index contrasts [G. Ouyan , S. G. Johnson].
Traditional optical fibers are limited to rather small refractive index difference between the core and the cladding (about 1.48:1.46). For photonic crystal fibers, this refractive index differences is significantly larger (1.00:1.46), and can be tailored to suit particular applications. It is this flexibility combined with the ability to vary the fiber geometry that enables the enhanced performance of the photonic crystal fibers. The PCFs can be designed to satisfy many specific purposes. For example, they can be single- mode over an extremely broad wavelength range, can support larger or smaller mode field diameters, can meet specific dispersion requirements, can increase or decrease nonlinearity, and can be highly birefringent for achieving improved polarization control.
For conventional optical fibers, the electromagnetic modes are guided by total internal reflection in the core region where the refractive index is raised by doping the base material. In PCFs, two distinct guiding mechanisms are possible: index guiding and band gap guiding. The guiding mechanism of index guiding
at of conventional fiber. It features a solid core surrounded by a regular array of microscopic holes extending along the fiber length. The solid core
can be viewed as a defect within the surrounding periodic structure formed by the regular array of air holes. The holey structure acts as the cladding to onfine the fundamental mode within the core of fiber, while allowing the t arises because e periodic holey structure creates an effective index difference between th
ved control ov
into the desired pattern, fusing the
c
higher order modes to leak out of the core. The confinemen th
e core and the surrounding material. On the other hand, the band gap guiding fibers, also termed as hollow core fibers, are constructed with a hollow core surrounded by a periodic structure of air-holes. The periodic structure generates a photonic band gap. When the light frequency is located within the band gap, the light can be confined in the core region and propagates along the fiber. In this study, the characteristics of holey fiber will be addressed.
Over the last seven years the PCFs have rapidly evolved from scientific curiosity to commercial products manufacturing and are sold by several companies. A central issue of PCFs from the early days to the present has been the reduction of optical loss, which initially was several hundred dB/km even for the simplest PCF design. Through the impro
er the homogeneity of the fiber structures and the application of highly purified silica as the base material, the loss has been brought down to a level of a few dB/km for the most important types of PCFs. The current world record is 0.37 dB/km [K. Tajima, 2003]. Thus, with respect to the optical loss, PCFs have undergone an evolution similar to that of standard fibers in the 1970s. Their application potentials have also increased accordingly. For some types of PCFs, the loss figures are still substantial and more work is definitely required. However, for many applications the optical loss has ceased to be a decisive barrier to the practical application of PCFs.
PCFs are most commonly fabricated by hand-stacking an array of doped or undoped silica capillary tubes or solid rods
stack into a preform , and then pulling the perform to a fiber at a temperature sufficiently low (~1900’C) to avoid the collapse of the holes.
The vast improvements of the fabrication process made in recent years have not only served to bring down the optical loss, but have also greatly increased the diversity of the fiber structures available to the designer.
Consequently, new PCF designs appear continuously, and it will probably take a few more years before the field can be said to have matured.
1.3 The Need of Full Vector Formulation
The scalar approximation of the wave equation is adequate for weakly guiding problems. But when the index difference between the cladding and core is large enough such as in PCFs, the x, y, z field components are no longer independent, and will be coupled together. Therefore, the scalar approximation is no longer valid. This idea can be explain clearly with the following mathematical description.
0
In the above vector wave equation, n=n(x, y) is the transverse dielectric profile, and k is the wave vector in free space. The double curl Ev
can be written as:
E
field of the j-th eigen mode of wave guide is express as:
z) By applying eq. (1.2) and (1.3) to eq. (1.1), the following three coupled full-vector equations can be obtained as follow
(1.6)
For a weakly guiding wave-guide, the index difference between the core
e.
and cladding is very small and hence the right hand side of eq. (1.4) to eq.
(1.6) can be neglected. These become the well known scalar Helmholtz equations. The filed components in x, y, z directions are all independent in this case.
For wave-guide with large index difference, such as photonic crystal fibers, the coupling among the three polarization field components through the boundaries should be taken into consideration. Therefore, the full-vector wave equation is demanded for calculating precise modal fields and propagation constants. Note that there is no TE or TM mode in this cas
Chapter 2 The Finite Element Method
2.1 The Finite Element Procedure
There are two approaches of finit on. One is the va
f the governing equation, and the solution corresponding to the equation should be the one which makes the varia functional to be
ze n the other
hand, the Galerkin’s method needs a set of test functions to perform the projection. For more details, the readers may refer to [J. F. Lee
The vector wave equation can be written as
e element formulati
riationl method, and the other is the Galerkin’s method (or the weighted residual method). For the variational method, one should first determine the functional o
tion of the ro. In this thesis, the variational method will be introduced. O
, 2002].
2.2 The Variational Method
(
[p][s])
ko2[q][s] 01
The functional of the vector wave equation eq. (2.1) is given as
*
When the whole area is divided into elements, F can be expressed as the summation of the integration over each element
o
the {U} and {V} are the vector edge interpolation functions, and {N} is the nodal vector interpolation functions listed in Table 1. The
) dal variable for each element respectively.
The variatio the edge and no
nδ of the functional F is given as F
∫
whereΓ is the outward boundary of the regionΩ , n is the outward unit normal vector. When φ is the solution ofδ =0, the following relations are F satisfied:
eq. (2.4) is the vector wave equation. This proves the solution of δ =0 is F also the solution of vector wave equation.
By applying eq. (2.3) to the functional F , taking the first variation of F, and settingδ =0. The following matrix equation can be obtained: F
Here
dxdy
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.
The propagation loss of two-ring air holes.