Band structure computations for dispersive photonic crystals

Department of Science and Technology

Institutionen för teknik och naturvetenskap

Linköpings universitet

Linköpings universitet

SE-601 74 Norrköping, Sweden

601 74 Norrköping

Examensarbete

LITH-ITN-ED-EX--07/015--SE

Band structure computations

for dispersive photonic

crystals

Fredrik Almén

Band structure computations

for dispersive photonic

crystals

Examensarbete utfört i ITN, Fysik och elektronik

vid Linköpings Tekniska Högskola, Campus

Norrköping

Fredrik Almén

Handledare Igor Zozoulenko

Examinator Igor Zozoulenko

Rapporttyp Report category Examensarbete B-uppsats C-uppsats D-uppsats _ ________________ Språk Language Svenska/Swedish Engelska/English _ ________________ Titel Title Författare Author Sammanfattning Abstract ISBN _____________________________________________________ ISRN _________________________________________________________________

Serietitel och serienummer ISSN

Title of series, numbering ___________________________________

Nyckelord

Keyword

Datum

Date

URL för elektronisk version

Avdelning, Institution

Division, Department

Institutionen för teknik och naturvetenskap Department of Science and Technology

2007-06-13

x

x

LITH-ITN-ED-EX--07/015--SE

Band structure computations for dispersive photonic crystals

Fredrik Almén

Photonic crystals are periodic structures that offers the possibility to control the propagation of light. The revised plane wave method has been implemented in order to compute band structures for photonic crystals. The main advantage of the revised plane wave method is that it can handle lossless dispersive materials. This can not be done with a conventional plane wave method. The computational challenge is comparable to the conventional plane wave method.

Band structures have been calculated for a square lattice of cylinders with different parameters. Both dispersive and non-dispersive materials have been studied as well as the influence of a surface roughness.

A small surface roughness does not affect the band structure, whereas larger inhomogeneities affect the higher bands by lowering their frequencies.

Detta dokument hålls tillgängligt på Internet – eller dess framtida ersättare –

under en längre tid från publiceringsdatum under förutsättning att inga

extra-ordinära omständigheter uppstår.

Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner,

skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för

ickekommersiell forskning och för undervisning. Överföring av upphovsrätten

vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av

dokumentet kräver upphovsmannens medgivande. För att garantera äktheten,

säkerheten och tillgängligheten finns det lösningar av teknisk och administrativ

art.

Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i

den omfattning som god sed kräver vid användning av dokumentet på ovan

beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan

form eller i sådant sammanhang som är kränkande för upphovsmannens litterära

eller konstnärliga anseende eller egenart.

För ytterligare information om Linköping University Electronic Press se

förlagets hemsida

http://www.ep.liu.se/

Copyright

The publishers will keep this document online on the Internet - or its possible

replacement - for a considerable time from the date of publication barring

exceptional circumstances.

The online availability of the document implies a permanent permission for

anyone to read, to download, to print out single copies for your own use and to

use it unchanged for any non-commercial research and educational purpose.

Subsequent transfers of copyright cannot revoke this permission. All other uses

of the document are conditional on the consent of the copyright owner. The

publisher has taken technical and administrative measures to assure authenticity,

security and accessibility.

According to intellectual property law the author has the right to be

mentioned when his/her work is accessed as described above and to be protected

against infringement.

For additional information about the Linköping University Electronic Press

and its procedures for publication and for assurance of document integrity,

please refer to its WWW home page:

http://www.ep.liu.se/

Band structure computations for dispersive photonic

crystals

Fredrik Almén

Supervisor: Aliaksandr Rahachou, Igor Zozoulenko Examiner: Igor Zozoulenko

Photonic crystals are periodic structures that oers the possibility to control the propagation of light.

The revised plane wave method has been implemented in order to compute band structures for photonic crystals. The main advantage of the revised plane wave method is that it can handle lossless dispersive materials. This can not be done with a conventional plane wave method. The computational challenge is comparable to the conventional plane wave method.

Band structures have been calculated for a square lattice of cylinders with dierent parameters. Both dispersive and non-dispersive materials have been studied as well as the inuence of a surface roughness.

A small surface roughness does not aect the band structure, whereas larger inhomogeneities aect the higher bands by lowering their frequencies.

Acknowledgments

I would like to thank Igor Zozoulenko for providing this interesting project and for being a great teacher. My life must have been a bit boring before electro-magnetism. A big thank you also goes to Aliaksandr Rahachou for helping me out in an invaluable way with coding, physics and other stu. I could not have made this project without your help and our discussions. In addition I would like to thank my computer for being loyal throughout all those hot CPU hours and not giving up on me.

This thesis concludes 19 uninterrupted years in the Swedish educational system. For that reason, I would to thank the somewhat square system itself, as well as all inspiring teachers (from Lindö, Djäkne, Ebersteinska and LiU). Some of you really contributed a lot to who I am today and also made me feel like the sky was the limit.

I use free software extensively and like it a lot. A thank you goes to RMS and the free software community in general, for spreading freedom and giving hope. Let's throw out those dirty licenses!

My friends and all people I know earn a thanks. I am sorry that I have not spend much time with you lately and that I often turn down social arrangements. I hope this will change in the future.

Finally I would like to thank my family for their innite love and support. Ad undas

Norrköping, June 25, 2007 Fredrik Almén

Contents

Abstract vi Acknowledgments vii 1 Introduction 1 1.1 Background . . . 1 1.2 Aim . . . 1 2 Photonic Crystals 3 2.1 History of photonic crystals . . . 3

2.2 Crystal structures in general . . . 4

2.2.1 Lattice and basis . . . 4

2.2.2 Symmetries . . . 5

2.2.3 Two- and three-dimensional lattices . . . 6

2.2.4 The reciprocal lattice . . . 6

2.2.5 The Brillouin zone . . . 8

2.3 Electromagnetism . . . 8

2.3.1 Basics of light . . . 9

2.3.2 Maxwell's equations . . . 10

2.3.3 Permittivity . . . 10

2.3.4 Permeability . . . 11

2.3.5 Maxwell's equations in dielectric media . . . 11

2.3.6 Dispersion . . . 14

3 Numerical computations 17 3.1 The plane wave expansion method . . . 17

3.2 The revised plane wave method . . . 18

3.2.1 Solving for dierent regions . . . 19

4 Results 23 4.1 Square lattice of alumina cylinders . . . 23

4.2 Square lattice of elliptic alumina structures . . . 24

4.3 Square lattice of polymer cylinders . . . 26

4.4 Square lattice of dispersive polymer cylinders . . . 28

4.5 Surface roughness . . . 29

5 Conclusions 35

Bibliography 37

A.2 Decopling of Maxwell's equations . . . 40 A.2.1 TM . . . 40 A.2.2 TE . . . 41 B Revised plane wave method 43 B.1 TE . . . 43 B.2 TM . . . 46

Chapter 1

Introduction

1.1 Background

During the stone, bronze and iron ages mankind learned to extract and make use of natural materials. Since that time there has been an enormous develop-ment. Today several materials, such as steel, concrete and plastics, are taken for granted.

A quite recent development has been the ability to control the electrical properties of certain materials. This was possible thanks to advances in semi-conductor and solid state physics. The resulting electronic revolution has left few unaected.

An advancement in the eld of photonics and optics has been the optical ber. The most useful application is perhaps the possibility to guide light over very long distances. This revolutionized telecommunications, today optical bers carry massive amounts of information all over the world.

