Document [original]

ANNEALING MECHANISMS OF POINT

DEFECTS IN SILICON CARBIDE

A Theoretical Investigation

Dissertation

zur Erlangung des akademischen Grades

Doktor der Naturwissenschaften (Dr. rer. nat. )

anerkannt von der Fakult¨at f¨ur Naturwissenschaften

der Universit¨at Paderborn

Dipl. Phys. Eva Rauls

Paderborn, 2003

Der Fakult¨at f¨ur Naturwissenschaften der Universit¨at Paderborn als Dissertation vorgelegt.

1. Gutachter: Prof. Th. Frauenheim

2. Gutachter: Prof. H. Overhof

Tag der Einreichung: 20. 05. 2003

Tag der m¨undlichen Pr¨ufung: 01. 07. 2003

Eva Rauls, Paderborn 2003

Contents

Introduction 1

Outline of this Work 2

1 Some Facts about Silicon Carbide 3

1.1 Historical Facts and Applications . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2 Doping of SiC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.3 Crystal Structure: the Polytypes of SiC . . . . . . . . . . . . . . . . . . . . 8

2 Theoretical Description of Defect Dynamics 13

2.1 Some General Aspects of Point Defects . . . . . . . . . . . . . . . . . . . . . 13

2.1.1 Point defects at thermodynamic equilibrium . . . . . . . . . . . . . . 13

2.1.2 The situation after an implantation process . . . . . . . . . . . . . . 14

2.1.3 Migration of defects . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.1.4 Jump probabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.1.5 Charge state effects on the migration of defects . . . . . . . . . . . . 16

2.2 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3 SCC–DFTB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.4 Modeling of Defects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.5 The Diffusion Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.5.1 The constrained relaxation technique . . . . . . . . . . . . . . . . . . 25

2.5.2 The activation relaxation technique (ART) . . . . . . . . . . . . . . 26

2.6 Lattice Vibrations and Free Energies . . . . . . . . . . . . . . . . . . . . . . 26

2.7 Applications and Test Calculations . . . . . . . . . . . . . . . . . . . . . . . 28

2.7.1 The vibrational spectrum of 3C-SiC bulk . . . . . . . . . . . . . . . 28

2.7.2 Heat capacity of diamond, silicon, and SiC . . . . . . . . . . . . . . 30

2.7.3 Absolute entropy in different methods . . . . . . . . . . . . . . . . . 31

2.7.4 Formation entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.7.5 The vacancy in diamond and silicon . . . . . . . . . . . . . . . . . . 37

ii CONTENTS

3 Vacancies and Interstitials 41

3.1 Vacancies ..................................... 41

3.1.1 Formation energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.1.2 Migration of vacancies . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.1.3 Vacancy – antisite pair formation . . . . . . . . . . . . . . . . . . . . 48

3.2 Interstitials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.3 Interstitial Recombination with Vacancies . . . . . . . . . . . . . . . . . . . 52

4 Aggregation of Antisites 55

4.1 The Antisite Pair CSi SiC............................ 55

4.1.1 Formation Energy of CSi SiC...................... 55

4.1.2 Properties of the Antisite Pair . . . . . . . . . . . . . . . . . . . . . 57

4.1.3 Creation of the CSi SiCPair in the perfect SiC-lattice . . . . . . . . . 59

4.1.4 Migration of Antisites in the ideal lattice . . . . . . . . . . . . . . . 60

4.1.5 Pair creation by vacancy migration . . . . . . . . . . . . . . . . . . . 60

4.2 Larger Aggregates of Antisites . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.2.1 Stability of various arrangements . . . . . . . . . . . . . . . . . . . . 62

4.2.2 A Vacancy’s Spiral Walk . . . . . . . . . . . . . . . . . . . . . . . . . 65

4.3 Contributions of the Vibrational Entropy . . . . . . . . . . . . . . . . . . . 69

4.4 Mobility of Isolated Antisites . . . . . . . . . . . . . . . . . . . . . . . . . . 72

4.4.1 The vacancy assisted mechanism . . . . . . . . . . . . . . . . . . . . 73

4.4.2 The SiC(CSi )2complex . . . . . . . . . . . . . . . . . . . . . . . . . 75

4.5 Influence of other Defects on Migration Processes . . . . . . . . . . . . . . . 77

4.6 Conclusions for the Formation of Antisite Aggregates . . . . . . . . . . . . 79

5 Nitrogen–related Defects 83

5.1 Nitrogen–related Pair Defects . . . . . . . . . . . . . . . . . . . . . . . . . . 84

5.1.1 Pair formation by aggregation of VSi and NC............. 85

5.1.2 Mobilizing NC— Creation of N-interstitials . . . . . . . . . . . . . . 86

5.1.3 Recombination of (NC)Cwith divacancies VCVSi . . . . . . . . . . 87

5.1.4 (NC)Cmeeting isolated vacancies . . . . . . . . . . . . . . . . . . . 89

5.1.5 The CSiNCpair as an alternative . . . . . . . . . . . . . . . . . . . 90

5.2 Dissociation or Aggregation . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

5.2.1 Dissociation of NCVSi pairs . . . . . . . . . . . . . . . . . . . . . . 91

5.2.2 VSi (NC)n-complexes . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

5.2.3 Formation of CSi(NC)ncomplexes . . . . . . . . . . . . . . . . . . . 94

5.3 Correlation of Activation Energies and Temperatures . . . . . . . . . . . . 95

5.3.1 Definition of an assignment . . . . . . . . . . . . . . . . . . . . . . . 96

5.3.2 Correlation of the calculated values with experimental findings . . . 97

CONTENTS iii

5.3.3 Entropy effects on nitrogen migration . . . . . . . . . . . . . . . . . 98

5.4 Implantation with Phosphorus . . . . . . . . . . . . . . . . . . . . . . . . . 99

6 Summary and Outlook 103

A Formation Energies 105

B Calculation of the Gibbs Free Enthalpy 107

C Basic Properties of Strain Fields 109

D Summary: Activation Energies 111

Acknowledgment 119

iv CONTENTS

List of Figures

1.1 The ones who found SiC: J. J. Berzelius and H. Moissan . . . . . . . . . . 3

1.2 SiC as Jewels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.3 Performance of SiC Schottky diodes . . . . . . . . . . . . . . . . . . . . . . 5

1.4 Polytypes of SiC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.1 Diffusion mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.2 Supercell and cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.3 3C-SiC: Vibrational density of states . . . . . . . . . . . . . . . . . . . . . . 28

2.4 4H-SiC: Vibrational density of states . . . . . . . . . . . . . . . . . . . . . . 29

2.5 Heat capacity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.6 Specific heat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.7 Absolute Entropy of SiC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.8 Embedded cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

2.9 Formation entropy of the vacancy in silicon and diamond . . . . . . . . . . 34

2.10 Sphere inscribed in a supercell . . . . . . . . . . . . . . . . . . . . . . . . . 36

2.11 Corrections to the formation entropy (Si, C). . . . . . . . . . . . . . . . . . 38

3.1 Transformation of the silicon vacancy to the CSi VCpair: geometry . . . . . 48

3.2 Transformation of the silicon vacancy to the CSi VCpair: energy curve . . . 49

3.3 The P6/P7 spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.4 Geometry of the (CC)Csplit-interstitial . . . . . . . . . . . . . . . . . . . . 51

3.5 Migration of (CC)Csplit-interstitials . . . . . . . . . . . . . . . . . . . . . . 51

3.6 Energy during molecular dynamics simulation at 1700 K: (CC)Cin 3C-SiC . 52

3.7 Molecular dynamics simulation: (CC)Cin 3C-SiC . . . . . . . . . . . . . . . 53

4.1 Geometry of CSi SiC............................... 57

4.2 Formation entropy of CSi SiC.......................... 57

4.3 Vibrational spectrum of CSi SiC........................ 58

4.4 The DI-spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.5 The Pandey mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.6 Energy change during the Pandey process . . . . . . . . . . . . . . . . . . . 60

vi LIST OF FIGURES

4.7 Antisite pair formation next to VC(left) and VC-CSi (right). . . . . . . . . 61

4.8 Energy diagram for antisite pair formation by vacancy migration . . . . . . 61

4.9 Possible orientations of two antisite pairs . . . . . . . . . . . . . . . . . . . 62

4.10 Number of ”wrong” bonds in the aggregate of nantisite pairs. . . . . . . . 63

4.11 Geometries of the most stable two-dimensional aggregates of antisites. . . . 63

4.12 Energies of 1D and 2D aggregates of antisites . . . . . . . . . . . . . . . . . 64

4.13 Geometries of the infinitely extended antisite aggregates . . . . . . . . . . . 65

4.14 Creation of antisite platelets by vacancy migration . . . . . . . . . . . . . . 66

4.15 Energy during the spiral walk, starting with VC............... 67

4.16 Energy during the spiral walk, starting with VCVSi ............. 68

4.17 Entropy-change at the saddle point during VCmigration . . . . . . . . . . . 69

4.18 Effect of entropy on VCmigration . . . . . . . . . . . . . . . . . . . . . . . . 70

4.19 Effect of entropy on the process VSi −→ VCCSi ................ 71

4.20 Influence of the entropy on the creation of antisite pairs . . . . . . . . . . . 72

4.21 Mechanism of vacancy assisted CSi movement . . . . . . . . . . . . . . . . . 73

4.22 Energy during the vacancy assisted CSi movement . . . . . . . . . . . . . . 74

4.23 Creation of the SiC(CSi )2and the SiC(CSi )3complex . . . . . . . . . . . . 75

4.24 Creation of the SiC(CSi )4............................ 76

4.25 The influence of other defects on VSi -migration next to CSi ......... 78

4.26 Potential wall and well of the inverted bilayer . . . . . . . . . . . . . . . . . 80

5.1 Recovery of the free carrier concentration . . . . . . . . . . . . . . . . . . . 83

5.2 Creation of a VSiNC-pair ............................ 85

5.3 Creation of (NC)Csplit-interstitials . . . . . . . . . . . . . . . . . . . . . . . 86

5.4 Rotation of the (NC)Csplit-interstitial. . . . . . . . . . . . . . . . . . . . . . 87

5.5 Recombination of (NC)Cwith divacancies . . . . . . . . . . . . . . . . . . . 88

5.6 Recombination of (NC)Cwith a divacancy . . . . . . . . . . . . . . . . . . . 88

5.7 Kick-out mechanism for (NC)Cand VSi .................... 89

5.8 Creation of CSiNCpairs ............................. 91

5.9 The VSi(NC)4complex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

5.10 The CSi(NC)ncomplex. ............................. 94

5.11 Creation of the CSi (NC)2complex . . . . . . . . . . . . . . . . . . . . . . . 95

5.12 Correlation of activation energies and temperatures . . . . . . . . . . . . . . 96

5.13 Influence of the vibrational entropy on the creation of (NC)C. . . . . . . . 99

5.14 Migration of phosphorus split-interstitials . . . . . . . . . . . . . . . . . . . 100

List of Tables

1.1 Properties of the three most common SiC-polytypes . . . . . . . . . . . . . 10

2.1 Phonon modes in 3C-SiC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.2 Data used for the calculation of the entropy of formation (Si, C). . . . . . . 37

3.1 Formation energies of vacancies in 3C-SiC . . . . . . . . . . . . . . . . . . . 42

3.2 Formation energies of vacancies in 4H-SiC . . . . . . . . . . . . . . . . . . . 42

3.3 Corrections to LDA formation energies of vacancies . . . . . . . . . . . . . . 45

3.4 Activation energies for sublattice migration of VCand VSi . . . . . . . . . . 46

3.5 Corrections to LDA activation energies of vacancies . . . . . . . . . . . . . . 47

4.1 Formation- and binding energies of antisites and antisite pairs in 4H-SiC. . 56

4.2 Formation- and binding energies of antisites and antisite pairs in 3C-SiC. . 56

5.1 Hyperfine interactions of the VSi NCpair, calculated within LMTO-ASA . . 85

5.2 Transition levels of VSi(NC)n.......................... 93

5.3 Transition levels of CSiNCand CSi(NC)2.................... 95

D.1 Activation energies (intrinsic defects) . . . . . . . . . . . . . . . . . . . . . . 111

D.2 Activation energies (dopants) . . . . . . . . . . . . . . . . . . . . . . . . . . 112

vii

Introduction

Silicon carbide (SiC) is a very promising material in semiconductor technology. Its hard-

ness, its high temperature stability and many other properties which will be described in

Chapter 1 give it a wider range of application than silicon. As a comparably new material,

there are, however, still many aspects of the material that are not fully understood and

which hinder its technological breakthrough in semiconductor industry. Understanding

and controlling the electronic properties of SiC is, therefore, the aim of current research in

this field.

Most electronic properties of semiconductors are determined by the kind and amount of

defects in the material. The surfaces of the crystal, interfaces to other materials or dislo-

cations in the bulk material influence the electronic properties of the material as well as

point defects of either intrinsic or extrinsic nature. The periodicity is broken by any kind

of defect, the lattice is locally distorted, electronic states may be changed or newly induced

as a consequence.

The origin of defects is manifold: the termination of a surface e. g. depends on the envi-

ronmental conditions during growth or e. g. oxidation, on the concentration of each species

and on the temperature. Dislocations can for example result from a lattice mismatch at

the interface with another material1.

Point defects are on the one hand unavoidably built in already during growth, on the other

hand, they are intentionally brought into the material in a special procedure in order to

tailor the electronic properties of the material.

For the fabrication of electronic devices, as e. g. field effect transistors, p-type and n-type

material is needed. This is created by doping with elements of the neighboring groups in

the periodical system. For SiC as a compound semiconductor of two group IV elements,

p-type material can be created by doping with group III elements, as boron or aluminum,

whereas n-type doping is achieved by implantation with group V elements, like nitrogen

or phosphorus.

But also intrinsic defects like vacancies or antisites affect the electronic properties of the

material, either in form of isolated defects or after complex formation with other intrinsic

defects or dopant atoms.

To understand how the electronic properties depend on the outer conditions and how tai-

loring of the electronic properties can be done most efficiently, defect dynamics have to be

investigated.

Experimental investigations alone can only supply incomplete information, and it is a

promising possibility to complement them with atomistic simulations. In many cases, such

simulations are, however, very costly, so that there were only few simulations in the past,

which could provide information that was useful to interpret experimental data. Due to lim-

ited resources very approximative simulation schemes had to be used, and, consequently, a

satisfying accuracy could often not be reached. More elaborate simulation methods, which

1This makes it, e. g., difficult to grow SiC on silicon substrate, where the ratio of 4:5 in the (110) planes

causes dislocations and thereby unavoidably electrically active interface states.

2INTRODUCTION AND OUTLINE

can provide the desired accuracy, are first now with increasing computing capacities getting

capable of supplying the results needed for the explanation of experimental observations.

This theoretical work deals with the behavior of selected point defects in SiC that can be

expected to be created during implantation and annealing processes commonly used for

the doping of SiC. Migration and aggregation of defects are covered with a special focus

on the understanding of the experimentally observed annealing phenomena. As a method

for atomistic simulations, a density-functional based tight-binding scheme has been used,

which is considerably more accurate than empirical tight-binding methods, but, especially

concerning the description of the electronic structure of defects, less accurate than first

principle methods. Its efficiency offers the possibility of a systematic investigation of the

defect dynamics, which in spite of today’s computing facilities is much too time consuming

for first principle methods. For more detailed investigations of selected defects, i. e. the

calculation of occupation levels or hyperfine interactions, various first principle methods

have been used, additionally.

Outline of this Work

In the first chapter, the investigated material, silicon carbide, is characterized. In a brief

historical sketch, the fabrication of SiC and with the improvement of material quality the

growing range of application is outlined. Silicon carbide’s characteristic property of crys-

tallizing in various crystal structures is discussed as well as current problems with doping

SiC in order to use it as a device material.

Subject of the second chapter are the possibilities to treat such problems theoretically.

Some general aspects of point defects, useful for the discussion in later chapters, are shown

up in the first section. In the following sections, the computational method, its theoret-

ical basis and algorithms for the calculation of diffusion paths or vibrational frequencies

are briefly described. The last section of this chapter is dedicated to the calculation of

entropies, needed for the description of high temperature processes. Various test calcula-

tions, also for silicon and diamond, are presented, here.

The third chapter summarizes the results of our calculations for formation and migration

energies of vacancies and interstitials in SiC. On the one hand, some reference calculations

allow a comparison of the results obtained within our method with literature data, and

thus an estimate of the expectable accuracy. On the other hand, the basic processes needed

in the following two chapters, i. e. the sublattice migration of vacancies and carbon split-

interstitials, are discussed. Antisites and the creation of pairs and larger aggregates of them

are investigated in the fourth chapter. Several mechanisms for aggregation of antisites to

one-, two- or three-dimensional structures and the influence of the vibrational entropy on

the activation energies for these processes are discussed.

The fifth chapter deals with the problem of n-type doping of SiC with nitrogen. Ex-

perimental observations concerning the carrier concentration in the implanted sample are

tentatively explained by a split-interstitial based migration mechanism of nitrogen and the

formation of inactive complexes. Also for the recovery of the saturated charge carriers we

discuss an atomistic model, and finally try to correlate our calculated activation energies

with annealing temperatures. A brief section points out the differences in the behavior of

phosphorus and nitrogen as n-type dopant.

Chapter 1

Some Facts about Silicon Carbide

1.1 Historical Facts and Applications

In recent years, silicon carbide (SiC) has become more and more important for the semi-

conductor industry. This field of application is, though, certainly not what people thought

of when the Swedish researcher J¨ons Jacob Berzelius (see Fig. 1.1) discovered SiC in 1824

while trying to synthesize diamond. In 1891, Edward Goodrich Acheson found a possi-

bility to synthesize SiC, melting carbon and aluminumsilicate. He named the material

”Carborundum”, believing it consisted of molten aluminum, which in mineralogy is called

Corundum. When he realised that it consisted of silicon and carbon, this name had already

been adopted into the language (and patented in 1893). Later on, Acheson used carbon

and sand to fabricate ”Carborundum” – a melting process that is still used today, though

in an improved way, and known as Acheson method.

However, this method is due to insufficient crystal quality not used for device fabrication,

but for emery, abrasives and polishing material. In metallurgy, silicon carbide is, further-

more, used for sintering or as a mean for desoxidation or alloys. In the field of ceramics,

SiC is widely used because of its hardness. Today, it is used as the basis for especially high

power ceramics where not only its hardness but also its high melting point is needed. Other

applications of SiC are certain kinds of gaskets, heating rods, overvoltage arresters, and last

but not least high temperature transistors and light emitting diodes (LEDs), just to name

two of the most common applications in the electronic industry. First achievements were

already made in the 1960s, and since then persistently new or improved semiconductor

devices based on extremely pure SiC have been presented by industry. For these applica-

tions, an extremely good material quality is required. For this aim, a sublimation process

Figure 1.1: The ones who found SiC: J. J. Berzelius and H. Moissan

4 CHAPTER 1. Some Facts about Silicon Carbide

Figure 1.2: The upper picture shows three moissanite crystals (0.12 ct) in nearly colorless till light yel-

lowish. The bottom picture shows (from left to right) a synthetic moissanite, a brilliant, and a synthetic

zirkonia.[1]

invented by J.A. Lely (1955) or a modification of it (Y.M.Tairov, 1977) is employed, which

results in monocrystallites with a very low defect density.

Siliconcarbide is not only a synthetical material, but occurs also (rarely) in nature. In

1904, Henry Moissan (see Fig. 1.1) found some small, dark, hexagonal crystallites next to

some small diamond crystallites when investigating meteorites in Canyon Diablo, Arizona.

He analyzed these crystallites as silicon carbide. As it requires extremely high tempera-

tures to be formed, SiC has only been found in meteoritic or volcanic stone, the so called

”Kimberlit”, in the USA, in the former CSFR, and in Sibiria. Mineralogists named the ma-

terial ”Moissanite” in honor of its discoverer. This name is still in use today, and leads us

to a field of application not named so far, where this name of SiC is commonly used: jewelry.

SiC has a Mohs hardness of 9.6, that is very close to that of diamond (10), and its density

is slightly lower than that of diamond. A high dispersion and refraction and a brilliance

only slightly less than in diamond make grinded SiC-stones look extremely similar to bril-

liants, see Fig. 1.2. Since the late 1990s the American company C3, Inc., North Carolina,

managed to fabricate achromatic moissanite of gemstone quality. Now that it is possible to

control the color, high quality grinded moissanites of 0.1 - 2.0 ct. in yellowish, greenish or

brownish colors as well as colorless are grown as monocrystals. The increasing perfection

of these crystals make it more and more difficult to prevent malpractice [1] in spite of many

cleverly thought-out characterization methods .

In spite of this wide range of application of silicon carbide, this work will focus on SiC

as a device material, only. While as a jewel, SiC is always compared to diamond1, in

this area, silicon is the material chosen most often for comparison. Today, for electronic

applications, silicon based devices are most widely-used. Due to intensive research, most

of the technologically relevant problems are already solved on this field. However, silicon

1The investigations of nitrogen–related complexes in Chapter 5 of this work benefit from earlier studies

on N in diamond, as well.

1.1. Historical Facts and Applications 5

is rather limited in its use for high-temperature and high-frequency applications. SiC is,

here, one of the most promising alternatives, because it has all the desired properties. Fur-

thermore, it is based on silicon – a very common material with yet widely developed and

understood technology. This is for compatibility reasons important for the construction of

Si/SiC-heterostructures.

page 2

S I L I C O N C A R B I D E S C H O T T K Y

D I O D E S S A T I S F Y M A R K E T

R E Q U I R E M E N T S

Silicon Carbide devices belong to the so-called

wide bandgap semiconductors,

so the voltage

range for Schottky diodes now can be extended to

more than 1000 V. This is possible because of

benefits specifically related to the SiC material:

• Low leakage currents are possible because the

metal semiconductor barrier is two times higher

than that of Si

• Very attractive, specific on-resistance com-

pared to Si and GaAs Schottky diodes is

possible because of tenfold breakdown field

strength (see Figure 1)

• High current densities are possible and allow

very small die sizes, because the thermal

conductivity is more three times larger than that

of Si (i. e. comparable to copper)

Figure 1 shows the minimum specific on-resistance

of Schottky diodes based on different semicon-

ductors versus the desired blocking voltage (only

drift region, substrate contribution to the resistivity

is neglected).

The ends of each line symbolize the usable voltage

range for the specific semiconductor.

S i C E N H A N C E S S C H O T T K Y D I O D E

B E N E F I T S T O A W I D E L Y E X P A N D E D

V O L T A G E R A N G E

SiC Schottky diodes offer a very low specific on-

resistance with high rated voltages. Figure 2

illustrates the typical forward and blocking

characteristics of 600 V SiC Schottky diodes up to

225 °C. Unlike Si and GaAs Schottky diodes, there

is only a moderate increase in leakage current with

increasing temperature. The area specific differen-

tial on-resistances of a 600 V SiC-Schottky diode

increases from about 0.9 mΩcm2 at room tem-

perature to 1.8 mΩcm2 at 150°C. This positive

temperature coefficient makes Schottky diodes well

suited to implementations in parallel—without the

risk of thermal runaway.

D Y N A M I C P E R F O R M A N C E O F S i C

S C H O T T K Y D I O D E S C L O S E R T H A N

E V E R T O I D E A L D I O D E

When switching off a Schottky diode, there is no

need to remove excess carriers from the n-region

as there is for pn diodes. Hence, there is no reverse

recovery current. Instead, only a displacement

current for charging the junction capacitance of the

diode can be observed. This is shown in Figure 3.

Up to very high frequencies, the current transient

depends solely on the external switching speed.

The charge transported by the displacement current

is very low compared to the reverse recovery

charge Qrr for pin diodes. Due to the different origin

of this charge, we have named it Capacitive Charge

Qc.

Figure 3 compares SiC Schottky Diodes to bench-

mark Si diodes: Si pin double diodes (2*300V serial

in one package) give better reverse current than

ultrafast Si pin benchmark diodes, but have much

higher forward voltage drop. The Capacitive Charge

Qcand the switching power losses of SiC Schottky

diodes are not merely ultra low compared to Silicon

Fig. 2: 6A/600 V SiC Schottky Diode Temperature

Dependence of Forward and Reverse Characteristics

Fig. 1: Comparison of On-Resistances and Blocking

Voltages

1,E-06

1,E-05

1,E-04

1,E-03

1,E-02

1,E-01

10 100 1000 10000

Blocking Voltage [V]

Specific on Resistance [Ω cm²]

SiC

GaAs

0 0.5 1 1.5 2 2.5 3 3.5 4

Forward Voltage [V]

Forward Current [A]

1,E-

0 100 200 300 400 500 600

Reverse Voltage [V]

Leakage Current [A]

T=25˚C

T=100˚C

T=150˚C

T=225˚C

SDP06S60

T=25

T=100

T=150

T=225

Figure 1.3: The diagrams in the left column show a comparison of the specific on-resistances and blocking

voltages of Si, SiC, and GaAs Schottky diodes. On the right hand side, the temperature dependence of

forward and reverse characteristics of a 6 A/600 V SiC Schottky diode is shown. Data taken from Ref. [2]

In 1983, S. Nishino et al. succeeded in growing SiC heteroepitaxially on silicon substrate.

Some years later, in 1987, the first SiC/Si heterobipolartransistor was built. In the fol-

lowing years, the growth procedure of bulk SiC was permanently improved, such that ever

larger monocrystals of device quality could be grown.

In the 1990s, the semiconductor industry concentrated strongly on the realization of unipo-

lar devices, especially Schottky diodes on the basis of SiC. About two years ago, Infineon

Technologies was among the first to fabricate SiC high power devices which are smaller and

cheaper than conventional power diodes on the basis of silicon or gallium arsenide (GaAs)

technology and, furthermore, do not require cooling devices or fans.

Additionally, SiC devices offer high switching frequencies with low losses and high reliabil-

ity. The Schottky barrier of SiC is higher than in the conventional device materials, the

breakdown field strength is by a factor of ten higher than in silicon, and the heat conduc-

tivity is comparable to that of copper, resulting in high current densities and low specific

on-resistances and leakage currents, compare Fig. 1.3. The upper left diagram in Fig. 1.3

6 CHAPTER 1. Some Facts about Silicon Carbide

shows the lower specific on-resistances compared to Si or GaAs at higher blocking voltages.

The diagrams on the right hand side illustrate typical forward and blocking characteristics

of a 600 V SiC Schottky diode at several temperatures. The leakage current increases much

less with increasing temperature than in Si or GaAs Schottky diodes. The bottom diagram

on the left hand side shows that the capacitive reverse charge Qcis nearly independent

from the temperature and dI/dt. A detailed discussion of these characteristics is beyond

the scope of this work. Here, these diagrams shall only illustrate the different behavior of

SiC devices and conventional Si or GaAs devices. Consequences for applications of these

devices are described in great technological detail in Ref. [2].

Schottky diodes based on silicon or GaAs reach their limits at about 200 V or 250 V,

with SiC, though, blocking voltages of 300 V up to 3500 V can be realized [2]. While the

maximum switching frequencies were with conventional materials restricted to less than

100 kHz, the new SiC-Schottky diodes make frequencies of more than 500 kHz possible.

Consequently, the number and size of passive components in integrated devices, i. e. coils

and capacitors, and thereby the costs of the total system, can be reduced.

All these properties make SiC an ideal material for the design of future electronic devices.

There are, as expected for a material that is still – compared to silicon – very new in

this field of application, unsolved problems. One problem to name here is the growth of

high quality crystals with low densities of micropipes and other extended defects. Another

problem is the doping of SiC for the creation of n-type and p-type material as is needed

for device fabrication. Because of the low diffusivities in SiC, doping has to be done by

ion implantation, which, though, damages the crystal lattice dramatically in the implanted

region. ”Repairing” the lattice is commonly done by post implantation annealing, but this

has to be done carefully, if the dopants shall not be removed again from their intended

sites. A lack of understanding of the processes in the annealing phase can be traced back

to a lack of understanding of the behavior of point defects, both intrinsic defects (created

as ”damage” by implantation) and dopant atoms. It is the aim of this work to bring some

light into this field and investigate some of the processes that are supposed to happen

directly after implantation.

In the remaining two sections of this chapter, some currently unsolved problems concerning

unidentified (most probably) intrinsic defects and n-type doping are discussed. Then the

structural principle of the various crystal structures of SiC and the consequences for the

investigations in this work are described.

1.2 Doping of SiC

Due to low defect diffusivities, SiC is commonly doped by ion implantation. Subsequently,

an annealing process reduces the lattice damage created and activates the dopants electri-

cally (i.e. promotes them to substitutional sites). The implantation/annealing procedure

has to be chosen in a way that unwelcome side effects, as e. g. the formation of electrically

inactive complexes, are kept at minimum. Therefore, the investigation of the annealing

behavior of intrinsic defects and dopant atoms is crucial for the optimization of the doping

procedure.

Early spectroscopic studies have revealed a multitude of radiation-induced intrinsic de-

fect centers, but their origin and structure are in many cases still unclear, hindering the

understanding of the dynamics of their creation and annihilation.

1.2. Doping of SiC 7

Intrinsic defects

Among the intrinsic point defects, the properties of vacancies [3, 4, 5, 6, 7, 8, 9, 10] and the

divacancy [11] have been already studied theoretically in great detail. The migration bar-

riers of vacancies and self-interstitials have also been investigated on the zero-temperature

potential energy surface [12, 13].

Experimentally, the annealing of vacancies has been studied since decades ago, based on

their EPR signals [14]. Such studies resulted in an annealing temperature of 200◦C for

carbon vacancies (VC) and of 750◦C for silicon vacancies (VSi ) in as-irradiated samples

[15, 16]. Based on a new assignment for VC, the annealing temperature of VChas been

modified to 500◦C [17, 18, 19, 20]. These values are valid, however, only in the presence of

self-interstitials, when the annihilation of vacancies can occur by recombination. In cases

where only out-diffusion of vacancies is relevant (i.e. in as-grown samples), VSi disappears

above 1400◦C and VCcan still be detected at temperatures as high as 1600◦C [21]. In

addition to single vacancies also divacancies VCVSi are common intrinsic defects in SiC.

Earlier works identified the P6/P7 signal (photoluminescence spectrum) with the diva-

cancy. Recent EPR-measurements and theoretical results could, however, show that this

signal is rather caused by an antisite-vacancy pair VCCSi , as described in more detail in

Chapter 3.

In contrast to vacancies, antisite defects received less attention although their formation

energies are even lower than those of vacancies [22, 23, 24, 25]. In fact, the carbon antisite

(CSi ) appears to have the lowest formation energy among all intrinsic point defects, but,

according to theoretical predictions, it is electrically, optically and magnetically inactive

and, therefore, experimentally almost invisible. In contrast, the calculations predict gap

levels related to the silicon antisite (SiC) in the vicinity of the valence band maximum

(VBM) and the possible existence of paramagnetic states [24].

Experimentally, the presence of SiCwas inferred in Au-irradiated samples based on chan-

neling results [26], and SiChas also been suggested to be the origin of the EI6 electron

paramagnetic resonance (EPR) center in electron-irradiated SiC [27]. The latter assign-

ment has, however, been challenged based on the theoretical calculation of the hyperfine

constants [13]. Still, it is inevitable that antisites are present after irradiation which un-

avoidably creates vacancies and self-interstitials. Namely, the recombination of the latter

during annealing can not only lead to the perfect lattice structure, but also to antisites,

compare Chapter 3. Thus isolated vacancies, interstitials, antisites, and their complexes

will coexist in the crystal in the early stages of annealing.

One of the oldest known irradiation induced intrinsic defect is the so-called DIcenter [28]

observed in 1972 with photoluminescence spectroscopy (mostly in irradiated material but

also in as grown samples quenched to room temperature). In Ref. [28], it was shown that

the center which survives annealing temperatures as high as 1700 ◦C has to be of intrinsic

nature, and the authors suggested the divacancy as atomistic model. Based on calculations

of the local vibrational modes, recent theoretical studies [29, 30] suggested, instead, the

antisite pair as the origin of DI. The creation of such antisite pairs and also larger antisite

complexes is discussed in Chapter 4.

n-type dopants in SiC

Implantation of SiC with nitrogen (N) is commonly used for the creation of n-type mate-

rial. In this work, the experimentally observed problem of the saturation of the free charge

8 CHAPTER 1. Some Facts about Silicon Carbide

carrier concentration at high nitrogen concentrations [31, 32] is investigated. According to

Hall- and DLTS-measurements, no concentrations higher than ≈1019 cm−3can be reached

for the free charge carriers by further nitrogen implantation. A post implantation treat-

ment consisting of further implantation with boron or aluminum and further annealing at

high temperatures can lead to a recovery of the carrier concentration [32, 33].

The processes connected with these observations, complex formation and migration mech-

anisms of the involved defects are subject of Chapter 5.

Phosphorus, which is used as alternative n-type dopant or also as a co-dopant together with

nitrogen in SiC, shows a completely different behavior than nitrogen [31]. In contrast to

nitrogen, our calculations show that phosphorus does not tend to form inactive complexes

(see Section 5.4). Thus higher concentrations of free charge carriers can be reached and

the post-implantation annealing process may be facilitated by using phosphorus instead of

nitrogen.

In Chapter 5, we show how the observation of the different behavior of nitrogen and phos-

phorus can be explained on the atomic scale. Based on investigations of possible mecha-

nisms of their creation, we also suggest some defect complexes that may cause the observed

EPR-signals, the interpretation of which is not yet clear.

The problem of interpreting signals measured with various methods and identifying them

with atomistic models raises the demand for a combined approach of experimental and

theoretical methods. Using and combining various methods is obviously mandatory to

improve the knowledge in defect physics.

On the experimental side, besides annealing studies, measurements of hyperfine constants

by electron paramagnetic resonance (EPR) and related methods, photoluminescence (PL)

spectroscopy, deep level transient spectroscopy (DLTS), and measurements of the Hall

constants give information about the nature and the properties of a defect center, see next

chapter.

A unique identification of the observed defect centers with atomistic models is, however,

most often not possible without augmenting and interpreting the experimental results by

theoretical investigations, which approach the identification problem from the other side

by starting with the atomistic model and deducing the properties which are measured by

various experimental methods.

Due to its ability of simulating large defect structures, the method used in this work

(SCC-DFTB) is suitable for a systematic investigation of migration and aggregation mech-

anisms. In order to obtain an accurate description of the electronic structure of selected

defect complexes, a first principle method should be chosen, additionally. To this aim, ad-

ditional calculations within the plane-wave based code of the Fritz-Haber-Institute (FHI)

[34], the Gaussian orbital based AIMPRO code (”Ab initio modelling programm”) [35],

and the Green’s function based LMTO-ASA code [36, 37] (”Linear muffin tin orbitals in

atomic spheres approximation”) were performed.

1.3 Crystal Structure: the Polytypes of SiC

In SiC, the badly needed understanding of defects is hampered by its characteristic property

