Loading...

First-Principles Prediction of Phase Stability of Electronic Materials

Master's Thesis 2013 162 Pages

Physics - Other

Excerpt

Contents

Zusammenfassung

I. Introduction

II. Theory

1. Basic Concepts: Density Functional Theory
1.1. Hohenberg-Kohn theorem
1.2. Kohn-Sham approach
1.3. Exchange correlation functionals
1.4. Van der Waals interactions in DFT

2. Phonons

3. Computational Aspects of DFT
3.1. FHI-aims
3.2. Phonopy

4. Ab Initio Atomistic Thermodynamics
4.1. Constrained Thermodynamic Equilibrium
4.2. Surface Free Energy
4.3. SiC Limits

III. Diamond Graphite Coexistence

5. Carbon
5.0.1. Graphene
5.0.2. Graphite
5.0.3. Diamond

6. Determination of the potential energy surface (PES) of Diamond and Graphite
6.1. Structure optimization
6.1.1. Minimal Energy for Fixed Volume Murnaghan Fit:
6.2. Temperature Dependent Energy Surfaces

7. Coexistence Line

IV. Surface Energies

8. Epitaxial Graphene Grown on Silicon Carbide
8.1. Silicon Carbide
8.1.1. SiC Polytypes
8.2. Surface Investigation
8.3. Case Study: C Terminated 1 × 1 Surface
8.4. (2 × 2) C Reconstruction
8.5. (3 × 3) Reconstruction
8.6. Graphene-Like Interface Structures

V. Conclusion

Acronyms

Bibliography

A. Covergence Tests
A.1. Free Atoms
A.2. Bulk Structures
A.2.1. k -grid

B. Bulk reference data

C. Reconstructions Models

Zusammenfassung

Mit der experimentellen Entdeckung der einzigartigen elektronischen Struktur von Graphen und der daraus folgenden hohen Elektronenmobilität beginnt eine neue Ära, deren Auswirkung nicht nur für die Physik und Nanoelektronik von Bedeutung sind. Im Gegensatz zu anderen Herstellungsverfahren hat die Graphitisierung der Polflächen von Siliziumkarbid (SiC) den wesentlichen Vorteil, dass es große Gra- phenlagen direkt auf einem Halbleiter bereitstellt. Vor diesem Hintergrund wird dieser Ansatz aktiv von vielen Forschungsgruppen wegen der Eignung für neuar- tige elektronische Geräte verfolgt. Da Graphen ein zweidimensionaler Kristall ist, muss seine Umgebung jedoch sorgfältig untersucht werden. Die höchste Ladungs- trägermobilität wurde bisher für freistehendes Graphen erreicht. Auf SiC ist es von besonderer Relevanz, die Interface-Struktur zu untersuchen und zu verstehen, wie diese auf die elektronischen Eigenschaften der Graphenschicht Einfluss nimmt. Für auf der Si-Seite von SiC gewachsenes Graphen ist die Schnittstelle bereits recht gut verstanden. Die erste Kohlenstoffschicht wechselwirkt stark mit dem Substrat und bildet kovalente Bindungen. Sie wirkt als elektronische Pufferschicht zwischen dem SiC und den darüberliegenden Graphenschichten, so dass die darüber liegen- de Kohlenstoffschicht niederenergetische elektronische Graphen-Eigenschaften auf- weist. Auf der Kohlenstoffseite des Substrats ist die Situation hingegen unklar. Di- verse Studien liefern scheinbar wiedersprüchliche Ergebnisse, die es zu vereinbaren gilt. Diese Arbeit zielt darauf ab, einen Überblick über die unterschiedlichen Konfi- gurationen auf dieser Oberfläche zu geben und Informationen über die Schnittstelle zwischen Graphen und SiC mit Hilfe von ab initio -Rechnungen zu erhalten. Hiermit lassen sich Strukturmodelle verifizieren und klare Aussagen über die Stabilität aller Konfigurationen vornehmen.

Zu Beginn wird der theoretische Hintergrund für das weitere Verständnis vermittelt. Darauf folgt eine ausführliche Untersuchung der Diamanten-Graphit-Koexistenz um das Stabilitätsverhalten im Kohlenstofflimit zu verstehen, welches entscheidend für das Wachstum von Graphen ist. Es wird der Beweis erbracht, dass trotz des in- dividuellen Charakters der Graphit- und Diamantbindungen, DFT und das Finite Displacement-Verfahren zahlreiche Eigenschaften beider Phasen in qualitativer und quantitativer Übereinstimmung mit experimentellen Daten liefern. Im Anschluss wird demonstriert, dass ab initio -Rechnungen auch in Hinblick auf die Vorher- sage der Stabilität von Oberflächenstrukturen ein mächtiges Werkzeug sind. Mit der vorgestellten Methodik wird die bekannte (2 × 2) C Struktur verifiziert und eine Vielzahl vorgeschlagener Modelle für die experimentell beobachtbare (3 × 3) Rekon- struktion falsifiziert. Indes wird ein von der Si Seite des Substrats adaptiertes (3 × 3) twist -Modell als die wahrscheinlich stabile Struktur vorgeschlagen. Die dargestell- ten Ergebnisse lassen dabei kaum Spielraum offen, dass es sich um eine andere als die genannte Struktur handelt. Zu guter Letzt wird gezeigt, dass sowohl kürz- lich experimentell beobachtete als auch von Si-Seite adaptierte, graphenähnliche Interface-Strukturen auf dem nackten Substrat instabil sind.