In the late 1980's, a new possibility to not only guide light but also to control it was introduced. This periodic structure that oers control of light propagation is called a photonic crystal. Photonic crystals consist of periodically arranged dielectric media. This arrangement mimics crystalline structures of semiconductor materials, as illustrated in Fig. 1.1.

Similar to semiconductors, the photonic crystal might also have a band gap. This photonic band gap forbids electromagnetic waves of certain energies to propagate through the material at some certain directions. If the gap extends into all directions it is called a complete band gap.

1.2 Aim

The aim of this work is to implement a revised plane wave method. This method makes it possible to calculate band structures of photonic crystals built with dis-persive materials, i.e. materials whose refractive index changes with frequency. The second aim is to see if defects such as surface roughness aect the band structure.

(a) Traditional multilayer lm

(b) 2-D photonic crystal

(c) 3-D photonic crystal

Chapter 2

Photonic Crystals

2.1 History of photonic crystals

The term photonic crystal is most often applied to articial high contrast struc-tures with two- or three-dimensional periodicity.

Photonic structures with one-dimensional periodicity is nothing new. Layers of thin lms can be found in many dierent applications. Some examples are antireection coating, lters and laser applications, such as the vertical cavity surface emitting laser (VCSEL). [1]

The history of photonic crystals can be traced back two independent propos-als by Eli Yablonovitch and Sajeev John respectively [2, 3]. Yablonovitch wanted to inhibit spontaneous emission in order to increase the eciency in lasers, whereas John wanted to create a light localization. Yablonovitch suggested that this could be achieved if a three-dimensional periodic dielectric structure had a photonic band gap overlapping the electronic band edge.

During the course of two years, 21 dierent structures were fabricated. Fi-nally a face-centered-cubic structure was identied to have a complete photonic band gap. [4]

Up till this time photonic band structures were calculated with a variety of methods. One of them used a scalar wave approximation. This approximation was traditionally employed for calculations of the electronic band structure. De-spite photons and electrons dier much in how they interact, their propagation is described by quite similar equations.

After some time the methods were adopted for photons and theorists intro-duced a fully vectorial plane wave method for Maxwell's equations [5, 6, 7, 8]. The results had a good agreement with the earlier experimental work, but in-dicated that the band gap was not complete. Only a pseudo gap existed. In addition, one of the papers [7] predicted that the diamond structure would have a complete band gap.

This nally resulted in that a successful photonic crystal with a complete band gap was fabricated [9]. This structure is a variant of the diamond structure and was named Yablonovite after its inventor.

Periodic photonic structures also exist in nature, both in ora and fauna, where they have a wide range of uses. They are for instance employed to cre-ate bright color and broad angle reectivity, reduced reectivity or increased

transparency.

An example of natural periodic photonic structures are light sensitive species of brittlestar. They use regular arrays of inorganic microstructures to collect light. This can act as a warning system for the presence of predators. Morpho rhetenor butteries have discrete multilayers of cuticle and air which produce their blue color and visibility of up to half a mile. Nipple arrays are used by some insects to reduce reectivity in their eyes in order to improve their visual system. The polychaete worm has setae which form a natural crystal ber with a partial photonic band gap. These are only a few examples of how nature uses periodic photonic structures. [10]

2.2 Crystal structures in general

Crystal structures have traditionally been studied in solid-state physics. Natural crystalline materials are made of atoms, ions or molecules.

Photonic crystals, on the other hand, are materials made out of dielectric media. Here, the periodicity has the scale of the wavelength of light. The general theories about periodic crystalline structures can however be applied to photonic crystals as well.

More detailed information about crystal structures can be found in a solid-state physics textbook [11], a general book about photonic crystals [12] or even a materials science book [13].

2.2.1 Lattice and basis

Crystals consist of identical building blocks placed in a regularly ordered and periodic pattern. A building block is called basis and can for instance be a single atom or a group of atoms. The regularly ordered and periodic pattern to which the basis is attached is called the lattice. Attaching the basis to the lattice results in a crystal structure, as illustrated in Fig. 2.1.

(a) Basis. (b) Lattice. (c) Crystal structure.

Figure 2.1: The basis, lattice and crystal structure. The basis (a) is attached to the lattice (b), thus forming the crystal structure (c).

A lattice is dened by lattice vectors. For a three-dimensional lattice they are represented as

2.2 Crystal structures in general 5 Here u1, u2, u3 represents arbitrary integers and a1, a2, a3 are the lattice

vectors. Fig. 2.2 illustrates lattice vectors in a two-dimensional lattice.

One important property is that the crystal should look the same for every integer multiple of these translation vectors.

(a) Lattice vectors. (b) Primitive lattice vec-tors.

Figure 2.2: Lattice vectors in an oblique lattice. The lattice vectors on the left (a) leave the crystal invariant under any translation. However, the lattice vectors on the right (b) are primitive, thus making it possible to translate any point in the lattice.

Lattice vectors which make it possible to translate any point in the lattice are called primitive lattice vectors. The region dened by the primitive lattice vectors is called a primitive lattice cell. The primitive cell is a minimum volume cell, containing all necessary information about the crystal structure. A crystal structure will form by repeatedly placing cells next to each other.

The primitive cell can be dened in other ways than by the primitive lattice vectors. Fig. 2.3 illustrates two primitive cells known as the Wigner-Seitz cell. The Wigner-Seitz cell is constructed in the following way: One rst draw lines from a chosen lattice point to all nearby lattice points. At the midpoint of these lines, one draw lines or planes which are normal to the original lines. The area enclosed by these normal lines or planes is the Wigner-Seitz cell.

(a) Hexagonal lattice. (b) Oblique lattice.

Figure 2.3: The Wigner-Seitz cell for two dierent lattices.

2.2.2 Symmetries

A crystal structure can exhibit several dierent symmetries. One has already been shown in the form of a lattice translation r = R + u1a1+ u2a2+ u3a3.

Discrete translational symmetry leaves the system invariant only for multi-ples of some basic step length. This basic step length is normally the primitive lattice vector, but other step lengths are used in special cases.

A system with continuous translational symmetry is invariant for all trans-lations in a certain direction. An example of continuous translational symmetry in all three directions is free space (Fig. 2.4a). An innite plane of glass stretch-ing out in the xy-plane (Fig. 2.4b) has continuous translational symmetry in two directions. A lattice of dielectric cylinders (Fig. 2.4c) is homogeneous in one direction and accordingly exhibits discrete translational symmetry in two directions.

(a) Free space. (b) An innite plane of glass. (c) Lattice of dielectric cylin-ders.

Figure 2.4: Examples of continuous and discrete translational symmetry. There also exist other symmetries which leave the system invariant. Ro-tational symmetry exists if a crystal is left invariant after a rotation. Mirror symmetry is specially interesting for two-dimensional crystals because it allows modes to be separated into two dierent polarizations.

2.2.3 Two- and three-dimensional lattices

Lattices can be categorized in dierent types. The lattice types dier geometri-cally and so does also the symmetry operations that leave the lattice invariant. There are ve types of two dimensional lattice types (shown in Fig. 2.5). Among these the square lattice is probably the simplest whereas the oblique lattice is the most general. Three-dimensional lattices are classied into seven dierent systems depending on the axes and angels. The dierent placing op-tions of atoms within each cell results in 14 dierent three-dimensional lattice types, summarized in table 2.1. The simple type has an atom attached to each lattice point. The body centered cubic type (bcc) resembles the simple type but has an additional atom in the center of the cell. The face centered cubic type (fcc) has one atom at each lattice point plus one at the center of each side or face of the cell.

2.2.4 The reciprocal lattice

