Master Thesis
Topological plasmons in periodically
modulated graphene
by
Roel Sluijter
Institute of Physics
Faculteit der Natuurwetenschappen en Informatica
University of Amsterdam
June 17, 2020
Title: Topological plasmons in periodically modulated graphene Author: Roel Sluijter
© June 2020
Superviser/Examiner: Dr. Erik van Heumen Second reviewer: Dr. Jasper van Wezel Daily supervisor: Stephan Bron
Conducted between: February 2019 - January 2020 Location: Amsterdam, the Netherlands
Abstract
The study presented in this thesis forms a preliminary theoretical investigation to model plasmonic behavior in 2D structures, with the aim to prepare for experiments on topologically protected edge states in graphene. The structures are created by designing artificial lattices with carefully chosen discrete symmetries, which are imposed on the material by either making holes in, or placing gold islands on, a sheet of graphene. A finite difference time domain method is used to solve the equations of the underlying problems, which is benchmarked to a known problem with non-trivial topological behavior, which consists of a hexagonal structured lattice created by holes in graphene in the presence of an applied uniform magnetic field. Subsequently, custom structures with specific symmetries contained in the unit cell are investigated to determine how band behavior can be tuned and possibly be used to create topological phases without the requirement of an external magnetic field. The symmetries of the investigated unit cells are recreated by using gold islands, a SiO2substrate and a h-BN protective layer to model an experimentally more relevant setup and asses its feasibility.
Abstract ii
Contents iii
Acknowledgements v
Introduction 1
1 From Maxwell’s equations to plasmons 4
1.1 Maxwell’s equations in matter . . . 4
1.2 The Drude-Lorentz model . . . 6
1.3 The wave equation and volume plasmons . . . 8
1.4 Surface plasmon polaritons (SPP’s) . . . 9
2 Modeling plasmonic lattices 13 2.1 Bloch’s theorem and the concept of band folding . . . 13
2.2 Calculating band structures: The FDTD approach . . . 18
2.3 Benchmarking Jin et al. . . 22
3 Tailoring the energy bands 31 3.1 Symmetry breaking: The creation of band gaps . . . 32
3.2 Substrates: Experimental necessities . . . 34
3.3 Custom designs: The results . . . 36
Conclusion and outlook 45
Contents iv
A Chern numbers 48
B Template scripts 54
First of all, I would like to express my gratitude towards Erik van Heumen, who has been my supervisor and advisor during every stage of this project. He especially helped me with some personal problems and in doing so gave me the confidence I needed to graduate. He also formed my biggest source for in depth physics discussions and was quick to arrange outside help when needed. He also guided me during the writing of this thesis and his many invaluable comments made it a much better final product.
I would also like to thank Stephan Bron, who I mostly worked with while doing experiments with s-SNOM and AFM-IR. Probably the most fun I had was while setting up the new CO2laser and the satisfaction of getting everything to work as envisioned. He was also one of the people I could always come to with questions, which I did on many occasions, leading to insightful discussions.
Furthermore, I would like to thank Jasper van Wezel and Jans Henke, for aiding me with their theoretical expertise regarding topological systems and most importantly gave me confidence and reassurance about being on the right track when I needed it.
I would also like to thank Xuanbo Feng, who became my go-to person whenever I had urgent theoretical questions. Apart from all his useful knowledge, he possesses a mind that handles problems in a more abstract way than mine, which helped me out many times by giving me a fresh point of view.
I would also like to thank Philippe Corboz, who was nice enough to allow me to use the institute’s cluster computer "Obelix" and advised me on how to execute tasks in an efficient manner. Without this resource, most results presented in this thesis would not have been obtainable during the time frame of this project.
Lastly, aside from the people I have already mentioned, I am very grateful to the rest of the quantum matter group who made my time spent there very enjoyable. Steef, George, Huaqian, Marc, Floris, Florian, Lewis, Anne, Mark, Linda, Ayako, Kumi and anyone I might have forgotten to mention, thank you all for the good times.
Introduction
During the past decade, the field of solid state physics has gone through an exciting and revolutionary period with the introduction of a new concept, previously viewed as purely mathematical, called "topology" applied to periodic quantum systems.
The most established theory used to explain phase transitions in matter is known as Landau theory [1]. In accor-dance with this theory, distinctive phases of matter can mostly be classified by the principle of spontaneous symmetry breaking. This occurs when the Hamiltonian of a system is invariant under a certain symmetry, but gets broken by the state it occupies. For instance, the energy ground states associated with electrons inside a ferromagnet all require spin-alignment, but differ by orientation. The states all have an equal probability to be occupied, but because the system hasto occupy a state, the rotation-symmetry will inherently be broken [2, 3]. Similarly, crystalline solids can be asso-ciated with the breaking of translation-symmetry and ordinary superconductors with the breaking of gauge-symmetry. Interestingly, this does not apply to newly discovered topological phases [4], which brings a whole new perspective to the field.
In the earlier stages of research concerning topology in combination with solid state physics, the focus was pri-marily directed towards electronic systems. This is not surprising, as the first experimental observation of topology in physics took place in a 2-dimensional electron gas at the interface of a metal and a semiconductor inside a uniform magnetic field. The experiment was conducted in 1980 by German physicist Klaus von Klitzing [5], who unexpectedly observed the phenomenon we now know as the integer quantum Hall effect, where the energy levels of electrons take on discrete quantized values. The effect could be explained with quantum mechanics and so called Landau levels, but no-one realized its association with topology. A big breakthrough came from Michael Berry, who theorized that a wavefunction undergoing an adiabatic evolution gains an accumulated phase that finds its origin in geometry. This phase is now known as Berry’s phase [6] and is directly related to the so called Chern number, which is the relevant topological invariant in this system. It would take over two decades more before the importance of topology really became appreciated, leading to a new characterization paradigm based on topological order. Consequently, creativity among theoretical physicists lead to all kinds of new predictions. Out of these, arguably the most prominent exper-imentally validated new discovery is called the topological insulator. It can be seen as a material that behaves as an
insulator in the bulk, but a conductor at the outer edges. This is because the electrons are "topologically protected" by the bulk-edge correspondence [4, 7, 8].
The above mentioned discoveries all involve fermionic electron systems. To what extent the concept of topology applies to other fields in physics is slowly becoming clear and stretches far beyond. One might say we are on the doorstep of all kinds of great new discoveries in this area. The research presented in this thesis aims to continue the voyage by examining topological effects on systems comprised of oscillating electromagnetic waves in the presence of matter. More specifically, when these waves couple to free electric charge inside a material, particles called plasmons are formed. These are "quasi" particles that behave as bosons. Analogous to electrons in topological insulators, a system could be created in which plasmons can travel along the edge of a material while being topologically protected from scattering on defects. This would give them significantly longer lifetimes.
Theoretically, systems exhibiting such states are not restricted to a specific dimensionality. For instance, the "edge" of a 3D system is a 2D surface encompassing it and for a 2D system it consists of the line surrounding it. As a matter of fact, experiments have already shown that topological insulators can be created successfully in 2D as well as 3D [9, 10]. In our case, practical experimental feasibility and promising theoretical research, among other reasons, have turned the focus on 2D plasmonic systems. Indirectly, this choice was influenced by the physicists Konstantin Novoselov and Andre Geim, because they were able to produce a single layer of graphite with a thickness of only one atom [11]. This material is now known as graphene. Due to its extreme thinness, the material effectively forms a 2 dimensional system. Unlike the setup of Klaus von Klitzing, who needed the interface between two materials to create an effective 2D electron gas, graphene has this property all by itself. Additionally, because the material has a strong light-matter interaction, electromagnetic waves are also strongly confined to the material [12]. This makes graphene a very suitable material for the creation of 2D plasmonic crystals. Currently, a number of similar 2D materials are also available, but graphene distinguishes itself by having a large tuneable carrier density and exceptionally long intrinsic relaxation times e.g. it is highly conductive [13].
The study presented in this thesis forms a preliminary investigation to aid lab experiments using a scattering scanning near field optical microscope (s-SNOM). This kind of apparatus can measure opto-electronic properties with high spatial resolution, thereby breaking the so called diffraction limit. It typically operates in the infrared regime. More specifically, the currently available lasers at our lab operate at 28-32 THz and 50-58 THz. In this respect, graphene is also a very suitable material, as the combination of its electrical and optical properties makes it a perfect platform to engineer systems in which the most interesting plasmonic behavior occurs in these frequency regions.
In condensed matter physics, material properties are mostly identified as emerging phenomena resulting from a macroscopic system consisting of periodically repeating structures in space. The fundamental mathematical tools to describe such a system are found in Bloch’s theorem. In his original paper, Swiss-born American physicist Felix Bloch describes how electrons, represented by probability density waves, inside a periodic potential can lead to the
Introduction 3
formation of electronic band structures [14]. Interestingly, his theorem can be used for any other system that consists of a periodic potential. This fact plays a prominent role in meta-material science and forms the basis of the research presented in this paper. Similar to electronic bands, electromagnetic waves inside a periodic potential may result in the formation of plasmonic band structures.
Periodic potentials are formed by repeating structures. These structures could exist of single molecules or even larger formations. The smallest portion of a crystal lattice is called a unit cell. In our case the unit cells are significantly larger than molecules and are more conveniently defined by their underlying dielectric geometry. Essentially, by constructing these structures in a clever way, we hope to create new topological phases.
The propagation of electromagnetic waves in the periodic potentials proposed in this thesis can not be solved analytically and demand numerical methods. Consequently, most of the research presented in this thesis will involve a specific computational technique known as the finite difference time domain (FDTD) method, also known as Yee’s method, named after the Chinese American applied mathematician Kane S. Yee [15].
The exponential increase over time in computational power predicted by Moore’s law, which was based on the assumption that the amount of transistors on a chip would double every 2 years, has been lagging for some time, as the atomic size limit is approaching quickly. In order to continue the trend, new technologies have to be found and incorporated. Possibly, topologically protected plasmons in graphene could act as a means to transmit information on chips. They support much higher frequencies and are far less lossy than conventional wires, thereby increasing computational speeds immensely [16]. Alternativaly, plasmons in graphene could be engineered to work as a two-qubit quantum logic gate, which could form the basis of a quantum computer [17].
The main content of this thesis is presented in the form of three chapters. The purpose of the first chapter is to provide theoretical background relevant to the subject and aims to provide a better understanding of the nature of plasmons. The second chapter focuses on the use of the FDTD method and its ability to obtain the band structures of plasmons in a sheet of graphene containing holes in a hexagonal pattern. To this end, results found by Jin et al. [13] were recalculated, as a means of benchmarking the method. The third chapter is devoted to the investigation of custom designs and exploring how the symmetries of a unit cell affect the band structures. Additionally, it provides simulations of more experimentally feasible structures that involve specific substrates and gold islands to create periodicity.
From Maxwell’s equations to plasmons
The key constituent of a plasmon, along with other specific characteristics, can be considered to be an oscillation in the electromagnetic field coupled with charge, most commonly electrons. Therefore, the starting point to investigate its origin, is sought in Maxwell’s equations in the presence of matter. More specifically, a solid contains a certain charge distribution that can interact with electromagnetic waves and alter its propagation.
1.1
Maxwell’s equations in matter
In order to build a framework in which electromagnetic waves can be described for macroscopic systems, composed of matter, Maxwell’s equations can be written in differential form as [18]:
∇ · D= ρext (1.1) ∇ · B= 0 (1.2) ∇ × E= −∂B ∂t (1.3) ∇ × H= Jext+ ∂D ∂t (1.4)
These quantities describe the electromagnetic field inside a material and in this respect, D represents the Electric displacement, E the electric field, H the magnetic field, B the magnetic induction or magnetic flux density, ρext the external charge density and Jextthe external current density.
These Fields can further be linked by considering the internal polarization P, the electric dipole moment per unit
From Maxwell’s equations to plasmons 5
volume, and magnetization M:
D= 0E+ P (1.5)
H= 1 µ0
B − M (1.6)
Where the constant 0is the electric permittivity in vacuum and µ0the magnetic permeability in vacuum. The polar-ization is caused by the alignment of microscopic dipoles with the magnetic field and is related to the internal charge density by ∇ · P= −ρ. To make everything complete, charge conservation implies that ∇ · J = −∂ρ/∂t, so the internal charge and current densities can be linked via:
J=∂P
∂t (1.7)
For all cases studied in this thesis, the response of the system to a perturbation can be considered linearly proportional. Therefore, one can remain in the regime of linear optics and define the following constitutive relations:
D= 0E (1.8)
B= µ0µH (1.9)
where is the so called relative permittivity (dielectric function) and µ the relative permeability. In the case of non-magnetic media, the value of µ will be equal or close to unity. Therefore, the focus is now turned to the dielectric function , which in all generality can be said to depend on time and location. Recognizing this fact, equation (1.8) can be written as:
D(r, t)= 0 Z dt0dr0(r − r0, t − t0)E(r0, t0) (1.10) J(r, t)= Z dt0dr0σ(r − r0, t − t0)E(r0, t0) (1.11)
It is informative to take the Fourier transform with respect to the wavevector K and frequency ω defined by:
Z
dtdrei(K·r−ωt) (1.12)
In this domain the convolution turns into a simple product and the constitutive relations are:
D(K, ω)= 0(K, ω)E(K, ω) (1.13)
Using equations 1.5 and 1.7 and recognizing that in the Fourier domain ∂t∂ → −iω, this leads to the fundamental relationship between the dielectric function and the conductivity:
(K, ω) = 1 +iσ(K, ω) 0ω
(1.15)
In case of large wavelengths, e.g. significantly larger than the characteristic dimensions of the system like the unit cell size, the dielectric response can be simplified to the spatially local limit via (K= 0, ω) = (ω). Both and σ can then in general be considered as complex functions of the form:
(ω) = 1(ω)+ i2(ω) (1.16)
σ(ω) = σ1(ω)+ iσ2(ω) (1.17)
For experimental practicality it is useful to relate the dielectric function of a material to the complex refractive index:
˜n(ω)= n(ω) + iκ(ω) (1.18)
It is defined through ˜n ≡ √ and accordingly κ is named the extinction coefficient and determines the optical absorption of electromagnetic waves propagating through the medium. This is especially useful for recovering the dielectric function from reflectivity studies, as this definition explicitly yields:
1= n2−κ2 (1.19) 2= 2nκ (1.20) n2= 1 2 + 1 2 q 2 1 + 2 2 (1.21) κ = 2 2n (1.22)
Lastly, another useful observable is the reflectivity of a (semi-infinite) solid, which is related to the dielectric function by: R(ω)= 1 − √(ω) 1+√(ω) 2 (1.23)
1.2
The Drude-Lorentz model
Now let’s take a closer look at the possible nature of the dielectric function. In 1900, German physicist Paul Drude modeled the conductivity of metals by imagining it to resemble a gas of free electrons against a background of fixed
From Maxwell’s equations to plasmons 7
positive ion cores [19]. For a single electron subjected to an external electric field in this "plasma sea", the associated equation of motion can be written as [18]:
m¨x(t)+ mγ˙x(t) = −eE (1.24)
Where γ is a damping factor associated with the velocity of the electron. For an oscillating time dependent electric field E(t)= E0e−iωt, the general solution is:
x(t)= x0e−iωt (1.25)
Which yields the relationship:
x(t)= e
m(ω2+ iωγ)E(t) (1.26)
The macroscopic polarization then equals:
P= −nex(t) = − ne 2
m(ω2+ iωγ)E(t) (1.27)
Remembering the relationship:
D= 0E = 0E+ P (1.28)
The dielectric function can be identified as:
(ω) = 1 − ω 2 p ω2+ iωγ (1.29) where ω2p= ne 2
m0 is called the plasma frequency. In 1905, dutch physicist Hendrik Antoon Lorentz extended this model by including a restoring force term, associated with a spatial deviation from the rest position of the electron. The equation of motion that applies in this case is:
m¨x(t)+ mγ˙x(t) + ω20x(t)= −eE (1.30)
Following the same procedure as before, this leads to the following solution:
(ω) = 1 + ω 2 p ω2 0−ω2− iωγ (1.31)
Which reduces to the Drude result for ω0 = 0. After quantum physics became more established, a more accurate description was formulated by replacing the single ω0 term by a sum of terms, each belonging to a specific atomic transition frequency with corresponding damping parameters γ belonging to the excited state’s lifetime. The final
result is a model in which the complete dielectric function can be described in the general form [20]: (ω) = ∞− ω2 p ω2+ iωγ+ X i σ2 i ω2 i −ω2− iωγ (1.32)
In which ∞ is the so called instantaneous dielectric, the effective response at "infinite" frequencies. This term is generally used to model a dielectric function in a specific frequency range, where left over tails from very high frequency oscillators can be approximated with a constant.
1.3
The wave equation and volume plasmons
To examine traveling-wave solutions of Maxwell’s equations, it is useful to look at the cross product terms with zero external current [18]:
∇ × E= −∂B
∂t (1.33)
∇ × H= ∂D
∂t (1.34)
Taking the cross-product on both sides of the first equation gives:
∇ × (∇ × E)= −∂∇ × B
∂t (1.35)
By using equation 1.34, the vector identity ∇ × (∇ × E)= ∇(∇·E)−∇2E and taking the Fourier transform, its frequency domain equivalence is identified as:
K(K · E) − K2E= −(K, ω)ω 2
c2E (1.36)
Where c = õ1
00 is the speed of light in vacuum. Solutions to this equation can be split up into transverse and longitudinal waves.
For longitudinal waves K(K · E)= K2E, which can only be accommodated for when:
(K, ω) = 0 (1.37)
In the other case, the transverse waves oscillate perpendicular to their propagation direction, so K · E= 0, which yields the following general dispersion relation:
K2= (K, ω)ω 2
c2 (1.38)
From Maxwell’s equations to plasmons 9
properties of electromagnetic waves in metals. Splitting this equation into its real and imaginary constituents gives:
r = 1 − ωp2τ2 1+ ω2τ2 (1.39) i= ωp2τ ω(1 + ω2τ2) (1.40)
where τ= 1/γ. In the frequency region with negligible damping, the product ωτ >> 1 and the dielectric function is predominantly real. Thus, for an undamped free electron plasma, the dielectric function can be written as:
1 −ω 2 p
ω2 (1.41)
Note that this approximation is only valid for negligible Drude damping and excluding any interband transitions, which could increase isignificantly. However, qualitatively it serves as a nice way to illustrate the classification of plasmonic waves.
Longitudinal waves require the dielectric function to be equal to zero, which is the case for ω= ωp. This solution represents a longitudinal charge density wave oscillating at the plasma frequency in a sea of electrons.
Using this same dielectric function in equation 1.38 gives the following dispersion relation for transverse waves:
ω2= ω2 p+ K
2
c2 (1.42)
Since the last term on the right hand side is positive, no solutions exist for ω < ωpand thus transverse waves don’t exist below the plasma frequency. For the ω > ωpregime, the plasma sea is predicted to maintain transverse waves propagating with a group velocity of dωdK < c, which are called volume plasmons.
1.4
Surface plasmon polaritons (SPP’s)
Besides the aforementioned volume plasmons in a single uniform metal, another kind of plasmonic wave phenomenon can be observed at the interface of a metal and a dielectric.
To illustrate such an excitation, we go back to equation 1.35, recall that B= µ0µD and setting µ = 1, which is valid for zero magnetization, and get the following equation [18]:
∇ × (∇ × E)= −µ0 ∂2D
∂t2 (1.43)
using the relation D= 0E, the equation takes on the following form:
∇(−1
E · ∇) − ∇2E= −µ00 ∂2E
∂t2 (1.44)
The wave equation for such a system can be simplified in case of negligible spatial variance of , and using the relation c= 1/õ00, to:
∇2E − c2
∂2E
∂t2 = 0 (1.45)
This equation has a similar magnetic counterpart and forms the most basic expression of electromagnetic wave theory. In all generality the electric field can be expressed with a harmonic time dependence E(r, t)= E(r)e−iωt, turning the equation to the form:
∇2E+ k20E (1.46)
where k0=ωc is the wave vector for a wave propagating in vacuum. In this form it is known as the Helmholtz equation [18].
The next step is to define the dielectric geometry of the problem. We are looking at a single flat interface between a metal and a dielectric. Lets say its surface spans the xy-plane and is located at z= 0, meaning = (z).
Figure 1.1: Surface plasmon polariton traveling in the x-direction. It is confined to the interface between two materials at z=0.
We are looking at waves traveling along the surface and to keep things simple choose the x-axis to point in the direction of propagation as depicted in figure 1.1. The field can then be described by E(x, y, z)= E(z)eiβx. Where β
From Maxwell’s equations to plasmons 11
is called the propagation constant of the traveling wave, which is complex in general and corresponds to the length of the wave-vector in the direction of propagation. The wave equation has now been reduced to:
∂2E(z) z2 + (k 2 0 − β 2 )E= 0 (1.47)
Before specifying the dielectric profile and solving the equation, another important aspect comes to light after using the curl equations 1.3 and 1.4 on above’s wave equation. Using the facts that for harmonic time dependence∂t∂ = −iω, for propagation in the x-direction∂x∂ = iβ and due to homogeneity in the y-direction ∂y∂ = 0, the following 6 equations remain:
∂Ey
∂z =−iωµ0Hx
∂Ex
∂z − iβEz= iωµ0Hy iβEy= iωµ0Hz ∂Hy
∂z =iω0Ex
∂Hx
∂z − iβHz= iω0Ey iβHy= iω0Ez (1.48) Out of these 6 coupled equations, two kinds of self-consisting solutions can be formed based on their polarization. The waves associated with the first kind are called transverse magnetic (TM). They consist of the components Ex, Ez and Hy, of which only the magnetic part of the in-plane components is transverse. The wave equation for TM-modes is: ∂2H y ∂z2 + (k 2 0 − β 2)H y= 0 (1.49)
The other two electric field components are then directly related by equations 1.48. In the same way one can construct a wave equation for transverse electric (TE) modes as:
∂2E y ∂z2 + (k 2 0 − β 2)E y= 0 (1.50)
The most simple dielectric profile containing surface plasmon polaritons (SPP’s) is that of an interface between a non-absorbing half-space with constant positive dielectric constant 1and a conducting half space with dielectric 2having a metallic character, meaning <(2) < 0. The solutions to equation 1.49 in the respective half-spaces are:
Hy(x, z)= A2eiβxe−k2z Ex(x, z)= iA2 1 ω01 k2eiβxe−k2z Ez(x, z)= −A1 β ω01 eiβxe−k2z z> 0 (1.51) Hy(x, z)= A1eiβxek1z Ex(x, z)= −iA1 1 ω02 k1eiβxek1z Ez(x, z)= −A1 β ω02 eiβxek1z z< 0 (1.52)
evanescently decaying in the z-direction. The fields have to obey continuity at the interface and therefore: A1= A2 k1 k2 = − 1 2 (1.53)
The expression for Hyfurther has to satisfy the wave equation 1.49, which yields:
k21= β2− k021 k22= β 2− k2
02 (1.54)
Using these equations, the dispersion relation at the interface of two half-spaces with opposite sign of the real part of their dielectrics, a requirement for surface confinement, is found to be [18]:
β = k0 r
12 1+ 2
(1.55)
Which is valid for complex and frequency dependent 2. To complete the analysis, it should be noted that a similar procedure for TE waves leads to the following continuity condition at the interface:
A1(k1+ k2)= 0 (1.56)
As the real parts of k1and k2are both positive, this condition can only be satisfied when A1= A2= 0, e.g. the system does allow any TE modes.
The above derivation serves as a qualitative way to give an understanding of the phenomenon of surface plasmon polaritons. In the remainder of this text, such excitations will be investigated in the presence of the 2d metallic material graphene.
Chapter 2
Modeling plasmonic lattices
In order to investigate topological plasmon excitations in a custom-designed structure of graphene, it is necessary to solve the underlying eigenvalue problem to construct a band diagram and calculate any other relevant quantities. This can be done by various numerical methods, all of which have their own advantages and disadvantages.
In our case, the purpose of the simulations is to investigate topological properties in a variety of different structures and materials. To accommodate for this, the choice was made to use a finite difference time domain (FDTD) method approach [21, 22], as it provides a lot of flexibility and can be applied to a broad range of set-ups. Besides its general versatility, the method also supplies a nice way to simulate a Drude-Lorentz response by directly calculating the polarization terms and their induced currents [18, 23]. The majority of calculations in this study is performed with an open source software package called Meep [24], which is an FDTD for computational electromagnetics. It was originally written as part of a graduate research at MIT and has since been maintained by Simpetus and a developer community at GitHub [25].
2.1
Bloch’s theorem and the concept of band folding
Before jumping into specifics of computational methods to solve the wave equation, a short intermezzo to explain the concepts of periodicity in a system and its consequences for the wave solutions is useful. When a system is described by a geometrically repeating pattern, which could be specified in terms of its potential or dielectric function, it is said to possess translation symmetry. This symmetry affects the basis states needed to describe the full set of solutions to the problem. Such a symmetry can be cast into mathematical form by considering a unit cell, which is the smallest possible geometric cell, along with a set of primitive lattice vectors. They contain the underlying symmetry and can be used to construct the full structure [26, 27]. In mathematical form, the symmetry of the dielectric function in three
dimensions can be described by:
(r + Rl)= (r) (2.1)
for any lattice vector:
Rl= l1a1+ l2a2+ l3a3 (2.2)
where a1, a2and a3represent the three primitive lattice vectors and l1, l2and l3are integer numbers. We now define the translation operator ˆTuby:
ˆ
Tuf(r)= f (r − u) (2.3)
Meaning that ˆTushifts the entire system by a displacement vector u. The fundamental electromagnetic wave equation can, similar to the time independent Schrödinger equation, be written as a linear hermitian eigenvalue problem. In fact, the two show profound similarities and their analogous relationship becomes apparent when comparing their underlying equations:
Quantum mechanics Electrodynamics
Field: Ψ(r, t) = Ψ(r)e−iEt/~ H(r, t)= H(r)e−iωt
Eigenvalue problem: HΨ = EΨˆ ΘH =ˆ ωc22H
Hermitian operator: Hˆ = −~2
2m+ V(r) Θ = ∇ × (ˆ
1 (r)∇×)
Table 2.1 source: Joannopoulos [27]
WhereΨ(r, t) is the quantum mechanical probability wave solution and H(r, t) the magnetic field wave solution. We are actually more interested in the electric field E(r, t), which can easily be recovered from H(r, t) by:
E(r, t)= i ω0(r)
∇ × H(r, t) (2.4)
However, because the relative permeability µ= 1, it has become mathematically more convenient to use H(r, t), as it directly leads to an hermitian eigenvalue problem.
The three translation operators that displace the entire system by a primitive lattice vector ˆTa1, ˆTa2and ˆTa3(or just ˆ
Tai), leave the system unchanged and therefore commute with each other and ˆΘ. Because ˆTai is a unitary operator, ˆ
T∗aiTˆai = ˆTaiTˆ ∗
ai = I, its eigenvalues must have unit norm and can therefore be written as e
−ikifor some phase factor k i. These conditions imply that for every solution of ˆΘ (or ˆH), apart from the band number n corresponding to increasing
Modeling plasmonic lattices 15
state energies, an extra label k can be added such that:
ˆ
ΘHnk(r)= ω2
n(k)
c2 Hnk(r) (2.5)
The label k is a vector who’s components correspond to the discrete translation symmetries of the system, such that:
Hnk(r+ Rl)= eik·rHnk(r) (2.6)
Since the complex exponent is a periodic function, the eigenvector with label kihas the same eigenvalue of the operator ˆ
Tai as an eigenvector with label kiplus any integer times 2π and therefore describe the same state. For this reason, the state with one dimensional discrete translation symmetry can be regarded as "living" on a unit circle line in k-space. For a system with translation symmetry in multiple dimensions it could be viewed as a torus in 2D and a 3-torus in 3D. These toroids can be regarded as topologies and the fact that they consist of closed spaces has big significance for the existence of special topological material properties.
To view Bloch’s theorem from a slightly different perspective, it is useful to construct a cell periodic function as:
unk(r)= e−ik·rHnk(r) (2.7)
In this way the phase dependent periodicity is taken out and the function obeys the regular periodic boundary condition:
unk(r+ R) = unk(r) (2.8)
With this definition, any Bloch wave can be regarded as a lattice periodic function unk(r) multiplied by a standing wave eik·r:
Hnk(r)= eik·runk(r) (2.9)
In this view it is convenient to construct a new space, that of k vectors. This space is called reciprocal space, because the eigenvalues of ˆTai have to be dimensionless, the components of k must cancel the units of r and therefore have a 1/length dimension. The space can be constructed by a set of primitive reciprocal lattice vectors bidefined as to be dual to ai, meaning ai· bj= 2πδi j. A reciprocal lattice vector then takes the form:
Due to the dual definition, the following condition is obeyed:
eiGm·Rl= 1 (2.11)
In this basis, any wavevector k can then be written as:
k= k1 2πb1+ k2 2πb2+ k3 2πb3 (2.12)
The initial eigenvalue problem, who’s solutions are now cast in the form of Bloch waves, was derived with the use of separation of variables and attributing a harmonic time dependence to the wave. The full solution thus has the form:
Hnk(r, t)= eik·r−ωtunk(r) (2.13)
It can be identified as a wave-packet, or as a traveling particle with a momentum P= ~k and an energy equal to ~ω. The relation between k and ω is called the dispersion relation, which is in general non-trivial and forms one of the most interesting properties of a system.
More specifically, in a periodic system the dispersion is said to undergo bandfolding. To demonstrate this phe-nomenon, let’s consider a dispersion relation of the following form:
ω ∝ p
|k| (2.14)
This square root relation also happens to be the dispersion of a plasmon at the surface of a uniform sheet of graphene in vacuum. For one dimension the dispersion relation is plotted in figure 2.1.
Modeling plasmonic lattices 17
The term homogeneous is used when the dielectric function of the material does not depend on the position. On a microscopic scale, graphene is actually formed by carbon atoms in a hexagonal pattern and contains periodic potential symmetry. However, in the regime where the wavelengths are much larger than the inter-atomic distance, which is true for the infrared spectrum, the macroscopic dielectric can still be considered homogeneous.
Now let’s pretend that the sheet of graphene is no longer homogeneous, but somehow acquires spatial periodicity with primitive unit cell size L e.g. (x+L) = (x). Then, according to Bloch’s theorem, each solution with wavenumber k has an identical state for k+ m2π/L, for any integer m. These identical solutions start to cross each other as depicted in figure 2.2.
Figure 2.2: Plot of all solutions to the wave equation for a system with periodic translation symmetry in one dimension. The bands with different origins in reciprocal space become "folded" onto each other.
Because the states at k and k+2πn
L are in fact identical, we can ignore all states outside −π
L < ω < π
L. Each region that holds a full set of solutions is known as a Brillouin zone as depicted in figure 2.2. Typically the first Brillouin zone forms the main focus and is therefore also known as just the Brillouin zone. In this example the symmetry was restricted to one dimension. In general, multiple symmetries in more than one dimension can be present. The energy bands of such a system are often presented by plotting the dispersion relation along lines of high symmetry, to include the maxima and minima of a band.
One important aspect of band structures not mentioned so far, is the possible existence of energy band gaps. In some cases, neighboring energy bands never "touch" each other, meaning there is an energy region in which no solutions exist. The search for a plasmonic band structure with such a gap and its implications will be of special interest in what follows.
2.2
Calculating band structures: The FDTD approach
The way in which an FDTD simulation generally works is by specifying a certain geometry inside a finite compu-tational volume and let it evolve in time according to the governing laws of physics. In the case of electromagnetic waves, these laws consist of Maxwell’s equations. Meep uses the following form of Maxwell’s equations to calculate the time evolution of the field [23]:
dB dt = −∇ × E − JB−σBB (2.15) B= µ(r)H (2.16) dD dt = ∇ × H − J − σDD (2.17) D= (r)E (2.18)
Besides the electrical current density J, the term JBcorresponds to a magnetic-charge current density, a fictional quan-tity that can be convenient for some simulation. The terms σBand σDrepresent frequency independent magnetic and electrical conductivities, respectively. A system with periodic translation-symmetry can be constructed by specifying a unit cell (or super cell) repeating itself an infinite amount of times according to a set of primitive basis vectors, which are defined to span the whole space. Unfortunately, such an infinite-sized arrangement can not be simulated in its en-tirety, as this would require an inexhaustible amount of computational resources. Luckily, the translation symmetries can be exploited in such a way that the full simulation cell is only slightly larger than a single unit cell.
The first step in solving the problem is to define the dielectric profile (r), which represents the geometry of the system. As an example, figure 2.3 shows the construction of an orthogonal patterned set of holes with a square unit cell repeating itself in two dimensions.
Modeling plasmonic lattices 19
Figure 2.3: Structure consisting of orthogonally patterned holes, constructed from a single unit cell (outlined in red) copied in two dimensions according to two primitive basis vectors.
A possible solution is a zero static field, but this is not very illuminating. The next step is to initiate some kind of field dynamics. The system needs to be kicked off by a source in such a way that the properties of interest can be examined. In the case of electromagnetic waves, this source is often conveniently specified by an oscillating point dipole with a Gaussian strength profile defined by:
E
↑= e
−iωte
−(t−to)22W2 (2.19)
Where E↑ denotes the electrical field at a fixed source location. Since the Fourier transform of a Gaussian is also a Gaussian, the temporal width W of the source is inversely proportional to the frequency width and thus by specifying this value, one can control the frequency range at which the systems modes are excited. The pulse also creates traveling wave packets, insuring that enough momentum is available to excite modes in the first Brillouin zone. It is usually best to position the source at a non-symmetric place as this increases the chance to excite all modes.
periodic, but after exciting the system with a source, the electromagnetic field will just spread out over all cells in a non-symmetric way. Luckily, the geometric symmetries also give rise to a specific periodicity of the field, as is constituted in Bloch’s theorem. Accordingly, the field in all neighboring cells can be defined to be the same as in the center cell, except for a phase difference eik·r. By fixing the wave vector k to a specific value, only the field in one cell has to be simulated, the fields in all other cells are the same except for a phase difference [23, 28].
The only thing that still has to be accounted for is the interference created by the fields of neighboring cells, which is the essential cause of dispersion and band folding in a periodic structure. In order to include these effects, the software adds an additional single pixel wide boarder along the edges of the cell and injects the values corresponding to the fields at the opposite side off the cell multiplied by the complex Bloch phase. This phase is only unambiguously defined for a fixed k-vector, so for each point in k-space, a separate simulation has to be performed. This ensues an extra set of boundary conditions to the full simulation cell, as illustrated in figure 2.4.
Figure 2.4: Full simulation cell example. It consists of a unit cell (marked by the thick red line) and a single pixel wide boundary (squares around the unit cell). Each small square represents a discretized point in space. The values at the boundary pixels are defined as the field at the opposite side of the cell multiplied by the corresponding Bloch phase, which is only uniquely defined for a fixed k. The values of two boundary pixels are given as an illustration. The subscripts of E and r correspond to the x and y coordinates respectively.
Modeling plasmonic lattices 21
This procedure forms the essence of calculating the dispersion relation using an FDTD method. The simulation is ran at a specified range of points in k-space, usually along the high symmetry lines. After the sources are turned of, the remaining electric field in each simulation consists of multiple bands with fixed k-vector:
Ek(r, t)= X
n
unk(r)eik·re−iωnt (2.20)
The data of each simulation is subsequently harmonically inversed, meaning that modes of the following form are extracted from the data [29]:
f(t)=X n
Ansin(Ωnt+ Φn)e−Γnt (Real valued field) (2.21)
f(t)=X n
anexp{−i(ωnt+ φn) − γnt} (Complex valued field) (2.22)
All coefficients represent real valued numbers. This is similar to a Fourier transform with complex frequency. Each extracted term in the sum then represents a solution belonging to a specific band and frequency. This way the whole band structure can be reconstructed.
If the only interest is in the band structure, analyzing one point in space is sufficient. This holds enough information to find the energy eigenvalues. However, the calculation of topological invariants, like the Chern number of a band, or the eigenvalues of symmetry operators, require the calculation of full spatial wave functions or eigenfields. They can be reconstructed by running a harmonic inversion at every simulated pixel of the unit cell.
The eigenfield is specified by the amplitude and phase of the field as a function of position. The amplitude represents the maximum value of the field, whereas the phase denotes the moment during the cycle of a harmonic oscillation. The phase is gauge dependent and without knowing how it is gauged, its value is only interesting relative to other points in the unit cell. For many purposes, like calculating Chern numbers, the eigenfields have to be integrated over k-space and it is important to know how the phases of the retrieved eigenfields using the FDTD method are gauged.
In case of a real valued field simulation, the harmonic inversion is performed on a list of real valued numbers representing the field at specific moments in time separated by equal time steps. The phase factor is then gauged by the first value in the list, t0, which represents the moment just after the sources are turned off. The phase value corresponds toΦnin equation 2.21.
With the exception of theΓ point, where eik·r = 1, the field has to be complex valued in order to account for the Bloch boundary conditions, because a real number doesn’t have a phase. Hence the initiating source is described by a complex wave using Euler’s formula (eix = cos(x) + isin(x)), facilitating in the possibility to relate the phase at every position to the source. Maxwell’s equations stay unmodified and the real and imaginary parts of the field only couple
via the boundary conditions [21]. The real part of the field inside the cell remains the physical representation of the field i.e.Φn= φn.
2.3
Benchmarking Jin et al.
In order to test the applicability of the FDTD approach to plasmonics, the first step is to reproduce results from Jin et al. [13] and thereby benchmarking the method. In this paper, a sheet of graphene containing holes in a hexagonal pattern is theoretically predicted to contain non-trivial topological bands when exposed to a large enough magnetic field. Their calculations are based on a two dimensional finite element method, which is a more efficient approach for this specific problem. However, it lacks the versatility to effectively model substrates and other modifications required for experimental purposes. If the results can be replicated correctly, it would suffice as a proof of concept and subsequently this allows for the investigation of all kinds of different geometries and conditions relevant to experiment. When a uniform sheet of graphene is altered by putting holes in them in a hexagonal pattern, this creates a spatially varying periodic potential. According to Bloch’s theorem, the wave solutions to such a problem can be written in the following form:
ψnk(r)= eik·runk(r) (2.23)
Where r is the position and unk(r) is a function that has the same periodicity as the underlying crystal. The eik·rterm represents a plane wave defined by a wave-vector k, as was explained previously.
In this case we consider a hexagonal pattern, which can be constructed from the unit cell and set of primitive lattice vectors as depicted in figure 2.5.
Figure 2.5: Left: Artistic impression of a sheet of graphene containing holes in a hexagonal pattern. Theoretically the solutions in k-space are only continuous for an infinitely repeating periodic potential and practically it should be considered large. Right: The unit cells are marked by the dashed lines, the arrows form a set of possible Bravais lattice basis vectors and the solid line connects the nearest neighbors, thereby forming a hexagon
Modeling plasmonic lattices 23
To identify the first Brillouin zone, the Wigner-Seitz cell in reciprocal space is identified by constructing the plane that contains all points that are closer to the center lattice point than any other lattice point. The first Brillouin zone also has a hexagonal shape. Due to the rotation symmetries, the first Brillouin zone can still be reduced. By identifying the high symmetry points,Γ, M and K, and drawing high symmetry lines by connecting these points, this marks an Irreducible Brillouin zone (IBZ). The resulting IBZ can be seen in figure 2.6.
Figure 2.6: The full hexagonal shape of the Brillouin zone is outlined in gray. By connecting the high symmetry points, marked by dots, the irreducible Brillouin zone is formed.
The Brillouin zone in this case is 2-dimensional. Along the high symmetry lines,Γ-M, M-K and K-Γ, the maxima and minima of the bands are located. Therefore, energy band diagram constructed by calculating the eigen energies at k-vectors on the edge of the IBZ will give information about band gaps.
The sofware package Meep, and most other FDTD’s, does not allow for a unit cell with non-orthogonal bravais lattice vectors [23]. A general work-around is to define a larger cell, a super cell, that can be used to construct the same lattice using orthogonal basis vectors. Figure 2.7 shows the geometry of the super cell used for the purpose of modeling the present system. Its basis vectors are orthogonal:
~b1= l √ 3 0 , ~b2= l 0 1
Figure 2.7: Depiction of the super cell used to describe the hexagonal system with orthogonal lattice vectors. The epicentres of the circular waves indicate the origin of the two sources. They are separated by a basis vector (blue arrow) that belongs to the hexagonal symmetry of the hole structure and their amplitudes differ by the appropriate Bloch phase.
The area of the supercell is exactly twice that of the unit cell, which means the supercell’s BZ in reciprocal space is exactly half the size of the unitcell’s BZ. This results in unwanted bandfolding, meaning that the final result will show additional modes that actually belong to different k-vectors. As explained in the previous paragraph, solving for each specific k-point is based on the premise that Bloch’s boundary conditions are correctly injected at the edges of a unitcell. To satisfy this condition in a supercell, the field created by the initiating source has to be duplicated a primitive lattice vector away with a phase difference of ei·k·R. Using this approach, the second source should theoretically cancel all unwanted modes. However, because the simulation grid has a finite size, some of these modes might still survive, but will have relatively low amplitude.
Besides the necessity of using a finite grid size, discretizing space also forces one to use a certain mesh shape. The implications are that the circular holes and spherically propagating sources have to be constructed from square meshes. The supposedly round edges are therefore a bit roughly modeled and can cause undesired scattering of the field. Effectively, this causes the "signal" to gain noise, which has to be dealt with during post processing. Suffice to say, increasing the grid resolution gives better results, but costs more computational resources.
Graphene consists of a single layer of carbon atoms and is in many cases regarded as a 2D material, whose thickness is mathematically described by a delta function. However, the electromagnetic field of a plasmon at the interface of a sheet of graphene suspended in vacuum, decays only evanescent with distance from the plane. Therefore, to adequately investigate its behavior, a 3D grid has to be used for the full simulation cell, with the implication that the thickness has to be represented by a finite mesh size (δ). This does not form a major concern as long as the plasmon wavelength
Modeling plasmonic lattices 25
(λp) remains much larger than the width e.g. δ λp. In actuality, the thickness of graphene is finite and of atomic size (≈3Å) anyway [30].
As mentioned, the boundary conditions at the edges of the unit cell, which follow from the translation symmetry of the in-plane directions of the sheet of graphene (xy-plane), are specified to be Bloch periodic.
Figure 2.8: Schematic overview of the geometry parallel to the z-direction.
In the direction orthogonal to the plane (z-direction), the size of the simulation cell has to be sufficiently large to account for the evanescent field, which is satisfied with approximately one wavelength. The boundary condi-tions in this direction are completely reflective, there-fore absorbers are placed to avoid unwanted reflections from the remaining evanescent field. Figure 2.8 shows the schematic geometry in the z-direction.
The sole dispersive material in the geometry consists of the sheet of graphene and its dielectric function is modeled by giving it appropriate Drude-Lorentz terms. The spectral region explored by Jin et.al. lies in the in-frared, as it supports a set of appropriate properties to establish topologically protected states. Coincidentally, many experimental techniques used for investigating lo-cal dielectric material properties, like AFM-IR and s-SNOM, operate in this regime. According to Jin et al. [13], the only significant remaining term is of Drude
form. This claim is supported by the graphene handbook [30], where the electric conductivity is approximated using the Kubo formalism as is illustrated in figure 2.9. A distinction is made between interband and intraband contributions, of which the latter leads to a Drude like character.
Figure 2.9: Graphs taken from graphene handbook chapter 19.2 [30], (a) showing the normalized conductivity of graphene, (b) the intraband contribution and (c) the interband contribution. Calculated with the Kubo formalism.
The specific value depends on the Fermi energy (Ef) which is proportional to the carrier density. This energy (denoted by µcin the figure) depends in the first place on production specifics, like type of substrate, and is dynamically tunable by applying a gate voltage. In this specific simulation Ef is assumed to be 0.2 eV, which means that for wavelengths above 4.0 µm (below 75 THz) the intraband term is dominant. As the electric conductivity of a 2d-material is specified as a surface conductivity, it needs to be divided by the thickness of a single pixel (δp) in the 3D simulations. The full expression for the dielectric function becomes:
(ω) = 1 − e2Ef δp0~2
1
ω(ω + iγ) (2.24)
In order to accurately extract the correct modes from the time-evolved data, the simulation time needs to be long enough to contain at least a couple of oscillation periods. Additionally, the modes are easier to distinguish if their decay time is long and accordingly the γ-term is set to zero. The effect of the The drude-loss term, estimated at γ/2π ∼ 1Thz, is that it "smears out" the gap region [13]. A typical distribution of the perpendicular component of the electric field (Ez) during the simulation is visualized in figure 2.10. The field appears to be strongest at the edges around the holes.
(a) Parallel to the sheet. (b) Perpendicular to the sheet. Figure 2.10: Visualization of a typical Ezdistribution during the simulation.
Modeling plasmonic lattices 27
In this system, only TM-modes are predicted to exist, therefore the choice is made to only concern ourselves with the z-component of the electric field. After the simulations are finished, the data is harmonically inverted. The thereby obtained mode spectrum still includes erroneous folded bands with low amplitude and spurious modes. These are filtered by demanding a sufficient amplitude and quality-factor (=π f /decay). The latter quantity, provided by harminv, is expected to be small for spurious modes. The band structure obtained in this manner, for a 400nm unit cell size and 200nm hole diameter, is shown next to the results found by Jin et al. in figure 2.11. The band diagram obtained by means of our alternate method, 3D FDTD, closely agrees with the results previously found with a 2D finite element method.
(a) 2D finite element method (Jin et.al.) (b) 3D FDTD method
Figure 2.11: Plasmonic band structure of a sheet of graphene with hexagonally patterned holes in vacuum without external magnetic field (B= 0T). Results obtained by Jin et al. [13] using a 2D finite element method (a) next to the results obtained by our 3D FDTD method (b).
Despite the many similarities, there remain some apparent inconsistencies that need to be addressed. First of all, the upper bands are missing a couple of dots, which is a consequence of the low amplitude filter. The data was only analyzed at a single point in space and the missing modes were probably close to a nodal point. This can be solved
by increasing the grid resolution and/or analyzing multiple points. Secondly, the lowest energy band is not completely resolved to k=0. This is a consequence of the simulated time being finite. For the lowest band, the energy approaches zero when k → 0. Consequently, the oscillation appears as time independent and cannot be resolved. If less then a single period is simulated, it becomes difficult to recover the mode. Thirdly, if band frequencies are closer together, they become harder to distinguish leading to a higher possibility of error, which can explain the multiple missing points of the fifth band. Lastly, the overall energies calculated by Jin.et.al. are slightly higher in value. This can be explained by the fact that due to limited computational resources, we had to settle with a resolution corresponding to a thickness of 4nm. This explanation is substantiated by plotting the energy of the first band at the M point as a function of resolution. As can be seen in figure 2.12 it indeed appears to converge to the correct value of 12.5 THz at higher resolution.
Figure 2.12: Energy of the first band at the M point plotted as a function of resolution (pixels/unit cell length), indicating convergence to 12.5 Thz.
The small deviations can all be dealt with by increasing the resolution, which should be done when more accuracy is required. However, for the current benchmark this accuracy is deemed sufficient.
The ambition is to eventually investigate non trivial topological properties and therefore our second benchmarking simulation involves bands with non-zero Chern numbers. Such states are said to arise when a homogeneous magnetic
Modeling plasmonic lattices 29
field perpendicular to the graphene sheet is applied to the former simulation, thereby breaking time-reversal symmetry. The effect of a magnetic field can be modeled by modifying the Drude-Lorentz equation. To show this, it is useful to look at the polarization equation used to simulate a Drude-Lorentz term [23]:
d2P dt2 + γ
dP dt + ω
2P= σω2E (2.25)
Where σ is a constant giving the strength of the resonance. The altered Drude-Lorentz model leads to the following form of the polarization equation:
d2P dt2 + γ dP dt − dP dt × b+ ω 2P= σω2E (2.26)
The added cross product term simulates the effect of the magnetic field and additionally breaks time reversal symmetry. For a Drude form, the ω2P term is omitted. The "bias vector" b points in the direction perpendicular to the rotation corresponding to the gyrotropic effect and its magnitude equals its angular frequency, which depends on the strength of the applied external magnetic field. The modification is applied to the simulation by specifying b/2π:
b 2π = 0 0 fc (2.27)
With fcthe cyclotron frequency given by:
fc= eB
2πm? (2.28)
Where B is the magnitude of the magnetic field, e the electron charge and m?the Drude mass, which is about 4% of the electron mass at Ef = 0.2 eV. The results for an applied magnetic field of 4 Tesla can be seen in figure 2.13.
(a) 2D finite element method (Jin et.al.) (b) 3D FDTD method
Figure 2.13: Plasmonic band structure of a sheet of graphene with hexagonally patterned holes in vacuum with an applied external magnetic field (B= 4T). Results obtained by Jin et al. [13] using a 2D finite element method (a) next to the results obtained by our 3D FDTD method (b).
Similar to the previous result, the band structures appear to be closely comparable, showing only minor differences. It should be noted that the fifth band in the right panel of figure 2.13 is difficult to discern. Even though this is most likely fixed by increasing the resolution, it does demonstrate a disadvantage of the FDTD method. Specifically, finding all solutions to the eigen-equation in this way, however likely, is never completely guaranteed.
Another missing element in our result is exposed by the red numbers added to each band in the diagram of Jin et al.. They denote the so called Chern number of the band, which forms a measure to the topological non-triviality of a state. The endeavor to recalculate these numbers is documented in appendix A.
Chapter 3
Tailoring the energy bands
Now that the method is sufficiently benchmarked for calculating plasmonic band structures in a single layer of graphene, we are in a position to search for potential topologically non-trivial configurations.
The results obtained by Jin et al. [13] exhibit very interesting topological phenomena, which, to the best of our knowledge, have not been observed experimentally at the time of writing. This is likely due to the theoretical nature of the article, which does not deal with some of the circumstances that are difficult to attain in experiment.
One of the most significant obstacles arises from the requirement to perform measurements in a large magnetic field, as measuring equipment is usually very delicate and vulnerable to large external fields. The reason for using it is twofold: It creates an energy gap and it breaks time reversal symmetry. Loosely speaking, to create a system that possesses topologically protected edge-states, you need bulk states with differing Chern numbers. When a photonic band has a degeneracy in one or more places, it will have a shared Chern number, known as a composite Chern number. Therefore, the only way to create bands with independent Chern numbers, is to break this degeneracy. Fortunately, there are other kinds of symmetry breaking tactics to create a band gap. In order to make the bands topologically non-trivial, however, the time-reversal symmetry still has to be broken with some kind of pseudo-magnetic field.
A second significant obstacle originates from the fact that a hypothetical single sheet of graphene suspended in a vacuüm is required, which is practically difficult to attain. A mono-layer of graphene is generally fixed on a substrate to give it sufficient structural stability. To insure experimental feasibility, alterations have to be made in this respect as well.
In short, it will be necessary to deviate from the structure proposed by Jin et al. [13] and in the remainder of this chapter, alternative ways to get topologically protected edge-states are examined.
3.1
Symmetry breaking: The creation of band gaps
The importance of symmetry in physics is often emphasized. This is not without reason, because by considering the symmetries of a system, one can induce many of its basic properties without solving the underlying problem in its entirety. For instance, the properties of translation symmetry lead to Bloch’s theorem, which allows for an elegant way of writing wave solutions as entities living in the Brillouin zone. In a similar way, other kinds of symmetries like rotational, mirror, time-reversal and inversion symmetry, can be used to make general predictions about the properties of a band structure, including the possible emergence of a band gap.
The exact relationship between the symmetries of a unit cell and the size of a band gap is unknown, but it is often possible to relate a degeneracy directly to a specific symmetry [31]. Because the fields of the modes belonging to different bands are solutions to a hermitian operator, they have to be orthogonal [32]. The inner product between two vector fields F(r) and G(r) is defined as [27]:
(F(r), G(r))= Z
d3rF∗(r) · G(r) (3.1)
We now return to the eigenvalue problem of Table 2.1 for the magnetic field H, with equations:
ˆ ΘH = ω2
c2H (3.2)
ˆ
Θ = ∇ × ((r)1 ∇×) (3.3)
By inserting the eigenvectors in the form of Bloch states:
Hk(r)= eik·ruk(r) (3.4)
it can be transformed into another hermitian eigenvalue problem [27]:
ˆ Θkuk= ω2 k c2uk (3.5) ˆ Θk= (ik + ∇) × 1 (r)(ik+ ∇)× (3.6)
For two normalized modes uk,1(r) and uk,2(r) with frequencies ωk,1and ωk,2respectively, hermiticity implies [27]:
ω2
k,1(uk,2(r), uk,1(r))= c2(uk,2(r), ˆΘkuk,1(r))= c2( ˆΘkuk,2(r), uk,1(r))= ω2k,2(uk,2(r), uk,1(r)) (3.7) → (ω2k,1−ω2k,2)(uk,1(r), uk,2(r))= 0 (3.8)
Tailoring the energy bands 33
When ωk,1 , ωk,2, this means the two states have to be orthogonal. When two modes share the same energy, ωk,1= ωk,2, the states are degenerate, but do not necessarily have to be orthogonal. However, because they are solutions to a linear operator, any linear combination of the two states also forms a valid solution with the same frequency. In this way, one can always construct two modes that are orthogonal and in all generality degenerate modes can be constructed to be orthogonal.
To interpret the nature of degenerate states, it is useful to look at the relation between their energies and field distributions. The energy of a mode is related to the fields by [32]:
E= 1 8π
Z 1
(r)|D|2+ |H|2dr (3.9)
Loosely speaking, the modes with their fields localized more in a region with higher dielectric constant typically have a lower energy. The existence of degenerate modes is thus closely related to symmetry, as eigenfields with the same energy can be distributed in such a way that their inner product, which involves an integral over the unit cell, is kept equal to zero. By the same token, breaking such a symmetry can lift a degeneracy and create a band gap.
To further identify the origins of band structure properties, the eigenvalues of the symmetry operators can be used to characterize the eigenstates. To illustrate this, we take a closer look at inversion symmetry. A system holds such a symmetry when an inversion center exists. When taken as the origin of a coordinate system, such a center has the property that every point (x,y,z) in the unit cell has an indistinguishable point at (-x,-y,-z). For a more mathematical description, we define an operator ˆI that inverts a vector a, such that ˆIa= −a. The operator ˆOIthat inverses a vector fieldcan then be defined as inverting both the vector field f as well as its argument [27]:
ˆ
OIf(r)= ˆIf( ˆIr) (3.10)
If a system has inversion symmetry, the result of applying the operator ˆΘk, is the same as first inverting the field, then applying ˆΘkand then undo the inversion. In other words, the two operators commute:
ˆ Θk= ˆO
−1
I ΘˆkOˆI (3.11)
→ [ ˆOI, ˆΘk]= ˆOIΘˆk− ˆΘkOˆI = 0 (3.12)
Consequently, for any mode of the system Ekthat forms a solution to the eigenvalue problem, the following relation has to be true: ˆ Θk( ˆOIEk)= ˆOI( ˆΘkEk)= ω2 k c2( ˆOIEk) (3.13)
differ by a constant. Applying the inversion operator twice should return the same state:
Ek= ˆOI( ˆOIEk)= α2Ek (3.14)
Therefore the value of α has to be either+1 or -1. This means that if the structure of a system contains an inversion symmetry, an extra label of+1 or -1 can be added to the state, thereby characterizing the eigenfunctions as being either even or uneven under inversion. This type of labeling will be used later to analyze the band structures presented in section 3.3.
3.2
Substrates: Experimental necessities
Figure 3.1: Modified structure viewed from the side (per-pendicular to the graphene plane). The periodicity is formed by the gold cylinders. A protective layer of hexagonal Boron Nitride (h-BN) is added on top of the graphene sheet. For structural stability, a substrate of sil-icon dioxide (SiO2) is used.
In order to reduce nano-metre scale strain, we use a sub-strate of silicon dioxide (SiO2), which allows for the pro-duction of high quality graphene. Furthermore, drilling holes in graphene is a precarious task and will most likely ruin it. To construct a periodic structure, lithographic methods are used to create circular gold islands on the graphene instead of holes. However, the photoresists used for lithography can also modify the electronic properties of graphene. For this reason, a single layer of hexa boron nitride (h-BN) is used to form a protective layer between the gold and the graphene. To bring about these changes to the simulations, the layout of the simulation cell is altered to the geometry depicted in figure 3.1.
To simulate the additional materials, we need the ap-propriate dielectric functions in the frequency range of in-terest. More specifically, as the purpose of the simulations is to prepare for experimental measurements, we need to make sure they are suitable for our CO2laser (27.85-32.37 Thz) and quantum cascade laser (QCL) (49.47-58.46 Thz). Because h-BN has a more or less frequency
indepen-dent dielectric function at the frequencies of interest, it can be approximated by a constant real value of 2.1 [33]. The optical response of gold in the IR region is dominated by intraband transitions and the dielectric response can
Tailoring the energy bands 35
be modeled very closely by the Drude model using ~ωp= 8.5 eV and τ = 1/γ = 14 fs [34]. For the purpose of mode extraction from the FDTD simulations, the value of γ is set to zero to provide better accuracy.
In the case of SiO2, the dielectric function shows a clear optical inter-band transition, recognized as a Lorentzian oscillator term, in the region between 800 cm−1and 1300 cm−1with its peak at 1080 cm−1[35]. A small oscillator with its resonance peak at 850 cm−1can also be identified. However, when fitting the data with two Lorentzian terms, the best fit consists of a main Lorentz peak at 1080 cm−1and a second peak at 54638 cm−1, which suggests the existence of one or more optical transitions outside the range of the plot. Alternatively, a fit with a single Lorentzian term results in the contribution of any high energy transitions to be captured inside the constant ∞. The former has a better overall fit, but the latter is more suitable for FDTD simulations as the real part is still matched closely, whereas the imaginary part goes to zero more quickly away from the main peak, which results in longer excitation lifetimes. Both fits can be seen in figure 3.2.
Figure 3.2: Dielectric fits of SiO2 using a single Lorentzian term (left) and two Lorentzian terms (right).
To check the usefulness of these fits, as well as qualitatively evaluate the difference of either using a single peak or double peak to model the dielectric response of SiO2, a simulation is performed for a homogeneous sheet of graphene on a SiO2substrate and compared to calculations by Fei et al. [36].
The system has continuous translation symmetry in the xy-plane and therefore band folding is not applicable, as the Brillouin zone becomes infinitely large. It is expected to show two bands, which are associated with plasmon excitations in graphene and phonon excitations in SiO2. More specifically, the individual dispersions of graphene plasmons and SiO2phonons cross each other. This can be seen in figure 3.3c, which shows the results of simulations with the individual materials. However, due to plasmon-phonon interactions around their resonance frequency, the bands will start forming hybrid plasmon-phonon modes and open a band gap. This can be seen in figures 3.3a and 3.3b, which show data from Fei et al. next to our own calculations. Due to the low band density, even excitations with relatively short lifetimes in the double peak SiO2model were not difficult to extract.