Part I. Introduction

A whole new era for physics as well as nanoelectronics (Geim and Novoselov, 2007; Geim, 2009) started with the experimental discovery of graphene’s unique electronic structure (Novoselov et al., 2005a; Zhang et al., 2005) and high electronic mobility (Du et al., 2008). In contrast to other fabrication techniques like mechanical exfoli- ation from a graphite crystal (Novoselov et al., 2005b) or catalytic decomposition of hydrocarbons on transition metals (Wintterlin and Bocquet, 2009), graphitization of the polar faces of silicon carbide (SiC) has the advantage of providing large- scale graphene sheets, directly on an insulating substrate (Berger et al., 2006a). The production of epitaxial graphene on SiC is therefore being actively pursued by many research groups because of the potential application of graphene for novel electronic devices (Hass et al., 2008b). Since graphene is an ultimate thin crystal, its environment has to be carefully considered. That the highest mobility so far has been achieved for suspended graphene proves this point (Zhang et al., 2005). As a consequence it is of particular importance to investigate the interface structure and to study how the substrate impacts the electronic properties of the graphene layer. For graphene grown on the so-called "Si side" of SiC, denoted as SiC(0001), the interface is quite well understood. The first carbon layer strongly interacts with the substrate and forms covalent bonds (Riedl and Starke, 2007; Emtsev et al., 2008a). The interface consists of a C-rich layer, possibly with defects (Qi et al., 2010), √ √ having 6 3 × 6 3-R30 1 periodicity with respect to the two-dimensional unit cell of the substrate surface (denoted 6 3 for short). It acts as an electronic buffer layer, between the SiC and overlying graphene layers, so that the a possible second graphitic plane on top of the buffer layer. exhibits graphene low-energy electronic properties (Mallet et al., 2007; Rutter et al., 2007; Brihuega et al., 2008; Varchon et al., 2007, 2008; Mattausch and Pankratov, 2007; Kim et al., 2008). For the C terminated surface (SiC(0001)), the situation is less well understood. Inconsistent results have been obtained in different studies: For ultrahigh vacuum- grown samples, most studies reveal coexisting native surface reconstructions of SiC (namely the SiC(3 × 3) (Hoster et al., 1997) and the SiC(2 × 2) C (Seubert et al., 2000b)), with graphene layers being weakly bound to the underlying SiC (Emtsev et al., 2008a; Hiebel et al., 2008; Starke and Riedl, 2009). The graphene-substrate interaction is much weaker than on the Si face. Graphene is almost freestanding on the SiC(3 × 3) reconstruction (Hiebel et al., 2009) and weakly interacting with the SiC(2 × 2) C reconstruction (Hiebel et al., 2008; Magaud et al., 2009). Monolayer samples on the C sied show a band structure that resembles that of freestanding graphene, but with an intrinsic doping in electrons (Emtsev et al., 2008a). On the other hand Hass, Feng et al. revealed an interface layer that was strongly bound to the SiC (i.e., similar to the Si-face situation) (Hass et al., 2007). The sample used in that latter study differed from those of the former in terms of its preparation (furnace grown vs vacuum prepared), so, Hass, Millán-Otoya, First, and Conrad suggested, it is possible that their interface structures could be different (Hass et al., 2008c). Srivastava, He, and Feenstra had resently demonstrated that the interface between graphene and SiC(0001) depends on the means of forming the graphene (Srivastava et al., 2012). For formation in vacuum, a 3 × 3 interface structure can be observed (Hass et al., 2008b; Emtsev et al., 2008a; Hiebel et al., 2008; Srivastava et al., 2012). When the graphene is formed in a Si-rich environment Srivastava √ √ et al. observe a carbon interface layer with 43 × 43- R ± 7 . 6 structure. This interface layer could be explained by a structure that is to some extent similar to √ √ the 6 3 × 6 3-R30 interface layer of the Si-face, but with the supercell for the carbon-face interface being rotated by only ± 7 . 6 , rather than 30 , relative to the substrat axes. With this interface layer lying directly on the C-face SiC substrate, it is found to form a buffer layer between the SiC and the graphene, similar to that modeled theoretically by Varchon et al. (Varchon et al., 2007). This thesis aims to give a comprehensive view over the various configurations either from literature or from well-known surface reconstructions and to obtain inform- ation on the interface between graphene and silicon carbide (SiC(0001)) using ab initio calculations. This can be used to verify the structure models and to make a clear statement about the stability of all configurations.

At the beginning the theoretical background for further understanding is conveyed. This is followed by a detailed examination of the diamond-graphite coexistence in order to understand the stability behavior in the carbon limit, which is crucial for the growth of graphene. It is shown that despite the distinct nature of the graph- ite and diamond bonds, DFT and the finite displacement method predict various properties of both phases in qualitative and quantitative agreement with the ex- perimental data. Subsequently it is demonstrated that ab initio calculations are a mighty tool in order to predict the stability of surface structures.