The periodicity of the crystal structure which leaves it invariant for certain translations has an important consequence. Namely that the local physical properties are invariant and periodic as well. An example of this is the atomic potential in a semiconductor or the dielectric function in a photonic crystal.

2.2 Crystal structures in general 7

(a) Square lattice. (b) Rectangular lattice. (c) Hexagonal lattice.

(d) Centered rectangular

lattice. (e) Oblique lattice.

Figure 2.5: Dierent types of two-dimensional lattices. System Lengths & angles Structure types Cubic a = b = c, simple, bcc, α = β = γ = 90◦ fcc Tetragonal a = b 6= c, simple, bcc α = β = γ = 90◦ Orthorhombic a 6= b 6= c simple, bcc, α = β = γ = 90◦ fcc, base-centered Rhombohedral a = b = c simple α = β = γ 6= 90◦ Hexagonal a = b 6= c simple α = β = 90◦, γ = 120◦ Monoclinic a 6= b 6= c simple, bcc α = β = 90◦6= γ Triclinic a 6= b 6= c simple α 6= β 6= γ 6= 90◦

Table 2.1: Three-dimensional lattice types.

Periodic functions are suitable for Fourier analysis, and can be expanded as x(r) =R xkei(k·r)dk

= x(r + R) =R xkei(k·r)ei(k·R)dk,

(2.2) where x(r) is an arbitrary periodic function and xk is the coecient for a plane

Only plane waves that fulll exp(i(k · R)) = 1 need to be included. This is thanks to the property that xk= xkexp(i(k · R)).

The wave vector of these plane waves, usually denoted by G, are the recip-rocal lattice vectors,

G = v1b1+ v2b2+ v3b3. (2.3)

The outcome is that the periodic function can be built as a summation over all reciprocal lattice vectors

x(r) =X

G

xGei(G·r). (2.4)

The reciprocal lattice vectors has the property bi· aj= 2πδij, which

corre-sponds to the conditions exp(i(G · r)) = 1 and G · r = n2π. It is usual that they are constructed in the following way:

b1= 2π a2× a3 a1· (a2× a3) b2= 2π a3× a1 a2· (a3× a1) b3= 2π a1× a2 a3· (a1× a2) . (2.5) If the translation vectors a1, a2and a3are primitive, then b1, b2and b3are

also primitive vectors of the reciprocal lattice.

2.2.5 The Brillouin zone

Just as in the real space lattice, a primitive cell can be dened in reciprocal space. The Wigner-Seitz cell in reciprocal space is called the Brillouin zone (see Fig. 2.6 and section 2.2.1 for construction of the Wigner-Seitz cell).

The Brillouin zone is important because it contains the complete set of wave vectors that can be diracted by the crystal. Every wave vector outside this primitive cell can be reduced back to the Brillouin zone, like k = k + G. The smallest region within the zone that is not related by symmetry is called the rst irreducible Brillouin zone.

2.3 Electromagnetism

What is light? Throughout history there have been numerous theories and answers to this question.

Up till the 17th century light was believed to be a white substance. Colors resulted when light was modied by medium. Sir Isaac Newton experimented with prisms and discovered that white light is a mixture of all colors.

The two dominating theories since that time have been that light is ei-ther particles or waves. The phenomenon that light produced could so far be explained by both theories. René Descartes and Newton considered light to consist of particles, Christiaan Huygens on the other hand proposed that light is a disturbance or pulse in the ether.

Newton's inuence made the particle theory dominating until the 19th cen-tury. Thomas Young's double slit experiment in 1802 could however not be explained if light was seen as particles. This together with the following work done by Jean Augustin Fresnel made the wave nature of light generally accepted. In 1861, James Clerk Maxwell made a modication to Ampère's law which would prove to be very important. Maxwell predicted that electromagnetic

2.3 Electromagnetism 9

(a) The reciprocal lattice. (b) The rst Brillouin zone.

Figure 2.6: The reciprocal lattice and the Brillouin zone for a two-dimensional square lattice. (a) The rst Brillouin zone is dened in the reciprocal lattice as a Wigner-Seitz cell. (b) Close-up of the rst Brillouin zone. The shaded area, enclosed by the high symmetry points Γ − X − M − Γ, is the rst irreducible Brillouin zone.

waves existed and that they would travel at the speed of light. This was veried in 1887 by Heinrich Hertz.

The outcome of Maxwell's work would not only be able to describe all elec-tromagnetic phenomena and show that light is an elecelec-tromagnetic wave, but also join the elds of optics and electromagnetism.

Today, with the concept of wave-particle dualism, light is considered both a wave and a particle (photon). [14]

2.3.1 Basics of light

When light enters a medium both its speed and wavelength changes such that v = c n λmed= v cλvac= λvac n . (2.6) Here c is the speed of light in vacuum, v is the speed of light in the medium and n is the refractive index. λvac and λmedis the wavelength in vacuum and

in the medium respectively. The frequency, f, is given by

fvac= c λvac , fmed= v λmed , (2.7) and in turn the energy is related to the frequency as

E = hf = ~ω, (2.8) where h is Planck's constant and ~ is the reduced Planck's constant (Dirac's constant). ω is the angular frequency and is related to the frequency as ω = 2πf. [15]

2.3.2 Maxwell's equations

Maxwell's equations describe the properties of electric and magnetic elds, and their interaction with matter. They can be expressed in both integral and dierential form, the latter reads as

∇ · E = ρ ε (2.9) ∇ · B = 0 (2.10) ∇ × E = −∂B ∂t (2.11) ∇ × H = J + ∂ ∂tD. (2.12) Here the symbols and their meanings are given in table 2.2. The rst equa-tion (2.9) is known as Gauss' law and describes how charges produce and electric eld. Equation (2.10) states that there are no magnetic charges, this is known as magnetic ux conservation. Faraday's law is expressed in equation (2.11). It can be interpreted as that a change of the magnetic eld produces an electric eld. Equation (2.12) is the Generalized Ampère's law (also known as Maxwell-Ampère's law). It states that a magnetic eld is produced by currents and a changing electric eld.

Symbol Meaning SI Unit E Electric eld intensity V /m D Electric displacement vector C/m2 H Magnetic eld A/m B Magnetic induction vector T J Current density A/m2

ρ Free electric charge density C/m2 ε Permittivity F/m

Table 2.2: Symbols in Maxwell's equations.

2.3.3 Permittivity

Provided that the material is linear and isotropic, the polarization vector, P , is P = ε0χeE , (2.13)

where ε0is the permittivity of free space and χeis the electric susceptibility. The

displacement vector is related to the electric eld and the polarization vector such that

D = ε0E + P = ε0(1 + χe)

| {z }

εr

E, (2.14) here εr is the relative permittivity of the material. This can be formulated in a

simpler form as

D = ε0εr(r)E = ε(r)E , (2.15)

such that D and E are related by the permittivity of the material. In this case a scalar dielectric function. [16]

2.3 Electromagnetism 11

2.3.4 Permeability

In a similar fashion as earlier, there is a relation between the magnetic induction vector and the magnetic eld. The magnetization vector, M, is given by

M = χmH . (2.16)

Here χm is the magnetic susceptibility of the material. Using this result, the

magnetic induction is equal to

B = µ0(H + M ) = µ0(1 + χm)

| {z }

µr

H , (2.17) where µ0 is the permeability of free space and µr is the relative permeability of

the material. The relation can be simplied into

B = µ0µr(r)H = µ(r)H , (2.18)

such that B and H are related by the permability of the material.

For most dielectric materials µr is close to unity such that B and H are

related by the permeability of free space. [16]

2.3.5 Maxwell's equations in dielectric media