of not crystallizing in a unique crystal structure. Unlike most common materials, e. g. its

constituents Silicon and Carbon, SiC can take various structures, the so-called ”polytypes”.

Many different polytypes are known (up to now more than 200), whereof the technologically

most important ones are the cubic 3C and the hexagonal 4H and 6H. This polytypism was

1.3. Crystal Structure: the Polytypes of SiC 9

already observed by H. Baumhauer in 1912, a detailed investigation of this effect was,

however, first made by B.G. Dubrovskii in 1971. All polytypes have the same local order

with fourfold coordinated, sp3-bonded Si- and C-atoms, but differ from each other by the

way, in which the SiC-bilayers, i. e. {111}(cubic polytype) or {0001}planes (hexagonal

polytypes), respectively, are stacked upon each other.

(quasi)cubic hexagonal

Figure 1.4: The most common polytypes of SiC (from left to right): 3C, 4H, 6H. The red line illustrates

the stacking order of the SiC-bilayers along the vertically shown [0001] direction. Inlet:Local structure of

(quasi)cubic (left) and hexagonal (right) sites in a binary lattice.

In the 3C-polytype, all the bilayers have the same orientation, resulting in a zinkblende

structure, stacking order ABC(ABC. . . ). The 3C–lattice is shown in the left part of

Fig. 1.4. In the figure, the [0001] direction2is shown vertically, the view-plane is of {11¯

20}

type. In the figure, two unit cells along the [0001] direction are shown, since the stacking

order leads to a periodicity of three bilayers. As a guide to the eye, the red line illustrates

the stacking sequence of the SiC-bilayers.

Before stacking the bilayers upon each other, they can as well be rotated by 60◦around the

[0001] axis. Performing this operation in defined ways leads to other polytypes with certain

stacking orders. The 4H-polytype, shown in the center of Fig. 1.4 has a periodicity of four

bilayers with the stacking sequence ABCB(ABCB. . . ). On the right, the 6H-polytype

2For better readability, the cubic [111] axis shall be denoted by the hexagonal notation [0001], as well,

which is equivalent but only differs in the choice of the crystallographic unit system.

10 CHAPTER 1. Some Facts about Silicon Carbide

Table 1.1: Properties of the three most common SiC-polytypes. The values are taken from Ref. [38].

The polytype 2H, wurtzite structure, does not exist but is listed for reasons of completeness. Defining a

hexagonality via the percentage of hexagonal lattice-sites, 3C and 2H are the limits of the scale, 4H lying

in the middle.

Polytype Stacking- Atoms per Latticeparameter Bandgap Hexagonality

(Ramsdell) sequence unit cell a[˚

A] c[˚

A] [eV] [%]

3C ABC 6 4.3596 4.3596 2.39 0

2H AB 4 3.0763 5.0480 3.33 100

4H ABCB 8 3.0730 10.0530 3.26 50

6H ABCACB 12 3.0806 15.1173 3.08 33

with a periodicity of six bilayers and the stacking sequence ABCACB(ABCACB. . . ) is

shown.

There exist other polytypes, e. g. the 15R with rhombohedral symmetry and the stacking

sequence ABCACBCABACABCB(AB. . . ), but all the polytypes do not differ much in

their energetical stability, and the main difference comes up already between 3C and 4H

or 6H.

Due to the rotation of a bilayer, the local order changes for the second neighbor surround-

ing of a lattice site. As a consequence, the lattice sites can be divided into (quasi-)cubic

(no rotation of the bilayer) and hexagonal (60◦rotation) sites. The local structure of the

different sites is shown in the small inlet in Fig. 1.4. While the 3C polytype consists of

cubic sites only, there are equal numbers of cubic and hexagonal sites in 4H-SiC. In 6H,

one third of the lattice sites is hexagonal, the other sites can be divided into two different

cubic sites, according to their next neighborhood.

For the investigation of surfaces of the material, especially the nonpolar (10¯

10)-surfaces, the

polytypes have to be treated separately, since completely different reconstruction mecha-

nisms can happen at the surfaces of the same direction in the different polytypes, as we have

investigated for the clean surfaces [39, 40] as well as for surfaces with adstructures. The

passivation of differently oriented surfaces with oxygen has been investigated in [41, 42].

For not so widely extended defects, as point defects and small complexes of them, the

calculated differences between the polytypes turn out to be very small in many cases.

Reflecting on the migration mechanisms of point defects in 4H or 6H, a change between

cubic and hexagonal sites is unavoidable for long-range diffusion processes. A defect that

migrates along the hexagonal axis of the crystal will move from a cubic to a hexagonal site

and so on (consider e. g. migration along the red line in Fig. 1.4). During migration inside

a (0001) plane, i. e. perpendicular to the hexagonal axis, however, migration leads either

from cubic to cubic or from hexagonal to hexagonal sites. Migration barriers vary slightly

because of this difference, but test calculations have shown that these energy differences

are usually small and often below the accuracy of the method.

Most of the investigations that are presented in this work were performed in 3C and 4H.

Variations concerning the symmetry of a small defect complex are most often as small as

the local deviation of 4H from C3vsymmetry. Where it turned out to give more insight or

where significant differences appeared, a distinction is made, otherwise the results are just

1.3. Crystal Structure: the Polytypes of SiC 11

described for one polytype.

This does, of course, not mean that for the description of point defects a distinction between

the polytypes is not needed in general. Attention has e. g. to be paid to the electronic

structure. Table 1.1 shows that the band gap in 3C-SiC is much smaller than in 4H and

6H. This can affect the number of existing charge states of a defect: A localized defectlevel

that is induced in the band gap by a defect in 4H or 6H can already lie in the conduction

band for the same defect in 3C, such that the number of existing charge states of a defect

can change with the polytype.

12 CHAPTER 1. Some Facts about Silicon Carbide

Chapter 2

Theoretical Description of Defect

Dynamics

How the problems revealed in the previous chapter can be treated on a theoretical basis

is subject to this chapter. Before the method of calculation used for atomistic simulations

is described, some general aspects of point defects are summarized as far as they play a

role for our investigations described in the following chapters. This is done on the one

hand with the object of showing up how much information about a problem can be gained

already before starting an atomistic simulation, and on the other hand to make clear which

effects cannot be covered by such simulations. This applies e. g. to the influence of high

defect concentrations on migration processes and charge state effects. Many general state-

ments, though, can be made material independent, and for a qualitative understanding, an

atomistic simulation of a specific material is not needed.

Next, a brief overview of the used computational method and its theoretical background,

density functional theory, is given. The last two sections of this chapter are dedicated to a

description of how migration paths, vibrational spectra and free energies can be calculated

with this method. The latter requires the calculation of entropies, for which extensive

test calculations are presented, which give a classification of the accuracy expected for the

method of calculation as compared to first principle calculations.

2.1 Some General Aspects of Point Defects

Although it is usually spoken about one defect, there is hardly ever only one defect in a

crystal, not even only one kind of defect. How many defects there are in a piece of material

is described by thermodynamics. This is, however, only true in thermodynamic equilib-

rium. As discussed in Section 2.1.2, an implantation process violates this assumption.

2.1.1 Point defects at thermodynamic equilibrium

In thermodynamic equilibrium, the concentration of the defects at a given temperature is

a function of their free energy of formation. For low defect concentrations, the different

kinds of defects can be treated independently, so that the number niof each type iof a

defect at the temperature T is determined by the free energy of formation ∆Gform

iof this

14 CHAPTER 2. Theoretical Description of Defect Dynamics

defect only [43]:

ni=ZiN·e−∆Gform

kBT.(2.1)

Here, Ndenotes the number of available sites in the lattice, while Ziaccounts for the

number of equivalent configurations the defect could take1. Realistic defect concentrations

are most often too high to justify a total neglect of interactions between the defects. Those

are of various origin. The fact that many defects exist in non-neutral charge states for a

given Fermi level results e. g. in a long-range Coulomb interaction between two defects of

unequal charges. This can support the formation of pairs of single defects with unequal

charges, or at least enhance a migration process of them towards each other.

Pair formation can also be caused by short-range interactions between two defects, if e. g.

the local lattice distortion can be reduced in this way. Considering such effects leads

to an expression for the number of pairs of defects that are present in thermodynamic

equilibrium. Since a detailed knowledge of the distribution of the various kinds of defects

is missing and is not even accessible, several assumptions have to be made. In any case,

expressions of the form of Eq. 2.1 can be derived for the number of vacancies, interstitials or

pairs of them, with the factor Ziadjusted to the described problem. For pairs of a vacancy

and an interstitial as an example, it becomes Zi=nvac·nint

N2. In case of pairs, ∆Gform

iis the

binding free energy between the constituents [43].

2.1.2 The situation after an implantation process

Ion implantation is used to create special point defects intentionally. In Chapter 5, the

n-type doping of SiC with nitrogen is discussed in more detail. Actually, the aim of the

implantation procedure is to bring N-atoms into the SiC-lattice, thereby forming substitu-

tional point defects on the carbon sublattice (NC). There they are known to act as shallow

donors, which are required as a source for free charge carriers.

It is, however, not possible to carry out this procedure without creating other intrinsic

defects at the same time. In addition to those defects that are present due to the equi-

librium concentrations under the given circumstances2, the bombardment process creates

besides the implanted species also vacancies and interstitials in equal amounts. Thus, in

the implanted region of the crystal, usually a rectangular profile, the defect concentra-

tion is substantially larger than the thermodynamic equilibrium concentration, as given by

Eq. 2.1.

Implantation is usually followed by an annealing phase. This phase is intended to supply the

energy needed to activate mechanisms for the disappearance of the unwanted by-products

of the implantation process, i. e. recombination or migration towards sinks at surfaces or

dislocations. At the same time the dopant atoms have to be promoted to their intended

positions and kept there if once built in.

In this post-implantation annealing process, the implanted sample is heated up either

gradually (isochronal) or in a procedure called RTA, rapid isothermal annealing, in which

the sample is rapidly heated up to a given temperature that is then kept constant for a

while, typically one half till one hour. The fact that during such an annealing phase also

new defect complexes can evolve can have unpleasant consequences, as discussed e. g. for

the process of nitrogen doping in Chapter 5.

1A substitutional defect with a D2dsymmetry can e. g. occur in three different configurations on a site

with originally Td-symmetry.

2that means defects that have already been created during growth

2.1. Some General Aspects of Point Defects 15

c1 d2

c2 d1

Figure 2.1: Principle ways how migration of point defects can take place in a binary material. The red

circles denote defect atoms, i. e. of intrinsic or extrinsic type, and are supposed to move along the green

arrows. See text for details.

2.1.3 Migration of defects

Regardless of the details of the annealing mechanism, the movement of atoms will con-

sist of one or a combination of several of the mechanisms that are schematically shown

in Fig. 2.13. Probably among the lowest energy migration mechanisms is the process la-

belled with (a), where an atom moves from one interstitial site to a neighboring one. Since

there are different kinds of interstitial sites in the polytypes of SiC with either hexagonal

or tetrahedral character, compare Fig. 1.4, the interstitial can move to an equivalent site

or to a different one. More complicated is the mechanism (b), where an interstitial kicks

out a substitutional atom to an interstitial site that in a next step may kick into another

substitutional site (kick-in/kick-out mechanism). The processes labelled with (c1) and (c2)

are based on split-interstitial mechanisms4. One of the two atoms sharing one lattice site

can either move to an equivalent position (c1) (or change to the other sublattice) or it can

move to an interstitial site (c2).

A substitutionally built in atom can migrate by exchange processes with either a direct (d1)

or a second (d2) neighbor. In the latter case only one sublattice is involved, while a direct

exchange mixes the sublattices. Another indirect exchange process is the basis of the ring

mechanism (e), in which several atoms simultaneously change places. Due to its rather

complicated structure requiring three atoms to move in a defined way, this mechanism is

not very likely to be among the most important mechanisms for the diffusion of a certain

defect.

The vacancy assisted mechanism (f), finally, consists of jumps of defect atoms into a va-

cancy. Like the exchange process this can happen on one sublattice exclusively or go along

with a change of the sublattice, depending on the capture radius of the vacancy (see also

Section 3.3). A directed long-range diffusion of a defect via vacancy jumps can then be

imagined as the three step sequence of a jump of the defect atom into the vacancy, the

dissociation of the vacancy from the defect and the arrival of a new vacancy. Alternatively

– especially in case of high binding energies between the defect atom and the vacancy – a

3The lattice shown here is of the cubic 3C-type, but of course the same holds for any crystal lattice.

4In Chapter 5 it will be shown that such split-interstitials play an important role for the mobility of

nitrogen in SiC.

16 CHAPTER 2. Theoretical Description of Defect Dynamics

ring-like movement of the vacancy around the defect atom is conceivable5.

2.1.4 Jump probabilities

It depends on the activation energy, which of the mechanisms discussed in the previous

section describes the migration of a certain defect, or at least is the main mechanism

responsible for a certain observation. Decisive are, i. e., the migration energy barriers for

the process and the probability that all necessary conditions for the process are fulfilled.

For the sublattice movement of an interstitial (Fig. 2.1 (a)) or a vacancy (Fig. 2.1 (f)) no

further assumptions are required, thus the jump probability is given by

K=ν·e−∆Gact

kBT(2.2)

with the height ∆Gact = ∆U−T·∆S+p·∆Vof the free energy migration barrier and

an ”effective frequency” νfor the jump [43].

The more complicated mechanisms require, additionally, the presence of other atoms or

vacancies, such that e. g. for the vacancy assisted mechanism the probability for a jump

is given by

K=ν·e−∆Gact

kBTe−∆Gform

vac

kBT.(2.3)

In this case, the jump probability depends on the free energy of formation of the assisting

vacancy.

For a migration process with equivalent start- and end-configurations, the rate constant

Kcan be further evaluated. The mainly on classical mechanics based rate theory leads to

the expression

ν=ν0·e

∆Smig

kB(2.4)

for the effective frequency ν[43]. It depends exponentially on the migration entropy ∆Smig,

i. e. the difference in entropy between the equilibrium- and the saddle point geometries.

ν0is the vibration frequency of the migrating atom in the equilibrium position.

For selected processes leading from a start-configuration to an equivalent end-configuration,

the jump probabilities and diffusivities have been calculated in Chapters 4 and 5.

2.1.5 Charge state effects on the migration of defects

The Fermi level determines, in which charge state a defect is stable. In thermal equilibrium,

it is constant, but especially in the implanted region, the Fermi level varies locally due to

an inhomogeneous distribution of the various defects, and consequently a defect can exist

in several charge states q1,q2with concentrations α1N,α2Ndetermined by the energetical

distance between the localized defect level EDand the Fermi level EF. The fraction α1of

defects in charge state q1compared to the fraction α2= 1 −α1of defects in charge state

q2is then α1N

α2N=α1

1−α1

=gD·e

ED−EF

kBT(2.5)

where the factor gDaccounts for the degeneracy of the defect level ED[43].

The migration enthalpy is changed accordingly, and as a consequence also the jump prob-

ability. Considering the coexistence of two charge states as above, the (total) jump prob-

ability can be written as the sum of the probabilities for a jump in each charge state:

Ktot =α1K1+α2K2=α1K1+ (1 −α1)K2(2.6)

5The vacancy assisted processes play an important role for the mobility of antisites in SiC, see Chapter 4.

2.2. Density Functional Theory 17

The effect of the lowering of an activation energy due to a change of a defect into another

charge state is known as ”ionization enhanced migration” and can in some cases have a

strong influence on the diffusion rate.

As in common literature in this work, migration of defects is only investigated for a non-

changing charge state during the process. Such charge state effects as well as the hindering

or enhancement of a migration process by the presence of other defects can not be ac-

counted for in a satisfying way, because such effects strongly vary with the environmental

conditions and do, therefore, not allow to deduce a general picture. Nevertheless, a super-

position of various competing and partly contrasting effects on the activation energies has

to be expected, such that the analysis of the pure processes might give a good description

of the average behavior of the investigated defect.

The considerations of the previous sections provide a qualitative understanding, especially

for general (i. e. not material dependent) information about the behavior of defects in

solids. Quantitative statements require, though, the use of more sophisticated methods.

In this work, an approximative method, based on density functional theory (see next

section), has been used for atomistic simulations. This method, the density functional

based tight binding (DFTB), has been described in great detail in various publications, see

Ref. [44, 45, 46], thus only a brief sketch of the idea shall be given here. The extensions

made in SCC-DFTB (selfconsistent charge- ) are illustrated in Ref. [47].

2.2 Density Functional Theory

In quantum mechanics, the properties of a crystal are usually described by a wavefunction

Ψ(R1, R2,...,RNnuc ;r1, r2,...,rNel ) that depends on the coordinates Riof Nnuc nuclei

and rjof Nel electrons. It is practically impossible to find the ground state of such a

system by minimizing the total energy of this many-body system , since a Schr¨odinger

equation for these 3Nnuc +3Nel coordinates would have to be solved6. The most well-known

approximation that is made to reduce this problem is the Born-Oppenheimer approximation

[48], where it is utilized that the motion of the nuclei due to their larger masses happens

on a substantially larger time-scale than the motion of the electrons. In other words, the

electrons follow the movement of the nuclei practically immediately. The motion of the

nuclei can therefore be decoupled from that of the electrons, and the key problem remains

to determine the electronic wavefunction for a given configuration of nuclei. In spite of this

reduction, for practical applications the number of parameters has to be reduced further

on7. One possibility to achieve this is to use the density functional theory (DFT).

Instead of using a wave function, the key quantities of the system are expressed in terms

of a particle density

n(r) = hΨ|X

δ(r−ri)|Ψi.(2.7)

The total energy of a quantummechanical system can then be written as a functional of

the electron density

E=E[n(r)] .(2.8)

6Consider that in a crystal a typical order of magnitude are ≈1023 atoms!

7For Nparticles the computing time grows > N5. Chosing twice the number of atoms in the model

leads therefore to a 32 times longer computing time.

18 CHAPTER 2. Theoretical Description of Defect Dynamics

That this is uniquely possible for a (non-degenerate) ground state of a system of electrons

has been shown by Hohenberg and Kohn [49].

For a system of Nelectrons, Ehas – in atomic units – the form

E[n] = T[n] + 1

2ZZ n(r1)n(r2)

|r1−r2|dr1dr2+Zn(r)Vext(r)dr,(2.9)

where T[n] is the kinetic energy of Ninteracting electrons, the second term is the Coulomb

interaction of all electrons, and the third term accounts for the external potential Vext for

the given configuration of the nuclei.

The kinetic energy T[n] is generally unknown, thus further approximations have to be

made. Starting with the kinetic energy of Nnon-interacting electrons

T0[n] =

i=1 ZΦ?

i(r)−∇2

2Φi(r)dr(2.10)

with

n(r) =

i=1 |Φi(r)|2,(2.11)

Kohn and Sham proposed to inherit this density for the description of interacting electrons

[50] and include the many-body effects in T[n] then by adding the so-called exchange–

correlation functional Exc[n]

T[n] = T0[n] + Exc[n].(2.12)

Most, but not all systems can be treated like this, but the validity of Eq. 2.11 is limited

to those states that can be expressed in a single Slater determinant. In cases where this

is not possible, other theories, e. g. Hartree-Fock methods, have to be used to obtain an

exact description of the electronic structure. This is, however, not always required, and at

least some properties of interest can also be obtained within the DFT description.

The wave functions Φi(r), the Kohn-Sham orbitals, can be obtained by the variation of

E[n] with respect to Φi(r)?and under the condition of constant particle number N. This

leads to the Kohn-Sham equations

−∇2

2+Veff [n(r)]Φi(r) = εiΦi(r) (2.13)

with the Lagrange parameters8εiof the constraint of constant Nand the effective potential

Veff =Vext(r) + Zn(r0)

|r−r0|dr0+δExc[n(r)]

δn(r).(2.14)

Equations 2.13 are coupled via the effective potential Veff and the exchange–correlation

potential Vxc =δExc[n(r)]/δn, which depend on all Kohn-Sham orbitals Φi(r), such that

a self–consistent solution is required.

The exchange correlation contribution to the total energy is generally unknown, but for

special cases of only small variations of this contribution, approximations based on a ho-

mogeneous electron gas can be used. Today, the most common approximation is the local

8Only if Koopman’s theorem is valid, i. e. if the removal of one electron does not lead to a change in the

Φiof the remaining N-1 electrons, the εican be interpreted as single particle energies. For large numbers

Nof electrons, it is probable that this presumption holds approximately. In general, this is, however, not

valid.

2.2. Density Functional Theory 19

density approximation (LDA)[51, 52, 53], which assumes that Exc can be expressed as a

functional the integral core of which only locally depends on the density n(r)

Exc[n] = Zn(r)Vxc[n(r)]dr≈Zn(r)VLDA

xc (n(r))dr.(2.15)

This approximation is well established and its limits as well as certain possible improve-

ments are known. As the probably most common example of disadvantages, the incorrect

description of the band gap in semiconductors has to be named. LDA usually underesti-

mates its width, which can lead to problems with the decision whether a certain charge

state of a defect exists or not. The Scissor operator [54] or the Baraff–Schl¨uter correction

[55] are examples of techniques to account for this failure.

Inserting the expressions for the exchange–correlation contribution and the kinetic energy

into Eq. 2.9 and considering the energy contribution of the nuclei, the final expression for

the total energy reads

Etot =

occ

njhΦj|ˆ

H|Φji− 1

2ZZ n(r1)n(r2)

|r1−r2|dr1dr2−Zn(r)Vxc[n(r)]dr(2.16)

+Exc[n(r)] + 1

Nnuc

I,J

ZI·ZJ

|RI−RJ|.

The terms in Eq. 2.16 have been slightly rearranged in order to obtain Etot in the form

needed in the following. The hamiltonian ˆ

Hin the first term contains the kinetic energy

part, the external potential, the Hartree energy and the exchange and correlation potential,

H=−1

2∇2+Vext +Zn(r0)

|r−r0|dr0+Vxc[n(r)] ,(2.17)

the second and third term in Eq. 2.16 correct for the double counting terms within this

hamiltonian. The fourth term of Eq. 2.16 is the exchange correlation energy, the last term

accounts for the nuclear repulsion, Eion−ion.

In principle, point defects, extended defects, surfaces, etc. could be described by first prin-

ciples within LDA-based DFT as presented here. For the description of larger defects, i. e.

not only extended defects as e. g. dislocations but also complexes of point defects, this

pure DFT is much too elaborate. Although the time scaling behavior could be reduced

from > N5in many-body theory [56] to ∼N3in LDA-DFT, the required resources exceed

most often the capacity of today’s advanced computing facilities. Molecular dynamical

simulations or as well static calculations of migration paths of defects – both procedures

which are in principle based on a large number of single total energy calculations –, are

much too time consuming in pure DFT.

Anyway, in many cases the good accuracy of this pure DFT is not needed, and more ap-

proximative methods can be used. Many different approaches have been made to develop

DFT-based methods with a good performance on the one hand and a low loss of accuracy

compared to pure DFT. A great variety of computer codes exists – the main difference

between them is the basis set used for the wave functions. Very common is the use of a

plane wave basis, as e. g. in the FHI-code [34], which has been used as a reference for some

calculations in this work. These methods are commonly considered to have a very high

accuracy, although, of course, also here typical drawbacks are known.

Considerably more approximative are tight-binding methods, which use a local basis set.

20 CHAPTER 2. Theoretical Description of Defect Dynamics

Most of these methods are, however, not based on density functional theory but are so-

called semi-empirical methods, i. e. depend on parameters obtained by fitting certain

quantities to experimental data. A clearly lower accuracy and often a bad transferability

of the parametrized description characterize many of these semi-empirical tight-binding

codes.

Both approaches have their advantages and disadvantages, and the method used in this

work (DFTB) tries to combine the accuracy of DFT-based methods and the high efficiency

of tight-binding approaches. Generally, the accuracy of this method, in which no empirical

parametrization is used, comes closer to that of first principle codes than that of con-

ventional empirical tight-binding methods. The basic ideas are described in the following

section, while for a detailed description, the reader is referred to Refs. [45, 46, 47].

2.3 SCC–DFTB

Starting with expression 2.16, the total energy used in DFTB can be derived, following an

idea of Foulkes and Haydock [57]: They suggested to expand the total energy Etot at a

reference density n0(r). Fluctuations δn(r) = n(r)−n0(r) are then treated as perturbation.

With this, Eq. 2.16 becomes

Etot =

occ

njhΦj|ˆ

H0|Φji− 1

2ZZ n0(r1)n0(r2)

|r1−r2|dr1dr2−Zn(r)Vxc[n0]dr+Exc[n0] + Eion−ion

2ZZ 1

|r1−r2|+∂2Exc

∂n(r1)∂n(r2)|n0δn(r1)δn(r2)dr1dr2(2.18)

with the ”bandstructure term” (first term in the first row) and all the other terms of

Eq. 2.16 depending only on the reference density n0(r), and a second order term (second

row) quadratically depending on the density fluctuations δn(r). Expressions linear in δn(r)

cancel exactly in the expansion. This equation can be rewritten as

Etot ≈

occ

njhΦj|ˆ

H0|Φji+Erep +1

µX

γµν∆qµ∆qν(2.19)

with the repulsive potential Erep and the second order term, depending on the Mulliken

charges. This term will be discussed later in this section.

In the standard formulation of the method, DFTB, only the first two terms were consid-

ered, while higher order corrections were neglected, and especially for homonuclear systems

a good description of various material properties was achieved. For systems consisting of

more than one species, and especially in case of strong charge transfer between the different

atoms (strong ionic bonding character), the extended version of the method, SCC-DFTB,

where ’SCC’ stands for ’self–consistent charge’ and the second order term is calculated,

yields clearly improved results 9.

The main distinction of the DFTB method from other (semi-empirical) tight-binding meth-

ods is the way how the hamiltonian for the bandstructure term and the repulsive interac-

tions are treated. In the original tight binding formulation of Slater and Koster [58], the

exact many-body hamiltonian operator is expressed with a parametrized hamilton matrix.

9In this work SCC-DFTB was used throughout.

2.3. SCC–DFTB 21

While the matrix elements were then obtained from fits to the experimentally determined

electronic structure of reference systems, in the DFTB method any empirical parameter-

ization is avoided. Instead, the hamiltonian and the overlap matrices are obtained from

preceding selfconsistent LDA-DFT calculations of suitable reference systems. The repul-

sive potential is then as a superposition of short-range repulsive two-particle potentials

constructed as the difference of the cohesive energy calculated within SCF-LDA and the

bandstructure energy calculated in the tight-binding approximation. The quality of the

obtained parameterization depends, of course, on the choice of the reference structures.

They should, therefore, be chosen carefully in a way to describe as many different proper-

ties of the material as well as possible. The advantage of this procedure over empirically

parametrized tight-binding methods is a substantially better transferability and results

closer to first principle calculations.

To evaluate the expression in Eq. 2.18, an initial density n0has to be known. In DFTB,

this initial or reference density is constructed as the superposition of atomic-like densities

nα

0centered at the atoms α:

n0=X

nα

0.(2.20)

The Kohn-Sham wavefunctions Φiare expanded in localized Slater-type orbitals ϕν:

Φi=X

cνi ϕν(r−Rα),(2.21)

where the ϕνare obtained from the solution of a modified Schr¨odinger equation:

[ˆ

T+Veff [nα

0] + r

r02

]ϕν(r) = ενϕν(r) (2.22)

with a compression radius r0depending on the covalent radius of the atom10. In numerous

calculations in various materials [59], the basis functions ϕνobtained with the compression

term ( r

r0)2have been shown to yield a better description of the energy and the geometrical

quantities, i. e. bond lengths and angles, than without this compression11[59].

The coefficients cνi in the expansion of the ϕiare obtained as solution of the secular

equation

cνi(H0

µν −ενSµν) = 0 (2.23)

which is obtained from the variation of the energy with respect to the cνi.H0

µν and Sµν

are the hamiltonian and overlap matrices: H0

µν =hϕµ|ˆ

H0|ϕνiand Sµν =hϕµ|ϕνi. The

diagonal elements of H0

µν are the energies of free atoms εatom

µ.

In the first versions of this method, the starting density n0was expressed as a superposi-

tion of the potentials of the (pseudo-)atoms. Later on, it turned out that a superposition

of the electron densities of these (pseudo-)atoms leads to a better description [45]. The

parameters used in this work were all created by the superposition of densities12. In any

case, the matrix elements H0

µν and Sµν depend only on the reference density n0and can be

tabulated with respect to the distance of two atoms. The repulsive potential constructed

as mentioned above can similarly be tabulated depending on this distance, so that the

10The standard value is r0≈1.85 ·rcovalent. The indices µand νin the summation are an abbreviation

for the quantum numbers n,l, and mthat appear in the expression for the Slater-type orbitals.

11This holds not only for solids but also for the description of molecules.

12For the parameterization of silicon, a wavefunction and a density compression radius r0of 6.7 and 3.3

was used, for carbon 7.0 and 2.7, for nitrogen 11.0 and 2.2, and for phosphorus 9.0 and 3.8 (in aBohr).

22 CHAPTER 2. Theoretical Description of Defect Dynamics

solution of Eq. 2.23 with these tabulated values yields the DFTB-energy of the system.

The second order term in Eq. 2.18 is calculated from the Mulliken charges ∆qα=qα−q0

with

qα=1

occ

niX

µ∈α

(c∗

µicνi Sµν +c∗

νicµi Sνµ).(2.24)

The function γαβ in Eq. 2.19 denotes the monopole approximation of an expansion of the

density fluctuations δn in this second order term in a series of radial and angular functions.

The analytic expression for γαβ is derived in Ref. [60]. In the limit of large atom distances

R−→ ∞,γαβ describes the Coulomb-interaction between two point charges, γαβ −→ 1

For one atom, i. e. R= 0, γαβ =γαα has the value of the Hubbard Parameter Uα(for

details see Refs. [45, 60]).

This construction of the second order correction necessitates a self-consistent treatment in

the energy calculation.

The interatomic forces for use in the conjugate gradient scheme or in molecular dynamic

simulations are obtained by taking the derivative of the total energy of Eq. 2.18 with re-

spect to the nuclear coordinates.

In the past, SCC-DFTB has been used successfully for the description of systems con-

taining up to approximately 700 atoms. The tabulation of the matrix elements make the

performance ideal for molecular dynamics simulations and diffusion processes. Such prob-

lems can hardly be treated by first-principle methods. The accuracy of SCC-DFTB has

to be classed between that of usual tight-binding and first-principle methods. Although

the description of geometries and total energies can compete with that of ab initio calcu-

lations, the electronic structure is clearly less accurate. It has therefore turned out to be

a good approach to combine several methods, such that the time-consuming investigation

of many different structures, molecular dynamics simulations, and diffusion mechanisms

are done within SCC-DFTB. For selected structures, ab initio calculations can then be

used to calculate the electronic structure, occupation levels, excitation energies, hyperfine

parameters etc., more accurately. This procedure has successfully been used in several of

our works, see e. g. Refs. [10, 61, 62, 63].

2.4 Modeling of Defects

What is desired from theory, is information about the properties of a single defect in an

otherwise perfect crystal. Although today’s computer facilities and the efficiency of the

methods generally applied to the modeling of defects are continually growing, the number

of atoms that these methods can treat is far from that in a real crystal. For the description

of an isolated defect, the simulation of an infinitely extended crystal would come closest

to realistic conditions. Except for the perfect crystal, this is, however, impossible within

usual atomistic simulations13. Only a bounded area around the defect can be described,

for which there are several possibilities.

The most common two possibilities are on the one hand the use of clusters and on the

other hand the use of periodic boundary conditions, compare Fig. 2.2.

In a cluster consisting of a part of the investigated material and a surface passivated

with (pseudo-)hydrogen atoms one isolated defect can be modeled. The passivation of the

13An exception is the elaborate Green’s functions approach [36, 37], where the infinite crystal can be

described with a Green’s function, so that a really isolated defect can then be modeled.

2.4. Modeling of Defects 23







































HH HH H H H HH HHH

HH H H H HH HHH



























  

!!!!