With the presented methodology, the well-known (2 × 2) C structure is verified and a variety of proposed models for the experimentally observable (3 × 3) reconstruction is proven to be wrong. However an adaption from the Si side of the substrate the (3 × 3) twist -model is proposed as the most likely stable structure. The results leave little scope that the wanted structure is different from the presented one.

Finally, it is shown that both, the recently experimentally observed 43 and the from the Si side adapted 6 3 graphene-like interface structures are unstable on the bare substrate.

Part II. Theory

Chapter 1. Basic Concepts: Density Functional Theory

The physical and chemical properties of any given system can be predicted by solving the corresponding Schrödinger equation. For the systems in this thesis the eigenvalue problem is reduced within the Born-Oppenheimer approximation (Born and Oppenheimer, 1927)1 to two "simpler" eigenvalue problems: the electronic and the nuclear one. The approach to calculate the energy E of the electronic Schrödinger equation

illustration not visible in this excerpt

for the Born Oppenheimer surface is based on Density Functional Theory (DFT), where He is the Hamilonian of the electronic system and ψ (ri; RI) the electronic one.

wave function. The key idea behind DFT is to characterize a given system via the observable electron density, not the wave functions. The advantage lies in the reduction of the dimensionality of the problem. Instead of dealing with a N -electron wave function with a dependence on 3 N spatial and N spin variables one has to deal with a density distribution depending on 3 coordinates and 2 spin components. The electronic energy of the system can then be expressed as a functional of the electron density E [ n (r)]. The minimization of this energy within a self consistent field (scf) calculation (figure 1.1) results in the ground state density and energy of the system for the "exact functional".2

1.1. Hohenberg-Kohn theorem

It is important that the electron density n (r) has to uniquely characterize the Hamilonian He to be able to determine all properties of the system of interest using the ground state density n 0. This was first proven by Walter Kohn and Pierre Hohenberg in 1964 (Hohenberg and Kohn, 1964) and is usually referred to as the first Hohenberg-Kohn Theorem. The second Hohenberg-Kohn theorem is based on a variational principle according to which the energy functional E [ n (r)], i.e. the total energy of the system as the functional of the system’s electron density n (r), gives the lowest energy only if the electron density is the true ground state density.

1.2. Kohn-Sham approach

The Hohenberg-Kohn theorem shows that it is possible to describe a system using the electron density and that the ground state density minimizes the total energy of the system. Still one requires a method to find n 0 or, in practice, a very good approximation to the true n 0. In 1965 Walter Kohn and Lu Jeu Sham pioneered for such a method by creating a decent mathematical basis (Kohn and Sham, 1965).3 The transformation of an interacting, inhomogeneous electron gas, moving in a potential v (r) into a system of noninteracting electrons, moving in an effective potential veff (r) is the fundamental idea. The new system can then be described by one-particle Schrödinger equations. Whereas the potential veff (r) depends on the electron density and as the electron density results from the solution of the Schrödinger equation, the whole system has to be solved using an iterative scf calculation.

The energy functional can be displayed by the following defining equation to a system of non-interacting electrons:

illustration not visible in this excerpt

with the electron density n (r) and the kinetic energy of a system of non interacting electrons TS [ n ]. EXC [ n ] contains the exchange and correlation energy as well as the difference between TS [ n ] and the actual kinetic energy of the system of interacting electrons. The exact form of EXC [ n ] is unknown, but several approximations exist, as shown in 1.3. Kohn and Sham have shown that this transformation to a system of non-interacting electrons in (1.2) can be performed, with respect to the first Kohn-Sham Theorem (Kohn and Sham, 1965). The ground state density of a non- interacting systems is given by:

illustration not visible in this excerpt

the one-particle Schrödinger equations, usually referred to as Kohn-Sham equations, with the one-electron wave functions as Kohn-Sham orbitals ψ j (r). For the effective potential is veff (r) = v + vH + vXC, with

illustration not visible in this excerpt

The density can then be calculated from the Kohn-Sham orbitals:

illustration not visible in this excerpt

however, veff thus depends on n (r), and n (r) depends on veff in Eq. 1.3. This is a chicken-egg problem. We need to find a density that generates the same potential that would generate this density. This is referred to as the self-consistency problem and here solved by an iterative approach, the so called scf cycle which is illustated in figure 1.1.

1.3. Exchange correlation functionals

DFT is an exact theory, in practical use the lack of knowledge of the shape of the exchange correlation (xc) functional represents a great interesting challenge. The local density approximation (LDA) is the simplest approach. In this case, one starts from a homogeneous electron gas and receives the approximation for the xc energy in the non-homogeneous case, by considering the system as a locally homogeneous electron gas at any point.

illustration not visible in this excerpt

illustration not visible in this excerpt

Figure 1.1.: Self consistent calculation using the Kohn-Sham approach

The xc energy per particle of the homogeneous electron gas ϵ hom XC canbedetermined by means of Monte Carlo simulations (Ceperley and Alder, 1980). This approx- imation allows for example the reliable determination of lattice constants; however cohesive energies are overestimated when using it (Haas et al., 2009). The more expedient approximations are the so-called generalized gradient approx- imations (GGA). By also including the gradient of the density in the determination of the xc energy one can obtain a potentially better approximation. This concept was first proposed by Langreth and Mehl in 1983 (Langreth and Mehl, 1983). The "simple" gradient expansion approximation (GEA) was already suggested by Kohn and Sham but does not work.4 The xc energy for generalized gradient approxima- tion (GGA) can be written as