In a dielectric material it can be assumed that there are no free charges and no currents, i.e. ρ = 0, J = 0. Adopting these assumptions and the relations introduced in section 2.3.3 and 2.3.4 results in

∇ · D = 0 (2.19) ∇ · B = 0 (2.20) ∇ × E = −∂B ∂t (2.21) ∇ × H = ∂D ∂t . (2.22) Time dependence

Electromagnetic elds are time dependent. The time dependence can be intro-duced as a harmonic mode times a complex exponential

H(r, t) = H(r)eiωt (2.23) E(r, t) = E(r)eiωt. (2.24) Electromagnetic waves

Maxwell's equations can be formulated as

∇2E = −εµω2E (2.25) ∇2_{H = −εµω}2_{H, (2.26)}

(For details, see appendix A.1). That is, two dierential equations, one for the electric eld and one for the magnetic eld.

It is very important to notice that these equations have the same form as the wave equation

∇2_{S = −}ω2

v2S. (2.27)

Making a comparison it is easy to notice that ω2

v2 → εµω

2_{, (2.28)}

which can be reduced to

v =√1

µε (2.29) Considering vacuum, ε and µ are equal to ε0and µ0respectively. This yields

c = √1 µ0ε0

= 3 · 108m

s , (2.30) which happens to be the speed of light. For dielectrics the corresponding case is v =√1 µε = c √ µrεr . (2.31) This implies that the velocity of electromagnetic waves decreases in materials, as seen in section 2.3.1.

Plane waves

Plane wave expansion is a convenient way to represent waves. A regular plane wave can be written as

s(r, t) = Aei(ωt+k·r), (2.32) or in a time-independent form as

s(r) = Aei(k·r). (2.33) A real wave vector will represent a propagating wave. A complex wave vector on the other hand represents a evanescent wave. This is illustrated in Fig. 2.7. TM and TE modes

For a two-dimensional case, Maxwell's equations can be decoupled into two sets of equations (see appendix A.2 for details). One set for transverse electric modes (TE) and one for transverse magnetic modes (TM).

TE modes have no electric eld in the direction of propagation, whereas TM modes have no magnetic eld in the direction of propagation.

That is, if the electromagnetic wave is traveling in the z-direction, then TM modes have Hx-, Hy- and Ez-components, whereas TE modes have Ex-, Ey

2.3 Electromagnetism 13

(a) Plane wave.

(b) Evanescent plane wave.

Figure 2.7: The propagating plane wave (a) has a real wave vector, whereas the evanescent wave (b) has a complex wave vector.

This makes it possible to simplify Maxwell's equations into three equations for each polarization. For TM the decoupled equations are

∂ ∂yEz(r) = −iωµ0Hx(r) (2.34) ∂ ∂xEz(r) = iωµ0Hy(r) (2.35) ∂ ∂xHy(r) − ∂

∂yHx(r) = iωε0εr(r, ω)Ez(r) , (2.36) and in a similar manner for TE they are

∂

∂yHz(r) = iωε0εr(r, ω)Ex(r) (2.37) − ∂ ∂xHz(r) = iωε0εr(r, ω)Ey(r) (2.38) ∂ ∂yEx(r) − ∂ ∂xEy(r) = iωµ0Hz(r) . (2.39)

2.3.6 Dispersion

If the refractive index in a material changes for dierent wavelengths, the mate-rial is said to be dispersive. The most famous or well known example is probably what happens in a prism. Each wavelength, or color, will have dierent refrac-tive index, thus resulting in a separation of colors. Fig. 2.8 illustrates how the dielectric function of silver and a polymer depends on the frequency.

(a) Silver.

(b) Representative polymer.

Figure 2.8: Dispersion in materials. Solid line (red) represents the real part of the dielectric function. Dashed line (blue) represents the imaginary part of the dielectric function.

In a non-absorbing material both the dielectric function and the refractive index are real. However, in a dispersive or absorbing material they will be complex such that

ε = ε0+ iε00 n = n + ik .˜ (2.40) Here the real part of the complex dielectric function represents diractive prop-erties whereas the imaginary part is related to losses.

2.3 Electromagnetism 15 The dielectric function has a direct relation to the index of refraction such

that

ε = ˜n2, (2.41) which yields

ε0 = n2− k2 _ε00_{= 2nk . (2.42)}

There also exists a relation between the real and imaginary parts of the dielectric function. In 1926-1927 Kramers and Kronig showed that the real part could be expressed as an integral of the imaginary part and vice versa. These kind of relations are today known as Kramers-Kronig relations.

The real part of the dielectric function can be obtained as ε0(ω) = 1 +2 πP Z ∞ 0 ε00_(ω0_)ω0 ω02_{− ω}2dω 0 _{, (2.43)}

whereas the imaginary part is obtained in a similar way ε00(ω) = −2ω π P Z ∞ 0 ε0_(ω0_{) − 1} ω02_{− ω}2 dω 0 _{. (2.44)}

Here P is the principle part of the Cauchy integral.

The reason behind this relation is that both dispersion and absorption orig-inate from the same process. This process is the excitation of dipoles.

If the dipoles can follow the eld there will be no absorption and no disper-sion. However, there might be frequency ranges where the dipoles will not be able to follow the eld entirely, but will naturally still try do so. This will aect the polarization vector (thus also the dielectric function) and the dipoles will absorb energy from the eld.

As a summary of the Kramers-Kronig relation it can be concluded that dispersion and absorption are physically linked together. If there is dispersion there will be absorption, and vice versa. [17, 18]

Chapter 3

Numerical computations

3.1 The plane wave expansion method

The plane wave expansion method (PWM) is the most popular technique for calculating photonic band structures. It results in an eigenvalue problem where eigenfrequencies are solved for a given wave vector. For a three-dimensional case the problem to be solved [12] is

∇ ×  ₁ ε(r)∇ × H(r)  =ω c 2 H(r) . (3.1) In a two dimensional case the decoupled set of equations for each polarization is formulated as an eigenvalue problem [19]. Starting with TM modes, the Hx

-and Hy-components are eliminated from the last equation in the decoupled set,

which yields 1 ε(r)  ∂2 ∂x2 + ∂2 ∂y2  Ez(r) = − ω2 c2Ez(r) . (3.2)

In a similar manner for TE modes, the Ex- and Ey-components are

elimi-nated leading to  ∂ ∂x  ₁ ε(r) ∂ ∂x  + ∂ ∂y 1 ε(r) ∂ ∂y  Hz(r) = − ω2 c2Hz(r) . (3.3)

As previously discussed in section 2.2.4, the properties for periodic structures make it possible to expand the dielectric function, ε(r), into a Fourier series of the form

ε(r) =X

G

κ(G)ei(G·r) (3.4) Because of the periodicity, the electric and magnetics elds are expanded according to Bloch's theorem as

H(r) =X G H(G)ei(k+G)·r (3.5) E(r) =X G E(G)ei(k+G)·r . (3.6)

This is called Bloch form and can interpreted as a plane wave modulated by a periodic function.

Inserting these expansions will yield the following eigenvalue problem for TM and TE respectively. X G00 κ(G − G00)|k + G00|2_E(G00_{) =}ω2 c2E(G) (3.7) X G00 κ(G − G00)(k + G00)(k + G)H(G00) = ω 2 c2H(G). (3.8)

The eigenvalue problem in equation (3.7) and (3.8) will be of innite size. It must therefore be truncated to include a nite number of plane waves, resulting in an eigenvalue problem with size N × N.

The band structure is computed by solving a frequency for a given wave vector. In addition, the eigenvalue problems will be hermitian, yielding real eigenvalues.

