by
Mohammad Amin Cheraghi Shirazi B.Sc., University of Tehran, 2011
A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of
MASTER OF APPLIED SCIENCE
in the Department of Electrical and Computer Engineering
© Mohammad Amin Cheraghi Shirazi, 2015 University of Victoria
All rights reserved. This thesis may not be reproduced in whole or in part, by photocopying or other means, without the permission of the author.
Beam Propagation Modelling of Whispering Gallery Microcavities
by
Mohammad Amin Cheraghi Shirazi B.Sc., University of Tehran, 2011
Supervisory Committee
Dr. Tao Lu, Supervisor
(Department of Electrical and Computer Engineering)
Dr. T. Aaron Gulliver, Departmental Member
Supervisory Committee
Dr. Tao Lu, Supervisor
(Department of Electrical and Computer Engineering)
Dr. T. Aaron Gulliver, Departmental Member
(Department of Electrical and Computer Engineering)
ABSTRACT
Whispering Gallery Mode (WGM) microcavities have a wide range of applications from fundamental physics researches to engineering applications due to their ultra high quality factor (Q). For example, an ultra-high Q WGM cavity can be used as an bio/nanosensor since a nano particle bound to the surface of the cavity will result in a resonance wavelength shift. In the last decade lots of research have been conducted on this topic, as a result, WGM biosensors are emerging as one of the mainstream senors.
This thesis presents an efficient beam propagation method (BPM) simulation tool to study the light propagation behaviour in WGM cavities. Using this tool, the perturbation of the cavity properties caused by a polystyrene nano bead attached to the surface of a WGM silica microcavity is investigated. Furthermore, we numerically verify a three times sensitivity enhancement by fabricating a nanohole at the surface of the WGM cavity sensor.
In addition, we study the open cavity structures, cavity-waveguide coupling, huge WGM cavities, and deformed microcavities radiation. Finally, the impact of fabri-cation inaccuracy on asymmetric WGM cavities is investigated in terms of quality factor degradation.
Contents
Supervisory Committee ii Abstract iii Table of Contents iv List of Figures vi Acknowledgements x Dedication xi 1 Introduction 1 1.1 General Introduction . . . 1 1.2 Thesis Outline . . . 3 2 Background 4 2.1 Quality factor . . . 4 2.2 Waveguide-Cavity Coupling . . . 6 2.3 Buildup Factor . . . 82.4 Finite Difference Scheme . . . 9
2.4.1 Uniform Finite Difference Method . . . 10
2.4.2 Non Uniform Finite Difference Method . . . 11
3 Finite Difference Mode Solver 13 3.1 Cartesian Mode Solver . . . 15
3.2 Cylindrical Mode Solver . . . 17
4 Finite Difference Beam Propagation Method 21 4.1 Two Dimensional Cartesian FD-BPM . . . 21
4.1.1 Formulation . . . 22
4.1.2 Transparent Boundary Condition . . . 26
4.2 Cylindrical FD-BPM . . . 29
4.2.1 2D FD-BPM . . . 30
4.2.2 Perfectly Matched Layer . . . 31
4.2.3 3D FD-BPM . . . 34
4.2.4 Characterization . . . 35
5 Applications 48 5.1 Nano Particle Sensing . . . 48
5.2 Enhanced Sensor Structure . . . 51
5.3 Waveguide-Cavity Coupling . . . 57
5.4 Huge Cavity Modelling . . . 60
5.5 Asymmetric Cavity . . . 61
6 Summary 64
List of Figures
Figure 2.1 An illustration of taper-cavity coupling. The dashed red box shows the coupling area as a 4-port transfer function. The 4 ports are the taper and cavity inputs and outputs. . . 6 Figure 2.2 The taper output transmission vs. the normalized coupling
pa-rameter K. . . 8 Figure 2.3 One dimensional finite difference scheme. . . 11 Figure 3.1 Illustration of the cartesian and cylindrical coordinate systems. 13 Figure 3.2 The first (red) and third (blue) order mode of a 8-µm silica slab
waveguide in the air (left axis). The dotted black line indicates the refractive index profile of the waveguide on the right axis. . 17 Figure 3.3 The first (blue), second (red), and third (green) order whispering
gallery modes of a silica ring resonator in the air (left axis). The major radius is 45µm and the width of the waveguide is 10µm. The dotted black line shows the refractive index profile of the waveguide on the right axis. . . 19 Figure 3.4 Four different whispering gallery modes of a silica micro toroid
in water in the central wavelength of 970nm. The one on the top left is the fundamental mode. The major radius of the toroid is 45µm and the minor diameter is 10µm. The black lines show the boundaries of the cavities. . . 20 Figure 4.1 Field intensity inside a 4-mm long, 8-µm core diameter bent fiber. 29 Figure 4.2 Reflectivity in dB vs. different σ0 values. σ0 = 2.5× 1016 is the
optimum value. . . 32 Figure 4.3 The field intensity of a straight waveguide, launching toward the
computation window edge, for three different σ0 values: 0 (top left) , 2.5× 1016 (bottom), and 1020 (top right). The white lines indicate the inner 2-µm PML edges. . . . 33
Figure 4.4 Discretization scheme in r-z plain. . . 34 Figure 4.5 Field intensity of a 45µm major radius and 10µm width silica
ring, in logarithm scale. The radiations is observable. The solid white lines are illustrating the ring borders and the dashed black lines are showing the PML layer boundaries. . . 36 Figure 4.6 Intensity distribution of a huge microdisk with the radius of 1mm. 37 Figure 4.7 (a) Quality factor vs. grid spacings in the r and ϕ direction. (b)
cross section view of (a) at ∆r = 3.1 nm and it’s relative error. (c) variations for the cross section at ∆ϕ = 0.0015 radians with the corresponding relative error. . . 39 Figure 4.8 Resonance wavelength of the ring resonator and its relevant error
vs. grid size in the radial direction. . . 40 Figure 4.9 The gaussian field pattern (red) after (a) 0, (b) 25, (c) 125 round
of propagation inside a 45-µm major radius and 10-µm channel width silica ring resonator. The dashed blue curve is the funda-mental mode profile. . . 42 Figure 4.10Magnitude of the overlap factor between the output profile and
mode profile vs. round trip number. The red curve illustrates the overlap factor departure from unity. . . 43 Figure 4.11Power at the output of a 45µm silica ring resonator after each
round of propagation normalized to the first round power for resonance (red) and for λres(1 + 2Q1 ) (green). . . 44 Figure 4.12Buildup factor and its relative error vs. grid spacing in the radial
direction. . . 45 Figure 4.13Quality factor of a micro toroid and its relative error vs. grid
spacing in the ϕ direction. R is the outer radius of the toroid (R = 50µm). . . . 46 Figure 4.14Resonance wavelength (a) and quality factor (b) of a micro toroid
Figure 5.1 r-ϕ plain cross section view of micro toroid electrical field pro-file (a) and field intensity in logarithm scale (b), perturbed by a 400-nm radius polystyrene bead. The black lines depict the cavity outer edge and the bead boundaries. (c) The resonance wavelength shift vs. the bead radius using 2D BPM (green), 3D BPM (red), and perturbation method (black). . . 49 Figure 5.2 Resonance wavelength shift arising from a 400-nm radius PS
bead and its relative error vs. grid spacing in the ϕ direction. . 50 Figure 5.3 Resonance wavelength shift arising from a 400-nm radius PS
bead and its relative error vs. grid spacing in the r and z direction. 51 Figure 5.4 (a) Field intensity (in logarithm scale) of a 1.15 mm radius disk
at r-ϕ cross section. The black line is the cavity boundaries. The
r-z cross section view of the field intensity in logarithm scale (b)
after the notch (ϕ = 5◦), (c) inside the notch (ϕ = 3◦), and (d) before the notch (ϕ = 1◦). . . 52 Figure 5.5 (a) Wavelength shift vs. the hole depth for a 150-nm radius
hole-bead. The dashed red curve is the cavity mode pattern. (b) Enhancement vs. hole-bead radius (blue) for a 475-nm depth hole. Black line is the ultimate enhancement for this cavity (i.e. 2.23). Red curve is the quality factor vs. the hole-bead radius (right axis) and the green line is the bare cavity Q. (c) Field intensity in a logarithm scale for a silica microtoroid with a
475-nm depth, 150-475-nm radius hole on it. The same radius PS bead
is placed inside the hole. . . 54 Figure 5.6 The ratio of resonance wavelength shift (∆λ) and cavity linewidth
(δλ) vs. the ratio of the hole depth to the bead-hole radius. . . 56 Figure 5.7 The field intensity in logarithm scale for (a) 300 nm (b) 625 nm
Figure 5.8 (a) The 1.2-µm diameter taper output transmission vs. the microsphere-taper gap, modeled by the FD-BPM (red) and mea-sured experimentally (black). (b) The normalized coupling pa-rameter K vs. the gap. The FD-BPM and experimental results are presented in black and red curves respectively. (c) The r-ϕ plain cross section view of the field intensity plot in logarithm scale. The black lines indicates the borders of the 67-µm diam-eter microsphere and 1.2-µm diamdiam-eter taper. The experimental results are extracted from [1]. . . 59 Figure 5.9 r-z plain view of the field intensity in a (a) 90-µm and a (b)
1-mm diameter toroid. The black line is the resonator edge. r-ϕ plain view of the field intensity in logarithm scale in a (c) 90-µm and a (d) 1-mm diameter toroid. (c) The purple line is showing the resonator outer edge and the black lines are showing the
2-µm PML. (d) The black line is showing the resonator outer edge
and the white lines are showing the 2-µm PML. . . . 61 Figure 5.10The field intensities at the second round of propagation for (a)
R0◦ = 45 µm and R90◦ = 37.5 µm as well as (b) R0◦ = 45 µm and R90◦ = 32.5 µm. (c) Quality factor vs. ellipticity. . . . 63
ACKNOWLEDGEMENTS
I would like to express my gratitude to my supervisor Dr. Tao Lu for his insightful guidance and helpful suggestions through my graduate studies. I would also like to thank my colleagues Niloofar Sadeghi, WenYan Yu, Serge Vincent, Xuan Du, and LeYuan Pan for their helpful discussion and assistance.
DEDICATION
To my family and Separ
Introduction
1.1
General Introduction
Trapping the light in a tiny piece of glass sounds unachievable and inapplicable but is feasible thanks to the advanced micro fabrication technologies. The power build up of the cavity allows us to have a very narrow band and powerful beam of light. In contrast to the traditional Fabry-Perot cavity [2], Whispering Gallery Mode (WGM) cavities trap the light based on total internal reflection (TIR). The WGMs, such as dielectric micro sphere, toroid, disk, or ring are featured for their smooth surface and low material absorption to light. Consequently, the light, when circulating along their surface, incurs little loss. When in resonance, the phase of the circulating light matches after each round trip. Therefore, the constructive interference of light wave from each round collectively establish a strong electro-magnetic field at the cavity surface whose intensity is many orders of magnitude higher than that from a single round. In addition, a particle binding on the cavity surface can be sensed because it will be exposed to a portion of circulating light. The fact that some amount of light propagate outside the cavity is a unique property of WGM cavities comparing with other optical resonators. That portion of light which pass through the particle will experience a longer optical path and this will cause a resonance wavelength shift. Monitoring this shift, one is able to sense nano particles as small as 12.5-nm in radius [3].
WGM microcavities attained the researchers’ attention on a variety of subjects ranging from biosensing, nonlinear optics, and laser physics, to fundamental physics such as cavity quantum electrodynamics [3–17]. As opposed to the rapid
experi-mental advances, numerical studies as in [18–23] have not progressed with the same pace. There are different methods to analyze light propagation in waveguides. These include the boundary element (BEM) [22–27], the finite element (FEM) [28], the finite-difference time-domain (FD-TD) [29], the mode matching [21, 30–34], and the free space radiation mode methods [35]. Some of first-principle methods such as FD-TD require that the entire computational domain be gridded and are computation-ally inefficient and inaccurate. The inefficiency will be even magnified in modelling multi-scale problems eg. when a nano scale particle is on the micron-scale cavity. Alternatively, one may use the eigenmode expansion techniques to improve the accu-racy and efficiency. However, such method still requires significant computation time when multiple eigen modes are involved in the studied structures.
On the other hand, Beam Propagation Method (BPM) is a mature and well-known simulation tool to model light propagation in different types of waveguides [36–44]. Comparing with other numerical techniques, BPM remains highly efficient without sacrificing substantial accuracy. Its speed and required memory scale linearly with cross section area [45] on the contrary with most of other techniques such as FDTD which scale proportional to the device volume. BPM is known to be efficient in modelling light propagation along both straight and curved waveguides. Also some eigenmode analyses have been done on WGM cavities using this method [46]. The conventional BPM is formulated based on the Fresnel approximation which is valid where light propagates close to the propagation axis [38–40], so it will fail to simulate WGM cavities because of their highly curved nature. So far to overcome this limitation for the bent structures, high-order algorithms known as wide-angle BPM [47] or the conformal mapping approach [48] are introduced. Here we offer to reformulate the BPM in cylindrical coordinate system to analyze the WGM cavities. In this thesis we formulate the scalar Finite Difference Beam Propagation Method (FD-BPM) to simulate the light propagation in different size cavities from micrometer to millimeter scale. The computed field distribution correctly includes the radiative field component, which a mode analysis technique would fail to simulate. On the contrary to other modelling methods, BPM does not need huge computation resource for large cavity simulations and to the best of our knowledge such a simulation by other numerical methods has yet to be reported. Spiral shaped delay lines, which recently have been proposed [49, 50], is a decent example of large WGM cavities.
In the end, some other properties of WGM cavities are studied as the applications of our FD-BPM tool. The WGM microcavities with a small size scatterer receive
notable attention due to their applications in unidirectional emission and sensing en-hancement. Lots of theoretical studies [51–53] have been done as well as numerical and experimental researches [54–56]. In the scope of sensitivity enhancement, a stronger local field is desirable. One way to achieve this goal is to tailor the plasmonics to pro-duce a higher intensity fields [57,58]. Hybrid WGM-plasmonic sensors are intropro-duced as an enhanced single particle sensor [59–61]. Notable quality factor degradation is the main disadvantage of WGM-plasmonic cavities. This would be even more signif-icant considering the residues of the deposition. Alternatively, a new non-plasmonic technique [62] has been proposed for enhancing the sensitivity of WGM single parti-cle sensors. Here we verified this structure by the FD-BPM tool. The results show that at least three times higher sensitivity is achievable with a negligible drop in the quality factor (Q> 108). Furthermore, Waveguide-cavity coupling and asymmetric cavities were modelled and the results are presented in the application chapter.
1.2
Thesis Outline
This thesis structure is as follows:
Chapter 2 provides principles of WGM cavities as well as the basics of finite difference scheme.
Chapter 3 illustrates the finite difference mode solver, starting from straight waveguide to WGM cavities.
Chapter 4 demonstrates the FD-BPM formulation in cartesian and cylindrical coordinate systems and for 2D and 3D cases. Two different boundary conditions are depicted in this chapter.
Chapter 5 illustrates some applications of WGM cavities using the FD-BPM tool. The new structured WGM cavity for enhanced nano particle sensing is also introduced in this chapter. Huge cavity simulation, waveguide-cavity coupling, nano particle detection, and asymmetric cavities are discussed in this chapter.
Chapter 2
Background
Whispering gallery mode cavities are at the frontier of research on a wide range of topics. The term ”whispering gallery” comes from the gallery in St. Paul’s Cathedral in London where a whisper can be heard from a long distance on the opposite side of the gallery. In fact, the acoustic waves were traveling in a circular path close to the wall which acts as an acoustic waveguide with low loss. Using the same phenomenon researchers started to investigate the optical whispering gallery resonators in 1980s [63]. In an optical whispering gallery cavity, any circular symmetric shaped dielectric can be used as the waveguide. It can be shaped as sphere, cylinder, disk, ring, or toroid. The light will propagate mostly inside the cavity close to the external surface. In the resonance condition the phenomenal fact is that the light will be confined in a small volume over a long period of time. Ultra high quality factor WGM cavities, which refer to Q> 108 [64], are achievable while two fabrication oriented criteria are met. First, the dielectric material should have low loss (< 0.2dB/km) and second, the external surface of the cavity should be atomically smooth. These kind of cavities are attractive for researchers working on a broad range of topics from biosensing [3, 6, 7] , nonlinear optics [65] , and laser physics [12, 15] , to fundamental physics such as cavity quantum electrodynamics [17].
2.1
Quality factor
Having low loss in an oscillation cycle is favorable for resonators. Quality factor (Q) is the quantitative parameter to characterize resonator loss. It is defined as 2π times the ratio of the energy stored in the resonator to the dissipated energy per cycle.
Equivalently for a WGM cavity we can define Q as
Q = 2πf0 ×
Estored
Ploss
= ω0τ (2.1)
where Estored and Ploss are the stored energy in the cavity and the lost power re-spectively. f0 is the resonance frequency, ω0 is the resonance angular frequency, and
τ is the photon life time in the cavity. Quality factor indicates the energy loss of
the cavity which is also related to the imaginary part of refractive index. Another way of defining the quality factor is based on imaginary part of mode number (m) or resonance wavelength (λres) [66, 67]
Q = Re[m]
2Im[m] =
Re[λres] 2Im[λres]
(2.2)
In numerical modelling, with the known input and output (after one round of prop-agation) power, one can calculate the quality factor by
Q = 2πm× Pi
Pi− Po
(2.3) in which Pi and Po are the input and output power respectively.
So far we have defined the intrinsic quality factor (Qint) which characterize the cavity losses. The total Q (Qtot) is defined as
1 Qtot = 1 Qint + 1 Qcoup (2.4)
Here, Qcoup is the coupling Q which shows the energy loss while coupling the light into the cavity and also from the cavity out to the coupler. In general, the intrinsic
Q is separable to four different terms as follows
1 Qint = 1 Qabs + 1 Qrad + 1 Qss + 1 Qcontam (2.5)
Qabs is absorption Q represents the waveguide material absorption. This is usually the dominant term to calculate the Qint. Knowing the attenuation coefficient (α) of the waveguide material, Qabs can be calculated by [68]
Qabs = 2πn
λresα
where λres is the resonance wavelength and n is the refractive index of the cavity. Radiation Q (Qrad) is related to the WGM radiation loss and is the dominant term for the smaller cavities (radius< 40µm). Qssis the scattering loss which represents the energy loss due to the cavity surface roughness. Surface absorption or contamination loss (Qcontam) arises from contaminants attached to the surface of the cavity.
2.2
Waveguide-Cavity Coupling
There are different ways to couple the light into a WGM cavity. Free space coupling is one of them [69]. It has low efficiency but still helpful in some certain applications. As the high circulating power is usually a desirable property, higher coupling efficiency is favorable. Evanescent field couplers are known to be more efficient and among them, taper-cavity coupling [1, 70] has the highest efficiency. Prism coupler [71, 72] and side polished fiber coupler [69, 73] are two other methods which work based on evanescent field coupling.
Figure 2.1: An illustration of taper-cavity coupling. The dashed red box shows the coupling area as a 4-port transfer function. The 4 ports are the taper and cavity inputs and outputs.
Here we describe some of the coupling parameters which will be used to charac-terize the taper-cavity coupling. More details and derivations can be found in [1, 68]. Consider the coupling system as a 4-port transfer function as illustrated in Fig.2.1.
The ratio of the field amplitude at the input and output of the taper is defined as reflection coefficient R. R is very close to unity due to the total internal reflection (TIR) at the surface of the taper. To simplify the equations, the reflection coefficient is assumed to be the same for the cavity. Another assumption is that the coupling to the cavity and to the taper is the same due to the symmetry (Fig.2.1). Coupling coefficient T is defined as the ratio of the field amplitudes after and before coupling. Neglecting the reflection loss at the coupling junction we have,
|T |2 +|R|2 = 1 (2.7)
R≃ 1 − T
2
2 ∼ 1 (2.8)
The coupling Q is defined as [68]
Qc=
Rω0τ0
2(1− R) (2.9)
where ω0 is the resonant angular frequency and τ0 is the circulation time for the field to propagate inside the cavity. It can be approximated to be
τ0 = 2πnr
c (2.10)
where r is the cavity radius, n is the refractive index, and c is the speed of light. Another parameter needs to be defined is the normalized coupling parameter K. It is defined as the ratio of the intrinsic Q to the coupling Q, K = Q0
Qc. Adding the
taper output field to the field coupled from the cavity to the taper, one is able to calculate the taper output transmission as
T ransmission = (1− K
1 + K)
2 (2.11)
The relation between the coupling factor K and the transmission is plotted in Fig. 2.2. Three different regions is illustrated on the figure which represent three states of coupling. When the taper is far from the cavity small amount of light will be coupled to the cavity. In this state, which is known as ”under coupled”, the coupling loss is smaller than the intrinsic loss or in other words Qc > Q0. That small amount of coupled light will accumulate in the cavity and then coupled back to the taper and interfere destructively with the taper output field causing a drop in the
trans-10−2 100 102 0 0.2 0.4 0.6 0.8 1 K Transmission Critical Coupling
Under Coupled Over Coupled
Figure 2.2: The taper output transmission vs. the normalized coupling parameter K.
mission. As the taper is moved closer to the cavity, coupling will get stronger and the transmission will get lower until at a certain distance it reaches to zero. In this state, which is called ”critical coupling”, the coupling loss equals the intrinsic loss and K = 1. Decreasing the taper-cavity distance, more light will couple to the cavity and in return a strong beam of light, which is bigger than the taper output field, will couple back to the taper. This will cause an increment in the taper transmission. In this region, which is called ”over coupled”, the coupling loss is bigger than the intrinsic loss, Qc< Q0.
2.3
Buildup Factor
A certain amount of power (Pin) is coupled to the cavity on each round of propagation. On the next round this power will add up to the circulating power, resulting in a bigger
circulating power (almost doubled, for high Q cavity and critical coupling situation). This continuously increases the circulating power until it saturates. Certain portion of circulating power will be lost in each round of propagation (mainly because of the absorption and radiation). The higher the circulating power, the higher the power loss per cycle. Saturation point will be reached when the power loss per cycle reaches the input power. Buildup factor (B) is defined as the ratio of the saturated power (Psat) to the input power (Pin).
Assuming the first field coupled to the cavity is E1 = ψ, the second round field (E2) will be the sum of propagated field and the new coupled field,
E2 = ψ + ψej(mr+jmi)2π (2.12) where mr and mi are the real and imaginary part of mode number respectively. The absolute value of E2 would be
|E2| = ψ(1 + e−mi2π) (2.13)
Therefor after n rounds of propagation we will have
|En| = ψ n−1 ∑ k=0 e−mi2πk (2.14) n → ∞ : |En| = ψ 1− e−mi2π (2.15)
For high quality cavities the imaginary part of mode number (mi) is a small number. To calculate the builup factor (B),neglecting the m2i in the Maclaurin series expansion of e−mi2π, we will have B = |E∞| 2 |E1|2 = 1 (2πmi)2 = ( Q πmr )2 (2.16)
2.4
Finite Difference Scheme
In this thesis we present an approach to solve a well-known differential equation, Helmh¨oltz equation. Finite difference method is a way of approximating the deriva-tives in order to solve the differential equations numerically. This method is known to be straight forward and easy to implement. In this section we will derive the first
order finite difference scheme both for uniform and nonuniform griding. Higher order finite difference scheme (five-point) will be presented as well.
2.4.1
Uniform Finite Difference Method
We want to find a first order approximation for the first and second derivatives of function f (x). Writing the Taylor Series expansion, we have
f (x0+ ∆x) = f (x0) + f′(x0) 1! ∆x + f′′(x0) 2! ∆x 2+f(3)(x0) 3! ∆x 3+ ... (2.17) Working f′(x0) out of Eq.2.17 we will have,
f′(x0) =
f (x0+ ∆x)− f(x0)
∆x −
f′′(x0)
2! ∆x− ... (2.18)
The first order approximation will be attained by neglecting the second term and on. In this case the error is in the order of step size (∆x). It is common to state this fact as
f′(x0) =
f (x0+ ∆x)− f(x0)
∆x + O(∆x) (2.19)
When approximating the derivatives in an equation using the finite difference method, the key point is to keep the same central points for all of the terms. For instance in Eq.2.19 the central point is x0+∆x2 . If it is needed to move that to x0 we will do the following fx′ 0+∆x2 (x0) = f (x0+ ∆x)− f(x0) ∆x (2.20) fx′ 0−∆x2 (x0) = f (x0)− f(x0− ∆x) ∆x (2.21) f′ x0+∆x2 (x0) + f ′ x0−∆x2 (x0) 2 = f ′ x0(x0) = f (x0 + ∆x)− f(x0− ∆x) 2∆x (2.22)
Now we derive the second derivative. From Eq.2.17 we have
f (x0+ ∆x) = f (x0) + f′(x0) 1! ∆x + f′′(x0) 2! ∆x 2+ f(3)(x0) 3! ∆x 3+ ... (2.23) f (x0− ∆x) = f(x0)− f′(x0) 1! ∆x + f′′(x0) 2! ∆x 2− f(3)(x0) 3! ∆x 3+ ... (2.24)
One can work out the second derivative of f (x0) as
f′′(x0) =
f (x0+ ∆x)− 2f(x0) + f (x0− ∆x)
∆x2 + O(∆x
2) (2.25)
For the second derivative approximation the order of error is ∆x2 and the central point is at x0. One can derive higher order finite difference scheme with the same approach. At the end of this chapter a five-point general finite difference scheme is presented.
2.4.2
Non Uniform Finite Difference Method
Smaller grid spacing is necessary when studying smaller features. For instance, if a 50nm particle attached to the surface of a cavity is to be modeled, 10nm grid spacing or smaller is needed. The minimum grid points is desirable in order to occupy less computation resources and have a faster simulation tool. Non uniform griding is the solution when we need more grids in a certain area and less grids in other places on the computation window.
Figure 2.3: One dimensional finite difference scheme.
Assume the non uniform grid spacing as in Fig.2.3. From Eq.2.17 we have
f (xi+1) = f (xi+ ∆xi) = f (xi) + f′(xi)∆xi+ 1 2f ′′(x i)∆x2i + ... (2.26) f (xi−1) = f (xi − ∆xi−1) = f (xi) + f′(xi)∆xi−1+ 1 2f ′′(x i)∆x2i−1+ ... (2.27) Dividing the first equation by ∆x2i and the second one by ∆x2i−1 and subtracting them, one is able to derive the first order derivative as
f′(xi) = ∆xi−1 ∆xi f (xi+1)− ( ∆xi−1 ∆xi − ∆xi ∆xi−1)f (xi)− ∆xi ∆xi−1f (xi−1) ∆xi+ ∆xi−1 + O(∆x) (2.28)
By dividing both sides of Eq.2.26 and Eq.2.27 by ∆xi and ∆xi−1 respectively and adding them together, the second derivative will be derived as
f′′(xi) =
∆xi−1f (xi+1)− (∆xi−1+ ∆xi)f (xi) + ∆xif (xi−1) 1
2∆xi∆xi−1(∆xi+ ∆xi−1)
+ O(∆x2) (2.29)
A general formula to calculate any order of finite difference scheme is reported in [74]. Here we just present the five-point finite difference approximations of the first and second derivatives. Assuming Fig. 2.3, we define α terms as
αj = xi+(j−3)− xi (2.30)
The first derivative would be
f′(xi) =− (α2α3α4+ α2α3α5+ α2α4α5+ α3α4α5)f (xi−2) (α1− α2)(α1− α3)(α1− α4)(α1− α5) (2.31) −(α1α3α4+ α1α3α5+ α1α4α5+ α3α4α5)f (xi−1) (α2− α1)(α2− α3)(α2− α4)(α2− α5) −(α1α2α4+ α1α2α5+ α1α4α5+ α2α4α5)f (xi) (α3− α1)(α3− α2)(α3− α4)(α3− α5) −(α1α2α3+ α1α2α5+ α1α3α5+ α2α3α5)f (xi+1) (α4− α1)(α4− α2)(α4− α3)(α4− α5) −(α1α2α3+ α1α2α4+ α1α3α4+ α2α3α4)f (xi+2) (α5− α1)(α5− α2)(α5− α3)(α5− α4) and the second derivative is derived as
f′′(xi) = + 2(α2α3+ α2α4+ α2α5+ α3α4+ α3α5+ α4α5)f (xi−2) (α1− α2)(α1− α3)(α1− α4)(α1− α5) (2.32) + 2(α1α3+ α1α4+ α1α5+ α3α4+ α3α5+ α4α5)f (xi−1) (α2− α1)(α2− α3)(α2− α4)(α2− α5) + 2(α1α2+ α1α4+ α1α5+ α2α4+ α2α5+ α4α5)f (xi) (α3− α1)(α3− α2)(α3− α4)(α3− α5) + 2(α1α2+ α1α3+ α1α5+ α2α3+ α2α5+ α3α5)f (xi−1) (α4− α1)(α4− α2)(α4− α3)(α4− α5) + 2(α1α2+ α1α3+ α1α4+ α2α3+ α2α4+ α3α4)f (xi+2) (α5− α1)(α5− α2)(α5− α3)(α5− α4)
Chapter 3
Finite Difference Mode Solver
For any certain waveguide, each pattern of electromagnetic wave which can repeat itself periodically while propagating through that waveguide, is a mode. The mode mainly depends on the structure of the waveguide; ie, the shape and size of the waveguide. It also depends on the wavelength of propagating wave. Each mode has a pattern of electromagnetic field and propagation properties which is defined by propagation constant. Finite Difference (FD) eigenmode solver is a well known approach to calculate the modes. In the present chapter we will derive the scalar FD mode solver in cartesian and cylindrical coordinate systems.
Figure 3.1: Illustration of the cartesian and cylindrical coordinate systems.
We start from the Maxwell’s equations to derive appropriate formulism for the electromagnetic waves. Later we will solve them numerically. The Maxwell’s equa-tions define the relaequa-tions between the electrical field ⃗E(in volts per meter), the
meters), and the magnetic flux density ⃗B(amperes per square meter) as [75] ⃗ ∇ × ⃗E = −∂ ⃗B ∂t, (3.1) ⃗ ∇ × ⃗H = ∂ ⃗D ∂t + ⃗J , (3.2) ⃗ ∇ · ⃗B = 0 (3.3) ⃗ ∇ · ⃗D = ρ (3.4)
where ⃗J is the free current density(in amperes per square meter) and is zero in
non-conductive materials and ρ is the free charge density which is also assumed to be zero in our case. Here is the relation between ⃗D, ⃗B, ⃗E, and ⃗H,
⃗ D = ε ⃗E
⃗
B = µ ⃗H
in which the permittivity ε and permeability µ are defined as
ε = ε0εr,
µ = µ0µr.
where ε0 and µ0 are the permittivity and permeability of the vacuum, and εr and µr are the relative permittivity and permeability of the material. The materials we are going to use are non-magnetic so µr= 1. For a monochromatic electromagnetic field at an angular frequency ω, we have, at a location defined by ⃗r, at an instance of time t,
⃗
E(⃗r, t) = ⃗F (⃗r)exp(jωt) (3.5)
So, the equations (3.1) and (3.2) would be
⃗
∇ × ⃗E = −jωµ0H,⃗ (3.6)
⃗
Applying a vectorial rotation operator ∇× to Eq.3.6, we will have
⃗
∇ × (⃗∇ × ⃗E) = −jωµ0∇ × ⃗⃗ H (3.8)
We consider the following relation for vectors [75]
∇ × (∇ × K) = ∇(∇.K) − ∇2K (3.9)
Substituting Eq.3.9 and Eq.3.7, in Eq.3.8,
∇(∇ · E) − ∇2E = k2
0εrE (3.10)
where k2
0 = ω2ε0µ0. By assuming the relative permittivity, εr, to be constant and
ρ = 0, we can find from Eq.3.4 that∇.E = 0. So the Eq.3.10 will reduce to Helmh¨oltz
equation
∇2E + k2
0εrE = 0 (3.11)
3.1
Cartesian Mode Solver
In a cartesian domain (x, y, z), assuming the separation of variables method, one of the solutions to Eq.3.11 is
E = ψ(x, y, z)ejβz (3.12)
where β is the propagation constant. Since the structure is z-invariant, we assumed that the envelope of the electrical field (ψ) is not a function of z (ψ(x, y)). So, we have
∂2E
∂z2 =−β
2E (3.13)
From the Helmh¨oltz equation we can find the mode equation by substituting Eq.3.13 into Eq.3.11 d2E dx2 + d2E dy2 + k 2 0n 2 (x, y)E = β2E (3.14)
where n2(x, y) = εr and n(x, y) is the waveguide refractive index profile on the trans-verse plain (x, y). For a slab waveguide (y-invariant) the mode equation will be even simpler by ignoring the derivatives with regards to y. So we will have
d2E dx2 + k
2 0n
2(x)E = β2E (3.15)
In Eq.3.15 the β2 is the eigenvalue for the following linear operator ˆ P = d 2 dx2 + k 2 0n 2 (x) (3.16)
These eigenvalues are the squares of guided modes propagation constants and the eigenvectors are the mode patterns. Now we discretize the x-axis by ∆x steps. We have x = p∆x, where p is an integer. Also we have Ep = E(p∆x) and np = n(p∆x). So, as derived in chapter 2, the first order finite difference scheme for the second order derivative in Eq.3.15 would be
d2E
dx2 = α 2
x(Ep−1− 2Ep+ Ep+1) (3.17) where αx = 1/∆x. And for the second term
k02n2(x)E = k20n2pEp (3.18) Now, we can find the matrix format for the operator ˆP . ˆPi,j is the ith row and jth column element of this matrix. From Eq.3.17 and Eq.3.18 we get
ˆ Pk,k =−2α2x+ k 2 0n 2 p (3.19) ˆ Pk−1,k = ˆPk,k+1= α2x (3.20) Calculating the eigenvalues of this matrix using MATLAB, one is able to find all the modes of a slab waveguide. The fundamental mode (first order mode) and the third order mode of an 8-µm silica slab waveguide is plotted in fig.3.2. The dotted black line is showing the refractive index profile on the right axis. The waveguide is in the air with the refractive index of unity.
−8 −6 −4 −2 0 2 4 6 8 −1 −0.5 0 0.5 1 x (µm)
Eelectrical field (a.u.)
−8 −6 −4 −2 0 2 4 6 80.5
1 1.5
Refractive index
Figure 3.2: The first (red) and third (blue) order mode of a 8-µm silica slab waveguide in the air (left axis). The dotted black line indicates the refractive index profile of the waveguide on the right axis.
3.2
Cylindrical Mode Solver
In this section the same derivation will be done in the cylindrical coordinate system (Fig.3.1) to find the whispering gallery mode of a cavity. The propagation direction is along the ϕ axis. The laplace operator for the cylindrical coordinates is
∇2E = 1 r ∂ ∂r(r ∂E ∂r) + 1 r2 ∂2E ∂ϕ2 + ∂2E ∂z2 (3.21)
By substituting Eq.3.21 into Eq.3.11 we get
r ∂ ∂r(r ∂E ∂r) + ∂2E ∂ϕ2 + r 2∂ 2E ∂z2 + r 2k2 0n2(z, r, ϕ)E = 0 (3.22)
where n(z, r, ϕ) is the refractive index of the structure. As we assume a ϕ-invariant structure to solve the modes, the refractive index profile and the envelope of the electrical field are not a function of ϕ. We assume the solutions to be
E(r, ϕ) = F (z, r)ejmϕ (3.23)
in which m is the mode number. The resonance condition will be satisfied when the real part of m is an integer value. Substituting this into Eq.3.22, we will find the mode equation ˆ KE = m2E (3.24) where ˆK is ˆ K = r2 ∂ 2 ∂r2 + r ∂ ∂r + r 2∂ 2E ∂z2 + r 2k2 0n2(z, r) (3.25) A z-invariant structure can be modelled in a 2D scheme. Also one can approximate the 3D structure into a 2D one by the Effective Index Method (EIM) [76]. In either cases, Eq.3.25 will shrink into
ˆ K = r2 ∂ 2 ∂r2 + r ∂ ∂r + r 2k2 0n 2(r) (3.26)
Solving the eigenvalues of this matrix numerically, one is able to find the eigenvalues which are the mode number (m) squares. The corresponding eigenvectors are the mode patterns. The last step is to employ the resonance condition criteria to obtain the resonance wavelength and the resonance mode number. In resonance condition the mode number would be an integer value; So one is able to find the resonance wavelength from,
λres = [m]
m λ0 (3.27)
where [m] is the absolute amount of m, and λ0 is the central wavelength. Fig.3.3 illustrates the three first whispering gallery modes of a 45-µm major radius silica ring resonator. Major radius is defined as the distance between the coordinate origin to the center of the waveguide channel (Fig.3.1). The width of the guiding channel is 10µm. The dotted black line indicates the refractive index profile of the waveguide
on the right axis. 38 40 42 44 46 48 50 52 −0.1 −0.05 0 0.05 0.1 0.15 radius (µm)
electrical field (a.u.)
38 40 42 44 46 48 50 52 0.6 0.8 1 1.2 1.4 1.6 refractive index
Figure 3.3: The first (blue), second (red), and third (green) order whispering gallery modes of a silica ring resonator in the air (left axis). The major radius is 45µm and the width of the waveguide is 10µm. The dotted black line shows the refractive index profile of the waveguide on the right axis.
To find the 3D structure whispering gallery modes, one is to discretize Eq.3.24 in r-z plain, with the ˆK defined in Eq.3.25. Fig.3.4 illustrates the four different
whispering gallery modes of a silica micro toroid in water in the central wavelength of 970nm. The red structure in Fig.3.1 is representing a toroid shape cavity. The major radius of the toroid is 45µm and the minor diameter is 10µm. The black lines are showing the boundaries of the cavities. The mode on the top left in Fig.3.4 is the fundamental mode. The calculated mode number for the fundamental mode is 457 and the resonance wavelength is 969.12nm.
Figure 3.4: Four different whispering gallery modes of a silica micro toroid in water in the central wavelength of 970nm. The one on the top left is the fundamental mode. The major radius of the toroid is 45µm and the minor diameter is 10µm. The black lines show the boundaries of the cavities.
Chapter 4
Finite Difference Beam
Propagation Method
This chapter will start with the derivation of the conventional FD-BPM formulation in the cartesian coordinate system. Transparent boundary condition will be added to the formulation to avoid the spurious reflections from the computation window edges. Then the derivations will be done in the cylindrical coordinate system to simulate the light propagation in the WGM cavities. For the later formulation, Perfectly Matched Layer (PML) will be utilized to minimize the reflections more efficiently. Numerical errors and stabilities will be studied as well as the convergence rates.
4.1
Two Dimensional Cartesian FD-BPM
In this section two dimensional scalar FD-BPM in cartesian coordinate system will be derived and used to model straight and slightly bent waveguide.
4.1.1
Formulation
The 2D wave equation can be derived from the Maxwell’s equations. Expanding Eq. 3.6 and Eq. 3.7, six set of Maxwell’s equations will be obtained as follows
∂Ez ∂y − ∂Ey ∂z =−jωµ0Hx, (4.1) ∂Ex ∂z − ∂Ez ∂x =−jωµ0Hy, (4.2) ∂Ey ∂x − ∂Ex ∂y =−jωµ0Hz, (4.3) ∂Hz ∂y − ∂Hy ∂z = jωε0εrEx, (4.4) ∂Hx ∂z − ∂Hz ∂x = jωε0εrEy, (4.5) ∂Hy ∂x − ∂Hx ∂y = jωε0εrEz, (4.6)
The direction of propagation is assumed to be along the z axis and due to the 2D modelling assumption, the structure is uniform in the y direction and consequently
∂/∂y = 0. Also, in the TE mode, the electric field is zero in the longitudinal direction(Ez = 0). Considering these, from Eq.4.6 we find, ∂Hy/∂x = 0. This means that Hy should be a constant in x direction. We know that the Hy in a long distance from the core should be equal to zero; So, Hy(x) has to be zero everywhere. By considering Hy = 0 and substituting this into Eq.4.4 we will have, Ex = 0. All in all, for the TE mode we have,
Ez = Ex = Hy = 0 (4.7)
Finding Hx and Hz as a function of Ey from Eq.4.1 and Eq.4.3 respectively and substituting them into Eq.4.5, we will have the Ey representation wave equation
∂2E y ∂z2 + ∂2E y ∂x2 + k 2 0εrEy = 0 (4.8)
where k0 is the free space wavenumber and is defined as k02 = ω2ε0µ0. The Hx can be derived by differentiating Eq.4.5 with respect to z and assuming that the variation of the relative permittivity, εr, along z (propagation axis) to be negligible(i.e
∂εr/∂z ≃ 0). So we have ∂2H x ∂z2 − ∂2H z ∂z∂x = jωε0εr ∂Ey ∂z (4.9)
Also, from Eq.3.3 we have (µr = 1),
∂Hx
∂x +
∂Hz
∂z = 0 (4.10)
Substituting this equation and Eq.4.1 in Eq.4.9, the Hx representation of the wave equation can be derived as follows
∂2Hx ∂z2 + ∂2Hx ∂x2 + k 2 0εrHx = 0 (4.11)
Using the same derivations for the TM mode, one will find
Hz = Hx = Ey = 0 (4.12)
As in TE mode, by doing the some derivation for TM mode one can extract the Ex and Hy representations of wave equation. These will be as follows,
∂2E x ∂z2 + ∂ ∂x( 1 εr ∂ ∂x(εrEx)) + k 2 0εrEx= 0 (4.13) ∂2H y ∂z2 + εr ∂ ∂x( 1 εr ∂Hy ∂x ) + k 2 0εrHy = 0 (4.14)
For a 2D waveguide and the TE mode, one can divide the principal field Ey(x, z) into two components, the envelope and the phase term,
Ey(x, z) = ψ(x, z)ejβz (4.15)
where β = nek0 is the propagation constant. ne is the effective refractive index or the reference refractive index. Free space wavenumber k0, is defined as 2π divided by the wavelength. The first and second derivatives of electrical field with respect to z
would be ∂Ey ∂z = ∂ψ ∂ze jnek0z+ jn ek0ψejnek0z (4.16) ∂2E y ∂z2 = ∂2ψ ∂z2e jnek0z+ 2jn ek0 ∂ψ ∂ze jnek0z − n2 ek 2 0ψe jnek0z (4.17)
Substituting Eq.4.15 in Eq.4.8, we have
∂2ψ ∂z2 + 2jnek0 ∂ψ ∂z − n 2 ek 2 0ψ + ∂2ψ ∂x2 + k 2 0n 2(x, z)ψ = 0 (4.18) Using the slowly varying envelop approximation (SVEA) i.e.
|∂2ψ
∂z2| ≪ 2nek0|
∂ψ
∂z| (4.19)
Eq.4.18 reduce to the Fresnel wave equation
j∂ψ ∂z = ˆHψ (4.20) where, ˆ H = 1 2nek0 [− ∂ 2 ∂x2 + (n 2 ek 2 0 − k 2 0n 2 (x, z))]
The SVEA holds as long as the refractive index variation along the propagation axis is small. Using the first order uniform finite difference(FD) scheme (chapter 2), one is able to solve Eq.4.20 numerically. Uniformly discretizing the x and z coordinates, we have
x = p∆x, (4.21)
z = l∆z, (4.22)
in which p and l are integers. Also, we will show the ψ(x, z) and n(x, z) as follows,
ψ(x, z) = ψlp, (4.23)
n(x, z) = nlp. (4.24)
the x axis of operator H. For the first term we have
∂2ψ
∂x2 =
ψp+1− 2ψp+ ψp−1
(∆x)2 (4.25)
And for the second term
(n2ek20− k02n2(x, z))ψ(x, z) = (ne2k20− k02n2p)ψp (4.26) We also will use these definitions
αx = 1 ∆x (4.27) αz = 1 ∆z (4.28)
Substituting Eq.4.25, Eq.4.26 and Eq.4.27 into Eq.4.20, we get
j∂ψ ∂z =
1 2nek0
[(n2ek02+ 2α2x− k20n2p)ψp− α2xψp+1− α2xψp−1] (4.29) The z axis discretization for the left hand side of Fresnel wave equation(Eq.4.20) is expressed as j∂ψ ∂z = j ψpl+1− ψpl ∆z = jαz(ψ l+1 p − ψ l p) (4.30)
As you can see in Eq.4.30, the difference center is the point l + 12 midway between l and l + 1. In order to have the same difference center in both side of Eq.34 we modify it as follows. So, the discretized form of Fresnel wave equation (Eq.4.20) will be
j(ψpl+1− ψpl) = 1 4nek0αz [(ne2k02+ 2α2x− k02(nlp)2)ψpl − α2 xψ l p+1− α 2 xψ l p−1] + 1 4nek0αz [(n2ek02+ 2α2x − k2 0(n l+1 p ) 2)ψl+1 p − α 2 xψ l+1 p+1− α 2 xψ l+1 p−1] (4.31)
A(p) = −α 2 x 4nek0αz , (4.32) B(p) = n 2 ek20 + 2α2x− k02(nl+1p )2 4nek0αz − j, (4.33) C(p) = −α 2 x 4nek0αz , (4.34) D(p) = 1 4nek0αz [−(4jαznek0+ n2ek 2 0 + 2α 2 x− k 2 0(n l p) 2)ψl p + α2xψlp+1+ α2xψpl−1] (4.35) Defining a tridiagonal matrix, T, in which the main diagonal is B(p), the diagonal+1 is C(p), and the diagonal−1 is A(p). In this case Eq.4.31 would be
T (p)ψl+1 = D(p). (4.36)
which is expanded in Eq.4.50.
4.1.2
Transparent Boundary Condition
Because the analysis window is not infinite, we have to consider the boundary con-ditions to compensate it. Transparent Boundary Condition(TBC), developed by Hadley [77], is a way of suppressing the reflections at boundaries. Actually, we con-sider two hypothetical nodes at p = 0 and p = M + 1(M is the number of discrete nodes in x direction) and find the effect of them on the field at p = 1 and p = M respectively. The wave function for the nodes p = 0, 1, 2 is
ψ0 = G(z)ejkx(x0−∆x) (4.37)
ψ1 = G(z)ejkxx0 (4.38)
ψ2 = G(z)ejkx(x0+∆x) (4.39)
where x0 is the location for the first node in the window and ψ0 is the field of hypo-thetical node just outside the window. Dividing Eq.4.38 by Eq.4.37 and Eq.4.39 by
Eq.4.38 ejkx(∆x)= ψ1 ψ0 (4.40) ejkx(∆x)= ψ2 ψ1 (4.41) From these two equations we will get
ψ0 = η1ψ1 (4.42) in which η1 = ψ1 ψ2 (4.43) By using the same approach we can find the ψM +1 as
ψM +1 = ηMψM (4.44) where ηM = ψM ψM−1 (4.45)
Now, we return to the Eq.4.36. As can be seen, by knowing the ψl at p = 1, 2, ..., M we are able to find the ψl+1. We are going to have the Eq.4.36 in the matrix form and then try to solve it numerically. To adopt the boundary condition in the BPM, the matrix element at p = 1 will be modified to B(1) according to
B′(1) = −α 2 x 4nek0 ηl1+ B(1) (4.46) where η1l = ψ l 1 ψl 2 (4.47)
By the same approach for p = M , instead of having C(M ), we will have a new definition for B(M ) B′(M ) = −α 2 x 4nek0 ηMl + B(M ) (4.48) in which ηlM = ψ l M ψl M−1 (4.49)
Finally, the matrix format for the Eq.4.36 is B′(1) C(1) A(2) B(2) C(2) . .. . .. ... . .. ... . .. C(M− 1) A(M ) B′(M ) ψ1l+1 ψ2l+1 ψ3l+1 .. . ψl+1M−1 ψMl+1 = D(1) D(2) D(3) .. . D(M − 1) D(M ) (4.50)
In Eq.4.50, by knowing the the refractive index in the whole structure, the coefficient matrix is known. So, in each step the input is ψl (D coefficient) and the output is
ψl+1. Repeating these steps, we will find the electrical field of any point in the 2D structure. This ψ is actually the envelope of the electrical field. In order to add the phase we multiply each ψl by ejnek0l∆z. As an example, the field intensity of an 8-µm
diameter bent fiber is plotted in Fig.4.1. The grid spacings in z and x direction are 4µm and 150nm respectively.
z (mm) x ( µ m) 0 1 2 3 4 −20 −15 −10 −5 0 5 10
Figure 4.1: Field intensity inside a 4-mm long, 8-µm core diameter bent fiber.
4.2
Cylindrical FD-BPM
In order to find the mode equation in cylindrical coordinate system, we will again start from the Helmh¨oltz equation (Eq. 3.11). In a WGM cavity, the electric field can be approximated as propagating in the form of
E(r, z, ϕ) = ψ(r, z, ϕ)ejmϕ (4.51)
where m is the mode number. The real part of m represents the azimuthal order of a cavity resonance mode while the imaginary part of that characterizes the attenuation of the field in the azimuthal direction. The first and second derivatives of electrical
field with respect to ϕ would be ∂E ∂ϕ = ∂ψ ∂ϕe jmϕ+ jmψejmϕ (4.52) ∂2E ∂ϕ2 = ∂2ψ ∂ϕ2e jmϕ + 2jm∂ψ ∂ϕe jmϕ− m2 ψejmϕ (4.53)
The mode equation can be derived by replacing Eq.4.51 into the Helmh¨oltz equation (Eq. 3.11). The BPM key assumption is that the field envelope varies slowly along the ϕ axis (slowly varing envelope approximation, SVEA). By this assumption we can ignore the second derivative with respect to the ϕ. The result is called the Fresnel wave equation, −2jm r2 ∂ψ ∂ϕ = ˆKψ (4.54) where ˆ K = ∂ 2 ∂r2 + ∂2 ∂z2 + 1 r ∂ ∂r + k 2 0n 2(r, z, ϕ)−m2 r2 (4.55)
4.2.1
2D FD-BPM
In the 2D case, where the structure is assumed to be z-invariant, Eq. 4.54 can be rewritten as −j ∂ ∂ϕψ = r2 2m ∂2ψ ∂r2 + r 2m ∂ψ ∂r + ( r2k2 0n2(r, ϕ) 2m − m 2)ψ (4.56)
Now we discretize the computation window uniformly so that the coordinates (rp, ϕl) of each grid (p, l) can be expressed as rp = p∆r, ϕl = l∆ϕ, and ψ(rp, ϕl) = ψp,l. One can evolve the field at ϕl from a previous azimuthal angle ϕl−1 according to
apψp−1,l+1+ (j∆ϕ1 + bp)ψp,l+1− cpψp+1,l+1 =−apψp−1,l+ (j∆ϕ1 − bp)ψp,l+ cpψp+1,l (4.57) where ap = p(18m−2p) bp = p2(2−∆r2k2 0n2p,l+1) 4m + m 4 cp = p(1+2p)8m (4.58)
Here ∆r and ∆ϕ are grid spacings along the ˆr and ˆϕ directions, as illustrated in
Fig.3.1. Also, np,l is the refractive index of the waveguide structure at each point. Rearranging Eq. (4.57) into a matrix form, we obtain
˜
Hψl+1 = ˜Dψl (4.59)
where ˜H and ˜D are two tridiagonal matrices. From the later equation, one is able
to find the full field profile from the excitation field at l = 0 in a step by step procedure [78].
4.2.2
Perfectly Matched Layer
Perfectly Matched Layer (PML) is another type of boundary condition which acts as an artificial absorption layer at the edges of computation window. Comparing with the TBC, PML is known to be more efficient in minimizing the spurious reflection at the boundaries. The PML is implemented by replacing the radial derivative with [79]
∂ ∂r → Σ ∂ ∂r (4.60) where Σ = 1 1 + jσ(r)ω (4.61) in which σ(r) is defined as σ(r) = { σ0(r− r0)2 Inside the PML 0 Elsewhere (4.62)
r0 is the inner radius of the lossy layer (PML layer). σ0 is a wavelength independent constant characterized by the light attenuation strength within the PML. The impor-tant point is to choose the right value for σ0. To find the optimized value, a beam of light is launched against the computation window wall through a straight waveguide and the reflection is measured for different values of σ0. The PML thickness is also important and here, we choose a 2µm one. The different values of reflectivity in dB versus σ0 are presented in Fig. 4.2.
1016 1017 1018 −50 −40 −30 −20 −10 0 σ0 Reflectivity (dB)
Figure 4.2: Reflectivity in dB vs. different σ0 values. σ0 = 2.5× 1016 is the optimum value.
Fig.4.3 show the field intensity for three different σ0 values, wherein the white lines are the inner PML edges. Fig.4.3 (a) illustrates an extreme case where the light reaches the outer edge of an undamped PML (i.e. σ0 = 0) and is reflected back. In the case of an overdamped PML (i.e. σ0 = 1020), as illustrated in Fig.4.3 (b), light is reflected at the PML’s inner edge. The optimum value of σ0 is found at 2.5× 1016 in Fig.4.3 (c), where a minimal reflectivity of -50 dB is attained. Note that over a wide span of σ0between 1016and 2×1017the reflectivity at the computation window edge is kept below -30 dB. The optimal value was then utilized for the remaining simulations given that σ0 was insensitive to, for example, the incident angle at the PML edge, the wavelength, and the waveguide structures inside the computation window [80].
(a) (b)
(c)
Figure 4.3: The field intensity of a straight waveguide, launching toward the compu-tation window edge, for three different σ0 values: 0 (top left) , 2.5× 1016 (bottom), and 1020 (top right). The white lines indicate the inner 2-µm PML edges.
4.2.3
3D FD-BPM
Adding the PML to the 3D Fresnel wave equation (Eq. 4.54), we obtain
−2jm r2 ∂ψ ∂ϕ = ˆKψ (4.63) where ˆ K = Σ2 ∂ 2 ∂r2 + Σ 2 ∂2 ∂z2 + ( Σ r − j ω ∂σ ∂rΣ 3)∂ ∂r − j ω ∂σ ∂zΣ 3 ∂ ∂z + k02n2(r, z, ϕ)−m 2 r2 (4.64)
In order to solve Eq. 4.63 numerically, we implement the uniform first order finite difference scheme (chapter 2) by discretizing the r-axis by ∆r steps and the z-axis by ∆z steps as shown in Fig. 4.4 .We have r = p∆r, z = l∆z, and ϕ = q∆ϕ, where p, l, and q are integer values. Therefore the 2D r-z plain would be transformed into a 1D vector of size p× l. Further we discretize the refractive index profile, the field, and the PML parameters(Σ and σ) so that we have nqp,l, ψp,lq , Σqp,l, and σqp,l.
Figure 4.4: Discretization scheme in r-z plain.
Rearranging the wave equation, the final discretized equation would be, ˜
Hψq+1 = ˜Dψq (4.65)
able to calculate the field distribution at ϕ = (q + 1)∆ϕ from the previous step field distribution at ϕ = q∆ϕ and so on. Therefor the propagation pattern can be obtained for any desired length of propagation. The diagonals of ˜H and ˜D matrices will be
H[−p] = r 2Σ2 ∆z2 + jr2 2ω∆z ∂σ ∂zΣ 3 (4.66) H[−1] = r 2Σ2 ∆r2 + jr2 2ω∆r ∂σ ∂rΣ 3− rΣ 2∆r (4.67) H[0] = −2r 2Σ2 ∆r2 − 2r2Σ2 ∆z2 + r 2k2 0n 2 p,l (4.68) H[+1] = r 2Σ2 ∆r2 − jr2 2ω∆r ∂σ ∂rΣ 3+ rΣ 2∆r (4.69) H[+p] = r 2Σ2 ∆z2 − jr2 2ω∆z ∂σ ∂zΣ 3 (4.70) D[−p] = r 2Σ2 ∆z2 + jr2 2ω∆z ∂σ ∂zΣ 3 (4.71) D[−1] = r 2Σ2 ∆r2 + jr2 2ω∆r ∂σ ∂rΣ 3− rΣ 2∆r (4.72) D[0] = −2r 2Σ2 ∆r2 − 2r2Σ2 ∆z2 + r 2k2 0n 2 p,l (4.73) D[+1] = r 2Σ2 ∆r2 − jr2 2ω∆r ∂σ ∂rΣ 3+ rΣ 2∆r (4.74) D[+p] = r 2Σ2 ∆z2 − jr2 2ω∆z ∂σ ∂zΣ 3 (4.75)
where H[0] and D[0] are the main diagonals of matrices ˜H and ˜D respectively. H[−p], H[−1], H[+1], and H[+p] are the diagonal minus p, minus one, plus one, and plus p of matrix ˜H respectively. And finally D[−p], D[−1], D[+1], and D[+p] are the
diagonal minus p, minus one, plus one, and plus p of matrix ˜D respectively.
4.2.4
Characterization
To characterize the 2D FD-BPM, we will first simulate a silica ring resonator immersed in water with the major radius of 45µm. The width of the waveguide is 10µm. The refractive index of the silica is 1.4508 + j(7.11× 10−12) [81, 82] at a wavelength of 970 nm and the surrounding water has a refractive index of 1.327+j(3.37×10−6) [83]. The intensity distribution for a quarter a round is plotted in Fig.4.5 in logarithm scale. The ring is excited with it’s fundamental mode. The window size in the r-direction
is 36µm with 1601 grid points. Also there are 750 grid points on the ϕ direction. Solid white lines indicate the two edges of the ring. As it can be seen, most of the power is propagating close to the outer edge. For this simulation a thick PML layer (5µm) is employed to make sure that the reflections are minimized. The boundaries of the PML is shown by the dashed black lines. The real part of mode number and the resonance wavelength calculated by the mode solver is 458 and 970.25nm respectively. Furthermore, since the intensity is plotted in the logarithm scale, the small portion of radiation is observable as well as the efficiency of the PML which absorbs the spurious reflections from the computation window outer edge.
Figure 4.5: Field intensity of a 45µm major radius and 10µm width silica ring, in logarithm scale. The radiations is observable. The solid white lines are illustrating the ring borders and the dashed black lines are showing the PML layer boundaries.
Next, a huge silica microdisk (1mm diameter) is modelled using the 2D FD-BPM. The field intensity for a portion of disk is provided in Fig.4.6. A supplementary video clip, showing the propagation in this cavity is available at [84]. Note that simulation
of the 1mm diameter microdisk in Fig.4.6 spanned 23 seconds on a desktop computer equipped with a 3.10 GHz AMD FX-8120 8-core processor and 32 GB of memory. In contrast, a simulation of such a large cavity has yet to be reported from other established methods, such as BEM, finite element method (FEM), or FDTD.
Figure 4.6: Intensity distribution of a huge microdisk with the radius of 1mm.
To characterize the 2D FD-BPM more quantitatively, quality factor is calculated for different grid spacings. As explained in chapter 2, one is able to work out the quality factor from
Q = 2πmr
P (ϕ = 0)
P (ϕ = 0)− P (ϕ = 2π) (4.76)
where P (ϕ) is the total power at the azimuthal angle ϕ and mr is the real part of mode number. Again the propagation of light in a 45µm major radius silica ring is modelled for half a round. This time a 20µm window size is chosen in the r direction.
Fig.4.7 (a) presents quality factor values as a function of radial and azimuthal grid spacing (∆r, ∆ϕ). As illustrated, the Q converges from 5.4×106 towards 4.99×106 by reducing the grid spacing from 200 nm to 3.1 nm in the r direction and 0.196 radians to 0.0015 radians in the ϕ direction. Fig.4.7 (b) and (c) are two cross section views of the Fig.4.7 (a). In Fig.4.7 (b) the quality factor and it’s relative error is plotted for different grid spacing in the azimuthal direction while the step size in the radial direction is kept at ∆r = 3.1nm. The relative error is obtained by choosing the Richardson extrapolated [85] value as the reference for the infinitesimal grid spacing. The black line is the least square fit to the error plot and is depicting a convergence rate of 0.9. Obviously, faster convergence rates is achievable by adopting higher order finite difference scheme (chapter 2). Note that the last two points are omitted for the line of best fit in Fig.4.7 (b). Similarly for the Fig.4.7 (c), the quality factor and it’s relative error is presented as a function of ∆r while ∆ϕ = 0.0015 radians. The line fitted to the error values (black line) shows a convergence rate of 2.8.
10−2 10−1 101 102 5 5.1 5.2 5.3 5.4 x 106 ∆φ (rad) ∆r (nm) Q (a) 10−3 10−2 10−1 4.9 4.95 5 5.05x 10 6 ∆φ(rad) Q 10−3 10−2 Relative Error relative error linear fit Quality factor (b) 10−2 10−1 5 5.1 5.2 5.3 5.4x 10 6 ∆r (µm) Q 10−6 10−4 10−2 100 Relative Error relative error linear fit Quality factor (c)
Figure 4.7: (a) Quality factor vs. grid spacings in the r and ϕ direction. (b) cross section view of (a) at ∆r = 3.1 nm and it’s relative error. (c) variations for the cross section at ∆ϕ = 0.0015 radians with the corresponding relative error.