illustration not visible in this excerpt

In this thesis we used the local density approximation based on Ceperley and Alder (Ceperley and Alder, 1980) as parametrized by Perdew and Wang 1992 (Perdew and Wang, 1992). For GGA we used the generalized gradient approximation of Perdew, Burke and Ernzerhof 1997 (Perdew et al., 1996) (PBE). Both functionals are implemented in the program package FHI-aims.

1.4. Van der Waals interactions in DFT

Within the DFT framework there are several methods to include long-range van der Waals interactions into approximate forms of Exc [ n ]. The scheme which is used in this thesis is considered closer is due to Tkatchenko and Scheffler (Tkatchenko and Scheffler, 2009).5 This is based on adding an empirical potential of the form C 6 R − 6 to the DFT energy with the interatomic distances R and the dispersion coefficients C 6. Hence, the correction is based on the optimized positions of the atomic nuclei and not directly on the electron density. The dispersion corrected total energy is then

illustration not visible in this excerpt

with the DFT energy EDFT obtained using a short-ranged xc functional and the empirical dispersion correction Edisp given by

illustration not visible in this excerpt

N denotes the number of atoms in the system, Cij 6 is the dispersion coefficient for interacting atoms i and j and Rij is the interatomic distance. To limit the correction to long ranged effects at large R a damping function given by

illustration not visible in this excerpt

at short distances (Tkatchenko et al., 2011). Certainly long-range electrodynamic response as well as non-additive many-body effects are not taken into account (Tkatchenko et al., 2012). Nevertheless the TS-vdW method can be used as ansatz for further development. By self-consistently including long-range screening effects in the effective atomic polarizabilities and using the coupled fluctuating dipole model (Tkatchenko et al., 2012) (CFDM) to include many-body effects the original method can be improved for more complex systems, e.g. large molecules, bulk solids and surfaces (Tkatchenko et al., 2012). A very few results in this thesis also show CFDM data for comparison. These results were calculated by Dr. Alberto Ambrosetti and applied to the framework constructed within the context of this thesis.

Chapter 2. Phonons

By itself, DFT just gives the ground-state, zero temperature density and total energy in the Born-Oppenheimer approximation for assumed classical nuclei. How- ever we can use the Born-Oppenheimer surface energy Etot [ r (j 1 , l 1) , ..., r (jn, lN)] as a starting point to bring back some of the temperature dependence of the full system, where r is the position of the j -th atom in the l -th unit cell, with n the number of atoms per cell and N the number of unit cells. The theory of crystal vibrations can be practically implemented by calculating the gradients of the Born-Oppenheimer energy surface hence the forces

illustration not visible in this excerpt

are obtained. With the acquired information the central quantity in the calculation of the phonon frequencies the force constant matrix

illustration not visible in this excerpt