A consequence of using this method is that the material has to be frequency independent. Dispersive materials are frequency dependent, which means that the frequency has to be known in order to determine the dielectric function. The frequency is unknown in the PWM because it is the output of the method.

3.2 The revised plane wave method

In order to avoid the limitations of the PWM, a revised plane wave method (RPWM) [20] can be used. The RPWM was introduced in early 2005. It was rst formulated for two-dimensional photonic crystals, but since then it has been expanded to three-dimensional cases as well [21].

The RPWM is basically a modied formulation of the regular PWM. The new formulation results in an eigenvalue problem where a wave vector compo-nent is solved for a given frequency. A main consequence, and advantage of this, is that dispersive materials can be used.

The expansions of the dielectric function and the harmonic modes are in-serted directly into the decoupled equations (see appendix B for details). This yields three equations for each polarization, for TM these are

(ky+ [Gy]) [ez] = −k0[hx] (3.9)

(kx+ [Gx]) [ez] = k0[hy] (3.10)

(ky+ [Gy]) [hx] − (kx+ [Gx]) [hy] = −k0[κ] [ez] , (3.11)

and similarly for TE

(ky+ [Gy]) [hz] = k0[κ] [ex] (3.12)

(kx+ [Gx]) [hz] = −k0[κ] [ey] (3.13)

(ky+ [Gy]) [ex] − (kx+ [Gx]) [ey] = k0[hz] . (3.14)

Using these equations two eigenvalue problems can be formulated for each polarization. The most important dierence between these two formulations is that the eigenvalue either will be the x- or y-component of the wave vector. For TM, the eigenvalue problem is achieved by using equation (3.9) or (3.10) to

3.2 The revised plane wave method 19 eliminate either [hx]or [hy]in equation (3.11). For TE this is done in a similar

way.

The resulting eigenvalue problems for TM are 1 k0   −k0[Gx] k02[κ] − (ky+ [Gy])2 k0 −k0[Gx]     [hy] [ez]  = kx   [hy] [ez]   (3.15) − 1 k0   k0[Gy] k20[κ] − (kx+ [Gx]) 2 k2 0 k0[Gy]     [hx] [ez]  = ky   [hx] [ez]   (3.16) The equations for TE have a similar form.

− 1 k0   k0[Gx] k02− (ky+ [Gy]) [κ]−1(ky+ [Gy]) k2 0[κ] k0[Gx]     [ey] [hz]  = kx   [ey] [hz]   (3.17) 1 k0   −k0[Gy] k20− (kx+ [Gx]) [κ]−1(kx+ Gx) k2 0[κ] −k0[Gy]     [ex] [hz]  = ky   [ex] [hz]  . (3.18) Here [ex,y,z] and [hx,y,z] are eigenvectors of size N and k0 = ω/c. [Gx]

and [Gy] are diagonal matrices with size N × N. The diagonal elements are

components of the reciprocal lattice vectors. [Gx]has form of a staircase while

[Gy]has a sawtooth-looking shape, as seen in Fig. 3.1.

[κ] contains the Fourier coecients of the dielectric function. They can be calculated eciently using the two-dimensional FFT or with

κ(G − G00) = 1 M X xi,yj ε (xi, yj) e−i{(Gx−G 00 x)xi+(Gy−G00y)yj}, (3.19)

where M is the size of the discretization of the cell.

In contrast to the PWM, the resulting eigenvalue problem has size 2N × 2N and is neither hermitian nor symmetric. This results in complex eigenvalues. As seen in section 2.3.5, a plane wave with a complex wave vector represents an evanescent wave. Thus, among the resulting 2N eigenvalues only those with a zero imaginary part represents propagating states.

The eigenvalue problems also reveal one important limitation, namely that the material has to be lossless. A lossy material would imply complex wave vectors, and there is no way to know what the imaginary part of the provided wave vector component should be. The consequence is that the material can be frequency dependent but not absorptive. This is contradictive to the nature of dispersion, as discussed earlier in chapter 2.3.6.

Despite this, the RPWM can be used for estimation of photonic crystals operated in frequency regions where absorptance of the material is low.

3.2.1 Solving for dierent regions

The approach to compute the band structure is quite straightforward in the ΓX and XM regions (see Fig. 2.6 for the dierent regions).

(a) X-components of reciprocal lattice vector.

(b) Y-components of reciprocal lattice vector.

Figure 3.1: Shape of matrices [Gx](a) and [Gy](b) containing components of

the reciprocal lattice vectors.

In the Γ−X region it is suitable to use equation (3.15) and (3.17) to solve kx

while ky is provided. In a similar manner for the X − M region, it is convenient

to use equation (3.16) and (3.18).

The computation of the MΓ region can be both dicult and computing intensive. The reason is that both the x- and y-component of the wave vector varies.

One solution is to solve a wave vector component for a given frequency and a provided wave vector component. Out of the resulting eigenvalues, only those with real part equal to the provided wave vector component should be considered (i.e. kx= ky). The eigenvalue problem will need to be solved for a number of

dierent provided wave vector components at each frequency point. The result is a more time consuming problem to solve.

Another solution is to make a coordinate change so that the new x-coordinate, x0, is in the MΓ direction (see Fig. 3.2). The cell needs to be rebuilt and the lattice vectors need to be adopted to reect the new cell. Thereafter the problem can be solved in the same way as in the ΓX region.

3.2 The revised plane wave method 21

(a) The original cell. (b) The cell after change of coor-dinates.

(c) Lattice vectors before coordinate

change. (d) Lattice vectors after coordinatechange.

Figure 3.2: Square lattice, before and after coordinate change.

original primitive cell. This will aect the reciprocal lattice as well, the resulting primitive cell in the reciprocal lattice will be twice as small.

Before the coordinate shift, the MΓ region extends from the origin to the upper right corner of the Brillouin zone. The solved wave vectors will terminate on this line such that 0 ≤ k ≤√2π/a.

In the coordinate shifted crystal structure, the original MΓ region is now located along the ΓX region. The solved wave vectors terminate on this line such that 0 ≤ k ≤√2π/2a. That is, the new region is twice as small, and so are also the solved wave vectors. A consequence of this is a folded band structure, containing double information. The resulting data will need post-processing in form of a band unfolding.

In this newly dened Brillouin zone, the bands will start at the Γ-point and then be reected at the X-point. The eigenvalues of this reected part need to be changed into G0

x− kx0.

The unfolding starts by selecting starting points from the data given by the computation of the X − M region. The selected bands are followed to the X-point and from now on the reected parts are processed.

Chapter 4

Results

The results were obtained by developing a computer program. All programming was done in Fortran 90. LAPACK[22] was utilized for solving of eigenvalue prob-lems and FFTW[23] was used to calculate the coecients of the dielectric func-tion. Plotting, band unfolding and other post-processing of the acquired band structure data was performed in Grace. Fields were visualized with Gnuplot.

4.1 Square lattice of alumina cylinders

The rst result (Fig. 4.1) show the computed band structure for a simple square lattice of dielectric cylinders embedded in air.

(a) TM. (b) TE.

Figure 4.1: Band structure for square lattice of dielectric cylinders. Each cylinder has dielectric function, εr = 8.9, and radius, r = 0.2a. This

result repeats the results shown in [12]. It also indicates that the method works as expected.

It is interesting to see that there is a complete band gap for TM modes, whereas there is only a partial band gap for TE modes.

In addition, the band structure for TM modes exhibits an ultra-at band in the ΓX region. To detect these ultra-at bands, the frequency step size must be suciently small.

