Fractal density of states in
graphene quasicrystals
THESIS
submitted in partial fulfillment of the requirements for the degree of
MASTER OF SCIENCE
in
PHYSICS
Author : Joana Fraxanet
Student ID : s2069105
Supervisor : Anton Akhmerov
Carlo Beenakker
Co-supervisors : Daniel Varjas
2ndcorrector : Yaroslav Blanter
Fractal density of states in
graphene quasicrystals
Joana Fraxanet
Instituut-Lorentz, Leiden University P.O. Box 9500, 2300 RA Leiden, The Netherlands
July 24, 2019
Abstract
Quasi-periodic Hamiltonians have a fractal density of states with a hierarchical structure of bands. One pathway to controllably creating quasi-periodic Hamiltonian is stacking layers of graphene twisted with
incommensurate angles. Nevertheless, due to the weakness of the inter-layer coupling, the fractality in these type of systems remains to be
observed. We investigate the fractality of graphene quasicrystals numerically. To achieve this, we develop an algorithm to compute the
fractal dimension of the density of states. Using this approach, we demonstrate that increasing the value of inter-layer coupling of dodecagonal graphene does generate fractal density of states. Finally, we
propose a four-layer system that utilizes the flat band appearing in the twisted bilayer graphene at the magic twist angles TBLG together with the dodecagonal stacking as a pathway towards observing fractal density
Contents
1 Introduction 1
2 Graphene systems 3
2.1 Introduction to graphene 3
2.2 Twisted bilayer graphene (TBLG) 7
2.2.1 Moir´e superlattices in TBLG 8
2.2.2 Tight-binding models for TBLG 10
2.3 Finite incommensurate systems 14
3 Quasicrystals 17
3.1 Almost periodic functions and quasi-periodic tilings 18 3.1.1 The canonical cut and projection method 19
3.2 Quasicrystals in condensed matter 21
3.2.1 Fractal electronic spectrum 22
vi CONTENTS
3.3 Analysis of graphene quasicrystals 25
4 Methods 29
4.1 The fractal dimension 29
4.1.1 Box-counting method for 1D curves 30
4.2 Simulations 32
4.2.1 Increasing the rotation angle by tuning the inter-layer
potential 32
4.2.2 The kernel polynomial method 35
5 Results 39
5.1 Analysis of 1D quasicrystals 39
5.1.1 Incommensurate chains 40
5.1.2 The Fibonacci chain 42
5.2 Graphene quasicrystals 45
5.2.1 Magic angle TBLG 47
5.2.2 Dodecagonal graphene 48
5.2.3 TBLTBLG 50
6 Future outlook and conclusions 55
7 Acknowledgements 59
vi
Chapter
1
Introduction
Quasicrystals are materials which are almost periodic and long range or-dered. They were first discovered in 1982 by Dan Shechtman [1], when he measured Bragg peaks corresponding to a 10-fold symmetric structure. Bragg peaks indicate the existence of long-range order, but the 10-fold symmetry is not allowed in periodic crystalline structures. At that time, these were thought to be the only existing structures in condensed matter. This discovery proved catalytic for a major effort into the search for novel structures containing non-crystallographic symmetries. Nowadays, not only new types of quasicrystals have been found in nature but techniques have also been developed to artificially synthesize them.
Recently, researchers were able to produce a new type of quasicrystals by stacking and twisting layers of graphene at incommensurate angles. This new method of obtaining quasicrystals opens up a new world of possibil-ities. By tuning parameters such as the number of layers and the rota-tion angles, we can obtain systems showing very diverse behaviour. Not long ago, exciting properties such as a correlated insulating phase [2] and unconventional superconductivity [3] were measured in a twisted bilayer system with a rotation angle of 1.08◦. Graphene quasicrystals are extrin-sic quaextrin-sicrystals, in the sense that they are built from periodic layers of graphene. It will become clear that the relevance of the quasi-periodicity
2 Introduction
in the system depends on the coupling between the layers.
Quasicrystals are interesting because their electronic spectrum has fea-tures that cannot be observed in conventional crystals. In general, qua-sicrystals can be understood as projections of higher-order periodic struc-tures. As a consequence, the indexing of the Bragg peaks is dense, and the reciprocal space looks intrinsically different to that of a periodic system. Many propagative and non-propagative modes coexist in a quasicrystal, enriching the electronic spectrum. This manifests in the density of states through the existence of many sharp peaks and gaps at different energies. In recent years, there has been an increasing amount of literature on these types of systems, but we are still lacking the necessary tools to de-scribe their electronic spectra. Our understanding of condensed matter systems heavily relies on the Bloch theorem, but since quasicrystals are not periodic it cannot be used to study them. Moreover, quasicrystals are too correlated to use the existing tools for disordered systems. Since a general theoretical approach has not yet been found, we will focus on de-tecting and analyzing the main signs of quasi-periodicity. In particular, we are going to focus on the fractal dimension of the electronic spectrum and use it as a classification metric of quasi-periodic systems. There are other interesting approaches that we will also take into account, such as studying the localization of the eigenstates due to the lack of periodicity.
The goal of this thesis is to design a mechanism to characterize the fractal dimension of the density of states of a system. We believe that the fractal dimension captures the strength of the quasi-periodicity of a system, and that we can use this to better understand quasicrystals. We will prove our method by applying it to two known one-dimensional qua-sicrystals: an incommensurate chain and the Fibonacci chain. Finally, we will turn to graphene quasicrystals. We will look into three different sys-tems. The first two, magic angle TBLG and dodecagonal graphene, have been widely studied [4][5]. We will show that unless the coupling between the layers is increased to unphysical values, these systems do not show a strong quasicrystallinity. That is why we propose a new system, twisted bilayer twisted bilayer of graphene (TBLTBLG), which has been specially designed to show strong signs of quasicrystallinity in its electronic spec-trum.
2
Chapter
2
Graphene systems
In this Chapter we will introduce the main properties of graphene and the models that we will use to describe its electronic structure. In Section 2.1, we will start by introducing monolayer graphene, focusing on the the tight-binding model for the non-interacting electron picture and empha-sizing the electronic properties in the low-energy range. In Section 2.2, we will study the possibility of stacking layers, focusing on the twisted bi-layer of graphene (TBLG), which is one of the central ideas of this thesis. We will describe different approaches to model its electronic structure in the case of commensurate rotation angles. Finally, in Section 2.3, we will present finite incommensurate systems. These systems are not periodic, which will bring us to Chapter 3.
2.1
Introduction to graphene
Graphene is a two-dimensional allotrope of carbon discovered in 2004 by Novoselov et al. [6]. Its atomic structure is characterized by two types of bonds which originate from the valence orbitals of carbon atoms (2s, 2px, 2py, 2pz). The σ bonds consist of (sp2) hybridized orbitals, and connect each carbon atom to its nearest neighbours. They are the responsible for the robustness of the two-dimensional lattice structure. The pz orbitals
4 Graphene systems
point outside of the two-dimensional plane and they form covalent bonds with carbon atoms nearby. These orbitals are called π orbitals and they correspond to one highly delocalized free electron [7]. These electrons are responsible for electronic transport in graphene.
Figure 2.1:(a) Lattice structure of graphene with two triangular sublattices A and B, primitive lattice vectors and primitive unit cell. (b) Brillouin zone of graphene with its symmetry points. (c) Graphene two-dimensional electronic band struc-ture and its projection through a one-dimensional path in the BZ (d).
The atoms in graphene are organized in a planar honeycomb lattice. The honeycomb lattice, as shown in Fig. 2.1 (a), is equivalent to a trian-gular lattice with two atoms in the basis. As a consequence, graphene has two sublattices. The distance between nearest neighbors is a0 = 1.42 nm and the vectors that connect the nearest neighbors are:
δ δδ1 = a0 2( √ 3, 1), δδδ2 = a0 2(− √ 3, 1), δδδ3 =a0(0, 1). (2.1) 4
2.1 Introduction to graphene 5
The primitive lattice vectors can be written as:
a1 =a(1, 0), a2 =a(1 2,
√
3
2 ), (2.2)
where a = √3a0. We can compute the reciprocal lattice vectors using
aibj =2πδij: b1 = 2π a (1, 0), b2 = 2π a ( 1 √ 3, 2 √ 3). (2.3)
The first Brillouin Zone is represented in Fig. 2.1(b) [8]. TheΓ point is the center of the Brillouin zone and the six corners can be classified into two inequivalent points, namely:
K= 2π 3a(1, 1 √ 3), K 0 = 2π 3a(1,− 1 √ 3). (2.4)
The K and K0points, as we will see later, play a central role in the study of the low energy excitations of graphene [9].
We are interested in the non-interacting electron picture. If we consider hopping only between nearest neighbors, the tight binding Hamiltonian for electrons in graphene is [10]:
ˆ
HTB,n.n. = −t
∑
σ,<i,j>
(a†iσbjσ+H.c.), (2.5)
where a†iσ(b†iσ)annihilates an electron with spin σ at site i in the A(B) sub-lattice. The brackets stand for nearest-neighbor hopping and H.c. is the Hermitian conjugate. Here we are considering the onsite energies to be zero.
In order to obtain the electronic structure of graphene we need to switch to k-space. Let us consider the Fourier transform of the electron operators:
ak,σ = √1 N
∑
n eikRAna nσ, bk,σ = √1 N∑
n eikRBna nσ. (2.6)Here, N is the total number of unit cells and n runs over all the unit cells of the crystal. The position of an atom in the n-th unit cell and sublattice α
6 Graphene systems
is denoted by Rα
n. If we plug in this expressions in the Hamiltonian from
eq. 2.5 we obtain: H =
∑
kσ Ψ† kσHˆkΨkσ, (2.7) where: ˆ Hk = −t 0 ∑3j=1exp(ikδδδj) ∑3 j=1exp(−ikδδδj) 0 ! , Ψkσ = (akσ, bkσ)T. (2.8) The single-electronic spectrum in Fig. 2.1 (c) can be obtained by diagonal-izing the Hamiltonian matrix from eq. 2.8. We obtain two eigenfunctions that depend on k: e(k1,2) = ±t exp (−ia0kx) 1+2exp 3ia0kx 2 cos √ 3a0ky 2 (2.9)The electronic spectrum consists of two bands, π and π∗, both degen-erate with respect to the spin σ. They correspond to a valence and a con-duction band respectively, as shown in Fig 2.1 (d). All the plots regarding lattice structures and electronic spectra in this thesis have been computed using Kwant [11]. At the corners of the Brillouin zone (K and K0 points) the energies e(K1,2,K)0 become zero. Therefore, the bands are degenerate,
creat-ing a Dirac point. Around the Dirac point the dispersion relation is linear (Dirac cones). For low energies, the dispersion relation near the K and K0 points can be approximated by:
E(K+q) = ±¯hvF q
q2
x+q2y, (2.10)
where|q| |K|and vF = 3a20t ≈106m/s is the Fermi velocity. The Fermi velocity corresponds to the slope of the bands [12] and does not depend on the energy or the momentum. In neutral graphene, the chemical potential is situated exactly at the Dirac point. Then, Graphene is a semimetal with unusual properties. The Dirac Hamiltonian describes masless fermions, and therefore many of the properties of QED appear in graphene, but at much lower speeds.
6
2.2 Twisted bilayer graphene (TBLG) 7
2.2
Twisted bilayer graphene (TBLG)
Stacking layers of graphene can change its electronic properties. In this Section we will focus on bilayer graphene, which can exist in three differ-ent forms [10]:
• The AA bilayer is the simplest stacking, but it is a metastable con-figuration. Each carbon atom of the second layer is placed above the corresponding atom in the first layer as shown in Fig. 2.2 (a). Its elec-tronic structure contains two Dirac points shifted in energy at each K and K0 point. AA-bilayer graphene is a metal even when is not doped.
• In the AB bilayer, or also called Bernal phase, half of the atoms of the second layer are placed above the centers of the hexagons of the first layer, and the other half are placed on top of other carbon atoms as shown in Fig. 2.2 (b). This is the most stable configuration. It can also be described as a twisted bilayer graphene with a rotation angle of 60◦. Its electronic structure contains one Fermi point at each K and K0 point, where two parabolic bands touch. A gap can easily be created by applying an electric field.
• In this thesis we will focus on the twisted bilayer of graphene (TBLG), which is also a stable configuration. We can consider either two AB or AA stacked layers, where the top layer is rotated by some angle θ with respect to the lower layer as in Fig. 2.2 (c). Its electronic struc-ture preserves the linear dispersion relation with well defined Dirac cones, but other electronic properties depend on the rotation angle and can vary considerably. No gap opens by applying an electric field.
We can use our knowledge of the single layer electronic structure to understand bilayer systems better. Even though the sp hybridization be-comes more complicated when we have more than one layer, the π band is still half-filled when the samples are not doped. The tight binding Hamil-tonian in eq. 2.5 can be easily extended by adding new hopping terms that couple the two layers. Then, we can study the electronic spectrum
8 Graphene systems
Figure 2.2: Three different stackings in bilayer graphene: AA-stacking (a), AB-stacking or Bernal phase (b) and twisted bilayer graphene (TBLG) with θ=30◦.
of TBLG at low energies in a similar way. Before introducing the tight binding model for twisted bilayer graphene we are going to classify these systems in terms of the rotation angle θ. This is important because depend-ing on the rotation angle, the tools to study the system will be completely different.
2.2.1
Moir´e superlattices in TBLG
The main feature of twisted bilayer graphene is the existence of a Moir´e pattern, which can be measured using STM techniques [13]. The Moir´e pattern exists for any rotation angle θ. The rotation angle is related to the measured Moir´e period L as follows:
L= a0 √ 3 2sin(θ/2), θ<30 ◦ (2.11)
The periodic repetition of a large unit cell, called the Moir´e unit cell, only happens for a discrete set of commensurate angles. Commensurate angles satisfy the following relation:
cos(θ) = 3m
2+3mr+n2/2
3m2+3mn+n2 , (2.12)
where m and n are co-prime positive integers [10]. The Moir´e unit cell 8
2.2 Twisted bilayer graphene (TBLG) 9
Figure 2.3:Moir´e pattern for TBLG with a rotation angle of 1.08◦. In this case, the Moir´e unit cell corresponds to one Moir´e period.
does not always coincide with the Moir´e pattern. One Moir´e cell might contain one or more Moir´e periods, as shown in Fig. 2.3.
To understand the meaning of the coefficients m and n we need to de-fine a set orthogonal lattice vectors for graphene:
t1 =a1, t2 =2a2−a1, (2.13) where a1,2are the usual primitive lattice vectors for the honeycomb lattice from eq. 2.2. One way to define a commensurate angle is to pick a rotation such that the site in (m, n) is moved to (−m, n) in the orthogonal basis. To simplify the computation, we can take m to be 1 and n to be an odd number. Then, the rotation angle is:
θ =2 arctan 1 √ 3n . (2.14)
By following this procedure, we make sure that the system has transla-tional symmetry. This symmetry is captured by the lattice vectors:
T1 =
nt2+mt1
2 , T2 =R60T1, (2.15)
where R60 denotes a rotation of 60◦. The Moir´e unit cell is always larger than the graphene unit cell. If we consider an incommensurate structure, the Moir´e unit cell becomes infinite and then the system is not periodic anymore. The number of atoms N in a Moir´e unit cell with coefficients m
10 Graphene systems
and n is computed as:
N =4|T1||T2|
|a1||a2|. (2.16)
Now we can to derive the corresponding Moir´e reciprocal lattice vectors using that Tib
(M)
j =2πδij[5]. It can be shown that they satisfy the follow-ing: b(1M) =b2−b02, b (M) 2 = (b 0 1+b 0 2) − (b1+b2), (2.17) where b01,2 are the reciprocal lattice vectors of the rotated layer. Since the Moir´e unit cell is larger than the graphene unit cell, the resulting Moir´e Brillouin zone will be smaller than the usual graphene Brillouin zone.
We can index the atoms in terms of their position in each of the layers and the sublattice they belong to. In the next Sections, we will refer to the following notation:
RnA = n1a1+n2a2+τττA, (2.18)
RBn = n1a1+n2a2+τττB, (2.19)
RAn0 = n1a01+n2a02+τττ0A, (2.20)
RBn0 = n1a01+n2a02+τττ0B, (2.21) where (n1, n2) are integers. Here, we denote the rotated primitive lattice vectors of the second layer as a01,2 = Rθa1,2. The basis of graphene is
de-noted by τττA, τττB and:
τ
ττ0A,B =RθτττA,B+dez, (2.22)
where d is the distance between the layers in the z axis.
We will see that incommensurate structures with small enough rotation angles can be approximated by a similar commensurate lattices with large Moir´e unit cells.
2.2.2
Tight-binding models for TBLG
We want a model that allows us to describe the low energy Moir´e bands of bilayer graphene for both commensurate and incommensurate systems. 10
2.2 Twisted bilayer graphene (TBLG) 11
Recent observations demonstrate that for a discrete number of small ro-tation angles, the magic angles, the bandstructure of TBLG shows certain interesting features that can lead to a correlated insulating phase [2] and unconventional superconductivity [3]. The largest magic angle is roughly 1.08◦ and it is the one we will study in this Section. If we look at the elec-tronic band structure of commensurate magic angle TBLG in Fig. 2.4, we can notice two relevant features. There are four nearly flat bands near zero energy, which are a consequence of the interaction of the four K points. We can also observe two gaps, one above EF and one below EF. They are cre-ated because of the flattening of the bands, and they isolate them from the rest of the bandstructure.
Figure 2.4:Electronic band structure and density of states of a magic angle TBLG. The electronic bandstructure (a) is computed for a commensurate system with m = 1 and n = 61, which corresponds to θ = 1.08◦. We can see the flat bands around−10 meV and the two gaps, which have a size of approximately 10 meV. The Fermi energy is shifted from zero as a result of the model that we used. The total density of states (b) also shows the flat bands. It is computed for a finite circular sample with radius 500 ˚A.
Several theoretical attempts to reproduce the features of magic angle TBLG have been published in the recent years. The model that will be used throughout this thesis is the continuum tight-binding model. It was the first model that was used to predict the magic angle TBLG in [14] and it is presented with some new implementations in several studies [5][15][16]. There are other examples of theoretical approaches in the literature. In [17]
12 Graphene systems
and [18], they propose an effective tight-binding model that captures the flat bands while also retaining all the symmetries of the system. In [19] they use DFT calculations to reproduce the flat bands. In [20], they study the system when adding electron-electron interactions using the Hubbard model. We are interested in the continuum tight-binding model because it is a flexible model that can be used for both commensurate and incommen-surate rotation angles and systems having several layers. This matches our purposes because we want to study different types of systems. Moreover, it is a simple, straight-forward model that can be easily implemented with the tools presented in Chapter 4 while it still captures the essence of the system.
We are going to focus on the implementation of the continuum tight-binding model presented by Lin et al. [5]. The parameters for this model are chosen such that the results are consistent with experimental data from monolayer graphene and AB stacked untwisted bilayer graphene. Let us consider a tight-binding Hamiltonian containing an intra-layer (Hk) and
an inter-layer (H⊥) term: H= Hk+H⊥ = −
∑
m,<i,j> ti,jm,m(c†m,icm,j+h.c.) −∑
m,i,j tm,mi,j +1(cm,i† cm+1,j+h.c.), (2.23) where cm,i(c†m,i) are the creation(annihilation) operators of an electron in the atomic site i of the layer m (m = 1, 2). In the intra-layer Hamiltonian, Hk, we consider only nearest neighbor hoppings. The value thatrepro-duces the correct Fermi velocity for monolayer graphene is tm,mi,j =Vppπ0 =
3.09 eV.
To describe the inter-layer coupling we will first consider two atoms that are on top of each other, separated by the distance between the layers d. We will set the hopping between them to be Vppσ0 . For the neighboring atoms, we will to compute the distance between the two sites projected onto one of the layers r =|r|. Then, the hopping scales as follows:
t(r) = Vppσ0 e−(
√ r2+d2
0−d)/λ d0
r2+d2. (2.24)
The parameters Vppσ0 = 3.09 eV and d ≈ 3.35 ˚A are extracted from obser-vations in AB stacked bilayer graphene. In twisted bilayer graphene, the 12
2.2 Twisted bilayer graphene (TBLG) 13
distance between the layers is spatially modulated, but we will approxi-mate it to a constant in this model. Finally the cutoff λ =0.27 ˚A makes sure that the hopping decreases properly when the distances are long. To check the validity of eq. 2.24, we can compute the hopping between neighbors of adjacent layers for AA-stacked bilayer graphene, for which r= √a
3. We obtain a value of 0.11 eV, which is agreement with the experimental data [21].
The low energy wavefunctions for this model can be spanned in the Bloch basis of pz orbitals on the two different layers 1 and 2.
|Ψ1,α(k)i = √1 N
∑
Rα n eikRαn|Rα ni, (2.25) Ψ2,α(k0) = √1 N∑
Rα0 n eik0Rα0n Rα0 n , (2.26)where N is the number of unit cells of the system, k and k0 are the Bloch wave vectors and |Rα
ni is the atomic orbital at site Rαn. The intra-layer Hamiltonian matrix elements for the first layer are given by:
Ψ1,α1(k1) Hk1|Ψ1,α2(k2)i = 3
∑
j=1 −Vppπ0 e−ik1δδδjδk1,k2, (2.27)where δδδj are the vectors connecting the nearest neighbors in monolayer graphene from eq. 2.1. The expression is equivalent for the second layer but considering the rotated momentum and the corresponding nearest-neighbor vectors. When the rotation is small, the inter-layer Hamiltonian matrix elements are:
hΨ2,α2(k1)|H⊥ Ψ1,α1(k2) = −
∑
G,G0 ˜t(k1+G)e−iGτττα1+iG 0 τττα2 δk2+G,k1+G0, (2.28) where G = m1b1+m2b2 and G0 = m01b10 +m02b02 are reciprocal lattice vectors. Here, ˜t denotes the Fourier transform of the inter-layer coupling. Since t(r)is isotropic, it can be expressed as follows:˜t(k) =2π Z ∞
0 rt
14 Graphene systems
Figure 2.5:Radial dependence of the hopping from eq. 2.24 and its Fourier trans-form in terms of the momentum. Both decay rapidly to zero.
where J0 is a Bessel function and k = |k|. The radial and momentum dependence for t(r)and ˜t(k)respectively are shown in Fig. 2.5.
The inter-layer interaction only happens when the Umklapp scattering conditions holds, i.e., k2+G = k1+G0 [4]. This can be derived from
the inter-layer matrix elements in eq. 2.28. Therefore, the space of wave functions associated to a momentum k0 is spanned by{|k, αi | k =k0+ G0−G}and {|k0, αi |k0 = k0+G−G0}, ∀Gand ∀G0. Since the Fourier transform of the hopping ˜t(k) decays for a large value of|k|, the relevant contribution for a given k consists only of a set of small reciprocal lattice vectors. In [5], they consider a set of 27 G-vectors for each k point. The idea that we need to keep in mind is that the K and K0 points from one layer only interact with other nearby K and K0 points from the other layer [4].
2.3
Finite incommensurate systems
The tight-binding model presented in Section 2.2 is valid both for com-mensurate and incomcom-mensurate structures. However, there are some no-tions, like the primitive unit cell or the Brillouin zone, that do not exist 14
2.3 Finite incommensurate systems 15
for incommensurate systems because of the lack of periodicity. Since the Bloch theorem does not hold anymore, the study of the electronic spec-trum is more complicated. Instead of focusing on the electronic band structure, we will rely on other tools such as the density of states and the spectral function.
For a finite sample, the density of states ρ(W)tells us about the distri-bution of the energies in the momentum space.
ρ(E) = 1 D D−1
∑
k=0 δ(E−Ek) = Tr(δ(E−H)). (2.30)Here, H is the Hamiltonian, Ekare its eigenvalues and D is the dimension of the Hilbert space. The density of states can be computed in several ways. One option is to diagonalize the Hamiltonian of the system, which can be an expensive computation for large systems. There are other more effective tools like the kernel polynomial method, which is explained in Chapter 4.
Sometimes, the density of states will not give us enough information about the system. Then, the spectral function is useful to study the en-ergy modes as as we move around k-space. The spectral function is the projection of the density of states at specific values of k:
ρk(E) = hk|δ(E−H)|ki. (2.31)
When the rotation is small, we can find a commensurate angle that is close enough to our incommensurate rotation. Then, we can compute the spec-tral function through the Brillouin zone of the commensurate system. The electronic band structure and the spectral function provide similar infor-mation about the system. The relevant difference is that we can only cap-ture the feacap-tures that are a consequence of the lack of periodicity using the spectral function.
In the following Sections, we will study the electronic spectrum of dif-ferent types of finite incommensurate structures. But first, we will intro-duce the notion of quasicrystals in the next Chapter.
Chapter
3
Quasicrystals
Quasicrystals are structures of atoms that are neither crystalline nor amor-phous. They present symmetries, but at the same time they are not peri-odic. These structures have been raising interest in the field of condensed matter physics because they offer a new range of possibilities that can not be found in conventional crystals and that are yet unexplored.
Figure 3.1: Stampfli tiling [22]. Extracted from https://geometricolor.wordpress.com.
Before the discovery of quasicrystals, quasi-periodic systems were al-ready known and studied by mathematicians. One example are the quasi-periodic tilings. A tiling is an arrangement of non-overlapping tiles that covers the plane. Two examples of famous tilings that are not periodic are the Ammann-Beenker tiling, which was discovered in 1961, and the
18 Quasicrystals
Penrose tiling, introduced by Roger Penrose in 1974 [23]. Quasi-periodic tilings can be used to model quasicrystals. In particular, we will see that the Stampfli tiling, shown in Fig 3.1, is directly related to one of the sys-tems that we want to study.
In this Chapter, we will first introduce the notions of almost periodic functions and quasi-periodic tilings, presenting two main methods that can be used to construct them. We will study the Fibonacci chain as an example and then we will go back to condensed matter theory to introduce the definition of quasicrystals and the systems that are relevant for this project.
3.1
Almost periodic functions and quasi-periodic
tilings
The aim of this Section is to introduce the main mathematical properties of quasi-periodic systems [23]. An almost periodic function is defined as a superposition of harmonics:
f(x) =
∑
ncneiknx. (3.1)
All periodic functions are included in the definition, but there are also other types of almost periodic functions. For example, the sum of two co-sine functions with periods that have an irrational ratio is a superposition of harmonics that is not periodic:
f(x) = cos(x) +cos(τ−1x). (3.2)
In this case, we are considering the golden ratio τ = (1+√5)/2. Al-most periodic functions are very repetitive, since they alAl-most overlap with copies of themselves that are shifted through the space. As a consequence, for any e, we can always find an almost period xe such that:
|f(x−xe) − f(x)| ≤ e. (3.3)
Quasi-periodic tilings are discrete versions of almost periodic functions in N dimensions. They are also highly repetitive and almost periodic. 18
3.1 Almost periodic functions and quasi-periodic tilings 19
Therefore, we can also find copies of the same pattern overlapping every-where. There is a mapping between a periodic tiling and a quasi-periodic lattice. All we need to do is to consider that the apendexes of the tiles are the atoms and picture the edges as the bonds between them.
The simplest method to obtain an almost periodic function is to add two or more periodic functions having an irrational ratio between their periods. The resulting structure is called an incommensurate system. This same idea can be applied to tilings. There are several ways to add pe-riodic tilings to obtain quasi-pepe-riodicity. In particular, we can add two-dimensional periodic tilings and rotate them by an incommensurate an-gle. This idea is equivalent to the graphene systems described in Section 2.2. As we will see, there are more elaborated ways to build quasi-periodic systems.
3.1.1
The canonical cut and projection method
In general, we can obtain almost periodic functions by considering a pe-riodic function in a higher dimension and then taking a projection of it. This also holds for quasi-periodic tilings and it is one of the main methods that is used to construct them. The canonical cut and projection method is described as follows [23]:
• First, we need to choose a D-dimensional hypercubic lattice living in a D-dimensional space E.
• Then, let us take a d-dimensional subspace of E, E0 ⊂ E. There are many subspaces in E. The choice of the subspace is often reffered to as the slope of the tiling.
• Now, we will select all the lattice points that are inside the infinite slice of E which is obtained by translating the unit cell of the hyper-cubic lattice along the subspace E0.
• Finally, we will project all the points into E0, obtaining a quasi-periodic lattice.
20 Quasicrystals
One example is the Fibonacci chain, which can be obtained using the canonical cut and projection method for D = 2 and d = 1. The slope of the subspace, which is the slope of the line cut in the two-dimensional cubic lattice, is the inverse of the golden ratio τ−1. One direct way to see that the Fibonacci chain is not periodic is to notice that the slope of the projection line is an irrational number. Since the lattice constant is 1, only the origin lattice point will fall exactly in the line, and that is why the sequence cannot be periodic. If the slope had been a rational number α =
p/q instead, the point with coordinates (p, q) would have fallen exactly in the line. We would have a sequence with period p+q. The Fibonacci chain, shown in Fig. 3.2, can be seen as a lattice with two possible bonds between the atoms that alternate in a quasi-periodic fashion.
Figure 3.2: Fibonacci chain obtained using the canonical cut and projection method. It can be represented as a one-dimensional chain of atoms with two types of bonds.
20
3.2 Quasicrystals in condensed matter 21
3.2
Quasicrystals in condensed matter
Quasicrystals in condensed matter were discovered in 1982 by Dan Shec-htman [1]. He was analyzing a sample of AlMn when he measured Bragg peaks corresponding to a 10-fold symmetric structure. Bragg peaks indi-cate the existence of long-range order, but the 10-fold symmetry is forbid-den in periodic crystalline structures. At that time, these were thought to be the only type of ordered structures in condensed matter. The crystal-lographic restriction theorem states that in a periodic lattice only a finite number of rotation symmetries are possible. In particular, only 2-fold, 3-fold, 4-fold and 6-fold symmetries are allowed. They found out that the AlMn sample had a configuration that was similar to the Penrose tiling and they classified it as a quasicrystal. Short after this discovery, new quasicrystals with 8-fold [24] and 12-fold [25] rotation symmetries were discovered. Other quasicrystalline structures were found in soft matter systems [26] or have been artificially engineered [27].
There is some debate regarding the definition of quasicrystals [28]. Some of the literature requires quasicrystals to have a rotational symmetry that is incompatible with periodicity, i.e. a non-crystallographic symme-try. If we follow this definition, there would be no one-dimensional qua-sicrystals. Moreover, any incommensurate structure obtained by stacking two or more periodic systems would not necessarily be considered a qua-sicrystal. The other option is to define quasicrystals as any aperiodic or-dered structure. Then, we define extrinsic quasicrystals as systems which are composed of two or more periodic structures. The systems that can-not be described in those terms will be called intrinsic quasicrystals. The main difference between these two types of crystals is that, in general, ex-trinsic quasicrystals can be described as a periodic system that has been slightly modified. This is because the underlying structure is still built on periodic subsystems. As we will see, the effectiveness of this approx-imation depends on the coupling between the subsystems. For intrinsic quasicrystals instead there is no possible connection to periodicity. These are the systems that contain non-crystallographic symmetries and a new approach must be used to study them.
22 Quasicrystals
reciprocal space. The density function of a periodic crystal in real space can be expressed in terms of the reciprocal lattice vectors G as follows [29]:
ρ(r) =
∑
G
ρ(G)eiGr. (3.4)
In a periodic crystal, any vector connecting sites in the lattice can be writ-ten as an integer linear combination of the reciprocal lattice vectors. In three dimensions:
G =hb1+kb2+lb3. (3.5)
Therefore, the Bragg peaks would be indexed by three integers, (h, k, l). The situation is different when we are considering quasicrystals (both ex-trinsic and inex-trinsic). Then, the number of linearly independent vectors that are required to span the reciprocal space exceeds the physical dimen-sion of the system. The first consequence of this fact is that the reciprocal lattice is not a discrete lattice anymore, but instead it consists of a dense set of points. This will cause the Bragg peaks to be indexed by a number of integers D which exceeds the number of spatial dimensions d. We can picture an abstract higher dimensional space to describe the positions of the Bragg peaks.
As a consequence of the reciprocal space being dense, the electronic spectrum of quasicrystals presents interesting features that are not found in periodic systems. We will focus on two signs of quasicrystallinity: frac-tality in the electronic spectrum and localization of the eigenstates.
3.2.1
Fractal electronic spectrum
The electronic modes of a quasicrystals are complicated. Since the recip-rocal space is dense, plane waves satisfy the diffraction condition for most of the k-vectors:
2G·k± |G|2 =0. (3.6)
Many different modes, both propagative and non-propagative, coexist in a quasicrystal, making the electronic spectrum a lot more richer. The re-arrangement of the energy bands leads to a spectrum containing many moderately flat bands at different energies, each band containing at least 22
3.2 Quasicrystals in condensed matter 23
one saddle point. This is translated to the density of states having many van Hove singularities and pseudogaps, which would be captured by the fractal dimension of the profile.
Figure 3.3:Spectral function of the tight-binding model of the Fibonacci chain. It clearly shows fractal behaviour. This result has been computed using the tools explained in Chapter 4. To better appreciate the features of the spectral function we have considered a color scale which is not linear and therefore it is not repre-sentative of the intensities.
Even though there are several theoretical approaches to describe the fractality of such electronic spectra, there is no clear idea of which is the information that can be extracted from it. Most of these approaches focus on intrinsic quasicrystals. Some examples are the study of the spectrum and the eigenstates of a Fibonacci chain [30],[27] or a Penrose photonic quasicrystal [31]. In [32] they study several intrinsic quasi-periodic chains and tilings, computing their multi-fractal spectra to show that they have critical eigenstates.
We believe that the quasicrystallinity is captured by the fractal dimen-sion of the density of states both in intrinsic and extrinsic quasicrystals. For extrinsic quasicrystals the signs of quasicrystallinity are more relevant when the coupling between the subsystems is greater. The goal of this the-sis is to show that the fractal dimension of the density of states is propor-tional to the quasicrystalline strength of system. In the following Chapters,
24 Quasicrystals
we will classify several systems in terms of the fractal dimension of their electronic spectra.
3.2.2
Localization of eigenstates
The quasicrystalline nature of a sample can also be confirmed by study-ing the localization of its eigenstates. When a system has translational symmetry, the eigenstates can be expressed according to the Bloch theo-rem and they are extended throughout the lattice, as shown in Fig 3.4 (a). On the contrary, if we consider a disordered system or an Anderson insu-lator, states are highly localized. Quasicrystals are believed to have low conductivity [33], which leads to the existence of some localized states as the one shown in 3.4 (b). Some studies for the one-dimensional and two-dimensional cases show that delocalized, localized and critical states coexist in quasicrystals [34], [35], [36]. The localization in quasicrystals is usually explained by the lack of periodicity and the existence of resonant states [4]. However, a good understanding of these type of states has not been found yet.
One way to characterize localization in a system is through the inverse participation ratio (IPR). The IPR for an eigenstate ψ is computed as fol-lows:
P−1(ψ) = ∑i|ψi|
4
∑i|ψ|2i
2, (3.7)
Where ∑i indicates the sum over all sites of the system [33]. The IPR is inversely proportional to the extension of a state. Systems with highly localized states will have a constant, large IPR. For metals we expect the IPR to decrease as O(1/N), where N is the number of states, as we increase the size of the sample [37]. We can detect quasicrystals because we expect them to show a slightly different behaviour.
Another measure of quasicrystallinity is given by the energy level statis-tics. The spacing ratios between the energy levels (the eigenvalues of the Hamiltonian) are given by:
rn =min(δn, δn+1)/max(δn, δn+1), (3.8) 24
3.3 Analysis of graphene quasicrystals 25
Figure 3.4: Local density of states of eigenstates of two samples of TBLG with radius . In (a) the rotation angle is commensurate and therefore the system is periodic. The eigenstate is extended and follows a periodic fashion. In (b) the rotation angle is incommensurate and therefore the system is not periodic. The state is seems to be slightly localized at the centre.
where δn is the level spacing between the (n+1)-th and the n-th energy eigenvalues. When the states of the system are delocalized, the energy lev-els repel from each other, leading to Gaussian orthogonal ensemble dis-tribution [33]. When the eigenstates of a system are localized, they are distributed randomly and therefore follow a Poisson distribution.
3.3
Analysis of graphene quasicrystals
In this Section we will connect our knowledge of quasicrystals and twisted bilayer of graphene. Incommensurate TBLG systems are extrinsic qua-sicrystals. They are interesting because they have already been artificially engineered and offer a great range of possibilities depending on the rota-tion angle and the number of layers. Moreover, monolayer graphene has a peculiar low-energy electronic spectrum that gives rise to interesting fea-tures in multi-layer systems. Intrinsic quasicrystals, on the contrary, are much more difficult to build or to find in nature, and their study demands
26 Quasicrystals
a completely new approach. We will use our knowledge of monolayer graphene to decide which type of systems are expected to show interest-ing features. The general idea is to look for systems in which the couplinterest-ing between the layers is large enough so that the quasi-periodicity becomes relevant.
Most of the research in twisted bilayer graphene has been focusing on the magic angle TBLG introduced in Section 2.2. After the discovery of correlated insulating [3] and unconventional superconductivity [2] the condensed matter community became really interested in these type of systems. But why is this system interesting when it comes to its quasicrys-tallinity?
At the low-energy range, we only need to focus on the Dirac cones. Since the magic angle is small ( 1.08◦), the K and K0points of one layer will be close to the K and K0of the other layer. Recall that the coupling between points in k-space is relevant only for nearby points because the Fourier transform of the hopping parameter decays with the distance. From this, we can conclude that in this case the coupling between the Dirac cones is strong. Nevertheless, turns out that since the rotation is small, we can always find a commensurate angle which is close enough to the magic an-gle. The system is effectively described by a commensurate TBLG, which is periodic. In this case, the Moir´e pattern governs the low-energy physics and since the Moir´e period is really large, it is difficult to observe quasi-periodic effects. For a really large system we would expect to start detect-ing some quasi-periodic effects, but otherwise we only see the commensu-rate approximation of the system.
Another system that has been studied in the recent years is the 30◦ TBLG (dodecagonal graphene). This system has raised special interest because it has a non-crystallographic symmetry combined with a mirror symmetry that exchanges the two layers: 12-fold + M. Samples of these type of systems have been artificially produced by some research groups [38], [39]. In [39], they show that if we project the two layers into one plane, the structure can be mapped into a Stampfli tiling (Fig. 3.1), which has a 12-fold rotational symmetry. This creates a direct connection between quasicrystals and quasi-periodic tilings, even if the TBLG is an extrinsic
26
3.3 Analysis of graphene quasicrystals 27
Figure 3.5: Lattice structures corresponding to magic angle TBLG (a), dodecago-nal graphene (b) and TBLTBLG (c). Note that magic angle TBLG has a 6-fold (+ mirror) symmetry, whereas dodecagonal graphene and TBLTBLG have a 12-fold (+ mirror) symmetry.
quasicrystal.
The angle of rotation in dodecagonal graphene is large. As a conse-quence, the K and K0 points of the two layers are far apart and the cou-pling between them is weak. This means that, at low energies, we will not be able to see interaction between them and instead we will observe the behaviour of the two layers separately. In this case the Moir´e period competes with the atomic length scale and therefore the approximation to a commensurate structure is not valid anymore.
In [4], P. Moon et al. propose a model that studies the electronic struc-ture of dodecagonal graphene through the inter-layer interaction between the Dirac cones. They show that that the quasicrystalline structure gives rise to peaks and pseudo gaps in the density of states at a range of ener-gies between−1 and−3 eV. They attribute these peaks to resonant states which appear because of the hybridization of the graphene Bloch states at the twelve K and K0 points situated at the corners of the two BZ. These Bloch states are degenerate in energy due to the symmetry.
Even though we detect some effects of quasi-periodicity, the electronic spectrum of this system is not fractal. Our goal is to increase the value of inter-layer coupling enough to compensate for the large rotation angle, so
28 Quasicrystals
that we can start measuring a fractal electronic spectrum. The downside is that the system then will not be realistic anymore, since we are considering a coupling that is much larger to that of real bilayer graphene.
Finally, we propose a new structure that has not been widely studied: twisted bilayer twisted bilayer graphene (TBLTBLG). This structure, intro-duced in [40], consists of four layers of graphene. Layers 1−2 and 3−4 are rotated by the magic angle 1.08◦ and the twist angle between layers 2−3 is 30◦. The exact sequence of rotation angles is 1.08◦−31.08◦−30◦ with respect to the first layer. This new system is interesting because it brings together the large coupling between nearby Dirac points from magic angle TBLG but at the same time it cannot be described using a Moir´e effective theory because the approximation is broken by the 30◦ ro-tation. For these reasons, we expect to find strong a large fractal dimension at a low-energy range, which is where flat bands will be.
In the following chapters we will compare and classify the three de-scribed systems by studying their density of states, but first, we need to introduce the methods that will be used to obtain the results.
28
Chapter
4
Methods
The goal of this Chapter is to introduce the methods and approximations that have been used to simulate the tight-binding systems under study. In Section 4.1, we will explain how to compute the fractal dimension of a one-dimensional profile. In Section 4.2, we will focus on the simulations, which are an essential part of this project. First, we will explain how to make the computations more effective by tuning the parameters of the system. Then, we will focus on the kernel polynomial method, which is the main tool that we have used to compute the electronic spectra.
4.1
The fractal dimension
There are several measures to determine if a surface is fractal-like. Most of these measures are based on the two main properties of fractals, which are roughness and self-similarity. We will focus on the Hausdorff dimension, which is a measure of how much does a curve fill up the space as we increase its resolution. This idea is equivalent to studying the complexity of the curve as we change the scale. Therefore, both self-similarity and roughness are captured in the Hausdorff dimension.
30 Methods
e-cover of X is defined as a finite or countable collection of balls Bi ∈ RD with diameter|Bi| ≤ e, such that they cover X [41]. Then, the δ−dimensional
Hausdorff measure of X is defined as: Hδ(X) = lim e→0inf ∞
∑
i=1 |Bi|δ : {Bi : i=1, 2...}is an e−cover of X . (4.1) There exists a unique positive value D such that Hδ(X) = ∞ if δ < Dand Hδ(X) = 0 if δ > D. D is the Hausdorff dimension of X. For a
d-dimensional smooth surface, the Hausdorff dimension will also be d. However, for a rough surface it will take values between d and d+1. The Hausdorff dimension of a pattern does not uniquely describe the pattern. There are several approaches to compute the Hausdorff dimension, such as the information dimension, the correlation dimension or the R´eni di-mension. In this work, we will use the box-counting dimension or also called Minkowski-Bouligand dimension.
4.1.1
Box-counting method for 1D curves
To measure the fractal dimension of a one-dimensional curve E using the box-counting algorithm, we will start by considering a sequence of lengths,
{ek}, such that limk→∞ek = 0. For each length ek, we will compute the number of two-dimensional boxes N(E, ek)with side ekthat are needed to cover E. Then, the fractal dimension is given by [42]:
D(X) = lim k→0
ln N(E, ek) ln(1/ek)
, (4.2)
which can be approximated by the slope of the line going through
(ln N(X, ek), ln(1/ek)). (4.3) In practice, we want to study the density of states of a system, which is a discrete curve X= {(x1, y1),(x2, y2), ...}. We have used a slightly different implementation of the algorithm, which consist of the following steps:
• Consider a sequence of resolutions R = {r1, r2...} such that rk is an integer for all k. Then, we can define X(rk)to be the curve at a res-olution rk, which means taking only rk equally spaced points from X.
30
4.1 The fractal dimension 31
• Then, ek is defined as L/rk, where L is the length of the curve. The parameter N(X, rk)will be computed as follows:
N(X, rk) =
∑
X0(rk)
, (4.4)
where X0(rk) is the derivative of the X(rk). In order to improve the computation, we will take the mean of N(X(rk), ek)for different sets of rkequally spaced points.
• The fractal dimension is 1−m where m is the slope of the line going through:
(ln N(X, rk), ln(L/rk)). (4.5)
Figure 4.1: Example of how to compute the box-counting dimension. This is the log-log plot of N(ek, DOS)as a function of ek for TBLG when the inter-layer
coupling between layers is V = 6. It will become clear later why we do not use the first data points to do the linear fit. D is the box-counting dimension.
The box-counting dimension does not always coincide with the Haus-dorff dimension of a curve [42]. Nevertheless, since our goal is to classify systems using their fractal dimension, the box-counting algorithm should be enough. We want our method to act equivalently in different systems, and this raises other problems that we need to address:
• The fractal dimension of the density of states should not depend on the size of the system. Having a larger system means having more
32 Methods
states in the density of states, which could increase the box-counting dimension.
• Having really large peaks in the density of states, which correspond to the van Hove singularities, can affect the computation of the frac-tal dimension.
To deal with these two problems, we will rescale and normalize the den-sity of states in the following way:
ρ0(E) = arctan ρ(E) N , (4.6)
where N is the total number of states. Applying this transformation to the density of states is interesting because it also allows us to easily detect when the resolution is too large and we are resolving individual peaks. When this happens, the function will converge to the constant 2N. On the contrary, this transformation means that we will not be able to differentiate between 1st and 2nd order van Hove singularities.
4.2
Simulations
The simulations for this project have been carried out using Python, and in particular, Kwant [11]. Kwant is a Python package for numerical quantum transport simulations which offers a straightforward way to work with tight-binding models. In this section we will comment on how to build the systems and effectively compute their properties.
4.2.1
Increasing the rotation angle by tuning the inter-layer
potential
In magic angle TBLG the contribution to the low-energy flat bands is around one electron per approximate Moir´e unit cell. This means that if we want to capture the behaviour at a low-energy range we would need approxi-mately 1000 data points in that range, which corresponds to a system with 32
4.2 Simulations 33
Figure 4.2:(a) Color map showing the Fermi velocity at the Dirac point as a func-tion of the inter-layer coupling V and the rotafunc-tion angle θ of TBLG. This is com-puted using the electronic band structure of the commensurate system. (b) Color map showing the size of the largest gap around the flat bands of TBLG as a func-tion of the inter-layer coupling and the rotafunc-tion angle θ. In general, the two gaps below and above the flat bands are almost symmetric. This is computed using the density of states of the system. The black dotted line shows the line with slope 2.5 and the red dotted lines show the combinations(θ, V) = (1.08◦, 0.39eV)and
(θ, V) = (3◦, 1eV). The behaviour below the black dotted line may be affected by
the shift of the Fermi energy.
1000 Moir´e unit cells. The number of Moir´e unit cells in a circular system of magic angle TBLG with radius 750 ˚A is around 104. Such a system has roughly 1000000 atoms and it is already quite heavy to simulate.
If we increase the rotation angle, the contribution to the flat bands also increases. This is because area of the Moir´e unit cell is smaller for larger angles. In order to increase the rotation angle without affecting the system, we also need tune the inter-layer coupling V. Recall that the correct value of the inter-layer coupling for graphene is 0.39 eV and the magic angle is 1.08◦. Fig. 4.2(a) shows the value of the Fermi velocity as a function of the rotation angle and the inter-layer coupling. If we now turn to Fig. 4.2(b), we can see the size of the gaps around the flat bands as a function of the rotation angle and the inter-layer coupling. By looking at the two results, it can be seen that if the ratio between θ and V is approximately 2.5 the behaviour of the magic angle TBLG is maintained. Up to a scale factor, the
34 Methods
Moir´e bands depend on the parameter α [14], where:
α= V
vFkθ
. (4.7)
Here, kθ = 2K sin θ
2 and vF is the Fermi velocity of monolayer graphene. In this case, θ is small and therefore kθ ∝ θ.
The electronic band structure for the system with a rotation angle θ =
3◦and inter-layer coupling V =1 eV is shown in Fig. 4.3. If we compare it with the electronic band structure of magic angle TBLG shown in 2.4, we see that the two are roughly equivalent. For this new combination of pa-rameters, a circular system with radius 750 contains 735 Moir´e unit cells, which is seven times greater than the value that we obtained for magic angle TBLG. For this reason, we have chosen to perform our simulations using this combination of parameters. Since we are not focusing on the magic angle properties but instead we are interested in the strong inter-layer coupling, this trick will allow us to better observe quasi-periodic ef-fects.
Figure 4.3:Electronic band structure and density of states of TBLG for V =1 eV and θ = 3◦. The electronic bandstructure (a) is computed for a commensurate system with m =1 and n = 23. We can see the flat bands around−60 meV and the two gaps, which have a size of approximately 25 meV. The Fermi energy is shifted from zero as a result of the model that we used. The total density of states (b) also shows the flat bands. It is computed for a finite circular sample with radius 500 ˚A.
34
4.2 Simulations 35
4.2.2
The kernel polynomial method
In order to obtain the electronic spectrum of a finite system we need to diagonalize its tight-binding Hamiltonian. The electronic spectrum is a measure of how the eigenvalues are distributed. In a system with N sites, the matrix that we need to diagonalize is N×N. Directly diagonalizing the Hamiltonian would mean, in general, that we need a memory of order N2and that the computation time scales as N3[43]. Since we are working with only nearest neighbors most of the entries of the Hamiltonian are be zero. Diagonalizing this type of matrices, which are called sparse matrices, is easier, but still we are interested in finding alternative approaches.
We will use the KPM method, which is already built in the library Kwant and it can be used to obtain the density of states and the expec-tation values of operators. The are two approximations inside the KPM implementation in Kwant. The first one is to use the stochastic trace esti-mation, which basically means expanding eq.2.30 on a small set of random vectors as follows:
ρA(E) ∼ 1
R
∑
r hr|Aδ(E−H)|ri, (4.8) where ρ is the density of states and A is any operator. The number of random vectors is R. Then, if we are interested in the local density of states we will switch the random vectors for a site i:ρiA(E) ∼ hi|Aδ(E−H)|ii, (4.9)
and if, instead, we want to compute the spectral function, we will use the plan waves corresponding to the vector k at which we want to project the density of states:
ρA(E) ∼ hk|Aδ(E−H)|ki. (4.10) The second approximation is using the Chebyshev expansion. In gen-eral, any function can be expressed as a sum of orthogonal polynomials [43]. In our case, we will be interested in the following expression:
f(x) = 1 π √ 1−x2 µ0+2 ∞
∑
n=1 µnTn(x) , (4.11)36 Methods
where Tn(x)are the Chebyshev polynomials of order n:
Tn(x) =cos(n arccos(x)), (4.12) and µn are the expansion coefficients (or also called moments):
µn =
Z −1
−1 f
(x)Tn(x)dx. (4.13)
Our goal is to find the best approximation to the density of states ρ(E)
by choosing the right the number of random vectors and the order of the Chebyshev expansion. Since we are truncating the infinite expansion of polynomials into a finite number, the result will show fluctuations. For this reason, the convolution of the function with a kernel is implemented in the algorithm. Then,
fKPM(x) = Z 1 −1π q 1−y2K N(x, y)f(y)dy. (4.14) There are many types of kernels that can be used for this purpose. The two main ones are the Jackson kernel and the Lorentz kernel. They have similar properties but the Lorentz kernel is more appropriate for our com-putations because we are interested in dumping the oscillations as much as possible. Otherwise, they would affect the computation of the fractal dimension. The effect of the Lorentz kernel is similar to making a convo-lution of our function with a Gaussian. This can lead to poor convergence at the boundaries, but it will not be a problem because we can get rid of them.
The resolution is directly related to number of moments used in the Chebyshev expansion as R = 1/M. The number of moments should not be too large because we do not want to resolve individual states, but we want it to be large enough to resolve the fractal behaviour. Usually, we will choose the resolution to be around:
R=10max(E) −min(E)
N , (4.15)
where (max E−min E) is the bandwidth and N is the number of states. 36
4.2 Simulations 37
Moreover, having a large number of random vectors makes the approxi-mation better and less noisy. If the system is large they are not needed.
Another thing to take into account is that the implemented KPM func-tion in Kwant outputs an array of energies that is not equally spaced, but it is the best fit to avoid unreal oscillations. The smoothing resolution for the Lorentz is around 4/M, where M is the number of moments. The max-imum width of the default energy intervals given by KPM is 2π/N. That is why, as we will see in the following Chapter, we will need to discard some resolutions used in the computation of the fractal dimension.
Chapter
5
Results
In this Chapter, we are going to present the final results of this thesis. The results are organized in two main Sections. In Section 5.1, we are going to study two one-dimensional quasicrystals: an incommensurate chain and the Fibonacci chain. In Section 5.2, we will turn to the graphene systems that were introduced in Section 3.3 and study how relevant is the quasi-periodicity in their electronic spectrum. We will focus on the box-counting dimension of the density of states of these systems. To complement our results, we will also look at the localization of the eigenstates.
5.1
Analysis of 1D quasicrystals
Our goal in this Section is to establish a reliable method to compute the fractal dimension of tight-binding systems. One-dimensional systems are interesting because, in general, they are simple and easier to model than higher order structures. Moreover, the Fibonacci chain is a well-known quasicrystal. This makes it a great starting point for our purposes. It is important to highlight that in this Section we are not looking for realistic systems, but instead modelling abstract systems that can help us under-stand better what causes the electronic spectrum to become fractal.
40 Results
Figure 5.1: Spectral functions of the incommensurate chain (a) and Fibonacci chain (b) as we increase the inter-layer coupling in (a) and the ratio between the hopping parameters in (b). The color maps correspond to V = 0, 1, 2 for the in-commensurate chain and K = 1, 0.5, 0 for the Fibonacci chain. Note that the color scale is chosen to enhance the visibility of the bands and it is not representative of the intensities.
5.1.1
Incommensurate chains
Let us consider a large sample (∼100000 sites) of two coupled one-dimensional chains with lattice constants 1 and the golden ration τ. We will take the intra-layer potential between nearest neighbors to be the inverse of the distance between them and we will define the inter-layer potential as fol-lows:
t(r) = Ve−r/e, (5.1)
where V is the inter-layer coupling and r is the distance between the two sites. Since we are not considering a fixed distance between the two layers, the sites can get arbitrarily close. That is why we have used a hopping that depends on the exponential of the distance and not directly the distance. The parameter e = 0.1 has been chosen to enhance the fractality of the spectrum. We will set a fixed number of three inter-layer neighbors for each site.
40
5.1 Analysis of 1D quasicrystals 41
Figure 5.2: Cumulative density of states and density of states for the incommen-surate chain when V = 0 (a)(b) and when V = 2 (c)(d). In (e) we see the box-counting dimension of the system as V increases from zero to two. We reach a maximum value around 1.45.
To analyze the fractal dimension of the system we will look at its evo-lution as we vary the inter-layer coupling. In in Fig. 5.1(a) we show the spectral function of the incommensurate chain as we increase the coupling between the two chains from V =0 to V = 2. When the coupling is zero, we see the two distinct dispersion relations with different periods. There is no quasi-periodicity because we are looking at isolated periodic systems. As the coupling is increased, resonant states appear in the system and it starts to look fractal. If we turn to the density of states, shown in Fig. 5.2, we can see that for zero coupling it is a smooth function that describes the sum of two distinct chains. When we increase the coupling to V = 2 it becomes a rough curve with many peaks and gaps. This behaviour can also be captured in the cumulative density of states, which looks like the Devil’s staircase or Cantor function [23]. The gaps are translated to flat
42 Results
segments and the peaks correspond to the steps. Nevertheless, the fractal dimension of the cumulative density of states would not give us informa-tion about the system.
In Fig. 5.2 (d) we show the evolution of the box-counting dimension as the inter-layer coupling is increased. The box-counting dimension in-creases monotonically until it reaches a value around 1.45, and then is stays more or less stable. From this we can deduce that this system reaches its maximum fractality when the inter-layer coupling is around 0.75 and then it maintains it as we increase the inter-layer coupling to 2. The initial box-counting dimension is 1.1 instead of 1, which would correspond to a smooth one-dimensional curve. We attribute this to noise coming from the KPM approximation. In Fig. 5.4, we show the log-log plot for the com-putation of the fractal dimension, in this case for V = 2. The linear fit characterizes the box-counting dimension.
We can complement our understanding of the system by having a look at the inverse participation ratio (IPR) (3.7) as a function of the radius of the sample. As mentioned previously, for a periodic system, the IPR should decrease as O(N). In Fig. 5.5(a) we show the IPR for different values of the inter-layer coupling V. When V = 0, there is aO(1/N) be-haviour. Nevertheless, this also seems to be the behaviour when the chains are coupled, and this could be a sign that the incommensurate chain is not a strong quasicrystal. Moreover, we have used a relatively small number of eigenstates (100) to compute the IPR. It might be that we are not com-puting them at the optimal part of the spectra.
5.1.2
The Fibonacci chain
Recall the Fibonacci chain described in Section 3.1. The distances between the sites on the chain are either 0.52 ˚A or 0.85 ˚A and they alternate in a quasi-periodic fashion. Let us picture it as a tight-binding one-dimensional chain with two types of hopping parameters. We will set the hopping pa-rameters between nearest neighbors to be equal to the distances between them and the length of the chain to be ∼ 100000 ˚A. When we force the two hopping parameters to be one, we obtain the Hamiltonian of a regu-42
5.1 Analysis of 1D quasicrystals 43
Figure 5.3: Cumulative density of states and density of states for the Fibonacci chain when the hopping parameters are set equal (a)(b) and when the hopping parameters correspond to the real values of the chain (c)(d). In (e) we see the box-counting dimension of the system as we increase the ratio between the hopping parameters. We reach a maximum value around 1.53.
lar periodic chain. We will study the fractal dimension of the system as we gradually increase the ratio between the hopping parameters until we get back to the ones that correspond to Fibonacci chain:
t(r) =r+K(1−r), (5.2) where r is the distance between the sites and K goes from 0 eV(Fibonacci chain) to 1 eV (regular periodic chain). In Fig. 5.3(a), we can see that for equal hopping parameters the density of states is smooth and the disper-sion relation corresponds to an ordinary one-dimendisper-sional chain. When we increase the ratio between the hopping parameters, the density of states becomes a rough self-similar function. Again, the cumulative density of states looks like the Devil’s staircase and agrees with the results from [27].