of the system can be constructed, where α and β are the Cartesian indices. The equation for the force constant matrix Φ αβ (jl, j ′illustration not visible in this excerpt

where F β (j ′ l ′; Δ r α (jl)) are the forces on atoms with a finite displacement Δ r α (jl) and usually F β (j ′ l ′) 0. This method of calculating the force constant matrix by explicitly displacing atoms is called the finite displacement method. The periodic boundary conditions used in DFT calculations can cause problems for the finite displacement method. Displacing one atom in a small unit cell creates forces on all the atoms in the same unit cell, but also on the periodic images of these atoms. Therefore large supercells are needed to accurately calculate the force constant matrix which is computationally demanding. The dynamical matrix is defined by

illustration not visible in this excerpt

where m is the atomic mass and q is the wave vector and the equation of motion is written as

illustration not visible in this excerpt

With the complete force constant matrix, the dynamical matrix and hence the frequencies ω (q, s) can be obtained at any particular wave vector q by the diagon- alization of the dynamical matrix. Phonons are quantizations of the normal modes of atomic oscillations in a crystal. The energies and entropies of vibrations, which can be calculated by integrating the phonon density of states, are a substantial contribution to finite temperature thermodynamics of crystals as for instance the Helmholz free energy defined by

illustration not visible in this excerpt

the entropy

illustration not visible in this excerpt

or the heat capacity of the crystal

illustration not visible in this excerpt

with ω (q, s) the eigenvectors, s band index, kB the Boltzman constant, T temper- ature. By performing DFT calculations (FHI-aims) and using a post-processing phonon analyzer (Phonopy) one can calculate the normal modes in the harmonic approximation of a crystal and thus calculate the phonon density of states (DOS). The first step to implement this theory in DFT is to relax the crystal and to bring it to its equilibrium geometry. This should first be done for the primitive cell of the crystal structure of interest. However, accurate phonon calculations require that (relatively) large supercells be used for calculations to ensure that the force constants between distant atoms go to zero.

Chapter 3. Computational Aspects of DFT

3.1. FHI-aims

The self-consistency problem which need to be solved was introduced in 1.2. The equation we want to solve is the Kohn-Sham equation (1.2). All calculations in this work were performed with FHI-aims program package, an abbreviation for "Fritz Haber Institute ab initio molecular simulation". It is an all-electron, full- potential code using numerically tabulated atom-centered basis functions to expand the Kohn-Sham orbitals ψ j (r). This allows high accuracy for the calculated total energies by using for each chemical element predefined hierarchical basis sets cal- culated on atom-centered grids of points. The total energy is the most important, experimental observable, property. From this quantity, one can obtain various prop- erties, e.g. equilibrium atomic structures, band structures, density of states and phonon dispersion curves. For more information about the construction strategy of the basis sets the reader is advised to read (Blum et al., 2009), however a briefly description is given by way of example in A.1. FHI-aims can compute the electronic structure and total energies, including relaxations and molecular dynamics, of both finite and periodic systems (Blum et al., 2009). Moreover, for instance, it allows to calculate core-level shifts and to simulate scanning tunneling microscopy (STM) images. Nearly linear scaling of the grid-based operations is guaranteed by the strictly localized orbitals, because different spatial regions are separated from one another. For instance, a system of polyalanine, i.e. a biopolymer with up to several hundred atoms, shows a nearly O (N) scaling with the system size N (Blum et al., 2009). Hence, solving the single-particle eigenvalue problem becomes the limiting part in the scalability for larger systems (Havu et al., 2009). This solver formally scales as O (N 3 ) with the system size.

3.2. Phonopy

This program can compute thermal properties from DFT ground state, zero temper- ature density and energy like described in chapter 2. Phonopy is an implementation of a post-process phonon analyzer, which calculates crystal phonon properties from input information calculated by external first-principles codes (Togo et al., 2008). It uses the Parlinski-Li-Kawazoe method (Parlinski et al., 1997) which is based on a supercell approach with the finite displacement method annotated in chapter 2. The calculation and symmetrization of force constants are executed by using singular-value decomposition (pseudo-inverse). The key of this method is the matrix formulations of equations, which leads to a coherent and flexible implementation (Parlinski et al., 1997).

Simple phonon calculations can be performed only at high-symmetry points in the Brillouin zone. With phonopy we can also perform phonon calculations at any point in the Brillouin zone using the supercell method and phonon interpolation. The size of the supercell is specified by the periodicity (the wave vector q) of the phonon displacements which are considered. In general phonopy calculates all the phonons corresponding to all the q -points belonging to a supercell mesh in reciprocal space. An interpolation scheme at general q -points with non-analytical term correction is implemented according to (Wang et al., 2010) (Pick et al., 1970; Giannozzi et al., 1991; Gonze and Lee, 1997; Wang et al., 2010).

Chapter 4. Ab Initio Atomistic Thermodynamics

The concept we will introduce in the following is a powerful tool to get the informa- tion of a surface. It allows to calculate the phase diagram of a surface system, i.e. a diagram which shows the surface free energy as a function of the chemical potential between the limits of the respective system. Figure 4.1 shows such a phase diagram as an example.

illustration not visible in this excerpt

Figure 4.1.: Surface free energy as a function of the chemical po- tential for the 3C-SiC(111) sur- face including the clean surface (blue), the (2 × 2) C reconstruc- tion (magenta) and the limit for the chemical potential for Si and C

The discussion of the shown phase diagram takes place in part 4 of this thesis. It establishes a connection between microscopic and mesoscopic regimes. In order to understand it, the framework on which this diagram is based as well as the indi- vidual components are explained below. As already explained in chapter 2 DFT is referred to as zero-temperature, zero-pressure-method. Ab initio calculations do not directly contain information about the thermodynamic behaviour of sys- tems under real conditions. The results of electronic structure calculations only correspond to T = 0 K and p = 0 atm for the electronic subsystem in an external potential of fixed, assumed classical nuclei. For this purpose ab initio thermody- namics combines DFT with classical thermodynamics (Reuter et al., 2004; Kaxiras et al., 1987; Scheffler and Dabrowski, 1988; Qian et al., 1988; Scheffler, 1988; Reuter and Scheffler, 2004, 2001). A brief overview how this is done within a "constrained" equilibrium is shown in the following.

4.1. Constrained Thermodynamic Equilibrium

If one investigates surfaces one most certainly wants to find out which surface struc- ture will be the thermodynamically stable surface structure under given conditions. In this case the key quantity is the Gibbs free energy. This thermodynamic po- tential allows us to identify the atomic structure with the lowest energy under real temperature and pressure conditions of a surrounding gas phase. The total energy of the system contributes to the internal energy U. The Gibbs free energy is given by

illustration not visible in this excerpt

The leading term Etot is the total energy of the system, which is directly obtained from DFT calculations, Fvib accounts for the vibrational free energy which provides a non-negligible contribution to the Gibbs free energy for the solid phase, T is the temperature, S the entropy, p the pressure, and V the volume. In general a reservoir of atoms to be in equilibrium with the surface is needes. A hypothetical or real gas phase can be such a reservoir. The system is characterized by the chemical potentials μ (T, p). In our case two independent gas phase reservoirs, Si an C, are assumed to be individually in equilibrium with the SiC surface or bulk. The thermodynamic stability of different surface reconstructions in this "constrained" equilibrium, with a gas reservoir containing molecular silicon and carbon as a function of tempertaure and particle pressure, can be investigated.

4.2. Surface Free Energy

The surface free energy represents the increase in internal energy if a solid is cleaved into two parts so that two surfaces appear, i.e. the energetic cost of this endothermic process (Zangwill, 2001). The internal energy of the bulk system Ebulk = Fvib + TS − pV + ∑ i N i μ i, where N i is the number of particles in the system, and μ i the chemical potential of the atomic species i. It is increased by a value proportional to the created surface area A. Therefore the internal energy of a cleaved solid material is given by

illustration not visible in this excerpt

The factor γ is defined as the surface excess free energy per unit area of the considered plane in a crystal (Zangwill, 2001). It is referred to as the surface free energy. Since the bulk is not exfoliating it is a positive quantity and by minimizing γ one will get the most stable surface.

The major issue in ab initio thermodynamics is that the system is assumed to be in thermodynamic equilibrium with its enviroment. Within this framework this system, as a multi-component system, is in equilibrium with different reservoirs which are again in equilibrium with the surface but do not necessarily have to be in equilibrium with each other. Using this concept of a "constrained" equilibrium to the surface being in equilibrium with a gas phase reservoir above this surface, the surface free energy expressed in terms of their chemical potentials μ i for system containing i components can be written as

illustration not visible in this excerpt

where Gsurf = Gtot − Gbulk, i.e. the Gibbs free energy of the cleaved surface without the Gibbs free energy contribution of the solid, Gbulk. As a first approximation differences between vibrational free energies Δ Fvib will be neglected in the following. This can be justified because Δ Fvib accounts to an assumed reservoir of C and Si which provides us with particles to put on or take away from our surface and therefore the contribution might be small. In the following we will consider the SiC system we are investigating. Another important quantity is the change of the Gibbs free energy the Gibbs free energy of formation, which is defined as

illustration not visible in this excerpt

For Δ GSiC (T, p) = 0 Si, C and SiC are in equilibrium, for positive values the system builds SiC and for negative ones Si and C precipitate.

4.3. SiC Limits

At first the limits of the chemical potential μ Si and μ C have to be defined. Under equilibrium conditions the permissible range of the chemical potentials is determ- ined by the elemental crystal phases of carbon and silicon. The maximum of the chemical potential μ C, max (μ C) =12 g C (T,p)canbeobtainedfromtheGibbsfree energy per atom (4.1). With the equilibrium condition of the Gibbs free energy of formation μ c + μ Si = gbulk SiC (T,p),onecanfind min (μ Si)= g S iC (T,p) 2 g C (T,p), which leads to the lower limit for the chemical potential μ Si. The maximum of the

chemical potential max (μ Si) =1 2 gS i (T, p) is given for pure silicon. The definition for the Gibbs free energy of formation in the case of Δ GSiC (T, p) < 0 gives us than the upper limit μ Si < gbulk SiC − 2 g C (T,p) Δ G SiC (T,p).

Finally one obtains for the limits of the chemical potential μ Si:

illustration not visible in this excerpt

Since all DFT ground state total energies energies, are simply zero temperature values in the Born-Oppenheimer approximation, G (T, p) = G (T = 0 K, p = 0 atm) = Etot. Therefore limits are written in terms of total energies

illustration not visible in this excerpt

For the purpose of determinating the formation energy, total energies had to be calculated for carbon and silicon. However the differences between the two carbon phases are very small. Covalent bonds like in diamond are well characterised by DFT, while inter-planar bonds like in graphite have to be treated with care. There- fore, the following part focuses on the thermodynamically legal carbon boundaries of the chemical potential in the phase diagram, for which an investigation of the pure bulk solids is needed.

Part III. Diamond Graphite Coexistence

In order to be able to make a statement about the stability of our investigated

surface reconstructions it is important to firstly investigate the carbon bulk phases. The phase diagram of carbon was controversial for many years; the higher regions are still poorly understood (Whittaker, 1977; WHITTAKER, 1978; Bundy, 1989; Steinbeck et al., 1985; Gustafson, 1985; Mitchell et al., 1986). The diagram is shown in Figure 4.2, as based on (Bundy, 1989; Steinbeck et al., 1985). There

illustration not visible in this excerpt

Figure 4.2.: Phase diagram of carbon (Zazula, 1997).

are at least two stable solid phases of carbon: from the structural point of view, graphite represents a crystalline hexagonal form, and diamond a tetrahedral form. Both forms can exist in the same quite wide range of thermodynamic conditions; however graphite can be transformed to diamond only at very high pressures above

10 GP a.1 Since there are two phases for the carbon-rich region which play a role in the epitaxial graphene growth, we want to examine the diamond-graphite coexistence at this point.

Chapter 5. Carbon

In solid state carbon atoms can arrange in several allotropes. Graphite and diamond are the best-known three-dimensional crystals but also for fewer dimensions allo- tropes exist, e.g. fullerenes, carbon nanotubes and graphene. The physical properties of carbon changes drastically with the allotrope that we consider, e.g. diamond is optically transparent, the hardest naturally-occurring material and has a low elect- ical conductivity, whereas graphite is optically opaque, soft enough to draw a line on a paper and an excellent conductor. 1 In this chapter the structure and electronic

illustration not visible in this excerpt

Figure 5.1.: Carbon allotropes, from left to right: diamond, graphite, fullerene, carbon nanotube and graphene. Picture taken from (Scarselli et al., 2012). properties of graphene, graphite and diamond are presented.

5.0.1. Graphene

In 1947 P.R. Wallace studied the band structure of the two-dimensional graphene crystal as a starting point to the investigation of graphite (Wallace, 1947). Back then it was only a theoretical object because due to thermodynamical instability two-dimensional crystals were not expected to exist. A detailed overview on the sub- ject is given in Neto et al. (2009). Graphene is an atomically flat two-dimensional crystal composed of sp 2 -hybridized carbon atoms. In sp 2 -hybridization, the 2 s orbital is mixed with two of the three available 2 p orbitals while the 2 pz orbital remains unchanged. In the graphene crystal three coplanar σ -bonds per carbon atom result from the overlapping of the sp 2 orbitals. The σ -band is responsible for graphene’s robust honeycomb structure. The unaffected p -orbital, which is per- pendicular to the honeycomb lattice, binds covalently between neighboring atoms to form a π -band. From the latter follow the electronic properties of graphene since only the π (π *) states are found in the vicinity of the Fermi level (Latil and Henrard, 2006).

5.0.2. Graphite

As shown in figure 5.1 graphite has a layered, planar structure. It can be seen as stacked layers of graphene. There are two known forms of graphite, alpha (hexagonal) and beta (rhombohedral), which differ in the graphene layers stacking but have very similar physical properties (Cousins, 2003). For the sake of simpli- city only alpha graphite of space group P 63 /mmc is investigated in this work. The hexagonal lattice with same atoms corresponds to a close-packing of spheres where the layers are packed in an ABA form, i.e. every other layer is the same. The unit cell is defined by the lattice vectors

illustration not visible in this excerpt

which are given by the two lattice parameters a and c. The cell has a four atomic basis with

atoms at illustration not visible in this excerpt

per atom by Vatom = 8 a 2 c.

5.0.3. Diamond

Diamond has a cubic crystal structure of tetrahedrally bonded carbon atoms in a covalent network lattice (sp 3 ). The lattice consists of two interpenetrating face- centered cubic (fcc) lattices. Each carbon atom is covalently bonded to four equi- valent neighboring atoms. It is a zincblende structure, belonging to space group

illustration not visible in this excerpt

The unit cell is defined by the lattice vectors(0 |a 2 | 2 ),(a 2 | 0 | 2 ) and(2 | 2 | 0

which are given by one lattice parameters a. The cell has a two atomic basis with

illustration not visible in this excerpt

atoms at (0 | 0 | 0) and 4 | 4 | 4 ) (Wyckoff, 1960; Wieferink et al., 2011). The unit cell volume is then given by Vuc =14 a 3 andthevolumeperatomby V a tom =8 a 3.

Chapter 6. Determination of the PES of Diamond and Graphite

6.1. Structure optimization

If one focuses on a new unknown system one needs to become familiar with its computational behaviour. Therefore the first step always comprises performing convergence tests. They are not only essential for an accurate description of the system but also necessary to accomplish further results in a reasonable amount of time and money. A brief glimpse is given in the appendix A. At this point only the Born-Oppenheimer surface is considered. The physically more relevant results with zero-point energy (ZPE) will be addressed further below. There are several ways to optimize structures with two lattice constants by minimizing the corresponding Born-Oppenheimer surface energy E (a, c). The following method has been employed in this thesis:

6.1.1. Minimal Energy for Fixed Volume Murnaghan Fit:

At first a fixed volume close to the experimental volume of the unit cell is picked. After that a number of lattice parameters a close to the experimental value are chosen. Corresponding to the fixed volume pairs (a, c) are obtained. For this pairs (a, c) the energy of the respective structure is calculated. If a polynomial fit is applied to the plotted energy over a the lattice parameter with the lowest energy can be determined and the energy of the corresponding structure calculated. This procedure is iterated for a couple of volumes in the vicinity of the experimentally obtained volume. Subsequently with the lowest energies for the fixed volumes one is able to plot the binding curve of this system, i.e. the energy over the volume. Continuum mechanics provides a few approaches for a plausible dependence of energy on volume, e.g. a harmonic solid, Murnaghan or Birch-Murnaghan equation of state (6.1).

For the latter one the energy per atom (E = Ecoh) is expressed as a function of the atomic volume (V = Vatom)

illustration not visible in this excerpt

The fitting parameters V 0 and E 0 are the equilibrium atomic volume and atomic energy, and B 0 is the bulk modulus and B ′ 0 itsderivativewithrespecttopressure. By fitting the results to the Birch-Murnaghan equation of state the volume with the lowest energy, as well as the bulk modulus B 0 of the structure is determined (Murnaghan, 1944; Birch, 1947; Scherz, 2005). If now the introduced procedure for a fixed volume to the volume with the lowest energy is applied we obtain the equilibrium lattice parameter (a, c) of the structure. For systems with only one lattice parameter a this method is reduced, because each fixed volume is corres- ponding to one lattice parameter a, i.e. only one calculation per iteration instead of several ones for a pair (a, c). Besides graphite this method was used to find the equilibrium lattice constant and total energy for diamond, bulk silicon and SiC.

Since the investigated structures are highly symmetric and no structural changes are expected, no atomic relaxation within the unit cell is needed. In the case of diamond the total energy was calculated for five values of the lattice constant a, a 0, a 0 ± 1% and a 0 ± 2% wheras a 0 is the lattice constant in equilib- rium. Using the Birch-Murnaghan equation of state the potential energy surface could then be fitted.

In contrast to diamond, graphite has a hexagonal crystal structure. The unit cell is defined by two lattice parameters a and c. For eight fixed volumes (7.9 Å3, 8.1 Å3, 8.3 Å3, 8.5 Å3, 8.7 Å3, 8.9 Å3, 9.1 Å3 and 9.3 Å3 ) the lattice parameter a was set to five values close to the equilibrium parameter a 0 (2.43 Å, 2.44 Å, 2.45 Å, 2.46 Å and 2.47 Å). The volume of the unit cell is given by V =1

illustration not visible in this excerpt

With that in mind one is able to calculate the lattice parameter c. Using these five isochoric geometries, a parabolic fit over the lattice parameter a then yields the minimum energy. The minimal energies plotted over the volume could then be fitted likewise to the diamond data with the Birch-Murnaghan equation (6.1). The Born-Oppenheimer PES for diamond and graphite calculated with the different xc functionals are presented in figure 6.1, 6.2 and 6.3.

For LDA the diamond phase is more stable than the graphite phase (see fig. 6.1), this behaviour is well know since 1984 (Yin and Cohen, 1984). The difference of the equilibrium energy Δ ELDA d − g ≈ 6 meV is quite small. In the case of PBE the energy surface of graphite is very flat (see fig. 6.2). To correctly describe the PES much more data points are need compared to LDA. The data set for eight fixed volumes is not enough. Here the difference of the equilibrium energy Δ EPBE d − g ≈ 1 eV is huge in comarison to LDA. Figure 6.3 shows the PBE result with vdW in the TS scheme introduced in section 1.4. The PES of graphite ist not flat anymore. The gradient of the graphite PES looks much more like the LDA gradient. However the beha- viour in terms of stability for the PBE+vdW TS scheme is correct. The difference of the equilibrium energy Δ EPBE + vdW (TS) d − g ≈ 60 meV is between the value for LDA

illustration not visible in this excerpt

Figure 6.1.: PES diamond and graphite calculated with LDA (tight species defaults, tier2, 24 × 24 × 24 k -grid in the case of diamond and 24 × 24 × 12 in the case of graphite). and pure PBE.

With the framework which was created in this work it was possible to calculate the PES for the PBE+vdW CFDM1 with a few specific calculations. The post processed results of the calculations performed by Dr. Alberto Ambrosetti are shown in figure 6.4 compared to the TS scheme. The comparing both methodes shows that the difference for the diamond PES Δ ETS − CFDM d ≈ 4 meV is smaller than the difference for the graphite PES Δ ETS − CFDM g ≈ 40 meV. The difference of the equilibrium energy Δ EPBE + vdW (CFDM) d − g ≈ 20 meV is reduced compared to the TS scheme. (i.e. including many-body effects to the original method) 1 Surface reconstructions are given in Woods’ notation (K. et al., 2003).

[...]


1 The electrons respond much faster to an external perturbation than a nucleus, due to the mass difference. Therefore the electrons can follow the movement of the nuclei instantaneously, thus they can be treated as moving in a constant field generated by fixed nuclei. This enables a separation of the Hamiltonian into a nuclear and an electronic part as well as a splitting of the total wave function into a nuclear and an electronic one.

2 In practice approximations are made, so the result is a good approximation to the true ground state density and energy of the system.

3 It was a formal proof that a generalization of Slater’s X-alpha

4 In comparison with GGA GEA produces an exchange correlation hole that is not normalized to 1.

5 The method in principle has a much longer history (Ahlrichs et al., 1977). illustration not visible in this excerpt + R 0 j isthesumoftheeffectiveatomicvanderWaals (vdW) radii of atoms i and j, the parameter d adjusts the steepness of fdamp and the scaling factor sR adapts fdamp to a particular DFT functional. Whereas early approach makes use of C 6 coefficients which are based on experimental values (Wu and Yang, 2002), Tkatchenko-Scheffler (TS) proposed a scheme which determines the required parameters from the DFT ground-state electron density (Tkatchenko and Scheffler, 2009), referred to as TS approach or scheme. The main idea of the TS method is to use the electron density for computing the relative rather than the absolute atomic dipolar polarizability and C 6 coefficients of an atom in a group of atoms based on highly accurate free atom values. The TS approach has several advantages; the method itself is easy to use with only the empirical parameter d and sR, the vdW energy for atoms in molecules and solids is computed using ground-state electron density with almost no additional cost beyond standard DFT calculations and hybridization effects on the polarizability are accurately treated

1 The intermediate metastable solid phases like carbyne and chaoite were discussed for several years (Whittaker, 1977). Evidence for a solid-to-solid phase transformation from diamond to a metallic state (Solid III), possibly cubic, is described at extreme high pressure (>100 GP a) (Gust, 1980).

Author

Share

Previous

Title: First-Principles Prediction of Phase Stability of Electronic Materials