Fig. 4.2 show the electric eld for TM modes. Fields are displayed for bands directly beneath and above the band gap at the X and M high-symmetry points. Each band is labeled with an integer, n. The numbering starts with the lowest band and continues in ascending order. In a similar manner, the notation Xn refers to the nth band at the X-point.

The eld for the rst band is concentrated within the dielectric cylinder. A low frequency band concentrates its energy within the high-ε regions. Because of this it is often called dielectric band.

The second band, located precisely above the band gap, has a eld mainly distributed outside the dielectric structure. A high frequency band tend to concentrate its energy in low-ε regions. For this reason it is called air band.

(a) X1. (b) X2.

(c) M1. (d) M2.

Figure 4.2: Electric eld for square lattice of dielectric cylinders (TM modes).

4.2 Square lattice of elliptic alumina structures

This result (Fig. 4.3) show the computed band structure for a square lattice of elliptic structures embedded in air.

The ratio between the semimajor and the semiminor axis is 1.5 and the radius of the semimajor axis is, r = 0.2a.

It is interesting to see that there are two complete band gaps for TM modes, whereas there is only a partial band gap for TE modes.

Fig. 4.4 show the electric eld for TM modes at the X-point.

The rst band, located directly beneath the rst band gap, has its eld concentrated within the dielectric structure.

4.2 Square lattice of elliptic alumina structures 25

(a) TM. (b) TE.

Figure 4.3: Band structure for square lattice of elliptic structures.

For the second band, located between the two band gaps, the eld is mainly distributed outside the dielectric structure. The band with highest energy, the third band, is located above the second band gap. In this case the eld is strongly concentrated outside the dielectric structure.

(a) X1. (b) X2.

(c) X3.

Figure 4.4: Electric eld at X-point, for square lattice of elliptic structures (TM modes).

Fig. 4.5 show the electric eld for TM modes at the M-point. In the same way as for the X-point, the low frequency band has its eld concentrated in

the dielectric structure. Bands with higher frequencies are on the other hand mainly distributed in the low-ε regions.

(a) M1. (b) M2.

(c) M3. (d) M4.

Figure 4.5: Electric eld at M-point, for square lattice of elliptic structures (TM modes).

4.3 Square lattice of polymer cylinders

In this case, the material in the cylinders has been exchanged to a polymer (εr = 2.25). The reason for using polymer materials is that they are cheaper

and that the fabrication process is simpler. The results, shown in Fig. 4.6, also include several dierent radii of the cylinder.

4.3 Square lattice of polymer cylinders 27

(a) TM, r = 0.15a. (b) TE, r = 0.15a.

(c) TM, r = 0.2a. (d) TE, r = 0.2a.

(e) TM, r = 0.3a. (f) TE, r = 0.3a.

(g) TM, r = 0.4a. (h) TE, r = 0.4a.

The band structures show that there is a pseudo band gap for TM modes. It is clear that the radius aects the band structure. The widest band gap seems to be achieved with r = 0.3a.

4.4 Square lattice of dispersive polymer cylinders

This section contains results for dispersive materials. It is not possible to achieve results for these materials with the original PWM.

The material used is a dispersive polymer with dielectric function as in Fig 2.8b. The results shown in Fig. 4.7 are for a regular square lattice of cylinders with dierent radii.

(a) TM, r = 0.3a (b) TE, r = 0.3a

(c) TM, r = 0.4a (d) TE, r = 0.4a

Figure 4.7: Band structure for square lattice of dispersive polymer cylinders.

When comparing these results with results for the non-dispersive polymer, it is clear that the frequency dependent dielectric function aects the band structure.

Compared to the ideal case, the bands have a lower frequency. From now on, this inuence is referred to as a compression.

This compression is specially visible for higher bands. The reason is that the dielectric function is larger for higher frequencies.

In addition, with dispersive materials, the coecients of the dielectric func-tion must be calculated for each frequency point.

4.5 Surface roughness 29

4.5 Surface roughness

A number of computations were done in order to see if a surface roughness could aect the band structure. This roughness was introduced as in Fig. 4.8.

(a) r = 0.2a, roughness 10nm. (b) r = 0.2a, roughness 20nm. (c) r = 0.2a, roughness 50nm. (d) r = 0.4a, roughness 10nm.

(e) r = 0.4a, roughness 20nm.

(f) r = 0.4a, roughness 50nm.

Figure 4.8: Cylinders with surface roughness.

Fig. 4.9 and Fig. 4.10 show the results for polymer and alumina respectively. The conclusion is that surface roughness does not aect the band structure. An explanation is that the magnitude of the roughness is very small compared to the wavelength. There are some very minor discrepancies from the ideal case, but it is not visible from these images.

(a) TM, roughness 10nm. (b) TE, roughness 10nm.

(c) TM, roughness 20nm. (d) TE, roughness 20nm.

(e) TM, roughness 50nm. (f) TE, roughness 50nm.

Figure 4.9: Band structure for square lattice of polymer cylinders with surface roughness. r = 0.4a, εr= 2.25.

4.5 Surface roughness 31

(a) TM, roughness 10nm. (b) TE, roughness 10nm.

(c) TM, roughness 20nm. (d) TE, roughness 20nm.

(e) TM, roughness 50nm. (f) TE, roughness 50nm.

Figure 4.10: Band structure for square lattice of cylinders with surface rough-ness. r = 0.2a, εr= 8.9.

In order to study aects of the correlation length of the inhomogeneities, band structures for two smooth structures, as in Fig. 4.11, were computed. The corresponding band structures are shown in Fig. 4.12.

(a) smooth 1 (b) smooth 2

Figure 4.11: Two smooth structures.

The band structure for smooth 1 appears to be very similar to the ideal cylinder (shown in Fig. 4.1). The dierence is a small compression of the band structure, which possibly originates from the larger volume of the rod. However, mode structures of this distorted cylinder is not sensitive to the exact shape of the rod.

To study this further, a heavily distorted structure smooth 2 is introduced. The band structure is still similar to smooth 1 and the ideal case, but the compression is more visible. It is interesting to see that the compression is larger for higher frequencies. An explanation is that the wavelength is closer to the size of the inhomogeneity for higher frequencies.

4.5 Surface roughness 33

(a) smooth1, TM (b) smooth1, TE

(c) smooth2, TM (d) smooth2, TE

Chapter 5

Conclusions

The RPWM is an alternative to the original PWM. An important limitation of the conventional PWM is that materials must be frequency independent. The main advantage of the RPWM is the ability to include frequency dependent materials. At the same time the computational problem is comparable to the PWM.

The RPWM has been used to calculate band structures for lossless dispersive materials but also to study aects of inhomogeneities. It can be concluded that small surface roughness, in principle, does not aect the band structure. Inhomogeneities with larger correlation length results in that the bands have a lower frequency compared to the ideal case. This change in frequency is larger when the wavelength is closer to the size of the inhomogeneity.

Possible future work in this eld could be investigations of dierent lattice types, structures and materials.

Bibliography

[1] Djafar K. Mynbaev and Lowell L. Scheiner. Fiber-Optic Communications Technology. Prentice Hall, 2001.

[2] Eli Yablonovitch. Inhibited Spontaneous Emission in Solid-State Physics and Electronics. Physical Review Letters, 58, 2059-2062, 1987

[3] Sajeev John. Strong Localization of Photons in Certain Disordered Dielec-tric Superlattices. Physical Review Letters, 58, 2486-2489, 1987