""""

####

$$$$

%%%%%

&&&&&

Figure 2.2: A defect modeled in a supercell (left) or in a (pseudo)hydrogen passivated cluster (right).

cluster surface is necessary to prevent the dangling bonds of the surface atoms. Since only

bulk material shall be simulated, the geometric positions and the charge transfer between

the outer atoms are modeled to be as they are in ideal bulk. This is achieved by defining

so-called pseudo-hydrogen atoms, which then have to be put in the right position to serve

their purpose. Nevertheless, the results usually depend on the symmetry of the cluster and

on the position of the defect14. Variations of the geometric relaxation and the energy of

the structure due to different distances of the defect to the passivated cluster-surface may

not be underestimated. Using a sufficiently large cluster can, of course, help to keep such

artificial effects small, but result in large numbers of atoms, also at the cluster-surface.

The other approach is to use periodic boundaries. This is usually achieved by using the

supercell concept, in which it is assumed that all margin atoms at one side of the cell

are neighbors to those atoms at the opposite side of the cell. Like this, interactions of

all atoms across the supercell boundaries are accounted for and no artificial surfaces as in

the cluster approach evolve. As a price, the risk of unphysical defect-defect interactions

between defects in ”neighboring cells” due to the created periodic images rises, since no

longer an isolated defect but rather a periodic grid of defects (with a ”lattice constant”

induced by the supercell vectors) is modeled, compare Fig. 2.2. Also here, using a large

supercell helps to reduce such influences till they are negligible. Because of the high com-

putational costs, especially first principle methods are limited to rather small numbers of

atoms, which consequently sets a limit to the size of defect complexes that can be described

in these methods.

The influences of periodic images in case of supercells or of passivated surfaces in case of

clusters can also be kept small by allowing only the atoms in a ”core region” around the

defect to relax while all outer atoms are kept fixed at their bulk positions. This, however,

can put unphysical constraints on the defect and its surrounding, thus a good compromise

has to be found in choosing such ”boundary conditions”.

In this work, all defect structures were modeled in supercells. The convergence has been

tested at several examples, and according to these results the supercells for further use

have been chosen. Calculations performed for the isolated vacancy in silicon have shown

14It is often difficult to position larger defects with an unsymmetric relaxation into a cluster, which

preferably should have spherical symmetry because of its surface.

24 CHAPTER 2. Theoretical Description of Defect Dynamics

a clearly better convergence of the results (geometry and formation energy) with the size

of the supercell in SCC-DFTB [61] than in the reference calculations with a plane wave

code [64]. Obviously, the long-range effects responsible for the interaction between periodic

images are suppressed more strongly in the tight-binding approach than in the plane wave

description (compare also the discussion in Chapter 3). This implies that the quality of

a result obtained in different methods but with the same boundary conditions cannot in

general be assumed to be similar. Independent from how elaborate the used method is,

there is always a need for test calculations and convergence checks. In Chapter 3 formation

energies and migration energies of vacancies (and pairs of them) in 3C-SiC and 4H-SiC have

been calculated in various supercells with varying boundary conditions, in Chapter 4 the

same has been done for antisites.

In numerical simulations, the integrals over periodic functions in real space (compare the

expressions in the previous sections) are most often evaluated by summation over their

Fourier transforms in the reciprocal space. A k-point sampling of the Brillouin zone is the

basis of this summation, which corresponds to an integration over the whole real space.

Utilizing symmetries the number of k-points by which a good description of the system

is achieved can be reduced substantially. The larger the structure in real space, the less

points are needed for the sampling in reciprocal space. Especially first principle methods

which are limited to smaller numbers of atoms in the simulations make, therefore, use of

larger k-point schemes together with small supercells. For large supercells, a small number

of k-points suffices to obtain a good description, so that very often only one k-point, namely

the Γ-point k= 0, is used15.

For selected structures the dependence on the k-point scheme has been checked for the

SCC-DFTB calculations, but differences were usually below 0.1 eV in energy and also the

geometrical differences were very small if only the Γ-point or an optimized special k-point

scheme was used. Tests were made for the common (2×2×2) Monkhorst-Pack scheme16

[65]. The plane wave and the AIMPRO reference calculations used either a small cell of 64

atoms with this special k-point scheme or a bigger cell with 128 or in selected cases also

216 atoms and only the Γ-point. In this work, the Γ-point was used throughout.

Unless otherwise stated, a (3×3×3) supercell, containing 216 atoms, was used for 3C-SiC,

and a (5×6×1) supercell with 240 atoms was used for 4H-SiC. Usually, as many atoms as

possible were allowed to relax.

2.5 The Diffusion Algorithm

When a defect complex has been calculated to be stable against dissociation (under certain

environmental conditions), the next question that arises is, how it can be created, or, how

it will behave during e. g. an annealing process. At which annealing temperature a defect

complex dissociates or a vacancy or an interstitial starts to migrate through the crystal

lattice is determined by the change in total energy of the structure during the migration

process. Thus, the initial and the final structure of the migration process are known, and

what is in demand is the saddle point geometry that separates them and whose energy

difference to the initial structure defines the activation energy of the process.

15The Γ-point has also some advantages compared to other k-points. Band dispersion and broken degen-

eracy of localized defect levels at k-points other than the Γ-point make it difficult to detect the position of

these levels relative to the valence band physically correct.

16Monkhorst and Pack have shown that the integration of periodic functions of a Bloch wave function

over the entire Brillouin zone can be reduced to an integration over ”special points” chosen in a special way

for each crystal lattice, so that the computational effort is reduced essentially [65].

2.5. The Diffusion Algorithm 25

Starting with the initial structure, the geometrical path of all atoms in the defective region

has to be determined in order to find the way that leads over the lowest energy barrier on

the potential energy surface, i. e. the saddle point geometry, towards the final structure.

While a minimum energy structure can in principle be obtained by relaxation of the system

using a simple conjugate gradient scheme based on the minimization of the forces, some

constraints have to be applied to the system for the calculation of migration paths.

Two different techniques were used in this work and are briefly described in the following

paragraphs.

2.5.1 The constrained relaxation technique

The following algorithm [66] is a possibility to calculate the whole migration path of an

atom from its position in the initial structure to its position in the final structure. It puts

constraints onto selected atoms that are supposed to undergo the largest changes during

the process, while all other atoms are allowed to relax freely.

1. Conjugate gradient relaxation of initial and final

structures (DFTB) ⇒rdiff (3N-dim)

2. Move mass center of selected atoms by ∆rdiff

3. Constrained conjugate gradient relaxation in planes

⊥rdiff (DFTB)

4. Calculate new rdiff to final structure

5. Repeat (2)-(4) until final structure is reached

6. Calculation of the vibrational spectrum of the saddle

point geometry

7. If required steepest descent relaxation from saddle

point geometry to equilibrium geometries

diff

rfinal rinitial





















The number of steps to take between the two equilibrium positions depends on the shape

of the potential energy surface. For a migration path along a rather flat energy surface,

less steps are required to find the saddle point geometry with the same accuracy as for a

migration path on an energy surface with large curvature in the range of the saddle point17.

The number of atoms which are constrained during the relaxation varies. In many cases

it is sufficient to constrain only one atom, e. g. for the sublattice migration of vacancies

the atom that changes its site with the vacancy. If more atoms are supposed to perform

a simultaneous movement, it can become necessary to constrain all these atoms. The

activation energies obtained like this can be higher due to these constraints, so that the

calculation possibly has to be refined in the region around the saddle point until the correct

saddle point geometry has been found. That a saddle point geometry has been reached can

be checked by calculating the vibrational spectrum of this structure. For a saddle point

structure, the spectrum will show one imaginary mode, belongig to the migrating atom

and indicating the direction of migration towards the two equilibrium structures.

From a saddle point geometry found with this technique, a steepest descent relaxation

without constraints is performed to obtain the whole migration path. With a very small

17This is one point where the efficiency of the method pays off. The large number of calculations along the

migration path with first principles, e. g. plane wave, codes is very time-consuming or requires parallelization

and the use of supercomputers, while SCC-DFTB-calculations can be performed on common workstations

in acceptable time.

26 CHAPTER 2. Theoretical Description of Defect Dynamics

movement of the constrained atoms in the positive and negative directions of the vibrational

mode belonging to the saddle point, relaxation will automatically lead to the equilibrium

geometries.

There are cases in which this technique is not very successful, and finding the saddle point

geometry is hindered by the applied constraints. In these cases, the technique described in

the following paragraph might help.

2.5.2 The activation relaxation technique (ART)

An alternative technique to determine the saddle point geometry is the activation relaxation

technique (ART) [67]. In contrast to the technique described above, no information about

the whole path is obtained, but a steepest descent relaxation is required afterwards.

Starting with an equilibrium geometry, the atoms are moved slightly in the direction of

the saddle point 18. A new force Gcan be defined with the 3N-dimensional force vector F

of the current structure and the 3N-dimensional unit vector ∆xpointing from the (initial)

equilibrium structure to the current configuration:

G=F−(1 + α)(F·∆x)∆x.(2.25)

Like in the usual conjugate gradient calculation, this redefined force Gis calculated and

the atoms moved according to it (instead of Fin the usual relaxation). The component

of Gthat is parallel to ∆xhas the opposite sign of F, so that the system is forced to

perform a valley up-hill movement in the energy landscape until the saddle point geometry

is reached. At this point, both forces Gand Fare zero. The positive constant αcan be

varied in order to accelerate the convergence behavior, usually α= 1 is a good choice.

After this activation phase, a relaxation phase yields the path from the saddle point struc-

ture into the equilibrium configuration. As described in the previous paragraph, a small

movement of selected atoms along both directions of the saddle point vibrational mode

and steepest descent relaxation leads the system to the equilibrium structures. In contrast

to the algorithm described afore, only the activation energy for the process is obtained

directly in the first step, while the migration path has to be calculated separately in a

second step.

2.6 Lattice Vibrations and Free Energies

For the calculation of formation energies of defects or activation energies of migration

processes, the energies of different structures have to be compared. This is most often done

by simply comparing the total energies Etot obtained from e. g. supercell calculations.

Spoken thermodynamically, however, the crystal (with the defect) is a great canonical

ensemble, and it would be thermodynamically correct to compare the Gibbs free enthalpies

G=U−T·S+p·V+X

µiNi(2.26)

of the structures. The last term accounts for the change in particle number Niwith the

chemical potential µi.19 In this work, calculations have been performed at constant

18Having a rough idea of the saddle point geometry is of advantage, since it can abbreviate the calculation

time massively to have a good initial configuration. Being too far from the saddle point geometry can in

some cases result in a loop around this configuration instead of converging towards the saddle point.

19In a diffusion process, in which the particle number remains constant, this term does not contribute,

but only e. g. for the calculation of the formation energy of a defect from perfect bulk. Since it does not

2.6. Lattice Vibrations and Free Energies 27

pressure and volume, so the pressure- and volume dependent term pV is not relevant, here.

The internal energy Uconsists of two parts:

U=Etot +Uvib (2.27)

The first term, the static part, is the total energy Etot, as obtained from the SCC-DFTB

calculation. The term Uvib is caused by lattice vibrations. Assuming a Planck-distribution

of harmonic oscillators, Uvib can be written as

Uvib =

i=1 



ωi

exp(



ωi/kBT)−1+1



ωi,(2.28)

following the Einstein model in statistical physics [71]. Here, ωiare the eigenfrequencies

obtained from the calculation of the vibrational spectrum of the defect. Tis the tempera-

ture.

The entropy Sin Eq.2.26 is usually neglected. This is justified in many cases by the

similarity of the structures that are compared, so that the change in entropy, ∆Sis, indeed,

very small. In some cases, as e. g. strong rearrangements of the lattice, it is, however,

important to consider the term T·S, especially at high temperatures (compare Section 4).

The entropy Sconsists of the configurational, the electron-hole-pair and the vibrational

entropy. The configurational term is given by the number of possible configurations in

which the defect can exist. The second term is usually negligible small and shall also not

be discussed further. The prevailing contributions come from the lattice vibrations and

can have a non-negligible influence on the energetics – depending on the temperature.

As described in detail in Appendix B, we can with

∂S

∂U V,N

T(2.29)

and Uvib from Eq. 2.28 derive an expression for the vibrational entropy Svib. The entropy

then takes the form

Svib =kB

i=1 (



ωi

kBTexp 



ωi

kBT−1−1

−ln 1−exp −



ωi

kBT).(2.30)

The frequencies ωiare the local vibrational modes of the defect structure. The vibrational

frequencies ωican be calculated within SCC-DFTB from the dynamical (Hessian) matrix

Dij =1

√mkml

∂2E

∂rik∂rjl

(2.31)

that is defined by the second derivatives of the energy, having the eigenvalues ω2

i. For

the calculation, each atom in the defective region is slightly displaced by a distance ±∆ri

in three directions. Based on the finite differences in the forces in these slightly different

geometries, the Hessian matrix can be built up. Its diagonalization yields the desired fre-

quencies.

change anything in the derivation and we mainly will apply this formalism to diffusion processes, we leave

this term in the following for clarity.

28 CHAPTER 2. Theoretical Description of Defect Dynamics

Unfortunately, the entropy is not easily accessible, neither theoretically20 nor experimen-

tally, making a comparison difficult. Only a few reference values, mostly force field calcu-

lations and partly of doubtful quality, are available, and the discrepancies between these

values are in some cases larger than the error due to a neglect of the entropy term. How-

ever, a classification of the quality of the calculated vibrational frequencies and of the

calculated entropies can be made for the vacancy in silicon. This requires the calculation

of the vibrational spectra of the perfect lattice and the lattice with the vacancy.

Before this will be done, it might be meaningful to investigate the quality of the calculated

vibrational spectrum of only one structure (instead of comparing two of them). First,

the vibrational spectrum of ideal bulk 3C-SiC will be discussed and the SCC-DFTB-result

compared to other methods, then the specific heat will be calculated for different materials,

for which experimental values are available in the literature.

2.7 Applications and Test Calculations

2.7.1 The vibrational spectrum of 3C-SiC bulk

At the example of a supercell of 3C-SiC containing 64 atoms, the vibrational spectrum

calculated within SCC-DFTB in Γ-point approximation has been compared to the spec-

tra obtained within AIMPRO [35, 68] and FHI[34, 69]. These spectra, broadened with

10 cm−1are shown in Fig. 2.3. The quantities discussed in the following sections were

calculated using the original frequencies as done with the SCC-DFTB-spectra. The gap

FHI

DFTB

3C−SiC,

64 atoms

AIMPRO

wave number [1/cm]

density of states, arbitrary units

−0.6

−0.4

−0.2

0.2

0.4

0.6

0.8

−0.8

1000

1.2

−1

−1.2 200 400 600 800 1000

Figure 2.3: The vibrational density of states for 3C-SiC, calculated using three different methods:

FHI(black line), AIMPRO(blue line), and SCC-DFTB(red line). For better clearness, the result of FHI has

been drawn to the positive ordinate, while for the other two curves the sign has been chosen opposite.

20The calculation of vibrational spectra is very time consuming, limiting ab initio methods to small

defective regions. Other, more approximative, methods often do not reach the accuracy in the description

of the forces, which is required for the calculation of vibrational frequencies.

2.7. Applications and Test Calculations 29

Table 2.1: Uppermost phonon modes below and lowest phonon modes above the frequency gap (3C-SiC).

The results of the different methods and the experimental values [70] are given. The values are given in

[meV], values in parentheses are in [cm−1].

Method Atoms ω1ω2Gap

SCC-DFTB 64 78.8 (635.5) 105.1 (847.7) 26.3 (212.1)

SCC-DFTB 216 78.7 (634.8) 106.4 (858.2) 27.7 (223.4)

AIMPRO 64 67.3 (542.8) 91.3 (736.4) 23.9 (192.8)

FHI 64 79.3 (639.6) 91.0 (734.0) 11.7 ( 94.4)

Experiment 76.5 (617.0) 91.0 (734.0) 14.5 (117.0)

between the acoustical and optical modes is described even better by the calculation with

the FHI-code[34], if a (2×2×2) Monckhorst-Pack k–point scheme [65] is used. Table 2.1

compares the results for the uppermost phonon modes below and lowest phonon modes

above the frequency gap are compared for the different methods and the experimental

data. A comparison of the states below the gap shows a very good agreement between the

SCC-DFTB-results and these FHI-results. Above all, the frequencies of the two highest

modes below this gap are at the correct position. The modes above the gap, i. e. the

optical modes have too high frequencies. The AIMPRO-results show a different behavior:

the optical modes are described correctly, while the frequencies of the modes below the

gap are too low. By a multiplication with a factor ≈0.9, the SCC-DFTB-frequencies can

nearly exactly be transformed into the AIMPRO-frequencies. A correction of the size of

the gap towards the experimental value can neither in SCC-DFTB nor in AIMPRO be

achieved by a simple scaling factor.

For 4H-SiC, the vibrational frequencies below the frequency gap are described as well

Wavenumber [1/cm]

4H−SiC

FHI

128 atoms

DFTB

−0.6

−0.4

−0.2

0.2

0.4

0.6

0.8

−0.8

1.2

−1.2 0 200 400 600 800 1000

density of states, arbitrary units

−1

Figure 2.4: The vibrational density of states for 4H-SiC, calculated in two different methods: FHI(black

line) and SCC-DFTB(red line). For better clearness, the result of FHI, has, again, been drawn to the

positive ordinate, while the SCC-DFTB result is mirrored to the negative ordinate.

30 CHAPTER 2. Theoretical Description of Defect Dynamics

as for 3C-SiC in SCC-DFTB compared to the FHI calculation, while, again, the higher

frequencies are overestimated, compare Fig. 2.4. For these calculations, a 4H-SiC super-

cell containing 128 atoms was used in both methods. In this case, also for the ab initio

calculation only the Γ-point was used. It can, therefore, not be attributed to the use of

the Γ-point approximation that the frequency gap is too large in the SCC-DFTB- and

AIMPRO-calculations. A reason for the better result of the FHI-calculation may be found

in the use of a plane-wave basis in FHI, in contrast to localized basis sets used in SCC-

DFTB and AIMPRO. The collective character of the vibrational frequencies might be

better described within this extended basis.

The qualitative agreement is sufficient for the purpose of this work, where only integrated

quantities are needed21. How these are affected by the deviations between the different

methods is subject to the following two sections.

2.7.2 Heat capacity of diamond, silicon, and SiC

The heat capacity Cvat constant volume can be calculated from Eq. 2.28 as the derivative

of Uwith respect to T:

Cv=∂Uvib

∂T V

i=1 





kB



ωi

kBTexp 



ωi

kBT

(exp 



ωi

kBT−1)2





(2.32)

For high temperatures, Cvapproaches 3NkB, in agreement with the rule of Dulong and

Petit22.

Silicon

3C−SiC

Diamond

Temperature T [K]

0 200 400 600 800 1000 1200 1400 1600 1800 2000

theoretical limit

Cv [J/(mol K)]

Figure 2.5: Temperature dependence of Cvfor diamond, silicon, and 3C-SiC. The dashed line denotes the

theoretical limit of 3NkB, which is reached by any solid for temperatures above the material dependent

Debye temperature ΘD.

21Further calculations have shown a correct description of localized modes due to defects in the lattice.

A nearly perfect agreement between AIMPRO and SCC-DFTB has e. g. been found for the modes induced

by the antisite pair, see Chapter 4.

22For very low temperatures, i. e. T < 0.1·ΘD, the description becomes wrong, and the Debye model

[71] has to be used. In this work, however, only the high temperature range is of interest

2.7. Applications and Test Calculations 31

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

0 200 400

0600

specific heat capacity C [J/(g K)]

800 1000 1200 1400 1600 1800 2000

Silicon

3C−SiC

Diamond

}exp. values

Temperature T [K]

Figure 2.6: Calculated specific heat capacity for diamond, silicon, and 3C-SiC compared to experimental

values at T=300 K.

In Fig. 2.5, Cvis plotted for diamond, silicon, and 3C-SiC in a temperature range between

zero and 2000 K. The shape and the limit of the curves and the differences between the

materials, as e. g. the Debye temperature ΘD(Si: 741 K (exp: 650 K), C: 2021 K (exp:

2230 K), SiC: 1456 K (exp: 1600 K)), are described correctly.

Calculating now the (mass dependent) specific heat c, we obtain c=0.450 J/(g·K) for di-

amond, c=0.699 J/(g·K) for silicon, and 0.563 J/(g·K) for 3C-SiC, which agrees well with

the experimental values of 0.502, 0.703, and 0.678 J/(g·K) as given in Ref. [72] for T=300

K. The whole curves are shown in Fig. 2.6.

Experimental results and the temperature dependence as predicted by the Einstein model

is described qualitatively correct by the vibrational spectrum calculated numerically within

SCC-DFTB.

2.7.3 Absolute entropy in different methods

Using the frequencies obtained by different computational methods as described above,

the vibrational entropy and the internal energy of vibration can be calculated. Fig. 2.7

shows a comparison of the temperature dependence of the absolute entropy, scaled by the

number of atoms in the supercell used. Since the vibrational spectrum calculated with

the FHI-code reproduces the experimental properties best (compare last section), it shall

be used as a reference curve. The slope of the AIMPRO curve is then slightly too large,

while that of the SCC-DFTB-curve, calculated in the 64-atom supercell, is by about the

same amount too flat. The SCC-DFTB-calculation in the 216-atom supercell, however,

reproduces the FHI-curve nearly exactly. The internal energy Udoes neither show strong

deviations between the different methods, except for the very low temperature range, which

is not described correctly in the underlying theoretical model23 and, furthermore, not of

interest in the applications.

23The description of the low temperature behavior of the thermal energy Uis described correctly by the

Debye Model, which in contrast to the Einstein model yields the experimentally observed Cv∼T3behavior

at low temperatures [71].

32 CHAPTER 2. Theoretical Description of Defect Dynamics

3C−SiC

AIMPRO

DFTB 64

0.1

0.2

0.3

0.4

0.5

00 200 400 600 800 1000

FHI

1200

Temperature T [K]

1400 1600 1800 2000

0.7

0.6

0.5

0.4

0.3

0.2

0.1

DFTB 216

S [meV]

U [eV]

Figure 2.7: Temperature dependence of the absolute entropy S(left ordinate) and the internal energy U

(rigt ordinate), scaled by the number of atoms in the used supercell.

Since the errors in the description of frequencies concern any structure calculated with one

method in a similar way, the variations of the quantity of interest between the different

methods, i. e. the difference in entropy ∆Sbetween two structures, can be expected to be

even smaller.

2.7.4 Formation entropy

The calculation of formation entropies of defects is by far more complicated and not as

straightforward as the calculation of formation energies. Calculating the entropy as de-

scribed before, a strong dependence on the cluster and supercell size is found. The reasons

for this dependence and how to correct for it is discussed considering the example of an

isolated vacancy in diamond and in silicon.

The periodic images of a defect strongly influence the formation entropy which is an un-

avoidable artifact due to the supercell approach, compare Section 2.4. This is certainly also

true for formation energies, but the effect is by far not as large as for the long-range char-

acter of formation entropies24. It has been found that, in order to simulate the behavior

of a defect in a nearly infinite crystal, modeling is done best by using an embedded cluster

scheme, i. e. dividing the fully relaxed supercell with the defect into a region around the

defect, the ”cluster”, in which all atoms are treated dynamically, and an outer region in

which the atoms are treated as static [73, 74, 75], compare Fig. 2.8.

Calculations with several different methods for the modeling of the vibrational behavior

have shown that the cluster size must be essentially smaller than the size of the supercell.

On the other hand the cluster size must be large enough to provide a correct description

of the defect. In Ref. [73], the authors used cells of about 5000 atoms to model an iso-

lated vacancy in copper, whereof only about 100 to 500 atoms were then included in the

calculation of the vibrational spectra. In Ref. [75], cells of up to 16384 atoms were used,

24Because of the change in the collective vibrational properties the defect causes the entropy is to a large

deal stored in the material outside the defect region.

2.7. Applications and Test Calculations 33

treating up to 1289 atoms dynamically25. Such large supercells are, of course, far beyond

the supercell sizes that can be treated within atomistic calculations, especially since the

diagonalization of the dynamical matrix does not only take a very long time but is also

associated with numerical problems.





supercell

core

cluster

Figure 2.8: The supercell used for modeling

the defect is divided into a spherical cluster

which contains the defect core and is treated

dynamically within the atomistic calculation

(dark violet), and an outer sphere (light blue)

which is treated by an continuum theoretical

approach.

Nevertheless, with the corrections proposed by

these authors, entropy calculations can also be per-

formed in our smaller supercells containing 216

atoms. The comparison with the calculation in a

supercell with 512 atoms for the vacancy in silicon

shows that a consistent picture can be obtained.

To calculate the formation entropy of a va-

cancy, the same calculations have to be per-

formed for the supercell with the vacancy and

for the perfect supercell using the same clus-

ter size. To correct for the different number of

atoms, the results for the perfect lattice have

to be multiplied by (N−1)/N, where Nis

the number of atoms in the perfect lattice clus-

ter.

The observed 1/N dependence [74] in the conver-

gence behavior of the so calculated entropy Sform

can be understood from elasticity theory by writ-

ing the total entropy as S=Score +Selastic. So far, only the first part, Score, has been

considered by the calculation described afore. Due to its long-range character, however, a

large part of the entropy, Selastic, is stored in the region outside the cluster. This elastic

part has, therefore, to be added to the values obtained by the atomistic calculation. It

can be divided into two ”correction terms” that can be calculated in elasticity theory. We

consider a spherical isotropic continuum of radius R. For a point defect with spherical

symmetry, as can be assumed for the vacancy, only radial strain occurs, and so we can get

with Eq. C.11 (see Appendix C) in spherical coordinates

∇2δ= 0 with δ=δ(r) (2.33)

for the dilatation δ. Since δ=∇·u, the displacement field has to be of the form

u(r) = A

r3+B·r(2.34)

with constants Aand Bto be determined by the boundary conditions.

The first (smaller) correction to the entropy can then be understood immediately: The

components of the strain tensor εij are defined as

εij =1

2∂ui

∂xj

+∂uj

∂xi,(2.35)

thus showing a 1/R3dependence for the u(r) of Eq. 2.34. As derived in detail in Ref. [76],

the R-dependence of the elastic constants transmits to the entropy, such that the entropy

25In these calculations Born-Mayer- and Morse-potentials were used for the entropy calculations.

34 CHAPTER 2. Theoretical Description of Defect Dynamics

form kB

[ ]

form

corrected

uncorrected

S][ B

216 atoms

Diamond

512, corrected

216, corrected

216, uncorrected

512, uncorrected

512 atoms

216 atoms

Silicon

−1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

−2

−1

−2

N/Nz

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

N/Nz

0.8 0.9 1

Figure 2.9: Formation entropy of the vacancy in diamond (upper diagram) calculated in a supercell

with NZ=216 atoms and in silicon (lower diagram) calculated in a supercell with NZ=216 atoms (yellow:

uncorrected, orange: corrected, red: linear fit) and with NZ=512 atoms (light blue: uncorrected, blue:

corrected, dark blue: linear fit), restricting the vibrations to Natoms in the shells surrounding the vacancy.

stored outside a sphere with radius Ris proportional to 1/R3. Since the number Nof

atoms is proportional to R3, the correction

∆S1=const.

N(2.36)

is obtained for the formation entropy [73] with a constant depending on the material.

Following the procedure proposed by Hatcher et al. [73], the entropy values obtained like

2.7. Applications and Test Calculations 35

this, require a further correction ∆S2that accounts for the image forces introduced by the

periodic images of the defect and the entropy stored in the up to now neglected region of

the crystal outside the defect cluster. The size of this correction can, again, be calculated

from the free energy [76], resulting in

∆S2=K·α·∆Vimage (2.37)

with the bulk modulus K, the thermal expansion coefficient αand the volume change

∆Vimage due to image forces. This term that can be split into two parts shall be motivated

in the following. The idea for this finite size correction comes again from continuum theory:

at the boundary of the supercell – for simplicity we regard the boundary of an outer sphere

of radius Raand assume isotropic dilatation – , the inevitably vanishing displacements of

the atoms (due to the boundary conditions) introduce a volume change

∆V=−∆Vtot Ri

Ra3

≈ −∆Vtot N

NZ(2.38)

at the inner sphere of radius Ri(which contains the dynamically treated cluster in our

calculations)[73, 76]. ∆Vtot is the relaxation volume of the defect in an infinitely extended

crystal. In reality however, we always have a finite crystal, terminated by the crystal sur-

faces. In supercell calculations, the simulated region is limited by the supercell boundaries.

The volume change is influenced by these constraints, and we will calculate what part of

the volume change ∆Vtot is caused by image forces. At the surface of the sphere (r=R)

representing the crystal no stress is left:

σrr(r=R) = 0 .(2.39)

Expressions for the strain εrr and dilatation δfollow from Eq. 2.34:

εrr =−2A

r3+B, δ =∇·u= 3B . (2.40)

With Eq. C.9 from Appendix C the constant Bcan be expressed in terms of A:

σrr(r=R) = 0 (2.41)

⇔2µεrr +λδ = 0

⇔B=4µ

2µ+ 3λ·A

R3.

The displacement field ubecomes then

u(r) = A1

r3+4µ

2µ+ 3λ·1

R3·r(2.42)

=us(r) + ud(r)

with a shear term usand a dilatation term ud[77].

Let the displacement at the radius rcore of the defect core be u(rcore) = ηrcore, then the

remaining constant Acan be written as

A=1

1 + 4µ r3

core

(2µ+3λ)R3·η r3

core ≈

rcore

R1ηr3

core .(2.43)

Now it can be seen that the dilatation

δ= 3B=12µη r3

core

(2µ+ 3λ)R3(2.44)

36 CHAPTER 2. Theoretical Description of Defect Dynamics

tends to zero for an infinite crystal, but this is not true for a finite crystal26. At the surface

of the outer sphere (r=R), the shear part of Eq. 2.42 causes the volume change

δVs= 4πR2·us(R) = 4πηr3

core ,(2.45)

i. e. the volume change the defect would cause in an infinitely extended crystal. The

dilatation term yields

δVd= 4πR2·ud(R) = 4π4µ

2µ+ 3ληr3

core .(2.46)

Using instead of the Lam´e constants µand λthe poisson ratio

ν=λ

2(λ+µ),(2.47)

we get

δVd= 4π2(1 −2ν)

1 + νηr3

core ,(2.48)

for the dilatation part and thus the total volume change becomes

δVtot = 4π3(1 −ν)

1 + νηr3

core .(2.49)

The dilatation part arises only due to image forces at the surface of the sphere. From

Eqns. 2.48 and 2.49 it follows that

δVd

δVtot

31−2ν

1−ν,(2.50)

setting the volume change due to image forces which is needed for the entropy correction

in relation to the total volume change [77].













VVK

’

supercell K

Ra/2

Figure 2.10: Sphere inscribed in a supercell

Up to now we used a sphere of radius Rin our considerations, but the atomistic calculations

are performed in cubic supercells with the same volume and edge length a. The postulation

of isotropy may not apply to the marginal regions of the supercell due to the influence of

neighboring cells. We, therefore, have to go over to a sphere with radius a/2 that is

inscribed in the supercell, see Fig. 2.10. Since the relation of these volumina is

=V0

Vsupercell

4π

3a

23

a3=π

6,(2.51)

26and especially not for the from the macroscopic view small supercells

2.7. Applications and Test Calculations 37

Table 2.2: Data used for the calculation of the entropy of formation of an isolated vacancy in diamond

and silicon. Experimental values from Ref. [78]

Diamond Silicon

B 536 GPa 100 GPa

ν0.1 0.29

α7.1·10−6K−12.6·10−6K−1

NZ216 216; 512

ncore 17 17

∆Vcore −1.330 ·10−30 m3−2.137 ·10−29 m3

Vtot 1.226 ·10−27 m30.431 resp. 1.024 ·10−26 m3

the smaller volume δVd·π

6has to be used in Eq. 2.37. Expressing the total volume change

by the volume change of the defect core

∆Vtot

∆Vcore ≈R3

core ≈NZ

ncore

,(2.52)

the final expression for the correction of the formation entropy becomes

∆S=const.

N+Kα ·2

(1 −2ν)

(1 −ν)

6−N

NZ·NZ

ncore ·∆Vcore .(2.53)

Due to the numerous assumptions and approximations, i. e. spherical symmetry of the de-

fect core, the cluster and the complete supercell in order to use the macroscopic continuum

theory for a furthermore rather small number of atoms as can be treated in our supercell

calculations, results can only be expected to be of qualitative accuracy. Nevertheless, the

results presented in the next section for isolated vacancies in diamond and silicon are quite

encouraging.

2.7.5 The vacancy in diamond and silicon

The procedure described above has been applied to the calculation of the formation en-

tropy of isolated vacancies in diamond and silicon, using the data given in Table 2.2. A

temperature of T=2000 K has been chosen, where the high temperature approximation is

valid and the formation entropy becomes temperature independent.

The results of these calculations are shown in Fig. 2.9. In the upper diagram, the values

of Sform of the vacancy in diamond obtained from the calculation without any corrections

are plotted over the relative cluster size N/NZ(black diamonds), the corrected values are

plotted in red. The straight red line has been fitted to these values, so that the desired

value for N=NZcan be extrapolated, resulting in Sform = 2.85kB. Variations of the cor-

rected values are less than ±0.3kB. For the isolated vacancy in diamond, the corrections

are, thus, very small.

The description of the vacancy in silicon turns out to be much more complicated. The di-

verging behavior of the uncorrected values for the entropy in the supercell with 216 atoms

38 CHAPTER 2. Theoretical Description of Defect Dynamics

216 atoms

Diamond

1and

(Nz=512)

∆S1

∆S2

∆S

∆

(Nz=216)

∆S1

∆S

(Nz=512)

(Nz=216)

and

1∆S2

∆S

216 atoms

512 atoms

Silicon

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

N/Nz

−1

−2

[kB]

−2

[kB]

−1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

N/Nz

Figure 2.11: Corrections to the formation entropy for the vacancy in diamond (upper diagram) and silicon

(bottom).

is found to be even worse in the larger supercell with 512 atoms, – an observation made

by e. g. Fern´andez et al., as well [74]. Taking a closer look at the corrections derived in

the previous section, this observation can be understood, because the second term of the

correction is proportional to ∆Vcore, which is by about one order of magnitude larger for

the vacancy in silicon than for the vacancy in diamond (see Table 2.2).

In Fig. 2.11, the two corrections to the formation entropy are plotted separately for the

vacancy in diamond (upper diagram) and in silicon (lower diagram). While in case of

diamond only a very small correction is obtained, the strong relaxation yields large image

forces and, accordingly, large corrections for the vacancy in silicon. The first correction

∆S1is only important for small cluster sizes N, where the largest part of the entropy is

stored in the region outside the cluster. Additionally, it is compensated in part by the

second correction ∆S2. This part of the correction is rather small for diamond, as well,

but it is of the same order of magnitude as the core entropy for the vacancy in silicon.

The total formation entropy of the vacancy in silicon, again obtained by extrapolation of

2.7. Applications and Test Calculations 39

the corrected values to N=NZ, results in 3.32 kB(216 atoms) and 3.33 kB(512 atoms),

compare Fig. 2.9. Variations are slightly larger as in case of diamond, but still in the range

of ±0.5kB.

The values obtained with the described corrections are of a proper order of magnitude, as

a comparison with the literature shows. In Refs. [43, 79], the entropy of formation (in the

high temperature approximation) is calculated for the unrelaxed vacancy in silicon, based

on a force constant model and including the nearest or next nearest neighbors of the va-

cancy, only. In Ref. [79] a Green’s function technique is used to investigate the influence of

the force constant model on the formation entropy. Although in most simple force constant

models a value of ≈1.7 kBcan be derived for an unrelaxed vacancy in any tetrahedrally

bonded material, the results obtained with different methods vary between 1.5 kBand up

to 8 kBand are highly sensitive to the force constants used. For the relaxed vacancy a

value of ≈3kBas is also obtained by first principle calculations [80], is commonly accepted.

Small local changes in the force constants can result in large changes in the formation en-

tropy, making a good description of the forces desirable.

The reason for the discrepancies between the formation entropies calculated within differ-

ent methods is to find in the fact that the vacancy in silicon is one of the most complicated

defects with respect to a converged energetical and electronic structure, compare also the

discussion in Chapter 3. Even for the formation energies convergence starts at supercell

sizes larger than 216 atoms, as has been extensively investigated in Ref. [61, 64]. Con-

vergence problems have to be ascribed to the extremely flat potential energy surface of

silicon, which in a high degree influences the quality of the dynamical matrix and thus the

vibrational frequencies.

Concluding this section, we have found a possibility to calculate the vibrational entropy of

defects. The absolute values calculated within SCC-DFTB are in rather good agreement

with ab initio results. Including certain corrections, we can obtain reasonable values for the

formation entropies of point defects, as shown at the example of isolated vacancies in silicon

and diamond27. The migration entropies discussed in Chapters 4 and 5 do not suffer from

the problems discussed in this section and can, therefore, be expected to be less sensitive

in their description. Firstly, forces in SiC can be calculated much more accurate than in

silicon, and, second, errors due to artificial constraints by the supercell boundaries or by

periodic images can be expected to be of similar magnitude in the compared structures,

whereas the comparison of a defective and a perfect supercell as required for the calculation

of formation entropies is something completely different. The most important quantity to

influence the calculated entropies is the accuracy of relaxation of the atoms surrounding

the defect, compare Ref. [81].

27The equations given can also be applied to other defect structures, as e. g. interstitials, with minor

changes. The treatment of dopant atoms requires further consideration. This is, however, not subject of

this work.

40 CHAPTER 2. Theoretical Description of Defect Dynamics

Chapter 3

Vacancies and Interstitials

The first section of this chapter gives a brief overview of our results for the formation

energies of silicon- and carbon–vacancies as well as next neighbor pairs of them. The

SCC-DFTB results are compared to literature data, in order to get an impression of the

expected accuracy of the energies calculated in various methods. Afterwards, the sublattice

migration of both vacancies is discussed, since it is not only of interest as a ”pure” process

but also plays a key role in the formation of aggregates of other defects, as is described

in Chapters 4 and 5. As a competing process to the sublattice migration, the formation

of vacancy–antisite pairs is discussed, reviewing the most important findings about the

VCCSi pair, which in turn appears in the discussion of the creation of antisite pairs (Chap-

ter 4) and the recovery of free charge carriers in nitrogen doped material (Chapter 5).

The second section of the chapter is dedicated to carbon split interstitials, and in the third

section, finally, the role silicon- and carbon vacancies play in the migration processes of

these interstitials is briefly described, clarifying the long-range effects defects can have on

activation energies in SiC.

The results of this chapter do not claim to be exhaustive, but only those matters that are

of interest for the following two chapters are pointed out, here.

3.1 Vacancies

As unavoidably in a large number created defects, vacancies are most important in the

situation directly after ion implantation: Vacancy assisted processes have been discussed

in Section 2.1.3, and in Chapter 4 they will turn out to be essential for the mobility of

antisites. The ion bombardement causes monovacancies on both sublattices as well as

pairs of a silicon- and a carbon vacancy, i. e. divacancies, and possibly also larger vacancy

aggregates.

3.1.1 Formation energies

A comparison of the formation energies of vacancies and divacancies in 3C- and 4H-SiC

is given in the following tables. Supercells of different sizes and the dependence of the

calculated formation energies on the number of atoms around the defect that were allowed

to relax during structure optimization have been compared in order to get an impression

of how our method (SCC-DFTB) behaves under varying outer conditions and under which

conditions results come closest to the available ab initio results. These literature values

42 CHAPTER 3. Vacancies and Interstitials

have been obtained with the Finnish plane wave based LDA-code ”FINGER” using a cubic

supercell with 128 atoms, which all were allowed to relax freely[11]. The results for 3C-SiC

are shown in Table 3.1.

The same calculations have been performed for the 4H-polytype of SiC, but here it has to

be distinguished between the two different sites – (quasi-)cubic or hexagonal – for each of

the vacancies, compare Fig. 1.4. All possible orientations have been calculated, and the

formation energies are shown in Table 3.2. The last four lines of Table 3.2 show the binding

energies of the divacancies on the different sites.

In a binary semiconductor, the formation energy of vacancies naturally depends on the

Table 3.1: Formation energies of vacancies in 3C-SiC, calculated within SCC-DFTB. The last column

shows literature values from Ref. [11]. Especially for the carbon vacancy and for the divacancy, the formation

energies vary rather strongly with the size of the supercell and the constraints in the calculation. In the

last line of the table the binding energy of the divacancy VCVSi , obtained by subtracting the third row

from the sum of the first two rows, is shown. In the columns denoted with ”2NN”, relaxation was limited

to nearest and next nearest neighbors of the defect. All values in [eV].

Defect 128 cell, 2NN 128 cell, all 216 cell, 2NN 216 cell, all 128 cell, all [11]

VSi 8.32 8.31 8.67 8.66 7.79

VC4.70 4.70 4.39 4.12 2.77

VCVSi 10.44 10.38 8.95 8.85 7.22

VCVSi 2.58 2.62 4.12 3.93 3.34

Table 3.2: Formation energies of vacancies in 4H-SiC. The differences between the hexagonal and cubic

sites are smaller than the variations due to different constraints in the calculation. These are smaller than

for 3C-SiC. The last four lines show, again, the binding energies of the divacancy in all possible orientations.

Sites 128 cell, 2NN 128 cell, all 240 cell, 2NN 240 cell, all Ref. [11]

VSi cub 8.81 8.78 8.57 8.54 8.37

VSi hex 8.82 8.79 8.58 8.55 8.26

VCcub 4.62 4.39 4.32 4.05 4.07

VChex 4.73 4.49 4.51 4.16 4.21

VCcubVSi cub 9.01 8.93 8.85 8.77 7.74

VCcubVSi hex 9.13 9.07 8.88 8.80 8.36

VChexVSi cub 9.12 9.03 9.01 8.90 8.34

VChexVSi hex 9.04 8.97 8.89 8.81 8.00

VCcubVSi cub 4.42 4.23 4.03 3.83 4.36

VCcubVSi hex 4.30 4.11 4.02 3.80 3.77

VChexVSi cub 4.42 4.24 4.07 3.80 3.90

VChexVSi hex 4.51 4.31 4.20 3.90 4.27

3.1. Vacancies 43

chemical potential of the constituents of the material, since the stoichiometry is changed

by the creation of a vacancy. To be able to compare our results with the available literature

data and not to destroy the clarity of the description needlessly, formation energies are only

given for one selected set of chemical potentials µSi and µC. Throughout this work, expres-

sion A.11 with ∆µ= 0 has been used to calculate formation energies (see Appendix A).

Since ∆µ= 0 means per definition (Eq. A.10) µSi =µC, this choice of ∆µ= 0 describes

stoichiometric conditions, which are also experimentally most common.

A comparison of the SCC-DFTB-results in Tables 3.1 and 3.2 for the second next neighbor

relaxed and fully relaxed structures shows for both supercell sizes only small variations, in

most cases ≈0.1 eV.

For defects that induce a stronger and/or more long-range lattice distortion in their neigh-

borhood, partly substantially larger differences can occur due to constrained relaxation,

compare the analogous discussion in Chapter 4 for antisites and antisite pairs. This must

be considered for the choice of the supercell so that the number of atoms that are allowed

to relax can be chosen large enough but artifacts due to periodic images of the defect are

kept small at the same time.

The energy differences between the calculations in the 128-atom supercell on the one hand

and the 240-atom supercell on the other are slightly larger, in most cases ≈0.2 eV. The

largest difference was found for the divacancy in 3C-SiC, where the 128-atom supercell is

obviously to small to describe the relaxation around the defect correctly. The fact that this

is not the case in 4H-SiC can be explained by the form of the 128-atom supercell, which has

a fcc-basis in the 3C- but a sc-basis in the 4H-polytype, leading to different defect-defect

distances across supercell boundaries.

The energy differences between the SCC-DFTB-results and the ab initio results are a little

larger than those between the different supercells within SCC-DFTB. This can have several

reasons. In principle, more accurate values should, of course, be expected from ab initio

calculations due to the less approximative description of the electronic structure in these

methods.

On the other hand, most ab initio calculations suffer from being limited to small super-

cells which can both lower or increase formation energies and migration barriers due to

artificial interactions between the periodic images of the defect. An example where large

variations in the formation energies occur if using different supercells and k-point sampling

schemes is the isolated vacancy in silicon1. In Ref. [64], Puska et al. present systematic

investigations on this matter, reporting a very slow and oscillating convergence behavior

of the formation energies and the geometrical relaxation of VSi in all charge states with

the supercell size. The principle behavior of the formation energy obtained by these plane

wave LDA-calculations could be confirmed by our DFTB-calculations for supercells of up

to 512 atoms [61, 82]. Convergence is, however, reached for much smaller supercells in our

calculations. This observation might be ascribed to a larger influence of periodic images

on the widely extended plane wave basis, while the tight-binding description obviously

suppresses such artificial effects better. This is also suggested by the induced relaxation

around the vacancy, which, in contrast to the plane wave calculations, does not change

1As already pointed out during the discussion of the calculation of formation entropies in the previous

chapter, silicon is due to its very flat potential energy surface an extremely critical case. The effects studied

at its example can, therefore, be expected to be considerably smaller in a material that is less critical in

this respect.

44 CHAPTER 3. Vacancies and Interstitials

substantially for supercells larger than 128 atoms in the DFTB-calculations [61, 82].

Furthermore, even different ab initio methods show variations of the same order of mag-

nitude. The low formation energy calculated for the carbon vacancy in Ref. [11] (with the

plane wave code FINGER) is e. g. not obtained with a similar plane wave code (FHI) [12],

the results of which (≈4 eV) are clearly closer to the SCC-DFTB results. A reason for

these deviations lies certainly in the free choice of certain parameters in the construction

of the pseudopotentials as well as other partly more technical quantities.

Having this in mind, one should rather set great store by a good consistent picture of the

results within the method concerning the absolute values. The qualitative results, how-

ever, can be compared with those of other methods. These compare quite well within

SCC-DFTB and the available ab initio data as the results for the vacancies in this chapter

as well as the results for antisites in Chapter 4 show. The 216-atom cell (3C) and the

240-atom cell (4H) with a large number of free relaxing atoms turn out to yield converged

results and come in most cases closest to the literature value (where also all atoms were

relaxed). These supercells are used in the following if nothing else is explicitly said.

Up to now, we have only challenged the accuracy of the absolute values for the formation

energies of vacancies and divacancies calculated with SCC-DFTB compared to results of

ab initio methods. The accuracy of SCC-DFTB has been discussed, while that of the

literature values has been accepted unquestioningly. However, further test calculations

with the Green’s functions based LMTO method have shown that also these first principle

values have to be handled with care, as explained in the following.

As discussed in the previous chapter, the local density approximation (LDA) in density

functional theory (DFT) has the disadvantage of an incorrect description of the energy

gap of semiconductors. Furthermore, it has been mentioned in Section 2.2 that there are

several techniques to account for this failure and correct for the size of the energy gap

and thereby also the position of localized levels in the gap. The energy gap obtained for

a material within LDA depends on the basis set used. The ab initio methods used for the

reference calculations, i. e. the FHI code and the FINGER code, use a plane wave basis

with s-, p-, and d-basis functions for the creation of the pseudopotentials. In SCC-DFTB,

a minimal basis of s- and p-orbitals is used, instead.

Since it is not possible without very large expenses in the plane wave codes or in SCC-

DFTB to change the basis set, the LMTO results shall illustrate the dependence of the gap

as obtained in pure LDA, i. e. without any corrections, on the basis set. For cubic SiC,

a minimal basis (sp) yields an energy gap of 2.91 eV, including the d-orbitals leads to a

considerably smaller gap of 1.29 eV, and further including the f-orbitals results in a gap of

1.63 eV. These values have to be compared to the experimental gap of 2.4 eV for 3C-SiC.

In SCC-DFTB, the minimal basis (sp) leads to a large band gap of 5.9 eV, thus clearly

further from the experimental value2. In spite of the obviously incorrect description of the

conduction band in SCC-DFTB with the sp-basis, it has also advantages for the description

of formation energies if the gap is too large, as will become clear in the following.

The effects of a Baraff-Schl¨uter correction [55], which is applied to the finished supercell

calculation, and a fully selfconsistent correction calculated in the Green’s function based

LMTO with the scissor operator technique [54] has been investigated at the example of

the silicon vacancy and the carbon vacancy in 3C-SiC.

In the Baraff-Schl¨uter correction, the calculated gap is scaled to the experimental gap,

and the localized defect levels in the gap are also shifted by this scaling factor, weighted

2This strong overestimation has not only to be ascribed to the minimal basis, but is typical for tight-

binding methods.

3.1. Vacancies 45

Table 3.3: The effect of two different corrections for the LDA-gap on the formation energies of VC,VSi ,

and VCVSi . ”B-S” shows the formation energy with the Baraff-Schl¨uter correction, ”Scissor” includes the

fully selfconsistent correction with the scissor operator, instead. For comparison, the last column shows

the SCC-DFTB results for the 216 atom supercell, allowing all atoms to relax, freely. Below the line the

binding energy is shown. All values in [eV].

Defect Eform [11] B-S (FHI) Scissor (LMTO) SCC-DFTB

VSi 7.79 8.26 8.40 8.66

VC2.77 3.57 3.78 4.12

VCVSi 7.22 – 8.47 8.85

VCVSi 3.34 – 3.71 3.93

with their overlap with the conduction band. In contrast to this correction, the scissor

operator allows to correct not only the localized gap levels, but the whole valence band

selfconsistently. In contrast to the Baraff-Schl¨uter correction of the gap levels only, resonant

levels in the valence band are accounted for by this method. Due to the empirical character

of this kind of corrections, most often no corrections are made in ab initio calculations. Only

rarely, the Baraff-Schl¨uter correction is used in ab initio methods, while the selfconsistency

and the neglected correction of the valence band is hoped not to have a significant influence,

which is, in fact, true in many cases3.

All literature values in Tables 3.1 and 3.2 are LDA-results without any corrections. Ta-

ble 3.3 shows the influence of the two different corrections on the formation energies of

VC, VSi and VCVSi . The third column shows the literature values from Ref. [11] plus

the Baraff-Schl¨uter correction calculated within the FHI code. Due to the overlap of the

occupied gap-levels with the conduction band, a correction of 0.47 eV has to be added to

the formation energy of the silicon vacancy. For the carbon vacancy, a larger value of 0.80

eV has been calculated. A similar calculation within the LMTO method, but correcting

the gap levels selfconsistently, yields very similar values. Adding the fully selfconsistent

correction of the gap levels and the whole valence band to the results of Ref. [11] results

in values even closer to the SCC-DFTB results, compare the fourth column. Especially

for the silicon vacancy, formation energies have become very similar. A larger difference

remains for the formation energy of the carbon vacancy, and consequently for the diva-

cancy. An FHI calculation yielded 7.75 eV for the formation energy of the silicon vacancy,

and 4.35 eV for the carbon vacancy (uncorrected values) [68]. Adding the Baraff-Schl¨uter

corrections to these results leads to 8.2 eV and 5.1 eV for VSi and VC, the full correction

results in 8.36 eV and 5.36 eV.

3.1.2 Migration of vacancies

In this section, the description of non-equilibrium structures is investigated at the example

of activation energies for the sublattice migration of isolated vacancies. In the binary lattice

of a compound semiconductor like SiC, a vacancy can migrate through the sublattice by

a second neighbor atom moving onto the vacant site. This sublattice migration process

(compare Section 2.2.3) is also the basis of many other migration processes: In Chapter

3In case of the saddle point structures for sublattice migration of VCand VSi the importance of the

correction of the whole valence band will become clear, see the next section.

46 CHAPTER 3. Vacancies and Interstitials

Table 3.4: Comparison of calculated activation energies (i. e. the energy difference between minimum and

saddle point geometry) for vacancies in 3C-SiC and 4H-SiC (only last column). For comparison, values

obtained by SCC-DFTB for the cubic (c) sites in 4H-SiC are also given. Besides the polytype, the number

of atoms in the supercell and the number of relaxed atoms is shown. All values in [eV].

Eactive Ref. [12, 13] Ref. [68] DFTB DFTB DFTB DFTB DFTB, 4H

216, 2NN 128, 2NN 128, 2NN 128, full 216, full 512, full 240, full (c)

VC2+ 5.2 6.4 5.8 5.8

VC3.5 4.0 5.1 4.6 4.8 4.8 4.7

VC2−4.4 3.9 4.1

VSi 2+ 5.3 4.3 4.4

VSi 3.4 3.6 4.9 3.9 3.9 4.1 4.1

VSi 2−2.4 4.3 3.5 3.6

4, it is shown that the mobility of vacancies is prerequisite for the mobility of isolated

antisites, and also for the formation of nitrogen complexes the mobility of vacancies is of

importance, see Chapter 5.

Unfortunately, the number of available ab initio calculations for comparison is very re-

stricted, which on the other hand suggests to use a more approximative method like SCC-

DFTB for a systematic series of investigations4.

The description of migration processes shows the same behavior as that of the formation

energies: in most cases, allowing more atoms to relax lowers the activation energies as

well as using larger supercells5. Though the SCC-DFTB-energies are throughout higher

than the available plane wave results, the same trend can be observed in both methods

for the different charge states: activation energies increase the more positively the defect

is charged.

A comparison of the SCC-DFTB-results of Tables 3.1- 3.4 draws a consistent picture of

both formation- and activation energies. The results obtained in the supercell containing

512 atoms confirms that the results in the smaller cell (216 atoms) are already converged.

The deviation of the activation energies of the sublattice migration of vacancies compared

to the plane wave results sets limits to the expectable accuracy of SCC-DFTB in the first

line, but as well of other, also ab initio methods. The attempt to correlate the calculated

activation energies with experimentally more direct accessible annealing temperatures (see

Chapter 5) assumes, however, first of all the consistency of the calculated data rather than

an absolute accuracy of less than 0.5 eV. Furthermore, the same considerations as in the

previous section suggest that the LDA results of Ref. [12, 13] and [68] require a correction

due to the underestimated band gap.

Including the Baraff-Schl¨uter corrections for the silicon vacancy and the saddle point struc-

4For more accurate values of the most important processes found with SCC-DFTB, an ab initio calcu-

lation can then follow. As pointed out in the previous chapter, the calculation of migration paths is rather

time consuming, for which reason many ab initio calculations are limited to those migrations processes,

where the saddle point geometry can be derived from e. g. symmetry arguments, so that only this structure

can then be relaxed under certain constraints to obtain the activation energy.

5Larger supercells can also result in higher formation energies, as for instance in the example of the

isolated vacancy in silicon [64]. A general statement about the variation of formation energies or migration

energies with the boundary conditions cannot be made.

3.1. Vacancies 47

Table 3.5: The effect of two different corrections for the LDA-gap on the activation energies of sublattice

migration of VCand VSi . ”B-S” shows the activation energy with the Baraff-Schl¨uter correction, ”Scissor”

includes the fully selfconsistent correction with the scissor operator, instead. For comparison, the last

column shows the SCC-DFTB activation energy (216 atom supercell, full relaxation). All values in [eV].

Defect Eactiv [68] B-S (FHI) Scissor (LMTO) SCC-DFTB

VSi 3.6 3.6 4.5 3.9

VC4.0 4.7 6.2 4.8

∆Eactiv 0.4 1.1 1.7 0.9

Defect Eactiv [12, 13] B-S (FHI) Scissor (LMTO) SCC-DFTB

VSi 3.4 3.4 4.3 3.9

VC3.5 4.2 5.7 4.8

∆Eactiv 0.1 0.8 1.4 0.9

ture does not change the results, since in both cases the correction amounts to 0.4 eV [68].

For the carbon vacancy, however, a large difference occurs, since several defect levels are

in the energy gap of the saddle point structure, but only one in case of the vacancy. The

Baraff-Schl¨uter corrections are 0.79 eV for the carbon vacancy and 1.44 eV for the saddle

point, leading to a correction of 0.65 eV for the energy barrier [68]. Adding these correc-

tions to the LDA results listed in the second and the third column yields 4.2 and 3.4 eV

for VCand VSi for the values of Ref. [12, 13] and 4.7 and 3.6 eV for the values of Ref. [68].

Tables 3.5 summarize the corrections applied to both the literature values of Ref. [12, 13]

and Ref. [68]. Including the full correction of the valence band and gap levels leads to

substantial changes, since large contributions to the saddle point energies arise from the

correction of the valence band, which are not covered by the simple Baraff-Schl¨uter correc-

tion, compare the third and fourth columns in Tables 3.5. Comparing now the SCC-DFTB

results in the last column with the corrected LDA values shows that at least part of the

required corrections is implicitely included in the SCC-DFTB energies, which lie between

the partly and the fully corrected LDA values. In contrast to the ab initio LDA calculations

with a considerably too small band gap, also localized levels that are close to the conduc-

tion band are clearly separated, and no correction due to the overlap with the conduction

band is required. Especially for the description of high negative charge states of defects,

less technical problems occur than in LDA calculations.

Regarding the energy differences also between different ab initio calculations with and

without the two corrections, the discussion of energy differences of less than 0.5 eV becomes

obviously meaningless. Considerably more important is the qualitative description, which

expresses itself e. g. in the difference of the activation energies for the sublattice migration

of carbon vacancies and silicon vacancies. The last lines in Tables 3.5 show this energy

difference. While especially the literature values of Ref. [12, 13] suggest an activation of

the sublattice migration of both vacancies at approximately the same temperature (since

the activation energy difference is negligibly small), the corrected LDA activation energies

and the SCC-DFTB activation energies suggest that sublattice migration of the carbon

vacancy should first be activated at several hundred degrees Kelvin higher temperatures6.

6A correlation between activation energies and temperatures will be made in Chapter 5

48 CHAPTER 3. Vacancies and Interstitials

VSi

SiSi

5.9%

−5.3%

−9%

−17%

−1%

2.8% 2.8%

−17%

−5.3%

2.8%

CCSi

Figure 3.1: By the migration of one of its carbon ligands towards the vacant site, the silicon vacancy can

transform into the CSi VCpair. Relaxation induced changes of the bond lengths are given with respect to

the ideal Si-C bond length 1.88 ˚

This finding could only recently be affirmed experimentally by EPR measurements on

annealing studies [21], in which the signal of the carbon vacancy in as-grown samples could

still be detected after that of the silicon vacancy had vanished. An exact value for the

annealing temperature of VC, i. e. the temperature that suffices to activate sublattice

migration, is not yet available.

3.1.3 Vacancy – antisite pair formation

A process that can be activated at much lower cost than the sublattice migration of a

silicon vacancy as described above is the movement of one of its carbon ligands into a

silicon vacancy. For the vacancy and the antisite on cubic lattice sites, the atomistic model

of this process is drawn schematically in Fig. 3.1. The relaxed geometry of the neutral

VCCSi pair has C1hsymmetry and is sketched at the right hand side of the figure, showing

a rather strong contraction of two of the C-C bonds around the carbon antisite, which itself

moves away from the vacancy, while the silicon ligands of the carbon vacancy relax inwards.

The initial and final structures as well as the saddle point geometry in 4H-SiC are also

shown in Fig. 3.2 together with the energy diagram that shows the change in energy during

the process. The resulting CSi VCpair is by 1.8 eV more stable than the silicon vacancy.

The carbon atom does not migrate along the direct connecting line between the vacant

silicon site and its initial site, but is slightly displaced from this line, leading to a slightly

asymmetric saddle point geometry. A comparatively low energy barrier of 1.7 eV has to be

overcome during this process, thus, according to Eq. 2.6, it should be expected to happen

much more frequently than sublattice migration (Eactiv ≈4 eV). At temperatures sufficient

to activate sublattice migration, the recombination barrier of 3.5 eV can be overcome,

such that the VCCSi pairs can transform back to silicon vacancies which then can move to

neighboring sites. With an activation energy of ≈4.6 eV, i. e. nearly the value obtained for

sublattice migration of VC, the dissociation of the carbon vacancy from the VCCSi pair is

clearly less likely than this back-transformation to the silicon vacancy.

Already in 1990, Itoh et al. observed VSi to anneal out at 750◦C [84]. An atomistic model

for the underlying diffusion mechanism was, however, missing. Since sublattice migration

of VSi requires much higher activation energies than the formation of VCCSi pairs which

would also cause a decrease in the signal intensity of VSi , we proposed the process VSi −→

VCCSi as an explanation for the first annealing step of VSi [9], while sublattice migration

should be expected to happen first at higher temperatures (compare the correlation made

3.1. Vacancies 49

1.7 eV

VSi

Saddle point Si

−0.5

0.5

−1.5

−1

1.5

−2

Figure 3.2: By the migration of one of its carbon ligands towards the vacant site, the silicon vacancy can

transform into the CSi VCpair. An activation energy of 1.7 eV is required for this migration process, and

the resulting pair defect is by 1.8 eV lower in energy. The saddle point structure shown above the energy

diagram is asymmetric, with the migrating carbon atom slightly displaced from the direct connecting line.

300 320 340 360 380 400

P7c

N + X

P7a

P7b

P6c

P6b

P6a



EPR signal [arb. units]

magnetic field [mT]

305 306 307

*) 29Si

13C

x 10 x 10



EPR signal [arb. units]

magnetic field [mT]

Figure 3.3: Two of the spectra that served for the identification of the P6/P7 spectrum with the VCCSipair

(taken from Ref. [10]). Left: X-band EPR spectrum of neutron irradiated and annealed 6H-SiC, measured

under illumination. Right: Hyperfine structure of the P6c low field EPR line. The intensity ration of the

hyperfine satellites to the central line corresponds to the interaction with 4-8 silicon nuclei and one carbon

nucleus.

in Chapter 5).

Based on these results, an identification of VCCSi with the so-called P6/P7 EPR spec-

trum, see the left part of Fig. 3.3, could later on be realized by using a combination of two

theoretical methods (SCC-DFTB and LMTO-ASA) with the experimental investigations

of neutron irradiated 6H-SiC by Lingner et al. [10, 83]. In earlier works, these P6/P7

spectra were attributed to nearest neighbor vacancy pairs [85], and in fact, the divacancy

50 CHAPTER 3. Vacancies and Interstitials

VCVSi can, according to our LMTO-ASA calculations, explain most of the experimen-

tal findings concerning the excitation scheme and the hyperfine splittings. The hyperfine

structure measured by Lingner indicates, however, the presence of one prominent carbon

nucleus, see the right diagram in Fig. 3.3. This observation and the rather high activation

energy for sublattice migration of VCand VSi rule out the divacancy as an atomistic model

for the P6/P7 spectra and limits the possible candidates to the antisite pair CSi SiCand

the VCCSi pair.

The mechanism for the creation of antisite pairs is expected to require higher temperatures

than 750◦C, see Chapter 4. Furthermore, the excitation energies of <0.2 eV calculated for

the antisite pair, deviate clearly from the ≈1 eV measured by MCDA (magnectic circu-

lar dichroism of the absorption). These optical transitions, the symmetry, the annealing

behavior and the observed hyperfine splittings can only be explained by the remaining

candidate, the VCCSi pair.

Using the optimized geometry of VCCSi obtained from SCC-DFTB calculations, the hy-

perfine splittings could be calculated within LMTO-ASA. The localization of the electronic

levels which give rise to the hyperfine splittings allowed us to determine the charge state

of the observed pair to be doubly positive7.

A detailed description of the properties of VCCSi can be found in Refs. [10, 83, 86]. The

identification of the proposed pair defect with an experimentally observed spectrum shows

first successful applications of SCC-DFTB to the description of point defects and mech-

anisms of their migration. In addition, the procedure underlines the importance of using

various methods, both experimental and theoretical, to obtain a complete picture of the

properties of defects and their annealing behavior.

The VCCSi pair itself turns out to play an important role in the creation of antisite– and

nitrogen–related complexes, which are subject to Chapters 4 and 5.

The simplicity of the above described mechanism for the silicon vacancy suggests that an

analogous mechanism should exist for the carbon vacancy, which is observed to anneal

out at even lower temperatures (≈200◦C) in irradiated material [14]. Calculations could,

however, show that the SiCVSi pair, which would result from the migration process, is

instable, with the silicon antisite recombining immediately with the silicon vacancy [9].

The process VC−→ SiCVSi does, accordingly not exist, as could also be confirmed by plane

wave calculations8[12].

3.2 Interstitials

During the implantation process, the implanted ions kick out atoms from their lattice sites

and move on as long as their energy suffices to do so, before they are finally built in.

Accordingly, interstitials are created in the same amount as vacancies, so that they are

supposed to play an important role in the annealing processes that are the central subject

of this work.

As carbon interstitials play the most important role in the mechanisms investigated in this

7The geometry of the doubly positive pair deviates less from C3v–symmetry than that of the neutral

pair shown in Fig. 3.1. Details can be found in Ref. [10, 83].

8The activation energies obtained with these calculations for the creation of the VCCSi pair (≈2 eV)

agrees quite well with our 1.7 eV.

3.2. Interstitials 51

work, we limit the description of interstitials to them. Silicon interstitials do not contribute

essentially to the process of especially nitrogen migration (see Chapter 5) which can be

attributed to their stability mainly in higher positive charge states [12], the strong lattice

deformation due to the size of the silicon atoms and thereby higher migration barriers. A

further reason is the attitude of nitrogen atoms to prefer the carbon site (see Chapter 5),

intimating a carbon related interstitial process.

Si −4.8 %

−26.6%

Figure 3.4: Geometry of the (CC)Csplit-

interstitial

The most stable carbon interstitials in SiC are

carbon split-interstitials, where two C–atoms

share one C-site, such that one of them and

its two Si-ligands are in a plane perpendic-

ular to the plane spun up by the other C

and its two ligands, see Fig. 3.4. The de-

fect has D2d-symmetry9with a highly con-

tracted C-C-bond (26.6 % compared to the

ideal SiC bond length) between the intersti-

tial atoms that is with 1.38 ˚

A even about 10

% shorter than the ideal bond length in dia-

mond. The silicon ligands are strongly pushed

outwards, so that for the four Si-C bonds

around the defect a contraction of only 4.8 % is

left.

This (CC)Csplit-interstitial can move through the lattice by the migration of one of its

C-atoms to the next C-site. In 3C-SiC, this migration requires ≈2.9 eV, while in 4H,

the activation energy depends on the direction, see Fig. 3.5. Here, the lowest activation

energy, ≈2.9 eV, is needed for a migration along the [0001]-direction, while a migration

along [10¯

10] is least favorable.

along [10-10]

along [11-20]

along [0001]

0reaction coordinate

Energy [eV]

Figure 3.5: Left: The migration of a (CC)Csplit-interstitial in 4H-SiC to a neighboring C-site along the

crystal directions [11¯

20], [10¯

10] and [0001]. Right: Energies along the migration path. The lowest activation

energy is needed for a migration along the [0001]-direction.

For the (CC)Csplit-interstitial in 3C-SiC, molecular dynamic (MD) simulations were per-

9This remains valid for all charge states less negative than 2-. In this case, the symmetry axis of the

defect is slightly inclined, resulting in a C1h symmetry. The energy difference, though, between this and

the D2d configuration is only very small.

52 CHAPTER 3. Vacancies and Interstitials

formed. The temperature simulation profile was chosen as a linear increase (200 K/fs) up

to a certain maximum temperature that was then kept constant for fifty femtoseconds,

followed by a linear decrease to zero temperature again.

0 100 200 300 400 500 600 700

simulation step

Etot

Ekin

Energy [eV]

Figure 3.6: The diagram shows the total energy

and the kinetic energy of the structure during the

MD simulation at T=1700 K. After 80 steps of lin-

ear increase, the maximum temperature has been

reached. The kinetic energy is determined by the

specified temperatures, the abrupt changes in the to-

tal energy indicate structural changes.

Already during this short simulation time,

some effects were observed: at a maximum

simulation temperature of T=1000 K, only

a rotation of the interstitial on its site could

be observed. At T=1500 K, the two carbon

atoms started to part with each other, but

the energy seemed to be not sufficient to dis-

rupt them definitely. A further increase to

a maximum temperature of T=1700 K was

enough to activate the sublattice migration

process described above. An additional ro-

tation of the defect at the new site was also

observed. Fig. 3.7 shows some snapshots

of this simulation. Below the structures in

Fig. 3.7, the total energy of the whole struc-

ture is shown over the simulation step. The

complete total energy curve as well as the

kinetic energy are shown in Fig. 3.6.

Starting from the (CC)Csplit-interstitial (upper left part of Fig. 3.7, defect atoms marked

with a red C), one of the C-atoms parts from the other (upper right, maximum of simulation

temperature reached). In the third step a bond to another C (labeled with a blue C) is

formed, with which in the fourth step a (CC)Csplit-interstitial forms. This (CC)Cturns

(several times) on its site (fifth picture), until the final structure (bottom right corner) is

reached (structure cooled down to T=0 K). Although the temperatures in MD simulations

have to be treated with care, the temperature found here agrees well with the temperature

range assigned to this process in Chapter 5.

3.3 Interstitial Recombination with Vacancies

Since especially low energy implantation creates a large number of close vacancy–interstitial

(Frenkel) pairs, recombination processes of such defects are as important for a complete

picture of what happens during annealing as long-range diffusion. For silicon- and carbon-

vacancies and carbon split-interstitials, the recombination mechanism depending on the

initial distance of the defects has been investigated. This leads to a definition of a capture

radius of the vacancies for this split-interstitial, which determines from which distance on

sublattice migration becomes the more important mechanism compared to direct recombi-

nation.

A carbon split-interstitial can be first neighbor of a silicon vacancy. The consequence is a

direct recombination without an energy barrier separating the structures.

If a carbon vacancy happens to be a second neighbor of a (CC)C, there are already

non-vanishing energy barriers to be overcome for the recombination. Depending on the

orientation of the (CC)Cto VC, energy barriers of 0.0, 0.15 or 0.34 eV were calculated with

3.3. Interstitial Recombination with Vacancies 53

Figure 3.7: Snapshots of a molecular dynamics simulation of the migration of a (CC)Cwith a maximum

temperature of T=1700 K. See text for details.

54 CHAPTER 3. Vacancies and Interstitials

energy gains of 11.1, 10.0 or 11.1 eV due to recombination.

For a VSi in a third neighbor distance, the highest barrier was found to be 0.91 eV (vari-

ations between different orientations are smaller than for the second neighbor VC), while

the energy gain is 12.2 eV.

Finally, a carbon vacancy in fourth neighbor distance results in a barrier of 0.98 eV, gaining

11.9 eV by recombination.

These results show the quite long-range influence of VCand VSi on (CC)Cwhich can be

understood from the large distortion the split-interstitial causes in its neighborhood. The

very small energies needed to activate the recombination together with the high energy

gains create a strong thermodynamical driving force for vacancy-interstitial recombination

which might explain the annealing stage at 200◦C for the carbon vacancy in as-irradiated

material. The energy barriers and gains obtained for the silicon vacancies show furthermore

that not only recombination to the perfect lattice, but even to carbon antisites is likely.

The sublattice migration mechanism for long-range diffusion will, according to these results,

become dominating for interstitial–vacancy distances larger than five neighbors at the

earliest. Depending on the defect densities and the uniformity of the defect distribution

after implantation, this mechanism may start after the recombination of close vacancy–

interstitial pairs has terminated, when vacancies and interstitials are left at larger distances,

only.

Chapter 4

Aggregation of Antisites

Antisites belong to the intrinsic defects with the lowest formation energies in SiC. Neverthe-

less, their existence has not yet been proved experimentally. The reported observation of

the silicon antisite with EPR and ENDOR measurements [87] could not be proved uniquely.

For the isolated carbon antisite, there is no hint from the experimental side, caused by its

electrical inactivity and therefore invisibility to common characterization methods.

Nevertheless, the results of Chapter 3 strongly suggest the existence of antisites after im-

plantation processes. During the post-implantation annealing phase, they may play an

important role in the formation of complexes of either intrinsic or extrinsic defects (com-

pare also Chapter 5). Therefore, their migration and aggregation behavior is worth to be

investigated.

Furthermore, recent photoluminescence measurements indicate a strong correlation be-

tween the local vibrational modes measured for the DIluminescence and those calculated

for the pair of a carbon and a silicon antisite. The conditions under which the DIlumines-

cence is measured also suggest the presence of silicon vacancies in the sample [88], which

agrees well with the mechanism we propose for the creation of antisites in the following

section.

Before we turn to the investigation of creation mechanisms for the antisite pair, its forma-

tion energy in 3C-SiC and 4H-SiC as well as its vibrational spectrum are discussed.

4.1 The Antisite Pair CSi SiC

4.1.1 Formation Energy of CSi SiC

In 4H-SiC, formation energies of the two types of single antisites CSi and SiCon the cubic

and hexagonal sites of the crystal lattice have been calculated. The formation energies

obtained with SCC-DFTB and ab initio results from Ref. [89] for the antisite pair in the

different possible orientations are listed in Table 4.1. Below the line, the binding energies

of the pairs are given. Various supercell sizes and numbers of relaxed atoms have been

tested. The analogous results for 3C-SiC are shown in Table 4.2.

The energy differences between the hexagonal and cubic sites of the 4H-lattice and the dif-

ferent orientations of the pairs are of the same order of magnitude as the energy differences

between the different supercells and restrictions. The same is true for a comparison of 3C-

SiC and 4H-SiC. In 3C-SiC, the comparison with the ab initio results from Ref. [7] shows

that especially the carbon antisite has a much lower formation energy within the SCC-

DFTB calculations than in the ab initio calculation. Consequently, also the CSi SiCpair

56 CHAPTER 4. Aggregation of Antisites

Table 4.1: Formation- and binding energies of antisites and antisite pairs in 4H-SiC. All values in [eV].

Sites 128 cell, 2NN 128 cell, all 4H-240, 2NN 4H-240, all

CSi cub 3.12 2.96 2.70 2.51

CSi hex 3.22 3.05 2.80 2.60

SiCcub 4.78 4.48 4.36 4.01

SiChex 4.68 4.36 4.65 4.01

SiCcub-CSi cub 5.19 4.95 4.76 4.49

SiCcub-CSi hex 5.13 4.86 4.71 4.42

SiChex-CSi cub 5.34 5.09 4.92 4.65

SiChex-CSi hex 5.23 5.01 4.80 4.54

SiCcub-CSi cub 2.71 2.48 2.30 2.03

SiCcub-CSi hex 2.87 2.67 2.45 2.19

SiChex-CSi cub 2.46 2.23 2.43 1.87

SiChex-CSi hex 2.67 2.40 2.65 2.07

Table 4.2: Formation- and binding energies of antisites and antisite pairs in 3C-SiC. All values in [eV].

Sites 128 cell, 2NN 128 cell, all 3C-216, 2NN 3C-216, all Ref. [7]

CSi 2.70 2.58 2.69 2.49 3.7

SiC4.47 4.25 4.42 4.06 4.6

SiC-CSi 4.93 4.72 4.85 4.59 5.8

SiC-CSi 2.24 2.11 2.26 1.96 2.5

is more than 1 eV lower in energy. In the binding energies, these deviations are lower. A

reason for these deviations may also lie in the relaxation constraints in the ab initio calcu-

lation, since especially around the SiCthe lattice is strongly deformed, and the relaxation

constraints in the SCC-DFTB calculations increases the formation energies by about 0.3

eV.

The similarity in the results lets us choose one polytype for further investigations: we chose

4H-SiC. The supercell was mainly the 240-atom cell, in particular for all calculations of

diffusion mechanisms. Total energy calculations, especially for the larger aggregates were

also performed in two larger supercells containing 384 or 480 atoms. Variations were found

to be in the same range as for the values in Table 4.1. In the last sections of this chapter,

for the discussion of larger onion-like aggregates, a cubic 3C-SiC supercell with 512 atoms

was used, additionally.

4.1. The Antisite Pair CSi SiC57

4.1.2 Properties of the Antisite Pair

SiC

−1 %

−12 %

+15 %

Figure 4.1: Geometry of the antisite

pair CSi SiC.

The antisite pair CSi SiChas C3vsymmetry in

the cubic polytype and has only slight devia-

tions from this symmetry in the 4H- and 6H-

polytypes, resulting in C1hsymmetry. The lat-

tice around the antisite pair is distorted similarly

as for the isolated antisites: The Si-Si-bonds at

the SiCare ≈15% (isolated SiC: 12 %) elongated

compared to the ideal Si-C-bond length, the C-

C-bonds at the CSi are contracted by ≈12% (iso-

lated CSi : 11 %). The Si-C-bond between

the two antisites is with ≈1% only slightly con-

tracted, thus relaxation results mainly in a shift

of the complete pair along its axis, compare

Fig. 4.1.

In order to calculate the internal energy and the entropy of formation of the antisite pair,

the vibrational spectra of the antisite pair and the ideal bulk have been calculated within

SCC-DFTB and AIMPRO, in a 64 atom supercell of 3C-SiC1. The antisites, its first and

CSi

41 atom vibr.

64 atom cell,

SCC−DFTB

AIMPRO

Temperature [K]

dU−TdS [eV]

600 800 1000 1200 1400 1600 1800 2000

−0.2

400

−0.3

−0.4

−0.5

−0.6

−0.7

−0.8

−0.9

−1 200

Figure 4.2: Using the vibrational spectrum calculated within SCC-DFTB and AIMPRO, the internal

energy and the entropy of formation have been calculated for the antisite pair in a 64 atom supercell of

3C-SiC. 41 atoms were allowed to vibrate.

second neighbors as well as the third neighbors of the silicon antisite were included in the

calculation of the vibrational spectrum. Including additionally the third neighbors of the

carbon antisite did not change the SCC-DFTB-result significantly2. The correction to the

total energy due to vibrations, ∆U−T·∆S, shows neither large deviations between the

two methods, see Fig. 4.2: for temperatures as high as T=2000 K, energy differences are

still less than 0.1 eV for an energy correction of ≈0.9 eV.

Taking a closer look at the vibrational spectra in Fig. 4.3, no significant dependence on

the supercell size (64 or 216 atoms) can be seen within the SCC-DFTB-calculations, but

1In SCC-DFTB, the same calculations were performed for a supercell containing 216 atoms, yielding

nearly exactly the same result.

2For AIMPRO this calculation was not available

58 CHAPTER 4. Aggregation of Antisites

DFTB, 64

DFTB, 216

AIMPRO, 64

−0.6

−0.4

−0.2

0.2

0.4

0.6

200 300 400 500 600 700 800 900 1000 1100

Intensity [arbitrary units]

Wavenumber [1/cm]

Figure 4.3: Positive y-axis: Vibrational spectrum calculated within SCC-DFTB in a 64-atom (red line)

and a 216-atom supercell (blue line). Negative y-axis: AIMPRO calculation in a 64-atom supercell. 41

atoms, that is two neighbor shells of the carbon antisite and three neighbor shells of the silicon antisite,

were allowed to vibrate. The green arrows indicate the localized modes.

substantial differences between the SCC-DFTB- and the AIMPRO results.

Although the phonon band gap is as usual overestimated in the SCC-DFTB calculations3

and the complete spectra look rather different in AIMPRO and SCC-DFTB (compare

Fig. 4.3), a comparison of the localized modes in the phonon band gap shows a very good

agreement between the two methods. Three localized modes are found in the phonon band

gap: a degenerate emode and an a1 mode are found at 626 cm−1and 635 cm−1induced by

stretching vibrations of the three Si-Si-bonds next to the SiC. The AIMPRO calculation

yields 627 cm−1and 641 cm−1for these modes. At 682 cm−1(AIMPRO 698 cm−1) a C-C

twisting mode is found, which is supposed to be forbidden for any PL transition[30].

Since the SCC-DFTB frequencies are much closer than the AIMPRO frequencies and,

Figure 4.4: Part of the DIphotoluminescence spectrum as measured in 1972 at two different temperatures

after annealing of an ion-bombarded sample at 1300◦C [28]. The localized mode is marked with an arrow.

with the same broadening, consequently overlap much stronger than these, the localization

3Compare the discussion in Chapter 2

4.1. The Antisite Pair CSi SiC59

can not be read out of the plotted spectrum as easily as in the AIMPRO spectrum. The

arrows in the diagram in Fig. 4.3 indicate the unbroadened frequency values, for which the

localization has been calculated by investigating the eigenvectors belonging to the respec-

tive modes.

The calculated localized modes suggest an identification of the antisite pair with the very

common DI-center, for which recent photo-luminescence measurements report two lines at

661.3 cm−1and 668.7 cm−1for 3C-SiC, which can be correlated with the calculated e-

and a1-modes. Probably the first photoluminescence spectra of the DI-center were already

measured by Choyke et al. in 1972 [28]. The localized vibrational mode can also be seen

in Fig. 4.4 which shows the spectra for two temperatures. Experimental and theoretical

details about this possible assignment can be found in Ref. [30].

Besides the possible technological relevance of the antisite pair, the comparison of the vibra-

tional spectra of the antisite pair has shown that localized modes and integrated quantities

like the vibrational entropy do not show significant deviations from the results obtained

with ab initio methods, although the total spectra look rather different. Variations between

these calculated values are in any case smaller than the variations that have to be con-

sidered for a comparison with experimental data and should therefore not be rated too high.

But how and under which conditions can such an antisite pair be created? In the fol-

lowing three sections all conceivable mechanisms are systematically investigated, starting

with exchange processes in the perfect lattice, and eventually considering vacancy assisted

mechanisms. An extension of the mechanisms to larger aggregates and an examination of

the influence of the vibrational entropy close the chapter.

4.1.3 Creation of the CSi SiCPair in the perfect SiC-lattice

In the perfect lattice an antisite pair can only be created by the exchange of a Si and a

C atom in a rotational movement of a Si-C pair. A simple rotation of a Si-C pair in a

(11¯

20) plane of the crystal is energetically very costly (≈12 eV). The energetically most

favorable mechanism we found is a concerted exchange mechanism as proposed by Pandey

for selfinterstitials in silicon [90].

reaction coordinate r

change of rotational direction

construction of the





r=x+y





Si Si

Figure 4.5: Geometry and movement of the atoms that

create the antisite pair. In the inlet, the construction of

the reaction coordinate ris shown (see text for details).

The Pandey process consists of a simul-

taneous movement of the pair with a

gradual change of the rotational plane

by 60◦. Rotation starts in a (11¯

20)

plane, then the rotation plane gradually

changes. After its rotation by 30◦, the

saddle point geometry is passed (here

the axis of the pair is in a (10¯

10) plane).

Finally, after another 30◦rotation of

the rotation plane the rotation of the

pair is finished in the next (11¯

20) plane

(60◦to the first plane). The movement

of the atoms is marked with arrows in

Fig. 4.5, the respective minimal energy

path is shown in Fig. 4.6. The reaction

coordinate rhas been constructed by

calculating the projection of the actual

positions of the diffusing C– and Si–atoms separately onto the connecting vectors from

60 CHAPTER 4. Aggregation of Antisites

their initial to final positions. In the inset in Fig. 4.5 these two projections are marked

with xfor the C-atom and yfor the Si-atom. The reaction coordinate is then obtained

by summing up these two contributions r=x+yfor each geometry during the diffusion

process.

Antisite pair

Ideal bulk

0 1 2 3

Energy [eV]

5 6 1.65

1.7

1.75

1.8

1.85

1.9

Reaction coordinate r

Bondlength [Å]

Figure 4.6: Energy difference between the actual geome-

try and the ideal crystal is shown during the process (left

ordinate); on the right ordinate the bondlength of the ro-

tating Si-C pair is shown.

At the saddle point geometry the dis-

tance between the moving C and Si

atoms reaches its minimum length with

a contraction of 11% compared to the

ideal bulk bond length, see Fig. 4.6.

The energy barrier between the ideal

bulk and the SiC-CSi defect for this con-

certed exchange mechanism has been

calculated to be 10.5 eV, i.e. ≈1 eV

lower than the barrier for the two atoms

rotating in a (11¯

20) plane. However,

this value of 10.5 eV is still very high,

indicating that direct antisite pair for-

mation is unlikely to happen in a per-

fect lattice.

4.1.4 Migration of Antisites in the ideal lattice

If isolated carbon- and silicon-antisites are present in the lattice, as can be expected due

to the recombination processes discussed in the previous chapter, pair formation might

happen by aggregation.

One of the antisites or both of them have to migrate through the lattice until they meet.

Assuming no other defects around, such a process is, again based on exchange processes,

but not with the direct neighbor as in the direct mechanism described before, but instead

with a second neighbor. These sublattice migration processes of SiCand CSi have been

investigated, but turn out to be energetically very costly: CSi migration can be activated

with 11.7 eV, for SiCmigration 11.6 eV are needed. The energies are even higher than for

the Pandey process, so that these mechanisms have to be ruled out, as well.

Consequently, the formation of antisites requires some support from other defects — the

most promising candidates are vacancies: The presence of silicon- and carbon–vacancies in

the material can facilitate the migration of other defects, since they can help to keep the

lattice distortion small in a migration process. This holds also for the creation of antisite

pairs [40]. Furthermore, their existence in a sufficient concentration is certain under the

assumed conditions.

4.1.5 Pair creation by vacancy migration

An exchange of a Si-C pair close to a vacancy requires considerably less energy than in

the perfect lattice. However, a vacancy can also help the formation of an antisite pair by

avoiding any costly exchange of atoms. For both vacancies, the activation energies of such

mechanisms have been found to be by 1 – 2 eV lower than for the respective exchange

4.1. The Antisite Pair CSi SiC61

Si Si

5.8 eV

EB4.7 eV

Figure 4.7: Antisite pair formation next to VC(left) and VC-CSi (right).

process.

The mechanism is as follows: Suppose we have a single carbon vacancy. Instead of ex-

changing its site with a neighboring carbon atom, one of the silicon ligands can as well

move into the vacancy, forming a silicon antisite. One of its carbon neighbors can take

the silicon site by following immediately. This mechanism is illustrated in the left part of

Fig. 4.7.

With the carbon vacancy having moved one site further on the carbon sublattice, an an-

tisite pair has thus been created. This process can be activated with an energy of 5.8 eV,

which is a large reduction compared to the direct exchange process in the perfect lattice.

In the ”2+” charge state, activation needs 6.4 eV, in ”2-” the activation energy is reduced

to 4.6 eV, confirming the trend to lower activation energies for higher negatively charged

structures as already discussed in Chapter 3.

E [eV]

CSiVC

VCSiC

CSi

CSiSiC

∆

1.7 eV

SiVC+

1.8 eV 2.9 eV

5.8 eV

−2

−1

Figure 4.8: Energy change during the vacancy migration process starting from either VC(blue curve) or

VSi (red curve). For better comparison, the energy of the initial structures has been chosen as reference

energy for each of the energy curves. In case of VSi , the mechanism consists of two steps: the formation of

VCCSi , followed by the CSi SiCcreation by the migration of VC.

Starting with a silicon vacancy VSi , instead, the first thing to happen is the formation of a

VC-CSi pair as described in Chapter 3. This process requires only about 1.7 eV activation

energy and leads to an energy gain of 1.8 eV compared to the silicon vacancy. In a next

step, an antisite pair can be created by migration of VCas described above, see the right

part of Fig. 4.7. The activation energy for this process is only 4.7 eV and the formation

energy of the antisite–pair from VC-CSi is only ≈2.0 eV [40]. The analogous process to

62 CHAPTER 4. Aggregation of Antisites

the VCmigration, i. e. a carbon atom moving into the silicon vacancy and the next silicon

atom following, is not possible, since the resulting structure of this process, an antisite

pair next to a silicon vacancy CSi SiC-VSi is not stable, compare Section 3.1.3. The silicon

antisite immediately recombines with the silicon vacancy.

The change in energy during the migration process is plotted in the diagram in Fig. 4.8

for both vacancies. For better comparison, both energy curves have been plotted with the

energy of their respective initial structures as reference. The red curve denoting the mech-

anism which starts from a silicon vacancy shows two barriers and a minimum in between.

For a fast process, no energy might dissipate in between, so that the overall barrier that

has to be overcome in order to create the CSi VC+CSi SiCcomplex is determined by the

highest of the two barriers, namely 2.9 eV. If the CSi VCpair already exists or the energy

can for other reasons dissipate at this step, the barrier of 4.7 eV determines the process.

In any case, assuming only the presence of vacancies in the crystal, this mechanism, which

includes a movement of a carbon vacancy, is the most favorable for the creation of antisite

pairs4.

4.2 Larger Aggregates of Antisites

With an activation energy not much higher than that calculated for the sublattice migra-

tion of vacancies, the creation of an antisite pair becomes imaginable as a high temperature

process. There is – at first sight – no reason, why aggregation of antisites should stop with

the formation of pairs, making the investigation of larger aggregates interesting.

The first question is, which geometrical arrangements are stable, which of them are most

favorable, then the question for their mechanism of creation has to be answered.

4.2.1 Stability of various arrangements

CSi

SiC

CSi

Si Si

Figure 4.9: Several orientations of two antisite

pairs: they can be oriented parallel to each other

(Left:), ring-like, as shown in the middle, or stacked

upon each other (right).

We start with an investigation of stoichiomet-

rically built up structures which can be imag-

ined to be aggregates of pairs of antisites.

For such aggregates consisting of nantisite

pairs, we have investigated the most stable

geometry. We distinguish three principally

different arrangements according to their di-

mensions. First, antisite pairs can be stacked

upon each other, forming a one-dimensional

chain of antisites. Second, the antisite pairs

can be arranged in a plane, e. g. parallel to

the polar direction of the crystal, resulting

in a two-dimensional planar defect. Finally,

three-dimensional aggregates can be formed.

4In Section 4.4, another vacancy assisted mechanism for the mobility of isolated antisites will be dis-

cussed. Here, a discussion of larger aggregates of antisites seems, though, to make more sense, since an

extension of the mechanism described in this section to the formation of larger aggregates suggests itself.

4.2. Larger Aggregates of Antisites 63

hexagonal

stacked

parallel

234

10 6 7

n antisite pairs

number of Si−Si and C−C bonds

Figure 4.10: Number of ”wrong” bonds in the aggregate of nantisite pairs.

The various possible orientations of two antisite pairs towards each other are shown in

Fig. 4.9. They can be arranged parallel as shown in the left part of Fig. 4.9, ring-like

(center part of Fig. 4.9) or stacked upon each other along a crystal direction (right part).

This distinction can be maintained for larger aggregates of more than two pairs, as well,

leading to the classification of one-, two- or three-dimensional aggregates5.

The number of ”non-SiC-bonds”, i. e. the number of (compared to the Si-C-bonds) short

C-C-bonds and long Si-Si-bonds induced by the antisites is a good quantity to compare the

different geometries. The resulting distortion is a measure for the formation energy and,

thus, for the stability of the structure. For two pairs, as shown in Fig. 4.9, six Si-Si- and

six C-C-bonds can be counted in the parallel arrangement, while only five of each non-SiC

bond are found in the two other structures. A comparison for aggregates of nantisite pairs

was done in Fig. 4.10, showing the number of non-SiC-bonds that a structure in each of

the three cases must have at least.

n=15

n=3 n=8

Figure 4.11: Geometries of the most stable two-

dimensional aggregates of antisites. Structures were cal-

culated in two different supercells, one of which was more

extended in the lateral direction (384 atom supercell), the

other more in the polar direction (480 atoms). The differ-

ences are found to be negligible.

In the parallel arrangement this num-

ber grows fastest with the number nof

antisite pairs. While there is no dif-

ference for two antisite pairs, the ring-

like (”hexagonal”) arrangement has the

smallest number of non-SiC-bonds for

three and more pairs. With three,

five,. . . antisite pairs, a whole ring can

be completed, leading to an especially

small number of non-SiC-bonds in this

arrangement.

The energetical behavior turns, in fact,

out to be analogous, showing that the

parallel arrangement is highest in en-

ergy, and even becomes instable for

large n. This can be understood, since

the strong expansion on the side of the

silicon antisites and the strong contrac-

5Non-stoichiometrical three-dimensional aggregates will be discussed later on

64 CHAPTER 4. Aggregation of Antisites

tion on the side of the carbon antisites can be balanced in no way. For the other two

arrangements, the energy

∆E(n) = E(n pairs)−nE(1 pair)

n(4.1)

is gained by bringing npairs together, see Fig. 4.12. The limiting cases are an infinite

chain, as shown here along the [0001] axis, in the one-dimensional case, and a completely

inverted bilayer in the two-dimensional case, see Fig. 4.13. In the latter, the bonds inside

the bilayer have the ideal Si-C bond length, and the longer Si-Si-bonds on one side are

compensated by the shorter C-C-bonds on the other side, causing only a small shift of the

complete bilayer along the [0001] direction.

∆E [eV]

Number of antisite pairs

−4

−3.5

−3

−2.5

−2

−1.5

−1

−0.5

0 5 10 15 20 25 30 35 40 45 50

Figure 4.12: One-dimensional (1D) and two-dimensional (2D) aggregates of antisites. The 1D structures

up to the infinite chain were calculated in a 480-atom supercell, while for the 2D structures the more

laterally extended 384-atom supercell was used, in order to be able to calculate some larger complexes. The

dashed lines indicate the limits of the infinite chain (1D) and the completely inverted bilayer (2D).

The energy gain ∆E(n) upon adding pairs of antisites to the aggregate, converges like

1/n to the limit of ∆E1D(∞) =-2.5 eV/pair (infinite chain) and ∆E2D(∞) =-3.6 eV/pair

(complete antisite layer), respectively.

Thus, especially the two-dimensional aggregation, where the main stabilizing factor is the

formation of ”inverted rings”, is energetically favorable, and mechanisms for the creation

of such aggregates should, therefore, be investigated.

For three-dimensional aggregates of antisites, we find the symmetrically, but non-stoichio-

metrically, built up ”onion-like” structures to be the most stable arrangements. Created

upon the alternating addition of single antisites of both types, these structures are centered

on one antisite and surrounded by one or more shells of antisites of alternating type. The

smallest ones, SiC(CSi )4and CSi (SiC)4have been found to be extremely stable: the

energy gains due to clustering are as high as 7.1 eV and 7.8 eV compared to the isolated

constituents. Substituting the next shell of atoms by antisites leads to complexes of the

form SiC(CSi )4(SiC)12 and CSi (SiC)4(CSi )12. Here, the energy gain compared to isolated

CSi and SiCantisites has been calculated in a 512 atom containing 3C-supercell to be 28.1

eV and 25.8 eV. In 4H-SiC, the formation energies of these structures might be slightly

4.2. Larger Aggregates of Antisites 65

Figure 4.13: Geometries of the infinitely extended chain of antisites (left) and the completely inverted

bilayer (right).

higher due to reduced symmetry.

A similar 1/r behavior of the energy gain, following the surface-to-volume ratio, can be

expected upon this kind of 3D–aggregation of individual antisites into concentric ”antiphase

shells”.

Now we turn to the question of how such large antisite aggregates can be created, starting

with the two-dimensional aggregates.

4.2.2 A Vacancy’s Spiral Walk for the Creation of Antisite Aggregates

In principle, it is conceivable to extend the vacancy migration mechanism that was de-

scribed for the pair in Section 4.1.4 to the creation of further pairs of antisites [91]. As is

indicated in Fig. 4.14, vacancy migration could go on along the dashed line by this mecha-

nism. The result of a seven-steps ”VC-walk” is shown for example on the right hand side

in Fig. 4.14. In principle, the whole bilayer could be inverted by this mechanism before it

stops.

The whole process and the actual energy value of the structure is shown in Fig. 4.15, start-

ing from an isolated carbon vacancy in the upper left picture and following the migration

process until the lower left corner is reached. Below the geometries of the minimum energy

structures, the part of the energy curve that has been passed by so far is shown. In each

figure, the atoms that will become antisites in the following step are marked with light

cyan circles, while already created antisites are marked with red circles. From the third

pair on, the migration barriers are noticeably lower, which seems to be caused by the ring

formation: In the fourth picture in Fig. 4.15, the first ring of antisites is complete, and

every second step one more ring of antisites is added. These ring structures are energeti-

cally very favorable, since they allow a very good balancing of different bondlengths — in

66 CHAPTER 4. Aggregation of Antisites

C14

VVC

Figure 4.14: Top view of a (0001) plane, creation of antisite pairs by migration of vacancies. Left:

In the first step, the atoms labelled with 1 are supposed to move along the dashed line towards VC.

After this step, the vacancy will be at the site of carbon atom 1. In principle, the mechanism can go

on, as is denoted with the dashed line, while perpetually pairs of antisites are created. Right: After

seven steps, the result looks like this: VChas moved along the dashed line up to the former site of

carbon atom 7. In contrast to the antisites on the edges of the resulting platelet, those antisites

created first are now coordinated as in the ideal lattice, except for the bond perpendicular to the

plane shown here.

the rings, there are only Si–C bonds, while the distortion due to longer Si–Si- or shorter

C–C–bonds is restricted to the direction of the c-axis.

In the last two steps of this mechanism, the effect of the barrier lowering might be super-

posed by the supercell constraints in the plane, caused by the periodic images of the defect,

which are in case of the large antisite aggregates already very close to each other.

Of course, there are not only carbon vacancies present, but also silicon vacancies and di-

vacancies. If starting from a silicon vacancy, instead, the mechanism is nearly the same.

The only difference is the first step, where the transformation to the VCCSi pair happens

before the carbon vacancy can follow the afore described procedure. The creation of the

VCCSi pair in the plane, where the spiral walk will take place one step later, is found to be

slightly more opportune for the activation energies of the following steps than the creation

of a pair perpendicular to this plane.

Slightly different looks the process if starting with a divacancy VCVSi . Although the

mechanism shown in Fig. 4.16, of a moving carbon vacancy is the same in this case as

in case of an isolated VC, ring formation does not occur during the first six steps. The

mechanism starts with the transformation of VCVSi to a VCCSi VCcomplex (which requires

≈2 eV to be activated), then one of the VCperforms the spiral walk centered around the

other VC(which in a larger aggregate has three C-ligands and is thus rather a VSi ). Not

before the seventh step, a ring-like structure has been created. However, the surrounded

vacancy seems to facilitate the migration process, resulting in distinctively lower migration

barriers compared to the mechanisms starting from isolated vacancies.

All the three processes described so far have one thing in common. Regarding again the

energy curve in Fig. 4.15, the problem of the recombination barriers becomes obvious. If

4.2. Larger Aggregates of Antisites 67

5.89 eV 5.95 eV

4.64 eV

4.68 eV

5.36 eV 5.54 eV 4.58 eV

Figure 4.15: Top view of a (0001) plane: stepwise creation of up to seven antisite pairs by migration

of a carbon vacancy. Below the geometries of the minimum energy structures, the part of the energy

curve that has been passed by so far is shown. In each figure, the atoms that will become antisites

in the following step are marked with light cyan circles, while already created antisites are marked

with red circles. The energy barriers are denoted above the arrows.

68 CHAPTER 4. Aggregation of Antisites

3.4

3.3 4.3

4.04.9

Figure 4.16: Top view of a (0001) plane: stepwise creation of antisite pairs by migration of a

carbon vacancy, but the origin is a divacancy VCVSi in this case. Again, the atoms that will

become antisites in the following step are marked with light cyan circles, while already created

antisites are marked with red circles. The energy barriers are denoted above the arrows.

4.3. Contributions of the Vibrational Entropy 69

the energy suffices to overcome the barrier in the direction of creation of antisite pairs, it is

certainly enough to overcome the barrier in the reverse direction which, in addition, leads

to a lower energy structure.

If starting with a silicon vacancy, the VCCSi pair is energetically favorable to be created,

and so is the first antisite pair. However, already for the second pair, the formation energy

is by ≈3 eV higher than that of VSi . Consequently, as shown in Fig. 4.20 (solid black

line), the energy barrier for recombination is distinctively lower than the barrier that has

to be overcome for the creation of the second pair, which requires about 6 eV, i. e. is much

too high. Furthermore, the recombination barriers are lower than the energy barriers cal-

culated for sublattice migration of VSi or VC. Thus, the dissociation of the vacancies after

some steps, which would leave a stable antisite platelet behind, does not seem likely.

4.3 Contributions of the Vibrational Entropy

During the migration process, the lattice has to be rearranged quite strongly, with four

Si-C–bonds being broken and two Si-Si- and two C-C–bonds coming up. Hence, the de-

scription of the energies by the total energies of the structures might be too inaccurate.

Thus, it might be necessary to consider the contributions from vibrational entropy. As

discussed in Section 2.7, their influence may in some cases not be negligible, especially

for processes involving substantial rearrangements of the lattice, which can give rise to

considerable changes in the vibrational frequencies.

−0.3

−0.25

−0.2

−0.15

−0.1

−0.05

0 200 400

−0.35 600

dU−TdS [eV]

800 1000 1200 1400 1600 1800 2000

Temperature [K]

all atoms

2 NB

Figure 4.17: Lowering of the sublattice migration

barrier of VC: the change in dU −T·dS at the

saddle point geometry of the process compared to

the minimum energy structure.

Although such a correction cannot be

expected to further stabilize the larger

aggregates essentially, it might have a

considerable influence on the barriers for

the above described processes on the one

hand and for vacancy migration on the

sublattice on the other hand, making the

dissociation of the vacancy more proba-

ble.

Furthermore, it can certainly not be ex-

pected that the entropical contributions

lower the activation energy needed for

the pair creation by an exchange process

in the ideal lattice (≈10 eV) sufficiently,

so that the process could be expected to

happen at usual temperatures6.

The migration of vacancies on their sub-

lattices and the vacancy assisted pro-

cesses for antisite creation, however, have to be re-investigated in this context.

6For a better readability, the difference between the energy barriers (calculated without consideration

of the entropical contribution) and the free energy barriers (including this term) is referred to as ”lowering

of the barrier”, although this is, of course, not a ”lowering” in the physical sense.

70 CHAPTER 4. Aggregation of Antisites

Example 1: Sublattice migration of vacancies

Considering the vibrational contributions dU −T·dS, the activation energies for the sub-

lattice migration of VCand VSi are lowered very similarly like shown in Fig. 4.17 for the

carbon vacancy. Including all atoms of the supercell yields the black line, restricting the

vibrations of nearest and next nearest neighbors results in the blue curve. At a typical

annealing temperature of T=1800 K, the activation energy for the sublattice migration of

VCis lowered by ∆E≈0.3 eV. The same holds for VSi . No qualitative changes result

from including the entropy in this case. The diagram in Fig. 4.18 shows how the energy

curve changes at this temperature (all atoms relaxed). Even at temperatures as high as

2000 K, the difference between the results with and without this constraint is only 0.05 eV

(compare Fig. 4.17), i. e. clearly less than the accuracy of the method.

CVC

TdS

Saddle point

Energy [eV]

Figure 4.18: Influence of the vibrational entropy on the activation energy for the sublattice migration

of the carbon vacancy.

Example 2: Transformation of VSi to the VCCSi pair

The transformation of the silicon vacancy into the VCCSi pair defect (compare Section

3.1.3) is affected similarly by the entropy contributions: At a temperature of T=1800

K, the activation energy (compared to the silicon vacancy as reference) is lowered by

∆U−T·∆S≈0.3 eV. The resulting pair defect, VCCSi , is lowered by approximately the

same amount compared to the silicon vacancy. Consequently, the recombination barrier

remains unchanged, and the VCCSi pair is further stabilized against recombination to VSi .

The free energy barrier for the process VSi −→ VC-CSi is, then, 1.4 eV, the pair VC-CSi is

by 2.1 eV lower in energy than the silicon vacancy.

4.3. Contributions of the Vibrational Entropy 71

VSi

Saddle point

VSi

TdS

−0.5

0.5

−1.5

−1

1.5

−2

Figure 4.19: Influence of the vibrational entropy on the transformation of the silicon vacancy into the

VCCSi pair.

Now, reconsidering the ”spiral walk” of the vacancy, i. e. the creation of pairs of antisites

along the way around the original vacancy site (see Section 4.2), we can calculate the en-

tropical contributions to all the equilibrium and saddle point structures. Because of the

stronger rearrangement of the lattice, these contributions can be expected to be larger than

those of the first two examples.

Example 3: The spiral walk

The first step, the transformation of VSi , has just been discussed in example 2. The acti-

vation energy for the next step of the process, that means the creation of the first antisite

pair VC-CSi −→ (VC-CSi ) + (CSi -SiC), is lowered to 3.7 eV. The final structure of this

step is then found to be ≈0.3 eV more stable than the silicon vacancy. The result is

shown in Fig. 4.20. For two different temperatures, T=1000 K (blue lines) and T=2000

K (red lines), the Gibbs free energies of minima and saddle point structures are shown in

the diagram together with the energies obtained without considering these terms (black

lines)7. Saddle point geometries are throughout slightly stronger affected than equilib-

rium structures, and the energy differences are as expected larger than those found for

the sublattice migration of vacancies. Nevertheless, no qualitative changes concerning the

energetical order of activation- and recombination barriers show up here, so that the inclu-

sion of entropy cannot change the conclusion that this mechanism will not lead to larger

aggregates of antisites but probably stop after the creation of (VC-CSi )+(CSi -SiC), which,

besides, has a slightly lower free energy than the isolated silicon vacancy at temperatures

above 1000 K.

At this point, it should be noted that a really substantial change in the results of this

section, i. e. a strong stabilization of the larger antisite aggregates, would mean that a

7The connecting lines are just guides for the eye

72 CHAPTER 4. Aggregation of Antisites

tot

VSi

VCCSi

VCCSi CSiSiC

+ 1

VCCSi CSiSiC

+ 2

G(T=2000 K)

G(T=1000 K)

−3

−2

−1

Energy [eV] (compared to Vsi)

Saddle

Figure 4.20: Influence of the entropy on the pair creation mechanism. For two temperatures,

T=1000 K and T=2000 K, the Gibbs free energy has been calculated for the equilibrium structures

and the saddle point geometries. Solid line: One (0001)-face of the supercell was kept fixed. Dashed

line: all atoms included in the vibrational spectrum

spontaneous inversion of the crystal lattice already during growth (which happens at suf-

ficiently high temperatures to activate the processes) became possible. A mechanism to

create larger aggregates of antisites has, therefore, to be more demanding, that means it

has to be dependent on preconditions that are not fulfilled during usual growth processes.

In order to find such a possibility, an alternative mechanism based on single antisite move-

ments is discussed in the next section.

4.4 Mobility of Isolated Antisites

We have seen in Chapter 3 that the recombination of interstitials with nearby vacancies is

among the first processes to happen during the annealing phase. Since this recombination

also leads to antisites, this suggests an aggregation mechanism of isolated antisites to an-

tisite pairs (or larger aggregates). However, we have seen in Section 4.1.3 that activation

energies for sublattice migration of antisites are much too high to be imaginable at realis-

tic temperatures. The ”spiral walk” of vacancies, compare Section 4.2.2, may be able to

explain pair formation, but not the formation of larger aggregates, thus vacancies alone

could neither help in this case.

The high vacancy concentrations may, however, assist the migration processes of isolated

antisites. The following sections deal with the role of vacancies in migration processes of

antisites.

4.4. Mobility of Isolated Antisites 73

4.4.1 The vacancy assisted mechanism

To activate a vacancy assisted mechanism for antisite migration, the vacancies must first

approach the antisites. Thus, sublattice migration of vacancies is essential in any case.

When a vacancy has reached a neighboring site to an antisite on either of the two sublat-

tices, an antisite migration mechanism can start. A silicon antisite SiCcan jump directly

into a neighboring carbon vacancy VC. With an energy of 4.1 eV the jump can be ac-

tivated. ∆Sis nearly zero for this jump, thus there is no lowering of the barrier due to

vibrations.

Si Si Si Si

CSi

3.0 eV

4.1 eV

CCC

SiSi Si

CCC

3.0 eV

SiSi

CCC

SiSi

CCCC

SiSi

CSi

4.2 eV

2.0 eV

Figure 4.21: Left: Two-step mechanism of a Si-vacancy assisted jump of the carbon antisite. The

first step requires 2.0 eV, the second 4.2 eV to be activated, but the intermediate structure VC-

CSi -VCis 2.2 eV lower in energy than the initial structure VSi -CSi .Right: Instead of dissociating,

VSi can as well migrate around CSi . The determining step is shown in the second picture, where

VSi has to move from one C-ligand of CSi to the other. The migration barrier of this step has nearly

the value as for usual sublattice migration of VSi , whereas the other two steps are less costly.

The carbon antisite CSi can perform a similar jump into a silicon vacancy VSi , requiring 4.9

eV activation energy. In this case, however, a two-step mechanism, as shown in Fig. 4.21, is

found to be much more favorable. The first step is a transformation of the silicon vacancy

to the VCCSi –complex, lowering the energy of the structure by about 2.2 eV (compare

the top and the middle structure in Fig. 4.21 and the energy curve in Fig. 4.22). The

activation energy for this transformation is only 2.0 eV. From this symmetric structure,

the first CSi can recombine with the VCto complete the jump (middle and bottom struc-

ture in Fig. 4.21). Assuming no energy dissipation at the intermediate stage, the overall

energy barrier which is important for the mechanism is 2.0 eV. Considering the entropical

contributions does even in this case not lead to a significant lowering of the barrier (for

T<2000K ∆U−T·∆Sis less than 0.1 eV).

74 CHAPTER 4. Aggregation of Antisites

VSi

VSi CSi

CSi

CSi CSi

Energy [eV]

−2

−1

Figure 4.22: In the two-step mechanism of the vacancy assisted CSi -migration an intermediate

structure with a lower energy than the initial and the final structure is crossed.

Up to now, only an exchange of a vacancy and an antisite has happened. To obtain a

directed diffusion process, we need the vacancy to either dissociate from the antisite —

then the antisite has to wait for the arrival of another vacancy — or to move around the

antisite, such that a further antisite jump can take place.

This latter mechanism is shown in the right part of Fig. 4.21 for the silicon vacancy moving

around the carbon antisite. It consists of three steps, whereof two are energetically equiv-

alent. First, the vacancy has to move to one of the other Si-sites surrounding the C-ligand

of the CSi (Eactiv = 3.0 eV). In the second step, migration has to continue to a Si-site next

to the other ligand of the CSi , which is a principally different process (Eactiv = 4.1 eV).

The third step is equivalent to the first one (Eactiv = 3.0 eV), but in the resulting structure

VSi and CSi have changed places – compare the top and bottom parts of Fig. 4.21.

The determining step is shown in the second picture, where VSi has to move from one C-

ligand of the carbon antisite to the other. The migration barrier of this step (Eactiv = 4.1

eV) has approximately the same value as that of usual sublattice migration of VSi (if no

antisite is present), whereas the other two steps are less costly in energy. This can be

understood considering the strong inward relaxation of the C-ligands around a CSi , which

destabilizes the bonds between the first and the second neighbors and thereby lowers the

energy that is necessary for moving the Si-atom to the vacant site. The entropical con-

tribution lowers the activation energy of the determining step again by about 0.2 eV at

T=1800 K.

The dissociation of the carbon vacancy from the silicon antisite can be activated with

4.8 eV. Considering the entropy, a free activation energy of only 4.4 eV is needed at

T=1800 K. Note that these values are slightly above those found for the sublattice migra-

tion of VC, thus determining the activation energy needed for the whole process.

Thus, vacancy assisted migration of both antisites is energetically comparable to plain

sublattice migration of vacancies. This result makes it conceivable that pairs of antisites

as well as possibly larger aggregates of antisites can be formed at conditions under which

4.4. Mobility of Isolated Antisites 75

SiSi Si

CSi

3.4 eV

0.8 eV

C CC

SiSi Si

CSi

C CC

CSi

VSi

C CC

SiSi

4.1 eV

Figure 4.23: Left:The SiC(CSi )2complex can be easily created from the VSi (CSi )2complex that is an

intermediate state during the vacancy assisted movement of CSi .Right: If a silicon vacancy approaches the

SiC(CSi )2complex up to the position shown in the upper structure, a barrier of 4.1 eV has to be overcome

to end up with a SiC(CSi )3complex. The energy is thereby lowered by 3.18 eV.

vacancy migration is observed.

4.4.2 The SiC(CSi )2complex

The two-step mechanism that has been discussed for the silicon vacancy assisted migration

of a carbon antisite leads in the first step to a very stable structure, the VC(CSi )2complex

(see the center structure of the left part of Fig. 4.21).

We consider the case that a carbon antisite via the vacancy assisted mechanism proposed

afore approaches a silicon antisite. If aggregation is supposed to happen this way, then

the situation shown on the left hand side of Fig. 4.23 has to occur in the last step8. As

indicated by the arrows in Fig. 4.23 (left), there are two possible movements of the silicon

antisite. A direct recombination of the SiCwith VC(green arrow) would require 3.4 eV, but

the mechanism indicated with the red arrows, in which the SiCpushes its silicon neighbor

into the vacancy and follows simultaneously, requires only 0.8 eV. The result is (in both

cases) the structure shown below, a SiC(CSi )2complex and in a second neighbor distance

a carbon vacancy that can dissociate in a next step.

Comparing the energy barriers for the second step of the vacancy assisted CSi migration

mechanism on the one hand and the formation of the SiC(CSi)2on the other hand, the latter

is clearly favored. That means, instead of the second step of the CSi migration mechanism

rather the addition of the silicon antisite to the CSi VCCSi complex will take place in this

situation. Furthermore, the structure resulting from the recombination process is by 2.7

eV lower in energy, while the second step of the vacancy assisted CSi migration ends up in

an energetically higher structure, instead.

Thus, the vacancy mediated mobility of a carbon antisite can lead to the formation of

larger complexes. It results, however, not in a simple pair formation, if a SiCis approached

8An intermediate structure with the assisting silicon vacancy next to the SiCwill probably not lead to

an antisite pair by the mechanism discussed here, but the silicon antisite will recombine with the vacancy,

leaving just a CSi and a VCin second neighbor distance behind

76 CHAPTER 4. Aggregation of Antisites

CSi

SiVC

CSi

VSi

CCC

CSi

VSi

CCC

6.4 eV 3.6 eV

Figure 4.24: Left: A ring mechanism for adding one more CSi to a SiC(CSi )3complex. Right: Another

mechanism which adds a fourth CSi to the complex. See text for details.

as described here, but the VSi transforms and two CSi are ”attached” to the SiC, while a

carbon vacancy remains in second neighbor distance and can dissociate later on.

The mechanism can proceed if sublattice migration of silicon vacancies is activated. If a

silicon vacancy approaches the complex obtained so far, consisting of two carbon antisites

and one silicon antisite (and possibly the carbon vacancy on the former site of the SiC,

if it is not yet dissociated), the process sketched on the right hand side of Fig. 4.23 can

take place. The silicon vacancy can move on as denoted by the arrow in the upper part

of Fig. 4.23 (right) with approximately the same activation energy (4.1 eV) as calculated

for undisturbed sublattice migration. The carbon ligand follows immediately without a

further barrier to the vacant site, such that a complex of a SiC, three CSi and one or

possibly two carbon vacancies is obtained. The total energy is lowered by 3.2 eV, showing

that the movement of the VSi in this direction is more likely than a movement onto another

Si-site, instead. The carbon vacancies can dissociate from the complex, and the growth of

onion-like antisite aggregates is settled.

When a silicon vacancy approaches the structure SiC(CSi )3just created, the structure

SiC(CSi )3VSi is obtained, which, in contrast to SiCVSi , is stable against spontaneous

recombination. However, only the very small energy barrier of 0.1 eV separates it from the

1.4 eV more favorable VC(CSi )3, after the silicon antisite has recombined with the silicon

vacancy.

A silicon vacancy that has approached the SiC(CSi )3complex up to a third neighbor site

of the SiCoffers two possibilities of adding a fourth CSi to the complex. The left part of

Fig. 4.24 describes a ring mechanism which involves the simultaneous movement of three

atoms, as denoted by the green arrows. In this mechanism, not only one CSi is created, but

one more pair of antisites9, which could be the starting point for a next shell of antisites

(for the creation of larger onions). The total energy of the two structures is nearly equal,

but they are separated by an energy barrier of 6.4 eV.

As an alternative mechanism, the process on the right hand side of Fig. 4.24 has been

investigated. Here, the silicon atom next to the SiCrecombines with the silicon vacancy,

9compare also the discussion of the spiral walk of vacancies in Section 4.2.2.

4.5. Influence of other Defects on Migration Processes 77

and its carbon neighbor moves to the CSi position. The carbon vacancy created in this

process can, again, dissociate in a next step. The total energy of the structure is reduced

by 3.3 eV by the process, and an energy barrier of 3.6 eV makes the addition of a fourth

carbon antisite and therewith the completion of the SiC(CSi )4quite probable, when the

temperatures are sufficient to activate migration of silicon vacancies (≈4 eV).

4.5 Influence of other Defects on Migration Processes

High defect concentrations due to high doping rates and/or severe lattice damage result in

short distances between the induced defects. So far, we only took into account ”isolated”

processes in an otherwise perfect crystal lattice. However, reality looks different, and

another defect that is present due to the created defect density or a mobile atom – maybe

involved in some different diffusion process – may appear in the vicinity of the regarded

process.

In this section, only a brief sketch of such influences shall be given, in order to rate possible

effects of the close environment of defects on the activation energies of the investigated

process. Considering a process with three defects, i. e. antisites or vacancies or even

nitrogen-atoms, involved, leads to a defect concentration of 1.2·1021 cm−3if using the

standard 4H-supercell with 240 atoms, as described in Chapter 2.

Especially during the first stage of annealing after the ion implantation process, the vacancy

concentration can be very high. For both VSi and VC, concentrations as high as 1021 cm−3

are reported [97]. Thus, the described arrangement with three defects in a supercell is in

the range of the defect concentration, and the assumptions agree with the experimental

conditions.

How the presence of a defect affects the height of the energy barriers and the direction

of a microscopic jump has been examined for the example of the movement of the silicon

vacancy next to a carbon antisite. As shown above, the mobility of the silicon vacancy

determines the efficiency of the vacancy assisted mechanism for CSi migration. Whether

migration is promoted or hindered depends on the defect, its distance from the place where

the migration process is supposed to happen, and the orientation towards the direction of

migration.

We have investigated two different directions for the migration of VSi : (a) to an equivalent

position with respect to the CSi (”direction (a)”, as would play a role in the process of

moving around CSi in order to obtain directed diffusion), or (b) one site further apart from

the CSi (”direction (b)”, as a first step of dissociation). For these two processes, we compare

the ”plain” process with those processes we obtain under the influence of another CSi , a

SiC, a VC, a VSi or a NC.

The geometrical arrangement of these defects and the direction of migration are indicated

in the six small pictures of Fig. 4.25 together with the energy curves along the diffusion

paths. The defect considered as a perturbation of the examined process is marked in ma-

genta. In the upper left picture, the ”plain” process, i. e. no other defect present, is shown.

The energetically more favorable direction in this case is direction (a). This is still true, if

another carbon antisite or a silicon vacancy are near (see the center left and right picture

in Fig. 4.25). The presence of another CSi lowers both migration barriers, and in case of

direction (a) the resulting (symmetrical) structure is also lower in energy (migration in the

other direction results in a structure equivalent to the initial structure). A silicon vacancy

leads to a lowering of the energy in both cases, but implies a transformation of one of the

silicon vacancies into a VCCSi pair, since two silicon vacancies are not stable on neighboring

78 CHAPTER 4. Aggregation of Antisites

Si Si

SiC

Si CC

Si C

2.91

1.74 VSi Si

2.81

Si CC

Si C

2.46

VSi Si

VSi

Si CC

Si C

2.68

3.34

VSi Si

steps

CSi

Si CC

Si C

3.78 2.61

VSi Si

Si CC

Si C

1.67 3.28

VSi Si

Si CC

Si C

4.05 3.04

Energy [eV]

steps

Energy [eV]

steps

Energy [eV]

steps

Energy [eV]

steps

Energy [eV]

steps

Energy [eV]

4 5 6 7 8 9 1021

−8

−6

−4

−2

12345678910

2.5

−0.5

0.5

1.5

2.5

3.5

1 2 3 4 5 6 7 8 9 10

1.5

−0.5

0.5

1.5

2.5

3.5

−1.5

4.5

1 2 3 4 5 6 7 8 9 10

0.5

−7

−6

−5

−4

−3

−2

−1

12345678910

−0.5

−3

−2

−1

12345678910

−1

Figure 4.25: The presence of other defects may essentially influence the energy barriers of a migra-

tion process. Shown above are two migration processes of VSi as might happen during a vacancy

assisted CSi -antisite jumps and their energy curves under the influence of CSi , SiC, VC, VSi or

NC(marked in magenta). The red (green) curve in the energy diagram belongs to the processes

marked with a red (green) arrow in the structure above the diagram.

4.6. Conclusions for the Formation of Antisite Aggregates 79

sites. A silicon antisite or a carbon vacancy, in contrast, promote migration along direction

(b). The formation of a divacancy leads to a stabilization of the structure shown in the

upper right picture in Fig. 4.25. The silicon antisite implies a recombination with the VSi ,

thus leaving only a VCbehind. A nitrogen atom NCpromotes VSi -migration along direction

(b), resulting in a VSi NCpair which might play an important role in the various annealing

steps of the silicon vacancy and will be discussed in detail in Chapter 5.

Thus, all these defects that can be expected to be created during the implantation process

or in the consecutive annealing process affect the migration barriers and in some case also

the resulting structures rather strongly, if they are in third or fourth next neighbor distance

to the migrating silicon vacancy.

For the process discussed in the previous section, this means that the formation of the

SiC(CSi )2complex from the initial structure, where a SiCapproaches a VC(CSi )2com-

plex, may be prevented by the presence of another defect close to the silicon antisite. A

silicon vacancy would for example rather lead to a recombination with the SiCthan to the

formation of the discussed complex.

4.6 Conclusions for the Formation of Antisite Aggregates

What have we achieved up to now? According to the energetical stability, the formation

of larger aggregates of antisites is not unlikely. The mechanisms described in the previous

sections can under the discussed assumptions lead to larger aggregates of antisites. The

vacancy assisted pair creation is only conceivable for complexes of one or two antisite pairs,

larger two-dimensional aggregates cannot be obtained like this. Onion-like structures can

result from the movement of single antisites. The results discussed in Section 4.4.2 suggest

that aggregation does not go stepwise, e. g. with single carbon antisites approaching a less

mobile silicon antisite one by one, but the SiC(CSi )2complex is formed at once. Mobility

of the silicon vacancy is the basis for further addition of carbon antisites to this complex

and consequently for ”onion-growth”.

First calculations have shown that the completely inverted bilayer, which was shown to

be very stable in Section 4.2.1, has some promising electronic properties. Fig. 4.26 shows

the potential of the structure calculated within SCC-DFTB and with the FHI code. At

the side of the C-C bonds, a potential well forms, while at the side of the Si-Si bonds a

potential wall is created. Some first one-dimensional ab initio calculations resulted in high

tunneling probabilities for free charge carriers through this quantum well/wall structure

[68]. Although two-dimensional antisite aggregation to complexes larger than two antisite

pairs can obviously not be expected to happen in a crystal after growth, it can possibly be

created already during growth.

With the help of atomic layer epitaxy (ALE) it is possible to grow monolayers of sili-

con and carbon in a controlled way on the substrate. ALE is a very new technique for

SiC but already established in other materials. In the 1990s it was for the first time

successfully applied to the growth of cubic and hexagonal SiC on silicon substrate. A de-

tailed description of the method and its application to various materials can be found in

Refs. [92, 93, 94, 95, 96]. The basic principle is the alternating supply of certain reactants

in a comparably low temperature range, typically 600-950◦C [92].

80 CHAPTER 4. Aggregation of Antisites

Energy [eV]

Atomic positions [Å]

[0001]

−10

−5

−10 −5 0 5 10

DFTB, 64 atoms

FHI

Figure 4.26: Potential wall and well of the inverted bilayer, calculated within SCC-DFTB and FHI.

These low temperatures lead naturally to a low growth rate, so that this method is not

suitable for bulk growth. Nevertheless, the low temperatures might be advantageous for

our aim: silicon and carbon monolayers can be grown either hetero- or homoepitactically,

so that during growth along [0001] the polarity of the material can be changed. If the

alternation of the reactants is stopped and the growth procedure, instead, repeated with

the same reactant a second time, a Si-Si or a C-C sequence is created. Like this, stress-free

sequences of silicon and carbon monolayers could be grown, including Si-Si-C-C sequences

like the inverted bilayer. This would offer a variety of possibilities to design two-dimensional

quantum devices.

One crucial point for the realization of such devices is the stability of the created sequences

of Si- and C-monolayers during the growth process. First molecular dynamics investigations

of an inverted bilayer at the surface of a crystal with and without hydrogen passivation

indicate that considerably higher temperatures than used in ALE should be required to

induce either a dissociation of the surface atoms or a re-ordering to a ”normal” Si-C bilayer

with the original polarity. Nevertheless, more detailed investigations are currently under

work both on the experimental side to optimize the growth technology by e. g. varying

the reactants and on the theoretical side, where we continue tunneling calculations and

molecular dynamics simulations.

Another result has been achieved in the previous sections: Having once found a mechanism

for the migration of a defect and having calculated the entropies for minimum energy and

saddle point geometry, we can give an estimate for the diffusivity Dof defects. With the

activation energy Eaof the process, the diffusivity can be expressed as

D=D0exp −Ea

kT ,(4.2)

4.6. Conclusions for the Formation of Antisite Aggregates 81

where D0is the diffusion coefficient

D0=d2ω

2p exp ∆S

k.(4.3)

Here, ∆Sis the entropy difference between the ground state and the saddle point structure,

dis the distance between the two equivalent minima (initial and final structure), and ω

is the attempt frequency of the jump, commonly approximated by the Debye-frequency

(ωD(SiC) = 1.6·1013 s−1).

In case of vacancy assisted diffusion, which will be discussed in the following, pis the

likelihood of finding a vacancy on a neighboring site, that is the vacancy concentration

divided by the number of sites in the sublattice. In case of direct diffusion, p= 1.

An important quantity for aggregate formation is also the diffusivity. An estimate for the

diffusivities of vacancies and antisites on the basis of the mechanism presented above can

be obtained using Eq. 4.2 and 4.3. For the migration of VCand VSi on their sublattices we

get a diffusivity of D(VC) = 2.49·10−15 cm2/s and D(VSi ) = 8.77·10−14 cm2/s at T=1800

K. Thus, the silicon vacancy has a slightly higher mobility than VC, which agrees with

the observation that in as-grown material, where no interstitials are available, VCsurvives

higher annealing temperatures than VSi [98]. Annealing at 1600◦C for one hour seems to

remove VC[21].

To calculate the diffusivity for the vacancy assisted migration of SiCand CSi , we need to

know the concentration of vacancies in the material. At the early stage of annealing directly

after implantation, concentrations as high as 1021 cm−3are observed for both VCand VSi .

Using this concentration, we obtain D(CSi) = 2.70 ·10−10 cm2/s and D(SiC)=4.09 ·10−16

cm2/s for the antisite migration. Thus, the diffusivity of the C-antisite is faster than

that of the Si-antisite or the vacancies, which has to be attributed not only to the lower

activation energy, but also to the larger change in entropy ∆Sat the saddle point. We get

the following ordering of diffusivities: D(CSi )>D(VSi )>D(VC)>D(SiC). Nevertheless,

since the mobility of VSi is required for the CSi motion, D(VSi ) is determining this process.

In case of SiC, D(SiC) is the slowest and thereby determining diffusivity.

Accordingly, the migration of CSi is much easier than SiCmigration. Consequently, the

Si-centered onion SiC(CSi )4should be more likely to be created than the inverse onion-

structure, the formation of which is hindered by the low mobility of SiC. A mechanism for

the stepwise creation of SiC(CSi )4has been proposed in the previous section. Furthermore,

other aggregates with e. g. equal numbers of antisites or an excess of CSi may be formed

by this process [91].

Naturally, the described mechanism of vacancy assisted migration of CSi and also SiCmay,

depending on the conditions, also apply to the formation of differently formed aggregates

of antisites. That the presence of other defects in the close neighborhood can have a

significant influence on the activation energy and also the result of a migration mechanism

has been shown in the previous section. Accordingly, the growth of onion-like aggregates

may be hindered by some other defect, leading to a differently formed, e. g. two-dimensional

antisite aggregate.

82 CHAPTER 4. Aggregation of Antisites

Chapter 5

Nitrogen related Defects in SiC

The construction of SiC-devices requires n-type and p-type doped crystals. As illustrated

in the introduction, nitrogen is the common dopant in SiC for creating n-type material,

because the nitrogen atoms are built into the SiC-lattice as shallow donors. As already

pointed out in Chapter 2, the low diffusivities of defects of nearly any type in SiC imply

that doping has to be done by ion implantation, which, however, has the disadvantage of

creating a massive lattice damage, up to the amorphization of crystal regions, if e. g. the

temperature during implantation is chosen too high [17, 104].

Unfortunately, doping with nitrogen does not lead to a higher concentration of free charge

carriers than ≈1019 cm−3, but saturates according to Hall- and DLTS-measurements at

this value [31]. A post implantation treatment consisting of further implantation with

hydrogen, helium, boron or aluminum and further annealing can lead to a reduction of the

free charge carrier compensation – with ≈30% recovery at ≈700◦C up to a full recovery at

annealing temperatures above ≈1050◦C [32], see Fig. 5.1.

Recovery after annealing [%]

Original free carrier concentration

Annealing Temperature [ C]

00 200 400 600 800 1000 1200

100

Figure 5.1: Recovery of the free carrier con-

centration versus annealing temperature as

measured by ˚

Aberg et al. [32, 33]. At ≈700◦C

30%, at ≈1050◦C 100% recovery are observed.

Since they observe a very efficient passivation,

Aberg et al. propose the formation of com-

plexes with elementary point defects, such as

vacancies or interstitials [32]. As the sili-

con vacancy unlike the carbon vacancy ap-

pears commonly negatively charged in n-type

material, and Coulomb attraction may en-

hance aggregation, they suggest the formation

of VSi NCpairs as reason for the passivation

of the nitrogen atoms and explain the observed

recovery process by the dissociation of these

pairs.

As will become clear later in this chapter, the de-

scribed behavior of the free charge carrier concen-

tration is closely related to the annealing stages of

the silicon vacancy. The low temperature anneal-

ing stage at 750◦C can be ascribed to either recombination with interstitials (Frenkel pairs)

or the transformation into a CSi VCpair, like discussed in detail in Chapter 3. Further weak

and strong annealing stages of mostly uncertain origin have been observed in capacitance

spectroscopy experiments [99], positron annihilation studies [17], and EPR measurements

84 CHAPTER 5. Nitrogen–related Defects

[21]. According to these experiments, a weak annealing stage is found at ≈1050◦C, while

the final strong annealing stage is found around 1400◦C [14, 17]. Also in these works, the

formation of nitrogen–related complexes, especially nitrogen-vacancy pairs, is proposed to

play an important role during the annealing processes.

The question, which nitrogen–related complexes are stable in SiC and for which of them

a formation process exists that has an activation energy that can be correlated to the ob-

served annealing temperatures is still unanswered in the literature. We will show in the

following sections that a consistent picture can be drawn on the basis of an interstitial

based diffusion mechanism.

5.1 Nitrogen–related Pair Defects

In principle, any defect that is present after the implantation process, could form a pair

defect with a nitrogen atom: other nitrogen atoms — leading to NCNSi pairs, vacancies —

leading to VSi NCor VCNSi pairs, or antisites — resulting in NCCSi or NSi SiCpairs. The

first two cases, nitrogen pairs and nitrogen–vacancy pairs are well known defects in dia-

mond [100, 101], and the similarity of both materials in many of their properties suggests

that similar defects can as well be expected in SiC. Experiments indicate that the isolated

nitrogen is usually built in on the carbon site in SiC as NC[102, 103], and the existence

of a NSi is rather uncertain and at least not very probable, as can also be affirmed by our

calculations, see Section 5.1.4.

The formation of NCNSi pairs becomes already doubtful by this, and additionally, the pair

is slightly less stable than its constituents NCand NSi .

Since the vacancy concentration in the implanted region can be expected to be about two

orders of magnitude higher than the concentration of the implanted species [104], it is very

probable that vacancies are involved in the aggregates. Nitrogen–vacancy pairs are, thus,

more likely than pure nitrogen aggregates. This is affirmed by the experimental results of

Kawasuso et al. as well as ˚

Aberg et al.[17, 32].

Based on their positron annihilation and DLTS measurements, they suggest the presence

of nitrogen-vacancy pairs – and, in fact, Vainer and Il’in ascribe their observation of the

so-called P12-center in N-doped 6H-SiC to a silicon vacancy–nitrogen pair [85]. Calcula-

tions of the hyperfine interactions of the VSi NCpair within the Green’s function based

Linear Muffin Tin Orbital (LMTO) method confirm this ascription, although presuming a

different symmetry of the electronic structure (C1hinstead of C3v), see Table 5.1. More

details of these calculations and the revised correlation of the EPR-signal with the defect

model can be found in Ref. [62].

But how can such pairs develop? We want to investigate the processes that could take

place directly after the ion implantation process. That means, to a great deal, nitrogen is

most likely built in substitutionally on carbon sites as NC, and vacancies, both isolated

or in form of divacancies, are present in the implanted and damaged region of the crystal.

Furthermore, carbon split-interstitials that have not yet recombined with close vacancies

can be expected to be found.

5.1. Nitrogen–related Pair Defects 85

Table 5.1: Hyperfine interactions of the VSi NCpair, calculated within LMTO-ASA

EPR Defect state 13C (-1,-1,1) 14N (1,1,1)

a b a b

P12 VSiNC2A0

1C1h 32.6 7.3 –1.6 0.26

Exp. 33.6 4.2 ±2.4 0.18

3.5 eV

6 eV

C C

4.1 eV

SiSi

12.5 eV

Figure 5.2: Creation of a VSiNC-pair by sublattice migration of either the silicon vacancy or the nitrogen

atom.

5.1.1 Pair formation by aggregation of VSi and NC

The most simple mechanism to imagine for the formation of VSi NCpairs is a sublattice

migration process of either of the constituents. To start the discussion of this mechanism

with the sublattice migration of the silicon vacancy towards the NC, we can recall the

results of Chapter 3: sublattice migration of VSi can be activated with an energy of ≈4.1

eV (in the otherwise perfect lattice). This energy would be needed for bringing together

distant VSi and NC. However, due to the high defect densities in the implanted region, it

is rather probable that VSi and NCare already quite close together, which does not leave

the energy barriers unaffected. For a third neighbor distance (as sketched in Fig. 5.2), the

migration barrier is reduced to 3.5 eV.

On the other hand, the nitrogen atom could migrate on the carbon sublattice until it meets

a silicon vacancy. Without the help of other defects, i. e. vacancies, this process is based on

very costly exchange processes between N- and C–atoms. In the otherwise perfect lattice,

such a process, similar to the Pandey-process1requires ≈12.5 eV to be activated and can,

therefore, be ruled out immediately. Regarding, again, a VSi and an NCin a third neighbor

distance, this activation energy is – due to the in this case only threefold coordinated C-

atom – reduced to 6 eV. The migration of the silicon vacancy is, thus, energetically always

more favorable, no matter how far the vacancy and the N-atom are separated.

Thus, if the external conditions, i. e. temperature and defect concentration, allow vacancy

1compare Chapter 4, where this mechanism is described for the formation of antisite pairs in the perfect

lattice

86 CHAPTER 5. Nitrogen–related Defects

migration, the pair might be created by this mechanism (with the nitrogen atom stay-

ing on its site). However, the energies needed to activate this mechanism are too high to

explain the experimental observation of the pair already before high temperature annealing.

A more complicated mechanism of aggregation with a lower activation energy seems to be

needed, therefore. Recalling our results for the migration of carbon split-interstitials, an

interstitial based mechanism might be an alternative. We have, therefore, examined the

interstitials that are present after implantation and their migration mechanisms.

5.1.2 Mobilizing NC— Creation of N-interstitials

Certainly present are carbon split-interstitials. Their migration can, as discussed in detail

in Chapter 3, be activated with 2.9 eV. If a (CC)Cmeets a carbon vacancy or a silicon

vacancy, it will recombine to either the ideal lattice or to a carbon antisite. The question

is now, what happens when a (CC)Csplit-interstitial meets an NC?

In principle, two things could happen: its migration mechanism could stop or change

C(NC)C

(CC)

Energy [eV]

2.0 eV

1.4 eV

3.5

2.5

1.5

0.5

Figure 5.3: A (CC)Csplit interstitial moves to the substitutionally built in NC, resulting in a (NC)Csplit-

interstitial.

direction. On the other hand, the C-atom could move to the site of NC, if this is energeti-

cally more favorable. Actually, the latter alternative is what happens: The C-atom of the

(CC)Ccan with an even lower energy be activated to move to the site of NCthan to carbon

site with a carbon atom, resulting in an (NC)Csplit-interstitial. In Fig. 5.3 the initial

and final structure of this process is shown above the energy along the migration path.

The newly created (NC)Csplit-interstitial is even lower in energy than the (CC)Csplit-

interstitial next to an NC, thus the formation of these split-interstitials is more likely than

the alternative of a change of the migration direction of the (CC)C.

The geometrical structure of (NC)Cis very similar to that of (CC)C, only with slightly

shorter bonds at the N-atom. As denoted in the energy diagram, the energy barrier is 2.0

5.1. Nitrogen–related Pair Defects 87

eV, while 1.4 eV are gained by this process.

N-atoms as ”usual” (instead of split-) interstitials have been found to be instable: without

an energy barrier, they move to a nearby C-site during the relaxation of the structures,

resulting in (NC)Csplit-interstitials. This has been found for all different interstitial sites in

4H- and 3C-SiC. According to this result, it becomes conceivable to find these (NC)Csplit-

interstitials besides substitutionally built in N-atoms directly after implantation, as well.

It is not surprising that these (NC)Csplit-interstitials can migrate in a similar way as the

(CC)Csplit-interstitials. The N-atom of the (NC)Ccan move a C-site further, requiring

≈2.5 eV, thus slightly lower than calculated for the (CC)Csplit-interstitial. The C-atom

of the (NC)Cwould need 3.4 eV to move to the next carbon site, so that the migration of

the N-atom is preferred and will maintain.

One further kind of movement is needed to render nitrogen migration possible through

the SiC-lattice: a change of the direction of migration might be required, e. g. when the

N-atom meets another defect that would otherwise stop the migration process. A rotation

Figure 5.4: Rotation of the (NC)Csplit-interstitial.

of the (NC)Csplit-interstitial on the C-site into another of the six possible orientations of

the interstitial (see Fig. 5.4) requires an activation energy of 0.8 eV – for the (CC)Csplit-

interstitial, a similar value of 0.6 eV was obtained. After the rotation, the sublattice

migration mechanism can be continued.

After all, we have found a way to mobilize NCand can now turn to our aim, to find a

mechanism for the nitrogen-vacancy pair formation.

5.1.3 Recombination of (NC)Cwith divacancies VCVSi

Provided sufficiently high annealing temperatures, nitrogen can be mobilized as described

above and move through the SiC-lattice in form of (NC)Csplit-interstitials. Due to the

high vacancy concentrations after implantation, these split-interstitials will probably very

soon come across a vacancy or a divacancy. Since it seems to be the most simple way

to create a nitrogen–vacancy pair, we will first discuss the case that an (NC)Cmeets a

divacancy VCVSi 2.

Two cases have to be distinguished: (NC)Ccan approach the divacancy either from the

side of the carbon vacancy VCor from the side of the silicon vacancy VSi . If an (NC)Chas

arrived at a second neighbor site of the carbon vacancy of a VCVSi (see Fig. 5.5 (left)), it

will fill up VC, resulting in a VSi NCpair. Only 0.2 eV are needed to activate this recombi-

nation process, and 7.4 eV are gained. On the other hand, (NC)Ccan approach VCVSi from

a C-site being third neighbor to the silicon vacancy of the VCVSi , as in Fig. 5.5 (right).

The N-atom does, however, not move through the divacancy in order to reach the vacant

2The presence of divacancies directly after implantation has been reported from various experimentalists,

and the strong binding energies (see Chapter 3) back up this observation.

88 CHAPTER 5. Nitrogen–related Defects

VSi

Figure 5.5: The N-atom of a (NC)Csplit-interstitial can fill up the carbon site (left) or the silicon site

(right) of a divacancy.

C-site, but stays on the Si-site. The resulting complex is a VCNSi pair. The recombination

process requires 1.0 eV to be activated, and 8.0 eV are gained. In Fig. 5.6, the change in

energy during both migration processes is shown. The ”inverse” pair, VCNSi , is only by

NSiV

E [eV]∆

(NC)C

VSi

1.0 eV

(NC)C

VC+

0.2 eV

Figure 5.6: The energy curves for the two possibilities a (NC)Ccan recombine with a divacancy.

0.8 eV less stable than the VSi NCpair. In the neutral charge state, a migration barrier

for the N-atom of 2.4 eV separates it from the VSi NCpair. In the negative charge state,

the VSi NCpair is by even 1.6 eV more stable, and the transformation of the inverse pair

to this structure costs only 2.1 eV.

However, in the positive charge state, the energetical order is reversed: in this case, the

VCNSi pair is by ≈0.1 eV lower in energy than the VSi NCpair, and the energy barrier

between the two pairs is ≈3.4 eV. According to these results, the nitrogen–vacancy pair

seems to be a bistable defect which, depending on its charge state, can appear in form of

a VSi NCpair (n-type) or a VCNSi pair (p-type). Since we are here generally focusing on

n-type material and the role of nitrogen in p-type material is only of secondary importance,

we can suppose that the concentration of VSi NCpairs prevails.

5.1. Nitrogen–related Pair Defects 89

VSi

Si N

Figure 5.7: If an (NC)Capproaches a silicon vacancy, it does not fill up the VSi but instead kicks out one

of the carbon ligands of the vacancy to an interstitial site ((CC)C).

5.1.4 (NC)Cmeeting isolated vacancies

Naturally, during its migration through the lattice an (NC)Csplit-interstitial can come

across not only divacancies but also isolated vacancies. Meeting a carbon vacancy, the N-

atom recombines with the vacancy as expected, resulting in the substitutional NC: From

a second neighbor distance, less than 0.1 eV are needed for this process, and the resulting

NCis 9.2 eV lower in energy. The migration mechanism is stopped, and the nitrogen atom

can only be re-mobilized, if a (CC)Ckicks it out into a split-interstitial geometry, as de-

scribed in Section 5.1.2.

The situation is substantially different for a silicon vacancy. One might expect a simple

recombination, like in case of the carbon vacancy, resulting in a substitutional NSi . Our

calculations show, however, that such a recombination does not take place, but one of the

carbon ligands of the vacancy is kicked out, instead. As shown in Fig. 5.7, the N-atom of

an (NC)Cthat has approached a silicon vacancy up to a third neighbor distance can kick

out one C-ligand with an activation energy of 2.9 eV, lowering the energy by 1.8 eV. The

N-atom is then three-fold coordinated built in on the C-site – a VSi NCpair results. The

C-atom has in this mechanism been moved to the next C-site where it forms a (CC)Csplit-

interstitial with the C-atom on this site. From there, it can dissociate in a next step and

migrate through the lattice as described before.

Thus, the kick-out mechanism with a ligand of a silicon vacancy is another possibility to

create a VSi NCpair. Though the activation energy needed for this process is higher than in

case of the recombination of (NC)Cwith divacancies, it is of the same order of magnitude

as the activation energy needed for the sublattice migration of carbon split-interstitials.

That means, at low annealing temperatures, a recombination process with divacancies will

be responsible for the creation of VSi NCpairs, while in cases where long-range diffusion is

expected due to larger defect-defect distances, and (CC)Csublattice migration as well as

(NC)Csublattice migration start to become important, the kick-out process can become

accountable for VSi NCpair creation, as well. Additionally, this kick-out mechanism can

contribute to the annealing of silicon vacancies in the respective temperature range3.

3A tentative correlation between the calculated activation energies and experimentally observed anneal-

ing temperatures is given in Section 5.3.

90 CHAPTER 5. Nitrogen–related Defects

5.1.5 The CSi NCpair as an alternative

The VSi NCpair and the inverse VCNSi pair are not the only stable nitrogen–related com-

plexes in SiC. Obviously, the probability of pair formation with vacancies is rather high

compared to the pairing with other defects, since under the assumed conditions the vacancy

concentration is extremely high. However, vacancies may also promote the formation of

other defect complexes without being directly involved in the resulting complex.

A very stable complex to be discussed in this context is the CSi NCpair. Stoichiometrically

equivalent to a NSi , it is by ≈3 eV lower in energy, and the exchange process of the N-atom

with the C-atom requires only ≈4 eV – much lower than the energies found for other ex-

change processes, compare the exchange of NCdefects with second neighbor carbon atoms,

which was discussed in the beginning of this chapter. Nevertheless, 4 eV are too high to

explain the existence of the CSi NCpair in a temperature range below the final annealing

stage of the silicon vacancy (activation energy for sublattice migration ≈4.1 eV), which

has been observed experimentally at ≈1400◦C.

There are, however, alternative ways to create the CSi NCpair with lower activation en-

ergies. Reconsidering the initial structure of the kick-out process between a (NC)Cand a

VSi in the previous section, we find an alternative final structure (compare top structures

in Fig. 5.8). The behavior of the N-atom remains unchanged, but instead of kicking the

carbon ligand to a (CC)Con a neighboring C-site, the N-atom kicks the C-atom throughout

the vacancy to the antisite position CSi . This mechanism requires only 2.0 eV, i. e. 0.9 eV

less than the kick-out process described before. Furthermore, a large energy gain is linked

with this process: the energy of the CSi NCpair is 10.6 eV lower than the (NC)Cnext to

the silicon vacancy.

That means, for the same initial situation, a (NC)Csplit-interstitial and a VSi in third

neighbor distance, these processes have to be regarded as competing processes with the

result, CSi NCor VSi NCand a (CC)C, depending on the disposable energy.

Alternatively, the CSi NCpair can be created as shown in Fig. 5.8 (center and bottom). Part

of the silicon vacancies may have transformed to VCCSi pairs4. An (NC)Csplit-interstitial

might therefore meet such a VCCSi pair rather than a VSi (see the center left structure in

Fig. 5.8). With 0.3 eV the recombination of the N-atom with the carbon vacancy can be

activated, gaining 8.8 eV.

As a third possibility, we considered the case of carbon split-interstitials (CC)Capproaching

a VSi NCpair, see the bottom structure in Fig. 5.8. A recombination of the C-atom of the

(CC)Cwith the silicon vacancy results in a carbon antisite, while the N-atom remains on

its site. With 1.9 eV, this process can be activated, i. e. a lower energy than needed for

the (CC)Csublattice migration, and 9.7 eV are gained.

Thus, we have found several ways to create CSi NCpairs at rather low energetical cost. The

large energy gains linked with these processes indicate the high stability of this pair defect

and make its creation yet more likely.

The role the discussed pair defects, VSi NCand CSi NC, play in the observed annealing

behavior of the silicon vacancy and the observed saturation of the free charge carrier

concentration in highly N-doped material and its recovery as reported by ˚

Aberg et al. is

the subject of the next section.

4This process requires 1.7 eV, compare Section 3.1.3.

5.2. Dissociation or Aggregation 91

()CSi

+Si

0.3 eV

()C

CSi

VNC

VCSi

()C+

SiNC

2.0 eV

1.9 eV

Figure 5.8: Three possible ways to create a CSi NCpair (right hand side). Top: (NC)Capproaches VSi .

Center: (NC)Capproaches VCCSi .Bottom: (CC)Capproaches VSi NCpair. The activation energies linked

with the mechanisms are indicated above the arrows.

5.2 Dissociation or Aggregation

5.2.1 Dissociation of NCVSi pairs

Aberg et al. have found that in N-doped 4H-SiC the reduction of the free charge carrier

concentration upon additional ion implantation with H, He, Al, and B is compensated

in part by post-implantation annealing above 700◦C [32, 33]. At a temperature around

1000◦C they obtain a full recovery of the free charge carriers (compare Fig. 5.1) and suggest

the dissociation of implantation-induced nitrogen-related complexes, probably VSi NCpairs,

as an explanation for this observation.

Complex formation with the newly implanted ions seems unlikely, because the concentra-

tions can not explain the observed efficiency of compensation. With the Debye frequency

ωDebye = 1.6·1013 s−1as attempt frequency ωattempt they derive a dissociation barrier ∆G

from the concentration of nitrogen donors

Nd(t) = N∞

d·1−exp −ωattempt ·exp(−∆G

kBTt) ,(5.1)

where N∞

dis the concentration after full doping recovery. With the measured concentra-

tions for 700◦C they obtain ∆G= 3.2 eV [33] as an estimation of the activation energy for

the dissociation.

If we now consider the mechanisms discussed in the previous sections for the creation of

NCVSi pairs, we find the lowest energy dissociation process for the silicon vacancy mov-

ing one Si-site further, namely 5.5 eV. As discussed in this context, sublattice migration

of N-atoms requires much higher energies, and due to the large energy gains in the re-

combination processes described for the (NC)Csplit-interstitials with divacancies, these

92 CHAPTER 5. Nitrogen–related Defects

mechanisms (in the opposite direction as described before) can be ruled out, as well.

These calculated activation energies can, thus, not explain the dissociation as proposed in

Ref. [33]. As Coulomb attraction may additionally hinder a dissociation into VSi −and

NC+, which further increases the discrepancy between our calculated 5.5 eV and the 3.2

eV proposed in Ref. [33].

An alternative explanation of the experimental findings may be found in further aggre-

gation of nitrogen. The formation of larger nitrogen–related complexes is, furthermore,

encouraged by the observed high ratio of compensation per vacancy.

5.2.2 Further Aggregation: the formation of VSi (NC)n-complexes

Provided high concentrations of nitrogen, the aggregation process does not have to stop

with the formation of VSi NCpairs. Especially at high temperatures, which allow the

activation of migration processes that require higher activation energies, further nitro-

gen atoms can be mobilized. Aggregation is no longer limited to the low energy recom-

bination processes of close nitrogen atoms and vacancies or divacancies, but long-range

diffusion becomes important, since sublattice migration of (NC)C(Eactiv=2.5 eV) and

(CC)C(Eactiv=2.9 eV) and thereby further creation of (NC)Cby the mechanism (CC)C+

NC−→ (NC)C(Eactiv=2.0 eV) can take place. We will, therefore, investigate what happens

when further (NC)Cor (CC)Capproach VSi NCpairs.

The migration of (CC)Cto VSi NCpairs leads to CSi NCpairs, as already discussed in Section

5.1.5. This CSi NCpair will be discussed in more detail in a later section.

Figure 5.9: The very stable VSi (NC)4com-

plex may be a final annealing product of the

silicon vacancy in n-type SiC.

If an (NC)Cmeets a VSi NCpair, it does

not recombine with the silicon vacancy of the

pair, which would result in an NSi NCpair.

Instead, performing a kick-out mechanism

as in case of the isolated VSi (see Sec-

tion 5.1.4), the N-atom can be built in on

the site of one of the remaining three C-

ligands of the silicon vacancy, resulting in

a VSi (NC)2complex and a (CC)Csplit-

interstitial that can dissociate in a next

step.

This mechanism can be imagined to happen also

for one or two more (NC)Capproaching the com-

plex: VSi (NC)ncomplexes can be created by a

subsequent addition of (NC)Csplit-interstitials to

VSi (NC)n−1complexes. Approaching the silicon

vacancy of a VSi (NC)n−1complex, one of the ligands can be kicked out, resulting in a

VSi (NC)ncomplex and a (CC)Csplit-interstitial. The activation energies for these pro-

cesses and the respective energy gains do not vary significantly with n: about 2.9 eV are

needed for activation, while the resulting complex is by ≈2 eV lower in energy.

With n= 4, a completely inactive complex has been created: the four nitrogen atoms of

5.2. Dissociation or Aggregation 93

Table 5.2: Transition levels of VSi (NC)ncalculated within LMTO and experimentally observed acceptor

levels.

Defect Transition LMTO Exp.

VSiNC–2/–3 –0.34 –0.35

VSi(NC)2–1/–2 –0.59 –0.60

VSi(NC)30/–1 –1.17 –1.10

VSi(NC)4inactive

VSi(NC)4(see Fig. 5.9) fully passivate the silicon vacancy. The formation of this complex

could explain the observed high compensation ratio at high N-concentrations.

At the first view, such a complex may seem ”exotic”, but it can be motivated from in-

vestigations of color centers in diamond. There, the respective VN4complex, known as

”B-center” is a very prominent defect center [105]. In SiC as well as in diamond this com-

plex turns out to be extremely stable: the binding energy has been calculated to be ≈10

eV. The complex has Td-symmetry, and the nitrogen atoms relax only slightly outwards

(4 %), see Fig. 5.9. Its electrical and optical inactivity and these geometrical facts will,

however, make its observation and thereby the proof of its existence difficult 5.

The stability of the VSi (NC)4complex is also reflected in the energy needed to activate

its dissociation. The process VSi (NC)4−→ VCVSi (NC)3+ (NC)Crequires as much as

9.9 eV to be activated, but since the energy of the resulting structure is by 8.9 eV higher,

recombination of the (NC)Cwith the VCwill take place immediately.

The lowest energy alternative to destroy the complex is by a (CC)Csplit-interstitial ap-

proaching it and kicking out one of the N-atoms: VSi (NC)4+ (CC)C−→ VSi (NC)3+

(NC)C. The activation energy of this mechanism is 4.9 eV, but still the resulting complex

is by ≈2.0 eV higher in energy. Dissociation is according to these results very unlikely, and

VSi

(NC)4is accordingly a good candidate for a final annealing product of the silicon vacancy.

The formation mechanism of the VSi(NC)4complex is based on the subsequent aggregation

of (NC)Cto the smaller VSi (NC)ncomplexes. It is not yet clear whether these complexes

have been observed experimentally, but there is strong evidence for this due to the results of

Ballandovich and Violina [99]. Using capacitance spectroscopy, they found three deep ac-

ceptor levels in 6H-SiC with energies at EC–0.35 eV, EC–0.6 eV, and EC–1.1 eV, thermally

stable up to 1050◦C [99]. Another experimental indication comes from the positron lifetime

spectroscopy experiments of Kawasuso et al. [17], according to which the disappearance of

these acceptor levels is correlated with a weak annealing stage for vacancy-related defects.

According to our LMTO-calculations, the energies for the uppermost acceptor levels of

VSi(NC)n, n=1,2,3 agree with the experimental data, compare Table 5.2.

If the assignment of the experimental findings to the VSi (NC)ncomplexes should be ten-

able, the observed vanishing of the acceptor levels has to be explained in our model, as

5In diamond, the B-center could neither be seen by positron annihilation, since a distinction between

an isolated vacancy and the complex could not be made.

94 CHAPTER 5. Nitrogen–related Defects

well. The dissociation of the complexes is rather improbable, since not only the VSiNCpair

and the VSi(NC)4complex but also the two complexes with two or three N-atoms have very

high binding energies, and the high dissociation barriers cannot be linked to the reported

temperature of ≈1050◦C.

The alternative explanation may, again, be found in further aggregation instead of disso-

ciation. An annealing product with a shallower donor level would explain the ”vanishing”

of the levels.

5.2.3 Formation of CSi (NC)ncomplexes

Reconsidering the results of Sections 3.2 and 5.1.2, the energy needed for the activa-

tion of (CC)Csublattice migration is only slightly higher than the activation energy of

(NC)Cmigration. Therefore, processes involving migration of (CC)Chave to be considered

at the same or slightly higher temperatures as processes based on (NC)Cmigration. That

means, it has to be investigated what happens when (CC)Csplit-interstitials meet VSi(NC)n

complexes.

For VSi(NC)4, we have already seen that it is energetically more opportune for the (CC)Cto

migrate to another C-site (leaving the total energy of the system unchanged) than to kick

out one N-atom (increasing the total energy by ≈2.0 eV) of the complex.

1.47 Å

Figure 5.10: The CSi (NC)ncomplex.

For the VSiNCpair, we have discussed the formation

of a CSi NCpair already in Section 5.1.5. Only 1.9

eV are needed for this recombination process that

results in the 9.7 eV lower CSi NCpair. The same

mechanism can take place for the VSi

(NC)2complex:

with 2.4 eV the recombination to a CSi(NC)2com-

plex can be activated, see Fig. 5.10. Here, 7.5 eV are

gained by the recombination. A strong off-center

relaxation of the antisite CSi produces together

with an inward relaxation of its two C-neighbors

strong C-C double-bonds (sp2-hybridization) that

are with ≈1.47 ˚

A shorter than in diamond (1.54

A).

This relaxation effect that stabilizes the CSi (NC)2

complex also causes the instability of the CSi(NC)3

complex, where only one C-C bond would remain.

Actually, the electronic properties of the CSiNCpair and the CSi(NC)2complex can explain

the observed vanishing of the acceptor levels. In Table 5.3, the occupation levels of these

two complexes, calculated within the LMTO method, are shown. The CSi NCpair turns

out to be a donor, similar to the isolated NC, CSi(NC)2has been calculated to be a double

donor. The +/0 level is for any polytype close to the conduction band minimum, similar

as the donor level of the isolated NC, which is even shallower. A trend to lower energies

can be seen beginning with the donor level of NCat ECB-0.05 eV, CSi NC: ECB-0.10 eV,

CSi(NC)2: ECB-0.3. . .0.4 eV.

Thus, the acceptor levels of VSi (NC)ncomplexes can turn into donor levels by the incor-

poration of a C-atom due to recombination with a (CC)C.

The CSi(NC)2complex can alternatively result from a (NC)Capproaching a VSi NCpair,

as shown in Fig. 5.11. This kick-out process, in which the N-atom kicks out one of the

C-ligands of the vacancy to the antisite position, requires 1.6 eV to be activated. The

5.3. Correlation of Activation Energies and Temperatures 95

Si Si

Figure 5.11: Creation of the CSi (NC)2complex when a (NC)Capproaches a VSi NCpair.

Table 5.3: Transition levels of CSi NCand CSi (NC)2in the most common polytypes, calculated within

LMTO.

Polytype Defect Transition LMTO

3C CSi NC+/0 2.33 eV

4H CSi NC+/0 3.10 eV

6H CSi NC+/0 2.95 eV

3C CSi (NC)2+/0 2.18 eV

3C CSi (NC)2++/+ 2.09 eV

4H CSi (NC)2+/0 2.66 eV

4H CSi (NC)2++/+ 2.58 eV

6H CSi (NC)2+/0 2.55 eV

6H CSi (NC)2++/+ 2.47 eV

energy gain of ≈9.3 eV affirms the stability of the complex.

5.3 Correlation of Activation Energies and Temperatures

So far, only activation energies have been discussed, but the experimentally accessible

quantity is rather the temperature at which a defect anneals out6. To get an idea of how

the calculated energies can be correlated with the observed annealing temperatures, we

make an attempt to interpret our data with a simple model.

6To deduce activation energies from experimental data by an Arrhenius plot is often not very accurate

and often technically difficult as it requires a large number of measurements in a wide temperature range.

Especially for high temperatures, results depend on the annealing time.

96 CHAPTER 5. Nitrogen–related Defects

5.3.1 Definition of an assignment

Assuming a Boltzmann distribution for the number of stable defects at a given temperature,

the expression

T=ln(ωD·τ)

∆Gact ≈ln(ωD·τ)

Eact (5.2)

can be derived [43], allowing a direct comparison of the annealing temperature Twith

the calculated activation energy Eact for a migration mechanism of a defect. The Debye

frequency ωD≈1.6·1013s−1is used as the attempt frequency for the jump, τdenotes the

lifetime of the defect under these conditions. As only little is known about this quantity,

we call a defect stable if it exists at least between one and one hundred seconds.

If there is a large contribution of the entropy at the saddle point compared to the minimum

[eV]

CSi C

VSi C

VCNSi

CSi C

( )2

CSi

VSi (NC)C

(from )

VSi

(from )

(NC)C

+VSi C

N(from )

VSiV

C+(NC) )

activ

(from

Temperature [K]

CSi

VC(NC)C

+(from )

VSiV

C(NC)

+(from )

(NC)Cmobile

(CC)Cmobile

VSi mobile

500 1500

1000 2000

Figure 5.12: Correlation of the calculated activation energies with a temperature at which the annealing

of the defect is observed. Due to uncertainties that are discussed in more detail in the text a range

of temperatures is assigned to an activation energy instead of one specific value. Sublattice migration

processes are denoted above the so obtained shaded region, other processes below.

energy structure, the free energy of migration ∆Gact =Eact −T·∆Shas to be used in

Eq. 5.2 instead of Eact. Since this term is temperature dependent, the expression for T

becomes

T=ln(ωD·τ)

1 + ln(ωD·τ)

kB∆S·Eact .(5.3)

The change in entropy ∆S=Ssaddlepoint −Sminimum is in most cases positive, so that the

additional factor is <1, and the activation energy Eact has to be assigned to a lower tem-

perature value.

5.3. Correlation of Activation Energies and Temperatures 97

While the entropical contributions lower the energy barriers and the correlated tempera-

tures, other influences can increase the barriers and therewith the temperatures. As shown

in Chapter 4 at the example of the migration of a silicon vacancy, the vicinity of other

defects can increase or lower the migration barriers. Furthermore, different charge states

can be the reason for an increase or decrease of the activation energies. Due to the inho-

mogeneity of a sample, the Fermi level may not be constant but locally varying to some

extent. This may cause the coexistence of various charge states of one defect, implying a

variation in the activation energies. Such effects cannot be known exactly and can there-

fore not be included in detail. Furthermore, the experimentally observed temperatures

can vary about ≈100 K due to technological reasons. Another uncertainty lies within the

calculated activation energies. A comparison of the activation energies calculated within

SCC-DFTB with available ab initio data (sublattice migration of vacancies, formation of

the VCCSi -pair) [12] suggests deviations up to ≈0.5 eV. As the assignment is special-

ized to our results of the silicon vacancy, these deviations may, however, be much smaller.

Nevertheless, within our calculations the obtained activation energies define at least upper

limits, if in spite of intensive search the saddle point configuration with the lowest energy

should have been missed. Since the main purpose is to give rather a qualitative overview

than an exact assignment of Tand Eact, these considerations lead us to the definition of a

range of temperatures that is correlated with Eact. This has been done in Fig. 5.12, where

the slope of the temperature–energy curve varies around an average value of 500 K ∼1

eV.

5.3.2 Correlation of the calculated values with experimental findings

The calculated values for the activation energies can now be compared with help of this

assignment. Above the shaded area in Fig. 5.12, the sublattice migration processes are

indicated. The silicon vacancy becomes mobile, i. e. long-range diffusion by sublattice

migration starts, at a temperature of 1600◦C [21], which agrees well with the calculated

≈4 eV for the activation energy: the above described assignment would suggest a corre-

sponding temperature of 2000 ±200 K.

The assignment of the activation energies for the mechanisms marked below the shaded

area has been made under the assumption of non-equilibrium conditions. Directly after

implantation, the defect- (especially vacancy-) concentrations are much higher than in ther-

mal equilibrium, creating a thermodynamic driving force that tends to bring the system

into equilibrium by certain annealing processes. We can, therefore, consider the presence

of all the defects involved in the mechanisms, in order to compare the activation energies,

only. Otherwise, in thermal equilibrium, the formation energies of e. g. vacancies for the

vacancy assisted processes have to be considered, additionally (compare Sections 2.1.1 and

2.1.2).

The transformation of the silicon vacancy to the CSi VCpair (activation energy ≈1.7 eV)

has to be assigned to a temperature of ≈850 K, which also agrees quite well with the

observed ≈750◦C [10].

In agreement with experimental observations is, furthermore, the vanishing of divacancies

below room temperature: the vanishing of the signal of the divacancy can be explained

by a direct recombination with (NC)Csplit-interstitials, resulting in either VSi NC- or

VCNSi -pairs, both processes are supposed to happen below 500 K. A direct recombination

of (NC)Csplit-interstitials with VCCSi pairs, which also may have been created directly by

the implantation process, can be explained in this temperature range. Consequently, the

98 CHAPTER 5. Nitrogen–related Defects

presence of both CSi NC- and VSi NC-pairs in as-implanted samples is conceivable.

The sublattice migration of (NC)Cand (CC)Csplit-interstitials can accordingly be expected

around 1250 K and 1450 K. Long-range diffusion and subsequent recombination of these

split-interstitials with other defects is a possible explanation for the observed annealing

stages at 1050◦C and 1400◦C.

The reduction of the free charge carrier compensation can be imagined like this: At lower

temperatures, mainly recombination processes, the activation of NCby close (CC)Cand

the migration of the created (NC)Cwill happen. Many nitrogen atoms will, however, get

immobile again due to the recombination with carbon vacancies. A further mobilization re-

quires new (CC)Cto approach the NCwhich further requires long-range diffusion of (CC)C.

If this is activated, further substitutionally built in nitrogen can be mobilized, and by the

mechanisms described in Section 5.2.3. CSi NCand CSi (NC)2complexes can be created

from various defects and complexes created before. The electronic properties of these com-

plexes can explain the recovery of the carrier concentration.

A competing process is, however, the formation of VSi (NC)ncomplexes which finally re-

sults in the inactive VSi (NC)4complex.

5.3.3 Entropy effects on nitrogen migration

The most important mechanisms, i. e. the creation of mobile nitrogen split-interstitials and

the sublattice migration processes, have been investigated in view of the effects of vibra-

tional entropy, as was done in case of vacancy migration and aggregation of antisite pairs

in Chapter 4. The influence of the entropy on the sublattice migration of substitutional

NCmight be rather strong, but the activation energy for this exchange process is much too

high to be brought into a realistic order of magnitude by the entropical contribution.

Nevertheless, the entropical contribution to the activation energies of the creation process

for nitrogen split-interstitials and the sublattice migration processes of the carbon- and ni-

trogen split-interstitials, has to be investigated. Our calculations show that the creation of

(NC)Csplit-interstitials is facilitated by ≈0.3 eV at a temperature of 1800 K, while (CC)Cis

stabilized by ≈0.06 eV compared to the resulting (NC)C. The change in energy along the

migration path associated with this mechanism is shown in Fig. 5.13. The activation en-

ergy of the process does not change considerably, even at these elevated temperatures.

The same calculations were performed for the sublattice migration processes of the two

split-interstitials. For the nitrogen split-interstitial (NC)C, calculations yield a negligible

contribution of T·∆S < 0.1 eV for temperatures lower than 2000 K. The activation of

(CC)Csublattice migration, though, is lowered by T·∆S≈0.7 eV at a temperature of

T= 1800 K. At a temperature of T= 1300 K, the barrier is lowered by ≈0.5 eV.

According to the stronger influence of vibrational entropy on the migration of carbon

split-interstitials than on that of nitrogen split-interstitials, the activation of (CC)Cand

(NC)Csublattice migration has to be expected to happen at approximately the same tem-

perature. Consequently, activation of nitrogen atoms that might e. g. have been trapped

by carbon vacancies on their way through the lattice is perpetually possible by (CC)C.

Furthermore, with the calculated entropies, we can now estimate the diffusivities of these

split-interstitials and compare them to those of the vacancies. At T= 1800 K, the cal-

culation of the diffusivities of the carbon and nitrogen split-interstitials yields D((CC)C)

5.4. Implantation with Phosphorus 99

C(NC)C

(CC)

TdS

Energy [eV]

3.5

1.5

2.5

0.5

Figure 5.13: Influence of the vibrational entropy on the creation of (NC)Csplit-interstitials.

= 2.67 ·10−9cm2/s and D((NC)C) = 4.64 ·10−10 cm2/s. For the sublattice migration

of carbon- and silicon-vacancies, we obtained values of D(VC) = 2.49 ·10−15 cm2/s and

D(VSi) = 8.77 ·10−14 cm2/s, as already discussed in Chapter 4. The split-interstitial mi-

gration is in both cases by about four orders of magnitude faster than vacancy migration,

so that the proposed mechanism is obviously very efficient.

5.4 Implantation with Phosphorus

An alternative to nitrogen implantation for the creation of n-type SiC-material is the

implantation with phosphorus (P) ions. With the implantation of phosphorus, a higher

electrical activation in the high concentration range (>1020 cm−3) can be reached [31, 106],

and a co-implantation of nitrogen and phosphorus was found to be advantageous. It is not

yet clear on the atomic scale why this saturation behavior of the free carrier concentration

is not (or possibly at higher concentrations) observed in P- or N-P co-implanted samples.

We have investigated migration mechanisms of the implanted phosphorus ions. For nitro-

gen, our calculations have shown a split-interstitial based mechanism to be very efficient.

We therefore start with the investigation of a similar interstitial based mechanism for

phosphorus. For this aim, we first have to investigate the stable interstitial structures of

phosphorus atoms in the SiC lattice.

Among all the possible interstitial configurations, phosphorus – like nitrogen – turns out to

prefer the split-interstitial position on the carbon site, (PC)C. But also as split-interstitial

on the silicon site, (PSi)Si , the phosphorus atom can be built in, although (PC)Cis favored

about 0.83 eV over for (PSi)Si . The geometries are similar to that of (NC)C, but the phos-

phorus atom induces a large lattice distortion, so that relaxation moves the defect axis to

a tilted position: in both structures, the P-atom is slightly moved out of the symmetry

100 CHAPTER 5. Nitrogen–related Defects

(PSi)

(PC)

C(PC)

Energy difference [eV]

0.5

1.5

2.5

Figure 5.14: Solid red line: Migration of phosphorus split-interstitials using both sublattices: (PC)C−→

(PSi)Si −→ (PC)Cnextsite.Dashed blue line: Direct migration (PC)C−→ (PC)Cnextsite requires a higher

activation energy. The atomic structure in the center shows the geometry of (PSi)Si , belonging to the two

step process.

plane.

The tetrahedrally coordinated interstitial positions are – in contrast to the case of nitrogen

– also stable, but about 7 eV higher than the split-interstitial configurations. Consequently,

a split-interstitial based migration mechanism for phosphorus, similar to that proposed for

nitrogen, can be expected to have the lowest energy of the possible interstitial based mech-

anisms.

A direct migration of (PC)Cfrom one carbon site to the next leads over a symmetric in-

terstitial structure, in which the phosphorus atom is bonded to both carbon atoms that

are part of the initial and the final (PC)Cstructure, and the closest silicon atoms. The

activation energy for this migration mechanism can be activated with 2.9 eV, compare the

dashed line in Fig. 5.14.

Because of the result that (PSi)Si and (PC)Care both stable interstitial configurations,

sublattice migration can also take place on both sublattices by first neighbor jumps with

the P-atom alternating from (PC)Cto (PSi)Si , see the solid line in Fig. 5.14, and further

on to the next carbon site. The jump from the carbon site to the silicon site ((PC)C−→

(PSi)Si ) can be activated with ≈2.2 eV, gaining 1.37 eV which are then needed to activate

the jump from the silicon site to the next carbon site ((PSi)Si −→ (PC)C).

That means, sublattice migration of phosphorus can be activated at slightly lower cost

than the analogous processes for nitrogen (2.5 eV) or carbon (2.9 eV). At elevated tem-

peratures, both described mechanisms (with and without a change of the sublattice) have

to be considered, thus we have found a very efficient diffusion mechanism for phosphorus.

Since the lattice is not at all perfect in the implanted region, the phosphorus atoms will

not very long migrate through the lattice before meeting another defect. As in case of

carbon split-interstitials in Chapter 3 and for nitrogen split-interstitials at the beginning

of this chapter, we have investigated the recombination of phosphorus split-interstitials

with vacancies on both sublattices for various dopant-vacancy distances.

In contrast to nitrogen, we have found that phosphorus fills up both silicon and car-

bon vacancies, and a process analogous to the kick-out mechanism described for nitrogen

(NC)Cinterstitials close to a silicon vacancy could not be found. Instead, simple recombi-

nation processes lead to large energy gains, and an at the first view surprisingly long-range

5.4. Implantation with Phosphorus 101

influence of both kinds of vacancies has been found. Since the energetical difference be-

tween the two phosphorus split-interstitials is rather small, recombination processes of both

(PC)Cand (PSi)Si with VCand VSi have to be considered.

A (PSi)Si split-interstitial being second neighbor to a silicon vacancy has been found to

be instable, as the phosphorus atom immediately and without any energy barrier moves

to the vacant site, resulting in a PSi . If instead a carbon vacancy is third neighbor to

a (PSi)Si , this structure has been calculated to be (meta)stable, but also in this case no

energy barrier has to be overcome, before the P can recombine with the vacancy to PC,

gaining 10.3 eV.

For (PC)Csplit-interstitials, results are similar: a (PC)Cas second neighbor to a VC

recombines with an activation energy of <0.3 eV and an energy gain of 8.7 eV to PC. If

a (PC)Chas approached a silicon vacancy up to a third neighbor distance, they recombine

spontaneously to PSi . Even if (PC)Chas approached a VSi up to a fifth neighbor distance,

only, the activation energy is still only ≈0.9 eV for recombination to PSi . The energy gain

due to this recombination amounts to 14.5 eV.

For a VCCSi pair with the carbon vacancy in a second neighbor distance to a (PC)C, an

activation energy of <1.0 eV has been calculated for the recombination to the PCCSi pair,

even this process accompanied by a large energy gain of 12.5 eV.

From these results, we can conclude that phosphorus can, in contrast to nitrogen, be built

in on both sublattices, as PCor as PSi .

The PCCSi pair has also to be considered important as a donor complex. Though it is about

2.9 eV higher in energy than the stoichiometrically equivalent PSi , its creation is likely by

recombination of a phosphorus interstitial with a VCCSi pair. Since the VSi SiCpair is insta-

ble against recombination to VC, an analogous mechanism for the creation of PSi SiCdoes

not exist. A creation starting from PCby a costly exchange process is unlikely, in particular

since the resulting PSi SiCpair is by 2.2 eV higher in energy.

Recombination of phosphorus split-interstitials with divacancies can lead to phosphorus–

vacancy pairs. These have already been observed by Veinger et al. [107] in 1986. For

nitrogen, two stable positions in a vacancy were found, NCVSi and NSi VC. For phospho-

rus, the situation is similar. Phosphorus prefers the silicon site of a divacancy by 1.1 eV

over the carbon site. In the PCVSi pair, the phosphorus is rather between the two vacancies

than on the carbon site (the three P-Si bonds are elongated by 22 %), leading to a small

barrier of 0.6 eV for the transformation into a PSi VCpair. In the negative charge state, the

PCVSi pair gets stabilized. It is now only 0.5 eV higher in energy than the PSi VCpair, and

the energy barrier separating the two pairs increases to 0.9 eV. The PSi VCpair is Jahn-

Teller instable and slightly more stable if the defect axis is tilted a few percent and two of

the silicon ligands of the vacancy relax to a larger distance between each other, resulting

in C1h-symmetry. The energy difference to the structure with C3v-symmetry is, however,

only ≈0.1 eV, so that at already rather low temperatures, the distortion will probably be

averaged and the pair will be observed in C3v-symmetry.

Its recombination behavior distinguishes phosphorus substantially from nitrogen. As it

simply fills up vacancies, the formation of larger P aggregates, like e. g. a VSi (PC)4or a

VC(PSi )4complex, is implausible. Only the existence of larger vacancy clusters in the im-

planted region could explain the creation of such complexes, but a compensation of the free

charge carriers as efficient as for nitrogen is not possible. If a P approaches a PSi VCpair,

it will fill up the vacancy, resulting in the PCPSi pair. For nitrogen–vacancy pairs, the

102 CHAPTER 5. Nitrogen–related Defects

formation of NSi NCpairs was hindered by the kick-out process or the also more favorable

formation of CSi (NC)2complexes. Furthermore, NSi NCwas calculated to have a slightly

higher energy than its constituents. On the contrary, the pair of two phosphorus atoms

has a binding energy of ≈2.4 eV over for isolated PCand PSi .

If phosphorus is used together with nitrogen in a co-doping procedure, there is, further-

more, the possibility of the formation of P-N pairs. The PCNSi pair has been calculated

to be by ≈5 eV higher in energy than the PSi NCpair. Both pairs can result from a re-

combination with phosphorus split-interstitials with the respective nitrogen–vacancy pairs.

Its lower energy and the more prevalent existence of the VSi NCpair over for the inverse

pair, supports the creation of PSi NCpairs. The PSi NCpair has C3v-symmetry, and a

binding energy of 2.5 eV compared to the isolated constituents PSi and NC. This binding

energy and the large energy gains due to recombination of both nitrogen and phosphorus

split-interstitials suggest that once created aggregates will not dissociate again.

On the other hand, according to the measurements of Laube et al. [106] the electronic be-

havior of the co-implanted sample is determined mainly by the phosphorus atoms, which

rather suggests the formation of separate complexes of P and N with intrinsic defects in

the first line.

The most important conclusions that can be drawn out of the results presented in this

section are that phosphorus as a n-type dopant behaves completely different from nitro-

gen. A very efficient diffusion mechanism based on split-interstitial jumps has been found

for long-range diffusion. The (for SiC) extremely low recombination barriers imply that

at high defect concentrations, as can be expected directly after the implantation process

(compare the discussion at the beginning of this chapter), the sublattice migration of phos-

phorus plays only a secondary role. Phosphorus split-interstitials will directly recombine

with vacancies in a large radius of about five atomic distances already at low annealing

temperatures. Once these processes have finished, lower defect concentrations demand

higher annealing temperatures to activate sublattice migration in order to promote the

phosphorus atoms to substitutional sites. These results indicate that doping with phos-

phorus has many advantages over for doping with nitrogen. Our finding that phosphorus

does not tend to form inactive aggregates suggests, moreover, that the sequence of im-

plantation steps with nitrogen and with phosphorus might make a difference. If already

inactive nitrogen complexes have been formed, phosphorus will not change the situation,

but if a large number of silicon vacancies are already filled up with phosphorus atoms,

the kick-out process of nitrogen split-interstitials and thereby the formation of inactive

VSi (NC)4complexes may be prohibited.

Chapter 6

Summary and Outlook

In this work, the selfconsistent charge density-functional based tight-binding method (SCC-

DFTB) was used to investigate some selected problems of point defects and aggregates

thereof in silicon carbide (SiC). Special importance was attached to mechanisms for the

formation and annealing of defect complexes starting from such defects that are supposed

to be created during ion implantation processes, both intrinsic and n-type dopants.

Antisites are not only among the intrinsic defects with the lowest formation energies but

also among those which received the least attention in the literature. On the experimental

side, this can be understood, since they are only hardly detectable by common methods.

We could, however, show that there is a strong indication that the common DIphoto-

luminescence center is related to the antisite pair CSi SiC. Among a variety of possible

mechanisms for the creation of an antisite pair, two vacancy assisted mechanisms have

been calculated to be energetically most favorable.

The stability of larger aggregates of antisites, especially a completely inverted bilayer as

the limit of two-dimensional antisite aggregation or the three-dimensional ”onions” of an-

tisites has been shown. For the creation of small onion-like aggregates centered around a

silicon antisite, a mechanism based on the vacancy assisted mobility of isolated antisites

has been described. A two-dimensional inversion of the lattice is, however, hindered by

recombination barriers that are smaller than the barriers for creation or for the dissociation

of the mediating vacancy from the complex, such that no complexes with more than two

antisite pairs will be created.

This result stays unchanged if taking into account the entropical contributions to for-

mation and activation energies. Although these contributions have been shown to be of

non-negligible order of magnitude at the high temperatures needed for the activation of

the discussed processes, recombination stays more likely for more than two pairs.

Thus, creation of a completely inverted bilayer which shows a promising quantum well

behavior and offers many possibilities for the design of quantum devices, has to be done

by a special epitaxial method during growth.

The annealing behavior of implanted nitrogen ions has been studied, resulting in a con-

sistent model by which the experimentally observed saturation behavior of the free charge

carrier concentration can be explained. With the help of carbon split-interstitials (CC)C,

nitrogen can be mobilized in form of nitrogen split-interstitials (NC)C. Calculations based

on the migration of these split-interstitials could confirm that nitrogen is preferably built in

on the carbon site by recombination with carbon vacancies. Instead of similarly filling up

103

104 CHAPTER 6. Summary and Outlook

silicon vacancies, resulting in NSi , the (NC)Cperform rather a kick-out mechanism which

leads to either VSi NCor CSi NCpairs. With VSi NCas the first step to the formation of

VSi(NC)4complexes which have been calculated to be electrically inactive, a very efficiently

passivating atomistic model for the observed saturation has been found. The formation of

CSi NCand CSi (NC)2complexes can cause the recovery of charge carriers upon further

implantation and annealing.

A similar, even more efficient, migration mechanism has been found for phosphorus as

alternative n-type dopant. According to our investigations of recombination processes of

phosphorus split-interstitials (PC)Cand (PSi)Si with vacancies, phosphorus can as well be

built in on silicon- as on carbon sites. Low activation energies and high energy gains go

along with these recombination processes, an no aggregation to inactive complexes could

be found, making P the preferred n-type dopant over for N.

In the near future, work on these matters has to be continued, especially in view of a

possible identification of the EPR-signals of nitrogen- and phosphorus-related complexes

with atomistic models – a field where there is today still much confusion in the literature.

For this purpose, it could be shown during this work that a combination of several meth-

ods with different advantages is a successful approach to obtain experimentally relevant

information about a topic.

Up to now, phosphorus and nitrogen have merely been treated separately. A more detailed

investigation of how P influences the tendency of nitrogen to form inactive complexes

during a real co-implantation process is required to find an efficient co-implantation and

annealing procedure for the creation of n-type material with high concentrations of free

charge carriers.

Furthermore, it has to be investigated both experimentally and theoretically, how the

formation of large antisite aggregates can be achieved. Some first calculations have shown

a promising behavior of such structures, especially the two-dimensional inverted bilayer,

as quantum well/wall structures with a large tunneling probability. If one could succeed in

deterministically growing silicon and carbon monolayers hetero- as well as homoepitaxially,

a large variety of quantum devices could be designed for special purposes.

Appendix A

Formation Energies

To determine the relative stability of defects, the formation energy is calculated. It de-

pends on the environmental conditions, i. e. the concentrations of all involved materials,

the temperature and the pressure. The temperature dependence does not only enter the

expression for the Gibbs free enthalpy with the term −T S (see chapter 2) but is also im-

plicitly included in the formation energy via the chemical potential. This dependence is,

though, only significant, if e. g. equilibrium of a solid with a gaseous phase is investigated

(compare the considerations of growth and oxidation of SiC surfaces in Ref. [41]). In the

processes discussed in this work, the implicit temperature dependence of the chemical po-

tentials is negligible as is the pressure dependence.

With the total energies Etot obtained from the supercell calculations, the formation energy

for any defect in SiC has the form

Eform =Etot −nSiµSi −nCµC−X

niµi.(A.1)

Here, nSi and nCdenote the concentrations of silicon- and carbon–atoms, nithe concen-

trations of all impurities iin the crystal. With µSi,µCand µi, the chemical potentials of

silicon, carbon and the impurities are taken into account.

In thermal equilibrium with the SiC crystal, the chemical potentials of silicon µSi and

carbon µCare connected by the relation

µSi +µC=µbulk

SiC ,(A.2)

with help of which one of the chemical potentials can be eliminated from the expression

for the formation energy.

With some further definitions, such as

∆µSi =µSi −µbulk

Si and ∆µC=µC−µbulk

C,(A.3)

denoting the deviation of the actual chemical potentials of silicon and carbon in the SiC-

crystal from their values in the most stable phase of the silicon- and the carbon-crystal,

and the formation enthalpy of SiC

∆HSiC

f=µbulk

SiC −µbulk

Si −µbulk

C,(A.4)

it follows that

∆µSi + ∆µC= ∆HSiC

f.(A.5)

105

106 CHAPTER A. Formation Energies

For intrinsic defects, Eq. A.1 can now be rewritten by substituting the chemical potential

of carbon µC:

Eform =Etot −nSiµSi −nCµC(A.6)

=Etot −nSiµSi −nC·(µbulk

SiC −µSi)

=Etot −(nSi −nC)·µSi −nCµbulk

SiC

=Etot −(nSi −nC)·(∆µSi +µbulk

Si )−nCµbulk

SiC .

Substituting instead the silicon chemical potential µSi leads to the equation:

Eform =Etot −(nC−nSi)·(∆µC+µbulk

C)−nSiµbulk

SiC .(A.7)

The addition of these two equations leads to a symmetric form of the expression:

Eform =Etot −1

2(nSi −nC)h·∆µSi +µbulk

Si −∆µC−µbulk

Ci(A.8)

−1

2(nSi +nC)·µbulk

SiC .

Using Eq. A.4, a new ∆µcan be defined as

∆µ= ∆µSi −1

2∆HSiC

f=−∆µC−1

2∆HSiC

f,(A.9)

leading to the relation

∆µSi −∆µC= 2∆µ . (A.10)

Substituting this difference in the second term of Eq. A.9 and reordering some terms leads

to the final expression for the formation energy:

Eform =Etot −1

2(nSi −nC)µbulk

Si −µbulk

C−1

2(nSi +nC)µbulk

SiC (A.11)

−(nSi −nC)∆µ .

Only the last term contains the energy’s dependence on the environmental conditions,

expressed in form of the chemical potentials of the silicon- and carbon–atoms, united in

the variable ∆µ. All other terms are independent from the environment and just correct

for the different numbers of atoms of either sort in the supercells used for the calculation

of different defects.

Varying ∆µmeans, thus, varying the C/Si ratio, as can experimentally be done during

growth. During growth the conditions can range from carbon–rich, i. e. µC=µbulk

C, to

silicon–rich, i. e. µSi =µbulk

Si . If stepping across these limits the silicon carbide phase is no

longer the most stable phase. If for example the silicon chemical potential becomes higher

than that of the silicon bulk, µSi ≥µbulk

Si , silicon will grow instead of SiC.

Using the definition of ∆µSi and ∆µC, it follows from Eq. A.9 that ∆µcan be varied in

the range 1

2∆HSiC

f≤∆µ≤ −1

2∆HSiC

f.(A.12)

If the crystal contains any impurities i, the formation energy has the same form as in

Eq. A.12, but the term Piniµi, as explained in Eq. A.1, has to be added. The chemical

potential and the range in which it can be varied has to be considered in each special case.

Appendix B

Calculation of the Gibbs Free

Enthalpy

In section 2.6, the calculation of the Gibbs free energy based on the calculation of vi-

brational spectra is discussed. Based on the Einstein model of harmonic oscillators, the

vibrational part of the internal energy

Uvib =

i=1 



ωi

exp(



ωi/kBT)−1+1



ωi(B.1)

was derived. From Eq. B.1, we can, following simple thermodynamics, derive an expression

for the vibrational entropy Svib. Assuming constant volume Vand number of particles N,

we have ∂S

∂U V,N

∂U

∂S V,N

T(S(U, V, N), V, N)=1

T(U, V, N)(B.2)

By solving Eq. B.1 for 1/T we obtain

i=1



ωi·ln Ui+1



ωi

Ui−1



ωi!=∂S

∂U =

i=1 ∂Si

∂Ui(B.3)

Integration of the expressions for the Siover U0

iranging from the zero point energy Ui=



ωi/2 to Uiyields

Si=kB



ωi(Ui+1



ωi)·ln(Ui+1



ωi)−(Ui−1



ωi)·ln(Ui−1



ωi)−



ωi·ln(



ωi),

(B.4)

and finally substituting Uvib from Eq. B.1 and summation over iresults in

Svib =kB

i=1 (



ωi

kBTexp 



ωi

kBT−1−1

−ln 1−exp −



ωi

kBT).(B.5)

For high temperatures T, the term



ωi/kBTbecomes small, and a linear approximation

can be made for the exponential function. This results in the high temperature expressions

Uvib =

i=1 kBT+1



ωi(B.6)

107

108 CHAPTER B. Calculation of the Gibbs Free Enthalpy

for the internal energy Uvib and

Svib =kB

i=1 1−ln 



ωi

kBT (B.7)

for the vibrational entropy. The formation entropy of a defect, Sdefect

vib −Sideal bulk

vib , becomes

in this limit temperature independent.

For the application to defect formation and migration in SiC, this high temperature ap-

proximation can, though, not be used, since it first becomes valid for temperatures clearly

above the Debye temperature of the material. The temperatures used in common annealing

processes, and thus the temperature range connected to most of the processes discussed in

this work, are, however, close to or below the Debye temperature of SiC1. Therefore, the

exact expressions of Equations B.1 and B.5 have been used throughout the work.

As the phonon spectrum should be rather continuous, a broadening of the calculated fre-

quencies may be senseful. Then, the summation over the 3Nfrequencies turns to an

integration. We can rewrite Eq. B.1 as

Uvib =Zωmax

g(ω)



exp(



ω/kBT)−1+1



ωdω (B.8)

with the partition function

g(ω) =

i=1

δ(ω−ωi).(B.9)

The δ-function in this expression can be substituted by a Gauß-broadening, which smoothens

it and achieves a better numerical stability:

g(ω) =

i=1

σ√2π·e−(ω−ωi)2

2σ2.(B.10)

The summation in Eq. B.5 can then be transformed into an integration similarly. In the

case of the applications discussed in this work, a sensitive broadening between 10 cm−1

and 20 cm−1did, though, not change the results for the entropy or the internal energy

notably.

1Literature values range from 1200◦C to 1600◦C [70]

Appendix C

Basic Properties of Strain Fields

In chapter 2 some formulas derived within elasticity theory are needed to calculate the

corrections to the formation entropy. Only these elementary equations are sketched in the

following.

The strained material can be described by a displacement field

~u(~r) = ~r 0−~r , (C.1)

where ~r are the positions in the unstrained, ~r 0in the strained lattice. The distance ~a of

two points in the material changes accordingly by

d~a =~u(~r +~a)−~u(~r),(C.2)

which can in components be linearly expanded to

dai=X

j∂ui

∂xjaj

j1

2∂ui

∂xj

+∂uj

∂xi+1

2∂ui

∂xj−∂uj

∂xi

(εij +πij).(C.3)

The πij only cause a rotation, no displacement, while the εij are the components of the

symmetrical strain tensor ε, which can be shown to have six independent components [77].

Considering the material to be built up of small cubes, the diagonal elements of εrepresent

the length changes of the edges while twice the off diagonal elements represent the change

in angle.

The volume change is, therefore, characterized by the dilatation

δ=δV

V=ε11 +ε22 +ε33 .(C.4)

With Hooke’s Law the stress tensor σis obtained from εas

σ=c·εwith σij =σji







σ11

σ22

σ33

σ23

σ31

σ12













c11 c12 c12 000

c12 c11 c12 000

c12 c12 c11 000

000c44 0 0

0000c44 0

00000c44













ε11

ε22

ε33

2ε23

2ε31

2ε12







(C.5)

109

110 CHAPTER C. Basic Properties of Strain Fields

The bulk modulus Kand the shear modulus γcan be obtained from the components of c

K=1

3(c11 + 2c12) and γ=c44 (C.6)

For an isotropic solid, the relation

c44 =1

2(c11 −c12) (C.7)

holds, furthermore. The Lam´e-constants λand µused in chapter 2 relate to the components

of clike

c11 =λ+ 2µ , c12 =λ , c44 =µ , (C.8)

so that the elements of the stress tensor σcan be written as

σij = 2µ εij +λ δ δij (C.9)

with the dilatation δand the Kronecker delta function δij.

A point defect inhomogeneously strains the crystal lattice. The stress σij changes over a

distance δxj=a(the edge of a small cube) by a(∂σij /∂xj), and, as in equilibrium the

total forces must vanish, we get

∂σij

∂xj

= 0 .(C.10)

Expressing εij in Eq. C.9 by the derivatives of the displacement field ~u(~r) as defined by

Eqns. C.3 yields

µ∇2~u + (λ+µ)∇(∇·~u) = 0 (C.11)

or, since ∇×(∇×~a)≡ ∇(∇~a)−∇2~a holds for any vector ~a,

(λ+ 2µ)∇(∇·~u)−µ∇×(∇×~u).(C.12)

Operating with ∇· on this equation yields (as ∇·(∇×~a)≡0)

∇2(∇·~u) = ∇2δ= 0 .(C.13)

Under the assumptions of isotropy the dilatation δhas, thus, to solve a Laplace equation,

for which the boundary conditions are given by the special problem, as e. g. a force free

surface of the surrounding material.

Appendix D

Summary: Activation Energies

To give a better overview and facilitate the comparison, the values calculated for the ac-

tivation energies and migration energies that are discussed in Chapters 3, 4 and 5 are

summarized in Table D.1 for processes of intrinsic defects only and in Table D.2 for those

processes that involve also nitrogen or phosphorus.

Values are sorted, so that the energy gain listed in the last column is always positive. The

recombination barriers, describing the inverse processes, can be obtained by adding the

energy barriers in the second and the energy gain in the fourth column.

Table D.1: Summary of the calculated activation energies for intrinsic defects.

Initial structure Energy barrier [eV] Final structure Energy gain [eV]

CSi 11.7 CSi next site –

SiC11.6 SiCnext site –

Si-C 10.5 CSi SiC4.5

VC5.8 CSi SiC+ VC3.5

VCCSi 4.7 CSi SiC+ VCCSi 2.0

SiC+ VC4.1 VC+ SiC–

CSi + VSi 2.0 VSi + CSi –

VC4.7 VCnext site –

VSi 4.1 VSi next site –

VSi 1.7 VCCSi 1.8

(CC)C2.9 (CC)Cnext site –

(CC)C0.6 (CC)Cturned –

In cases where varying activation energies were obtained for the different lattice sites or

orientations of the pair defects, only an average value has been given in the table. Some-

times the nomenclature does not clearly describe the mechanism. In that case, the reader

is referred to the respective chapters for more information.

111

112 CHAPTER D. Summary: Activation Energies

Table D.2: Summary of the calculated activation energies for processes that involve dopant atoms.

Initial structure Energy barrier [eV] Final structure Energy gain [eV]

VSi + NC3.5 VSi NC2.0

(CC)C+ NC2.0 (NC)C1.4

(CC)C+ CSi NC1.8 CSi (NC)C1.9

(CC)C+ VSi NC1.9 CSi NC9.7

(CC)C+ VSi (NC)22.4 CSi (NC)27.5

(NC)C2.5 (NC)Cnext site –

(NC)C0.8 (NC)Cturned –

(NC)C+ VSi 2.0 CSi NC10.6

(NC)C+ VCCSi 0.3 CSi NC8.8

(NC)C+ VSi NC1.6 CSi (NC)29.3

(NC)C+ VSi 2.9 VSi NC+ (CC)C1.8

(NC)C+ VSi (NC)32.9 VSi (NC)4+ (CC)C2.0

(NC)C+ VCVSi 0.2 VSi NC7.4

(NC)C+ VSi VC1.0 VCNSi 8.0

(NC)C+ VCVSi (NC)31.0 VSi (NC)48.9

VCNSi 2.5 VSi NC0.8

(PC)C2.9 (PC)C–

(PSi)Si 1.4 (PC)C0.8

(PC)C+ VC2nd neighbor <0.3 PC8.7

(PC)C+ VSi 5nd neighbor <0.9 PSi 14.5

(PC)C+ VCCSi 2nd neighbor <1.0 PCCSi 12.5

Bibliography

[1] BeyArs GmbH, St. Gilgen, Austria, http://www.beyars.com.

[2] Infineon Technologies AG, M¨unchen, http://www.infineon.com, ”thinQ!TM SiC Schot-

tky Diodes”.

[3] T. Wimbauer, B. Meyer, A. Hofstaetter, A. Scharmann, and H. Overhof, Phys. Rev.

B 56, 7384 (1997).

[4] A. Zywietz, J. Furthm¨uller, and F. Bechstedt, Phys. Status Solidi B 210, 13 (1998).

[5] A. Zywietz, J. Furthm¨uller, and F. Bechstedt, Phys. Rev. B 59, 15 166 (1999).

[6] A. Zywietz, J. Furthm¨uller, and F. Bechstedt, Phys. Rev. B 62 (11), 6854 (2000).

[7] L. Torpo, R. M. Nieminen, K. E. Laasonen, and S. P¨oykk¨o, Appl. Phys. Lett. 74, 221

(1999).

[8] P. De´ak, J. Mir´o, A. Gali, L. Udvardi, H. Overhof, N. T. Son, App. Phys. Lett. 75,

2103-2105 (1999 ).

[9] E. Rauls, Th. Lingner, Z. Hajnal, S. Greulich–Weber, Th. Frauenheim, J.-M. Spaeth,

phys. stat. sol. (b), 217/2 (2000), p. R1.

[10] Th. Lingner, U. Gerstmann, E. Rauls, S. Greulich-Weber, J.-M. Spaeth, H. Overhof,

and Th. Frauenheim, Phys. Rev. B, 64, 245212, 2001.

[11] L. Torpo, T.E.M.Staab, and R. Nieminen, Phys. Rev. B 65, 085202 (2002).

[12] A. Mattausch, M. Bockstedte and O. Pankratov, Mat. Sci. Forum 353-356, 323

(2001).

[13] M. Bockstedte, M. Heid, A. Mattausch, and O. Pankratov, Mater. Sci. Forum 389-393,

???(Tsukuba) (2002).

[14] H. Itoh, A. Uedono, T. Ohshima, Y. Aoki, M. Yoshikawa, I. Nashiyama, S. Tanigawa,

H. Okumura, and S. Yoshida, Appl. Phys. A 65, 315 (1997).

[15] Recently it has been realized that the assignment of the early EPR centers to vacancies

was in error.

[16] B. Aradi, A. Gali, P. De´ak, J. E. Lowther, N. T. Son, E. Janz´en, and W. J. Choyke,

Phys. Rev. B, 6 3, 245202 (2001).

[17] A. Kawasuso, H. Itoh, and S. Okada, J. Appl. Phys. 80, 5639 (1996).

[18] N. T. Son, E. S¨orman, W. M. Chen, C. Hallin, O. Kordina, B. Monemar, and E.

Janz´en, Phys. Rev. B 55 (5), 2863 (1997).

113

114 BIBLIOGRAPHY

[19] E. S¨orman, N. T. Son, W. M. Chen, O. Kordina, C. Hallin, and E. Janz´en, Phys. Rev.

B61 (4), 2613 (2000).

[20] M. Wagner, B. Magnusson, W. M. Chen, E. Janz´en, E. S¨orman, C. Hallin, and J. L.

Lindstr¨om, Phys. Rev. B 62 (24), 16555 (2000).

[21] N. T. Son, Mat Sci. Forum, Proceedings of ECSCRM2002, in print.

[22] C. Wang, J. Bernholc, R. F. Davis, Phys. Rev. B 38, 12752 (1988).

[23] L. Torpo, S. P¨oykk¨o, R. M. Nieminen, Phys. Rev. B 57 (1998), p. 6243.

[24] L. Torpo, R. M. Nieminen, Mat. Sci. Eng. B 61-62 (1999), p. 593.

[25] E. Rauls, Z. Hajnal, A. Gali, P. De´ak, and Th. Frauenheim, Mat. Sci. Forum, Vol.

353-356, 435-439 (2001).

[26] F. Gao, W. J. Weber, W. Jian, Phys. Rev. B 63, 214106 (2001).

[27] N. T. Son, P. N. Hai, E. Janz´en, Phys. Rev. Lett. 87, 045502 (2001).

[28] L. Patrick and W. J. Choyke, Phys. Rev. B 5, 3252 (1972).

[29] T.A.G. Eberlein, C.J. Fall, R. Jones, P.R. Briddon, and S. ¨

Oberg, Phys. Rev. B 65,

184108 (2002).

[30] A. Gali, P. De´ak, E. Rauls, N. T. Son, I. G. Ivanov, F. H. C. Carlsson, E.Janz´en, and

W. J. Choyke, Phys. Rev. B, 67, 155203 (2003).

[31] M. Laube, F. Schmid, G. Pensl, and G. Wagner, Mat. Sci. Forum 389-393, 791 (2001).

[32] D. ˚

Aberg, A. Hall´en, P. Pellegrino, and B. G. Svensson, Appl. Phys. Lett. 78, 2908

(2001).

[33] D. ˚

Aberg, A. Hall´en, P. Pellegrino, and B. G. Svensson, Appl. Surf. Sci. 184, 263

(2001).

[34] M.Bockstedte, A. Kley, J. Neugebauer and M. Scheffler, Comput. Phys. Commun.

107 (1997), p. 187.

[35] R. Jones, Philos. Trans. R. Soc. Ser. A 341, 351 (1992).

[36] O. Gunnarsson, O. Jepsen, and O. K. Andersen, Phys. Rev. B 27, 7144 (1983).

[37] U. Gerstmann, H. Overhof, Physica B 308, 561 (2001).

[38] G. L. Harris (Hrsg.), Properties of Silicon Carbide, EMIS Datareview Series 13, IN-

SPEC (1995).

[39] E. Rauls, R. Gutierrez, J. Elsner, and Th. Frauenheim, Sol. State Comm. 111, 459

(1999).

[40] E. Rauls, Z. Hajnal, A. Gali, P. De´ak, and Th. Frauenheim, Mat. Sci. Forum, Vol.

353-356, 435-439 (2000).

[41] E. Rauls, Diplomarbeit, Universit¨at Paderborn, (1999).

[42] E. Rauls, Z. Hajnal, P. De´ak, and Th. Frauenheim, Phys. Rev. B 64, 245323 (2001).

BIBLIOGRAPHY 115

[43] M. Lannoo and J. C. Bourgoin, Solid State Commun. 32, 913 (1979).

[44] D. Porezag, Th. Frauenheim, T. K¨ohler, G. Seifert, R. Kaschner, Phys. Rev. B 52,

(1995) p. 11837.

[45] Th. Frauenheim, G. Seifert, M. Elstner, Z. Hajnal, G. Jungnickel, D. Porezag, S.

Suhai, R. Scholz, phys. stat. sol. (b), 217, 41 (2000).

[46] M. Sternberg, Dissertation, Universit¨at Paderborn, (2001).

[47] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai,

G. Seifert, Phys. Rev. B 58 (1998) p. 7260.

[48] K. Ohno, K. Esfarjani, and Y. Kawazoe, Computational Materials Science, Springer

Verlag, Berlin (1999).

[49] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1965).

[50] W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965).

[51] D.M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).

[52] J.P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 (1986).

[53] J.P. Perdew, K. Burke, and M. Enzerhof, Phys. Rev. Lett. 77, 3865 (1996).

[54] F. Beeler, O.K.Andersen, and M. Scheffler, Phys. Rev. B 41, 1603 (1990).

[55] G.A. Baraff, E.O.Kane, and M. Schl¨uter, Phys. Rev. B 21, 5662 (1980).

[56] F. Jensen, Introduction to Computatinal Chemistry, John Wiley & Sons, Chichester,

New York (2001).

[57] W.M.C. Foulkes and R. Haydock, Phys. Rev. B 39, 12520 (1989).

[58] J.C. Slater and G.F. Koster, Phys. Rev. B 94, 1498 (1954).

[59] H. Eschrig, The Optimized LCAO Method and Electronic Structure of Extended Sys-

tems, Akademie-Verlag, Berlin (1988).

[60] J. Elsner, Dissertation, Universit¨at Paderborn, (1999).

[61] U. Gerstmann, E. Rauls, H. Overhof, and Th. Frauenheim, Phys. Rev. B 65, 195201

(2002).

[62] U. Gerstmann, E. Rauls, Th. Frauenheim, and H. Overhof, Proceedings of EC-

SCRM2002, Physica B, in print, (Oct. 2002).

[63] U. Gerstmann, E. Rauls, Th. Frauenheim, and H. Overhof, Phys. Rev. B 67, in press,

(2003).

[64] M. J. Puska, S. P¨oykk¨o, M. Pesola, and R. M. Nieminen, Phys. Rev. B 58, 1318

(1998).

[65] H. J. Monkhorst and J. K. Pack, Phys. Rev. B 13, 5188 (1976).

[66] M. Kaukonen, T. Frauenheim et al., Phys. Rev. B 57, 9965 (1998).

[67] G. T. Barkema and N. Mousseau, Phys. Rev. Lett. 77, (1996) p. 4358.

116 BIBLIOGRAPHY

[68] A. Gali, private communication (2003).

[69] A. Mattausch, private communication (2003).

[70] Landolt-B¨ornstein, New Series III/41 A1α, 2001 and references therein.

[71] Ch. Kittel Physik der W¨arme, R. Oldenbourg Verlag, (1973).

[72] H. Kuchling, Taschenbuch der Physik, (1986).

[73] R.D. Hatcher, R. Zeller, and P.H. Dederichs, Phys. Rev. B 19, 5083 (1979).

[74] J. R. Fern´andez, A.M. Monti, and R. C. Pasianot, Phys. Stat. Sol.(b) 219, 245 (2000).

[75] Y. Mishin, M.R. Sorensen, and A. F. Voter, Phil. Mag. A 81, 2591 (2001).

[76] L. D. Landau, E. M. Lifschitz, Lehrbuch der theoretischen Physik, Band VII, Elas-

tizit¨atstheorie, Akademie-Verlag Berlin (1975).

[77] C. P. Flynn, Point defects and diffusion, Clarendon Press Oxford (1972).

[78] Cox, J. D., Wagman, D. D., and Medvedev, V. A., CODATA Key Values for Ther-

modynamics, Hemisphere Publishing Corp., New York, (1989).

[79] M. Lannoo and G. Allan, Phys. Rev. B 25 (6), 4089 (1982).

[80] M. Scheffler, J. Dabrowski, Phil. Mag. (1988).

[81] D. Morgan, A. van de Walle, G. Ceder, J.D. Althoff, and D. de Fontaine, Modelling

Simul. Mater. Sci. Eng. 8, 295 (2000).

[82] U. Gerstmann, E. Rauls, H. Overhof, Th. Frauenheim, Physica B, 308-310, 497-501,

2001.

[83] Th. Lingner, S. Greulich-Weber, J.-M. Spaeth, U. Gerstmann, E. Rauls, and H. Over-

hof, Physica B 308, 625 (2001).

[84] H. Itoh, M. Yoshikawa, I. Nashiyama, S. Misawa, H. Okumura, and S. Yoshida, IEEE

Trans. Nucl. Sci. 37, 1732 (1990).

[85] V. S. Vainer and V. A. Il’in, Sov. Phys. Solid State 23, 2126 (1981)

[86] Th. Lingner, Dissertation, Universit¨at Paderborn, (2001).

[87] N. T. Son, P. N. Hai, M. Wagner, W. M. Chen, A. Ellison, C. Hallin, B. Monemar, E.

Janz´en, Semicond. Sci. Technol.14 , 1141 (1999).

[88] N. T. Son, private communication (2003).

[89] L. Torpo, M. Marlo, T. E. M. Staab, and R.M. Nieminen, J. Phys. C 13, 6203 (2001).

[90] K. C. Pandey, Phys. Rev. Lett. 57 (18) (1986), p. 2287.

[91] E. Rauls, A. Gali, P. De´ak, and Th. Frauenheim, in Proceedings Vol., Mat. Sci. Forum

(ECSCRM2002), in print.

[92] J. J. Sumakeris, L. B. Rowland, R. S. Kern, S. Tanaka, and R. F. Davis. Thin Solid

Films, 225:219, 1993.

BIBLIOGRAPHY 117

[93] D. G. J. Sutherland, L. J. Terminello, J. A. Carlisle, I. Jimenez, Himpsel F. J., K. M.

Baines, D. K. Shuh, and W. M. Tong. J. Appl. Phys., 82:3567, 1997.

[94] T. Fuyuki, T. Yoshinobu, and H. Matsunami. Thin Solid Films, 225:225, 1993.

[95] H. Nagasawa and Y. Yamaguchi. Thin Solid Films, 225:230, 1993.

[96] S. Hara, T. Meguro, Y. Aoyagi, M. Kawai, S. Misawa, E. Sakuma, and S. Yoshida.

Thin Solid Films, 225:240, 1993.

[97] T. Lohner, Nguyen Quoc Khanh, private communication.

[98] A. Ellison, B. Magnusson, C. Hemmingsson, W. Magnusson, T. Iakimov, L. Storasta,

A. Henry, N. Henelius, and E. Janz´en, Mat. Res. Soc. Symp. 640, (2001), p. H1.2

[99] V. S. Ballandovich and G. N. Violina, Cryst. Lattice Defects Amorphous Mater. 13,

189 (1987).

[100] G. Davies, S.C. Lawson, A. T. Collins, A. Mainwood, and S. J. Sharp, Phys. Rev. B

46, 13157 (1992).

[101] A. Mainwood, Phys. Rev. B 49, 7934 (1994).

[102] S. Greulich-Weber, phys. stat. sol. (a) 162, 95 (1997).

[103] S. Greulich-Weber, phys. stat. sol. (b) 210, 415 (1998).

[104] M. Laube, Erlangen, private communication (2002).

[105] T. Evans, The properties of natural ans synthetic diamond, ed. J.E. Field (Academic

Press: London), 259 (1992).

[106] M. Laube, F. Schmid, G. Pensl, G. Wagner, M. Linnarsson, M. Maier, J. Appl. Phys.

92, 549 (2002).

[107] A. I. Veinger, A. G. Zabrodskii, G. A. Lomakina, and E. N. Mokhov, Sov. Phys.

Solid State 28, 917 (1986).

118 BIBLIOGRAPHY

Acknowledgment

First of all I would like to thank my thesis advisor Prof. Th. Frauenheim for his support

and for allowing me to pursue my own research interests. For several advises and for

writing one of the certificates for this work, I am grateful to Prof. H. Overhof.

For his interest in my work on antisites and also on surface structures, which are, though,

not included in this thesis, I want to express special thanks to Prof. Peter De´ak, TU

Budapest, — also for his invitations to Budapest which I followed a few times. Our

intensive discussions during his regular stays at Paderborn as well as during my stays at

Budapest always helped keeping an eye on the correlation of our theoretical studies to

experiment.

He was also the one who initiated my stay in the group of E. Janz´en at Link¨oping, Swe-

den, during last spring. For interesting discussions about e. g. the annealing behavior of

vacancies or the origin of the DIcenter during this time and afterwards, I would like to

thank Nguyen Tien Son, especially for giving me some new insight on these topics from

the point of view of an experimentalist.

For providing me with reference calculations with the FHI-code or AIMPRO I would like to

thank Alexander Mattausch, Universit¨at Erlangen, who calculated the vibrational spectra

cited in this work with the FHI-code, and ´

Ad´am Gali, TU Budapest, who did the same

with AIMPRO and performed all other calculations of selected defect structures with the

FHI-code.

For several articles we wrote together and for sharing his experience in defect physics with

me, I am especially grateful to Uwe Gerstmann. He has also performed all LMTO-ASA

calculations cited in this work, which e. g. helped to understand the origin of the differences

between the formation and migration energies of vacancies calculated within SCC-DFTB

and ab initio LDA methods.

Furthermore, I want to thank all members of the group of Th. Frauenheim and all other

”inhabitants” of the N3-floor that have not been named yet. This includes Zolt´an Hajnal,

who provided me with the parameters for nitrogen and phosphorus, and the system ad-

ministration, in the first time in person of Michael Sternberg, and in the last two years in

person of Christof K¨ohler and Peter K¨onig, for keeping our machines in order and, thus,

providing us with the necessary computing facilities.

Last but not least, my warmest thanks to my mother and to Bobby for their patience

during the last years!

119