[4] Eli Yablonovitch and T. J. Gmitter. Photonic Band Structure: The Face-Centered-Cubic Case. Physical Review Letters, 63, 1950-1953, 1989 [5] K. M. Leung and Y. F. Liu. Full Vector Wave Calculation of Photonic

Band Structures in Face-Centered-Cubic Dielectric Media. Physical Review Letters, 65, 2646-2649, 1990

[6] Ze Zhang and Sashi Satpathy. Electromagnetic Wave Propagation in Pe-riodic Structures: Bloch Wave Solution of Maxwell's Equations. Physical Review Letters, 65, 2650-2653, 1990

[7] K. M. Ho, C. T. Chan and C. M. Soukoulis. Existence of a Photonic Band Gap in Periodic Dielectric Structures. Physical Review Letters, 65, 3152-3155, 1990

[8] M. Plihal, A. Shambrook, A. A. Maradudin and Ping Sheng. Two-dimensional photonic band structures. Optics Communications, 80, 199-204, 1991

[9] Eli Yablonovitch and T. J. Gmitter. Photonic Band Structure: The Face-Centered-Cubic Case Employing Nonspherical Atoms. Physical Review Let-ters, 67, 2295-2298, 1991

[10] Pete Vukusic and J. Roy Sambles. Photonic structures in biology. Nature, vol. 424, 852-855 & 680, 14 August 2003

[11] Charles Kittel. Introduction to Solid State Physics. John Wiley & Sons, eight edition, 2005.

[12] John D. Joannopoulos, Robert D. Meade and Joshua N. Winn. Photonic crystals: molding the ow of light. Princeton University Press, 1995. [13] James F. Shackelford. Introduction to Materials Science for Engineers.

[14] Harris Benson. University Physics. John Wiley & Sons, revised edition, 1995.

[15] Igor Zozoulenko. Lecture notes on photonics. LiU, 2006.

[16] Igor Zozoulenko. Lecture notes on electromagnetic eld theory. LiU, 2004. [17] George B. Arfken and Hans J. Weber. Mathematical models for physicists.

Academic Press, fourth edition, 1995.

[18] Hans Arwin. Thin Film Optics. Lecture notes. LiU, fth edition, 2000. [19] Kazuaki Sakoda. Optical properties of photonic crystals. Springer, 2001. [20] Shouyan Shi, Caihua Chen and Dennis W. Prather. Revised plane wave

method for dispersive material and its application to band structure calcu-lations of photonic crystal slabs. Applied Physics Letters, 86, 043104, 2005 [21] C S Feng, L M Mei, L Z Cai, X L Yang, S S Wei and P Li. A plane-wave-based approach for complex photonic band structure and its applications to semi-innite and nite systems. Journal of Physics D: Applied Physics, 39, 4316-4323, 2006

[22] E. Andersson, Z. Bai, C. Bischof, et al. LAPACK Users' Guide. Society for Industrial and Applied Mathematics, third edition, 1999

[23] Matteo Frigo and Steven G. Johnson. The Design and Implementation of FFTW3. Proceedings of the IEEE 93 (2), 216-231, 2005.

Appendix A

Maxwell's equations

A.1 Electromagnetic waves

Starting from Maxwell's equations

∇×E(r, t) + µ∂H

∂t = 0 (A.1) ∇×H(r, t) − ε∂E

∂t = 0. (A.2) Time dependance is introduced and the derivatives are evaluated

∇×E(r)eiωt_{+ µ}∂ ∂tHe iωt_{= 0 (A.3)} ∇×H(r)eiωt_{− ε}∂ ∂tEe iωt_{= 0, (A.4)}

with the result

∇×E + iωµH = 0 (A.5) ∇×H − iωεE = 0. (A.6) Taking the curl and using the property

∇ × ∇ × V = ∇ (∇ · V ) − ∇2_{V (A.7)} yields ∇ (∇ · E) | {z } =0 −∇2_{E + iωµ ∇ × H} | {z } iωεE = 0 (A.8) ∇ (∇ · H) | {z } =0 −∇2H − iωε ∇ × E | {z } −iωµH = 0. (A.9) Using properties of Gauss' law and magnetic ux conservation give the nal result

∇2_{E = −εµω}2_{E (A.10)}

A.2 Decopling of Maxwell's equations

∇ · H(r, t) = 0 ∇ × H(r, t) = ε0ε(r)_∂t∂ E(r, t)

∇ · ε0ε(r)E(r, t) = 0 ∇ × E(r, t) = −µ0_∂t∂H(r, t)

(A.12)

A.2.1 TM

Beginning with Faraday's law and the generalized Ampère's law, their left-hand side are ∇×E(r, t) = ˆi∂ ∂yEz− ∂ ∂zEy  + ˆj ∂ ∂zEx− ∂ ∂xEz  + ˆk ∂ ∂xEy− ∂ ∂yEx  (A.13) ∇×H(r, t) = ˆi∂ ∂yHz− ∂ ∂zHy  + ˆj ∂ ∂zHx− ∂ ∂xHz  + ˆk ∂ ∂xHy− ∂ ∂yHx  . (A.14) TM only have Hx- ,Hy- and Ez-components. This makes it possible to simplify

these two equations into ∂ ∂yEz(r, t) = −µ0 ∂ ∂tHx(r, t) (A.15) − ∂ ∂xEz(r, t) = −µ0 ∂ ∂tHy(r, t) (A.16) ∂ ∂xHy(r, t) − ∂ ∂yHx(r, t) = ε0ε(r) ∂ ∂tEz(r, t). (A.17) Then, time dependence is introduced

∂ ∂y  Ez(r)e−iωt  = −µ0 ∂ ∂t  Hx(r)e−iωt  (A.18) ∂ ∂x  Ez(r)e−iωt  = µ0 ∂ ∂t  Hy(r)e−iωt  (A.19) ∂ ∂x  Hy(r)e−iωt  − ∂ ∂y  Hx(r)e−iωt  = ε0ε(r) ∂ ∂t  Ez(r)e−iωt  , (A.20) and the derivative d/dt are evaluated. The exponential terms cancel resulting in the decoupled equations for transverse magnetic modes:

∂

∂yEz(r) = iωµ0Hx(r) (A.21) ∂

∂xEz(r) = −iωµ0Hy(r) (A.22) ∂

∂xHy(r) − ∂

A.2 Decopling of Maxwell's equations 41

A.2.2 TE

As with TM, begin with Faraday's law and the generalized Ampère's law. TE have Ex- ,Ey- and Hz-components, which yields

∂ ∂yHz(r, t) = ε0ε(r) ∂ ∂tEx(r, t) (A.24) ∂ ∂xHz(r, t) = −ε0ε(r) ∂ ∂tEy(r, t) (A.25) ∂ ∂xEy(r, t) − ∂ ∂yEx(r, t) = −µ0 ∂ ∂tHz(r, t). (A.26) The time dependence is introduced

∂ ∂y  Hz(r)e−iωt  = ε0ε(r) ∂ ∂t  Ex(r)e−iωt  (A.27) ∂ ∂x  Hz(r)e−iωt  = −ε0ε(r) ∂ ∂t  Ey(r)e−iωt  (A.28) ∂ ∂x  Ey(r)e−iωt  − ∂ ∂y  Ex(r)e−iωt  = −µ0 ∂ ∂t  Hz(r)e−iωt  (A.29) and the d/dt derivatives are evaluated. The exponential terms cancels resulting in the decoupled equations for TE:

∂

∂yHz(r) = −iωε0ε(r)Ex(r) (A.30) ∂

∂xHz(r) = iωε0ε(r)Ey(r) (A.31) ∂

∂xEy(r) − ∂

Appendix B

Revised plane wave method

Expansions

Expansion used for the revised plane wave method are H(r) = 1 Z0 X G H(G)ei(k+G)·r (B.1) E(r) =X G E(G)ei(k+G)·r (B.2) εr(r, ω) = X G κ(G)ei(G·r) (B.3) Here, Z0 is the impedance of free space,

Z0= r µ0 ε0 = µ0c = 1 ε0c . (B.4)

B.1 TE

Beginning with maxwell's decoupled equations for TE ∂

∂yHz(r) = iωε0εr(r, ω)Ex(r) (B.5) − ∂ ∂xHz(r) = iωε0εr(r, ω)Ey(r) (B.6) ∂ ∂yEx(r) − ∂ ∂xEy(r) = iωµ0Hz(r) (B.7) First equation

Inserting expansions into equation (B.5) and evaluating ∂/∂y yields 1

Z0

X

G

i(ky+ Gy)Hz(G)ei(k+G)·r

= iωε0 X G0 κ(G0)ei(G0·r)X G00 Ex(G00)ei(k+G 00_)·r . (B.8)

Here i cancels, 1

Z0 and ε0 are simplied into

ω c. X G (ky+ Gy)Hz(G)ei(k+G)·r =ω c X G0 κ(G0)ei(G0·r)X G00 Ex(G00)ei(k+G 00_)·r (B.9)

Using substitution G = G0_{+ G}00 _yields

X G (ky+ Gy)Hz(G)ei(k+G)·r = ω c X G X G00

κ(G − G00)ei(G·r)e−i(G00·r)Ex(G00)ei(k+G

00_)·r

. (B.10) Some of the exponential terms on the right hand side vanishes

X G (ky+ Gy)Hz(G)ei(k+G)·r =ω c X G X G00 κ(G − G00)Ex(G00)ei(k+G)·r, (B.11)

as well as two summations and the exponential terms. The result is (ky+ Gy)Hz(G) = ω c X G00 κ(G − G00)Ex(G00) (B.12) Second equation

Equation (B.5) is quite similar to equation (B.6). Inserting expansions into equation (B.6), evaluating ∂/∂x and performing similar operations yield

(kx+ Gx)Hz(G) = − ω c X G00 κ(G − G00)Ey(G00). (B.13) Third equation

Inserting expansions into equation (B.7) and and evaluating derivatives results in iX G [(ky+ Gy)Ex− (kx+ Gx)Ey] ei(k+G)·r = iωµ0 1 Z0 X G0 Hz(G0)ei(k+G 0_)·r , (B.14)

which simplies into

(ky+ Gy)Ex− (kx+ Gx)Ey=

ω

B.1 TE 45 Eigenvalue problem

After inserting the expansions we have the following equations, (ky+ Gy)Hz(G) = ω c X G00 κ(G − G00)Ex(G00) (B.16) (kx+ Gx)Hz(G) = − ω c X G00 κ(G − G00)Ey(G00) (B.17) (ky+ Gy)Ex− (kx+ Gx)Ey= ω cHz(G) , (B.18) which can be formulated as

(ky+ [Gy]) [hz] = k0[κ] [ex] (B.19)

(kx+ [Gx]) [hz] = −k0[κ] [ey] (B.20)

(ky+ [Gy]) [ex] − (kx+ [Gx]) [ey] = k0[hz] (B.21)

First rearrange equation (B.20) as [ey] = −

1 k0

[κ]−1(kx+ [Gx]) [hz] , (B.22)

then insert it into equation (B.21). This results in ky[ex] + [Gy] [ex] + (kx+ [Gx])  1 k0 [κ]−1(kx+ [Gx]) [hz]  = k0[hz] . (B.23)

ky[ex]is separated and the right hand side multiplied with k0/k0 yielding

ky[ex] = k₀2 k0 [hz] − k0 k0 [Gy] [ex] − 1 k0 (kx+ [Gx]) [κ]−1(kx+ [Gx]) [hz] , (B.24) and further as ky[ex] = 1 k0  −k0[Gy] [ex] + k20[hz] − (kx+ [Gx]) [κ] −1 (kx+ [Gx]) [hz]  . (B.25) Now, equation (B.19) can be rearranged into

ky[hz] = k0[κ] [ex] − [Gy] [hz] , (B.26) which is equal to ky[hz] = 1 k0 k02[κ] [ex] − k0[Gy] [hz]
. (B.27) At this state, equation (B.25) and (B.27) are formulated as an eigenvalue problem 1 k0   −k0[Gy] k02− (kx+ [Gx]) [κ]−1(kx+ Gx) k2 0[κ] −k0[Gy]     [ex] [hz]  = ky   [ex] [hz]   (B.28) An eigenvalue problem with kxas eigenvalue is achieved in a similar manner.

Equation (B.19) is rearranged and inserted into equation (B.21) followed by resembling operations.

B.2 TM

Maxwell equations

Starting from Maxwell's decoupled equations for transverse magnetic modes, ∂ ∂yEz(r) = −iωµ0Hx(r) (B.29) ∂ ∂xEz(r) = iωµ0Hy(r) (B.30) ∂ ∂xHy(r) − ∂

∂yHx(r) = iωε0εr(r, ω)Ez(r). (B.31) The expansions are inserted into these equations.

First equation

Inserting expansions into equation (B.29) and evaluating ∂/∂y yields X G (ky+ Gy)Ez(G)ei(k+G)·r = −ωµ0 1 Z0 X G0 Hx(G0)ei(k+G 0_)·r , (B.32) which is simplied into

(ky+ Gy)Ez(G) = −

ω

cHx(G). (B.33) Second equation

Equation (B.30) does not dier much from the rst. Inserting expansions yields (kx+ Gx)Ez(G) =

ω

cHy(G). (B.34) Third equation

Inserting expansions into equation (B.31) yields X G [(kx+ Gx)Hy(G) − (ky+ Gy)Hx(G)] ei(k+G)·r =ω c X G0 κ(G0)ei(G0·r)X G00 Ez(G00)ei(k+G 00_)·r (B.35)

Using the substitution G = G0_{+ G}00 _{results in}

(ky+ Gy)Hx− (kx+ Gx)Hy= − ω c X G00 κ(G − G00)Ez(G00). (B.36)

Rewriting as an eigenvalue problem

At this point the three equations are formulated as

(ky+ [Gy]) [ez] = −k0[hx] (B.37)

(kx+ [Gx]) [ez] = k0[hy] (B.38)

B.2 TM 47 Equation (B.38) is rearranged as [hy] = (kx+ [Gx]) [ez] k0 , (B.40) and then inserted into equation (B.39). This yields

ky[hx] + [Gy] [hx] − (kx+ [Gx])

 (kx+ [Gx]) [ez]

k0



= −k0[κ] [ez] , (B.41)

which can be written as − 1 k0 h k0[Gy] [hx] +  k2₀[κ] − (kx+ [Gx])2  [ez] i = ky[hx] . (B.42)

Now, equation (B.37) is rearranged into

− k0[hx] − [Gy] [ez] = ky[ez] , (B.43)

and then one step further resulting in − 1

k0

k2

0[hx] + k0[Gy] [ez] = ky[ez] . (B.44)

This result is combined with equation (B.42) into an eigenvalue problem − 1 k0   k0[Gy] k20[κ] − (kx+ [Gx]) 2 k02 k0[Gy]     [hx] [ez]  = ky   [hx] [ez]   . (B.45) An eigenvalue problem with kxas eigenvalue is achieved in a similar manner.

Equation (B.37) is rearranged and inserted into equation (B.39) followed by resembling operations.