ABSTRACT
Title of dissertation: DFT AND RELATED MODELING OF
POST-SILICON VALENCE 4 MATERIALS:
SiC AND Ge
Christopher Darmody, 2020
Dissertation directed by: Dr. Neil Goldsman
Department of Electrical and Computer Engineering
Though silicon (Si) is in many ways the material of choice for many electronic
applications due in part to its mature processing technology, its intrinsic properties
are not always suited for every challenge. Specialized high power and high temper-
ature devices benefit from using semiconductors with a larger band-gap and higher
thermal conductivity such as silicon carbide (SiC). Additionally, the 1.1eV bandgap
of Si makes it unable to effectively absorb infrared photons so a material with a
smaller bandgap, like germanium (Ge), is more suited to the task.
Currently SiC power transistors are commercially available but suffer from
poor channel mobility due to interface roughness which limits their performance.
To predict the maximum theoretically achievable mobility for different crystallo-
graphic interfaces I developed a novel technique for extracting an atomic-roughness
scattering rate from an arbitrary atomic surface. The term atomic-roughness here
means an interface purely due to the variation of atom species and position without
the presence of a crystallographic miscut due to epitaxial growth considerations.
I used Density Functional Theory (DFT) to obtain a perturbation potential from
which I can calculate a scattering rate. This scattering rate can then be used in a
Monte Carlo simulation to predict mobility for a given field configuration.
In addition to SiC?s low channel mobility, SiC p-type dopant species also ex-
hibit an abnormally large ionization energy compared to its n-type dopants and
to the primary dopants in many other semiconductors. This fact can cause is-
sues such as unexpectedly high resistance regions at lower operating temperatures
- causing the need to dope at significantly higher concentration. To characterize
the incomplete ionization fraction p/NA, I first gathered nearly all existing pub-
lished data on the ionization energy of aluminum (Al) in 4H-SiC and created an
empirical concentration-dependent model of this function. Then I put together a
physics-based model of the entire acceptor and valence band system and used my
concentration-dependent ionization energy as an input to predict p/NA. I verify my
physics-based model result against a separate experimental dataset derived from
nearly-exhaustive literature measurements of Hall mobility and resistivity. Finally,
I transform fully temperature-dependent result of p/NA from a complex numerical
computation to a more easily implementable parameterized function with the use
of a genetic algorithm.
The remaining part of my work was performed on Germanium which has
interesting application in short-wave infrared imaging due its 0.66eV indirect and
0.85 eV direct bandgaps, which corresponds closely to the peak illumination of
the ?night glow? at 0.75 eV. Optical devices greatly benefit from direct gap band
structures to increase photon absorption and emission efficiency. Though Ge is an
indirect gap material, it can be alloyed with a direct gap material, namely tin (Sn),
to transition it to a direct gap material at a certain molar fraction. Through DFT
calculations I investigate the nature of this transition and determine theoretically
the minimum molar fraction needed to achieve a direct bandgap.
DFT AND RELATED MODELING OF
POST-SILICON VALENCE 4 MATERIALS:
SiC AND Ge
by
Christopher Darmody
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2020
Advisory Committee:
Professor Neil Goldsman, Chair/Advisor
Professor Aristos Christou
Professor Kevin Daniels
Professor Romel Gomez
Professor Agisilaos Iliadis
? Copyright by
Christopher Darmody
2020
Dedication
To my loving parents, James and Jacquie.
ii
Acknowledgments
I would first and foremost like to express my eternal gratitude to my advisor Dr.
Neil Goldsman for his years of support, encouragement, and insightful discussions.
Without your guidance I would not be where I am today.
Thank you to my thesis defense committee: Dr. Kevin Daniels, Dr. Romel
Gomez, Dr. Agis Iliadis, and Dr. Aristos Christou for taking time away from your
busy schedules to participate in this important step in my academic life. I would
also like to thank Dr. John Melngailis for introducing me to the world of device
physics which has become the basis of my graduate studies.
I am also grateful to have had such a great group of colleges and friends who
have made my graduate school experience such a joy - Ittai Baum, Yumeng Cui,
Xiyang Xiao, Alex Mazzoni, Franklin Nouketcha, Dev Ettisserry, Peiwen He, and
Eric Krokos.
Lastly, I would like to thank my family for their loving support. Thanks Mom
and Dad for always believing in me.
iii
Table of Contents
List of Tables vii
List of Figures viii
List of Abbreviations xiii
1 Introduction 1
1.1 Post Silicon Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Overview of Group IV(a) Semiconductors . . . . . . . . . . . . . . . . 4
1.2.1 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.2 Bonding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Silicon Carbide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.1 Background and Applications . . . . . . . . . . . . . . . . . . 10
1.3.2 Crystal Structure and Properties . . . . . . . . . . . . . . . . 12
1.3.3 Issues in SiC Power Devices . . . . . . . . . . . . . . . . . . . 17
1.4 Germanium Photonics: Benefits and Challenges . . . . . . . . . . . . 17
1.5 Research Accomplishments . . . . . . . . . . . . . . . . . . . . . . . . 19
1.6 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2 DFT and Electron Transport Modeling and Theoretical Framework 23
2.1 Many-Body Schrodinger Equation . . . . . . . . . . . . . . . . . . . . 23
2.1.1 Born-Oppenheimer Approximation . . . . . . . . . . . . . . . 24
2.2 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.1 Origins of DFT . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.2 Kohn-Sham Method . . . . . . . . . . . . . . . . . . . . . . . 28
2.3 DFT Calculation Results . . . . . . . . . . . . . . . . . . . . . . . . . 31
3 Atomic-Level Electron Transport 34
3.1 Perturbation Scattering Theory . . . . . . . . . . . . . . . . . . . . . 34
3.2 Scattering and Transport Model . . . . . . . . . . . . . . . . . . . . . 40
3.3 Scattering Rate Calculations . . . . . . . . . . . . . . . . . . . . . . . 45
3.3.1 Acoustic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3.2 Non-Polar Optical Scattering . . . . . . . . . . . . . . . . . . 47
iv
3.3.3 Polar Optical Scattering . . . . . . . . . . . . . . . . . . . . . 48
3.4 Novel Interface Atomic Roughness Scattering . . . . . . . . . . . . . . 49
3.4.1 Existing Surface Roughness Model . . . . . . . . . . . . . . . 49
3.4.2 My Perturbation and DFT-Based Model: Calculation and
Scattering Rate Results . . . . . . . . . . . . . . . . . . . . . 50
3.4.3 Comparison of Intrinsic and Extrinsic Mobilities . . . . . . . . 56
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4 Incomplete Ionization in p-Type SiC 67
4.1 Background and Introduction . . . . . . . . . . . . . . . . . . . . . . 67
4.2 Standard Incomplete Ionization Formulation . . . . . . . . . . . . . . 71
4.3 Doping-Dependent Ionization Energy . . . . . . . . . . . . . . . . . . 83
4.4 Acceptor Density of States . . . . . . . . . . . . . . . . . . . . . . . . 86
4.5 Disorder Effects and Band Tailing . . . . . . . . . . . . . . . . . . . . 88
4.6 Theoretical Incomplete Ionization Model . . . . . . . . . . . . . . . . 91
4.7 Experimental Incomplete Ionization Verification . . . . . . . . . . . . 97
4.8 Temperature Dependence . . . . . . . . . . . . . . . . . . . . . . . . 106
4.9 Results and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5 Transport Model Validation and Genetic Algorithm Application 112
5.1 The Monte Carlo Method . . . . . . . . . . . . . . . . . . . . . . . . 112
5.2 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . 119
5.3 Genetic Algorithm For Parameter Extraction . . . . . . . . . . . . . . 121
5.3.1 Method Overview . . . . . . . . . . . . . . . . . . . . . . . . . 121
5.3.2 Breeding Process Details . . . . . . . . . . . . . . . . . . . . . 126
5.3.3 Genetic Algorithm Applied to Doping and Temperature-Dependent
Hall Mobility . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
6 Germanium Modeling 135
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
6.2 DFT-Based Analysis of Ge . . . . . . . . . . . . . . . . . . . . . . . . 139
6.3 DFT-Based Analysis of GeSn Alloys . . . . . . . . . . . . . . . . . . 145
A DFT Evolution and History 150
A.1 Hartree-Fock Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
A.2 Thomas-Fermi Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
A.3 Hohenberg-Kohn Theory . . . . . . . . . . . . . . . . . . . . . . . . . 153
B Semiconductor Physics 156
B.1 Low Doping Limit of Hole Concentration . . . . . . . . . . . . . . . . 156
B.1.1 Conductivity Mobility Parameterization . . . . . . . . . . . . 158
B.2 Modified Fermi-Dirac for Dopant States . . . . . . . . . . . . . . . . 159
B.3 Valence Band Density of States . . . . . . . . . . . . . . . . . . . . . 165
B.3.1 Dopant DoS Spreading . . . . . . . . . . . . . . . . . . . . . . 166
B.4 Compensated Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 168
v
C Standard Component Mobility Formulations 170
C.1 Bulk Mobility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
C.2 Surface Phonon Mobility . . . . . . . . . . . . . . . . . . . . . . . . . 170
C.3 Combined Empirical Surface Roughness Mobility . . . . . . . . . . . 171
C.4 Coulomb Mobility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172
Bibliography 174
vi
List of Tables
1.1 Material Properties of Column IV . . . . . . . . . . . . . . . . . . . . 6
1.2 SiC/Si Electrical Properties . . . . . . . . . . . . . . . . . . . . . . . 14
2.1 SiC DFT Bandgap Comparison to Experiment . . . . . . . . . . . . . 33
2.2 DFT Effective Mass Comparison to Experiment . . . . . . . . . . . . 33
3.1 Phonon and SiC Monte Carlo Properties . . . . . . . . . . . . . . . . 45
3.2 Inelastic Acoustic Phonon Integration Limits . . . . . . . . . . . . . . 47
3.3 Important component mobilities at the SiC/SiO2 interface for T=300K 57
4.1 Parameters characterizing Al-doped 4H-SiC . . . . . . . . . . . . . . 100
5.1 Optimized parameters for empirical Hall mobility in Al-doped 4H-SiC
obtained using a genetic algorithm . . . . . . . . . . . . . . . . . . . 134
6.1 Ge DFT Bandgap Comparison to Experiment . . . . . . . . . . . . . 144
C.1 Extrinsic Coulomb Scattering Mobility Parameters [1, 2] . . . . . . . . 173
vii
List of Figures
1.1 Plot of popular semiconductor elements and compounds relating atomic
bond length to energy bandgap. The scatter data shows a general
trend of smaller bond distances correlate to larger bandgaps. (Figure
from [3]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Metalloid section of the periodic table where most of the important
semiconductor elements reside. Elements in each column have valence
3, 4, and 5 from left to right, and increase in radius from top to bottom. 5
1.3 Cartoon visualization of the four sp3 orbital lobes. . . . . . . . . . . . 7
1.4 Actual form of sp3 tetrahedron (level set of the wavefunction) showing
overlaid orbitals. Yellow and blue colors denote opposite signs the
wavefunction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.5 Molecular wavefunctions for a linear chain of 6 atoms. Degenerate en-
ergies split as interatomic distance decreases due to the interaction of
neighboring atomic orbitals. Light and dark circles represent positive
and negative orbital lobes respectively. By Tem5psu (Own work) [CC
BY-SA 4.0 (http://creativecommons.org/licenses/by-sa/4.0)], via Wiki-
media Commons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 Splitting and occupancy of 2s (3s,4s) and 2p (3p,4p) states of car-
bon (silicon, germanium) as a function of inter-atomic spacing. The
distance between the top of the filled band and the bottom of the
unfilled band shows formation of various bandgap values at experi-
mental inter-atomic distances. Inset images depict spatial visualiza-
tions of the kinds of orbitals found in each band. Orbitals shown:
2s orbital (bottom right), 2p orbital (top right), overlaid hybrid sp3
orbitals (center), sp3 sigma bonding orbital (bottom left), and sigma?
antibonding orbital (top left). Color represents sign of wavefunction. [4] 11
1.7 SiC polytypes depicting only key atoms forming the bonding chain
in the stacking direction. Primitive cells are shown for the hexagonal
polytypes and a non-primitive hexagonal cell is used for 3C. Polytypes
differ only in the stacking order and periodicity of the close-packed
atomic planes. For the cubic variant, this is the {111} plane and for
the hexagonal variants, the {0001} plane. (Figure from [5]). . . . . . 13
viii
1.8 Wigner-Seitz cell of hexagonal real-space lattice structure and rele-
vant crystal planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.9 4H-SiC lattice depicting the stacking arrangement and crystal faces . 16
2.1 Theoretical band structure of 4H-SiC calculated with DFT. . . . . . . 32
2.2 Constant energy ellipses for the lowest conduction band minimum at
the M point, calculated using DFT. . . . . . . . . . . . . . . . . . . . 33
3.1 Classification of scattering mechanisms in semiconductors. Taken
from [6] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 Illustration of initial vector k scattering to a different state k? via the
absorption of a phonon with wavevector q. Note: ? is in general not
in the same plane as ?k, which is simply used to define the position
of k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.3 Brillouin Zone of a hexagonal lattice. . . . . . . . . . . . . . . . . . . 43
3.4 Simplified band illustration for 4H-SiC with associated energies from
literature [7] (Figure from [8]). . . . . . . . . . . . . . . . . . . . . . . 44
3.5 4H-SiC/SiO2 interface supercell created by [9], used to extract calcu-
lated potential. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.6 Interface potential taken at the semiconductor-insulator interface.
The blue (well) areas depicted are within the oxygen atoms. . . . . . 61
3.7 Interface perturbation potential taken at the semiconductor-insulator
interface. The blue (well) areas are the oxygen atoms, and the yel-
low (peak) areas are where the carbon sites are located if the bulk
continued beyond the interface. . . . . . . . . . . . . . . . . . . . . . 62
3.8 Real-space band extrema cartoon depicting the well structure which
forms at the SiC-SiO2 interface when a channel is formed. The sub-
bands are enumerated 1, 2, 3... and exist as states found only within
the well near the surface. On the right, the triangular approximation
to the well is made and the corresponding approximate wavefunctions
are shown which increase in number of nodes with energy. . . . . . . 63
3.9 Calculated transition rates for electrons due to atomic-roughness scat-
tering with DFT-extracted interface perturbation potential. Initial
state has 1eV of kinetic energy and is in the second subband of a
triangular well with surface field of 1MV/cm. Scattering within the
second subband is shown in red, scattering to the third subband in
orange, and to the first subband in blue. . . . . . . . . . . . . . . . . 64
3.10 2D atomic roughness scattering rates for electrons in the 1st (blue),
2nd (red), and 3rd (orange) subbands. The x axis shows total energy
of the electron, so each curve begins at the corresponding subband
energy for the given field. Insets show the relative shape of the trian-
gular well that each plot represents. Bottom right plot contains both
energy and field variation. . . . . . . . . . . . . . . . . . . . . . . . . 65
3.11 Field-dependent scattering within the 1st subband for electrons with
various kinetic energies. . . . . . . . . . . . . . . . . . . . . . . . . . 66
ix
4.1 Doping levels in 4H-SiC. . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.2 Doping levels in Ge, Si, and GaAs. . . . . . . . . . . . . . . . . . . . 70
4.3 Hole concentration calculated using the standard incomplete ioniza-
tion model for various NA doping. Dashed lines show the quadratic
solution which ignores ni and the solid lines show the cubic solution
which includes it. Calculated for 4H-SiC with ?EA = 200meV. . . . . 77
4.4 p/NA calculated for different doping concentrations as a function of
temperature for 4H-SiC with ?EA = 200meV. . . . . . . . . . . . . . 78
4.5 Coulombic potential of impurity atom core. the solid line is the bare
hydrogenic 1/r potential and the dashed line is the screened potential
under the effects of mobile carriers redistributing themselves. . . . . . 79
4.6 Local arrangement of atoms leading to inequivalent lattice sites in
4H-SiC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.7 Inequivalent sites (h=hexagonal, k=cubic) labeled in SiC polytypes.
The black circles represent carbon atoms and the white circles repre-
sent Si atoms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.8 Doping density-dependent fit to experimental acceptor ionization en-
ergy for 4H-SiC doped with Al. The blue line is the Pearson-Bardeen
model and the black line is the logistical parameterization. . . . . . . 85
4.9 P-doped semiconductor density of states diagram. Left: Shows the
position-dependent variation of the local band structure (negative of
potential). Right: Dashed blue lines show the standard valence band
density of states and the rectangular acceptor band density of states.
After the inclusion of disorder effects, the resulting density of states
with band tailing is shown with solid blue lines. Note: Image is not
to scale - depending on the doping concentration, the doping band
will be significantly closer to and may even overlap the valence band
density of states tail. The doping band width B as well as its height
will also vary considerably depending on the doping density. . . . . . 92
4.10 Evolution of the density of states under various conditions: Row
a) Acceptor density of states starts as a single energy level at low
doping, spreads due to wavefunction overlap at higher concentrations,
then smears due to disorder Row b) Valence band density of states
smearing due to disorder effects row c) Combined picture of the
valence band and acceptor density of states . . . . . . . . . . . . . . 93
4.11 Hall mobility fit to experimental measurements of 4H-SiC doped with
Al . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.12 Conductivity mobility derived from experimental measurements of
resistivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.13 Resistivity fit to experimental measurements of 4H-SiC doped with
Al. Dashed black line comes from using the predicted resistivity from
Heera [10] using our values for mobility and ionization energy. . . . . 103
x
4.14 Incomplete ionization ratio p/NA for 4H-SiC doped with Al at T =
300K. Experimental data fit (?Cond(NA)/?Hall(NA)) is shown in solid
blue. Our theoretical model, calculated with Equation 4.51 is shown
in dashed red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.15 Incomplete ionization ratio p/NA calculated using our Model at ele-
vated temperatures T = 300K ? 800K (solid line) compared to our
empirical parameterization (dashed). . . . . . . . . . . . . . . . . . . 109
4.16 Hole concentration in 4H-SiC doped with Al calculated using our
Model at elevated temperatures T = 300K ? 800K. . . . . . . . . . . 110
5.1 Flowchart of Monte Carlo Algorithm for a single field configuration
using 4 random numbers r1, r2, r3, and r4 . . . . . . . . . . . . . . . . 114
5.2 Scattering rates for electrons with different total energy. Rates as-
sociated with phonon emission are shown as dashed lines and solid
lines are used for absorption. The 2D Atomic Roughness scattering
rates out of the first three subbands are also shown for a surface field
of 2MV/cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
5.3 Average electron velocity (in the direction of the applied lateral field)
plotted against field strength. Monte Carlo simulation results are
solid lines and data-fitted curve from [11] is shown with the dashed
line. Curves shown are for different temperatures: Blue 300K, Red
600K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.4 Plot of average electron energy for given lateral field strength. Zero
energy corresponds to electrons at the conduction band minimum.
Only M-Valley scattering is considered in this simulation. . . . . . . . 119
5.5 Distributions of wavevector component in the direction of applied
lateral field (top) and electron energy (bottom) for various lateral
field magnitudes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
5.6 absolute field mobility for electrons plotted against applied lateral
field, taken as ?(F )/F . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
5.7 Differential electron field mobility as a function of field, extracted
from Monte Carlo simulation results. . . . . . . . . . . . . . . . . . . 121
5.8 Flowchart of the basic Genetic Algorithm process. . . . . . . . . . . . 123
5.9 Operations performed to create a new generation of individuals in a
genetic algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.10 Random binary bit array used to generate an offspring from two par-
ent individuals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.11 Random mutation applied after the crossover to finally create the
offspring. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.12 Example bimodal probability density function which could be used
to pick random gene mixing fractions. . . . . . . . . . . . . . . . . . . 129
5.13 Child created using linear combinations of its parents? genes. . . . . . 129
5.14 Fitness across generations for different runs of the algorithm. . . . . . 131
5.15 Genetic Algorithm fit of the doping and temperature-dependent Hall
mobility of Al-doped 4H-SiC. . . . . . . . . . . . . . . . . . . . . . . 133
xi
6.1 Absorption coefficient for various semiconductor materials. (Plot
from [12]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
6.2 Atomic supercells of Ge1?xSnx used in DFT calculations for 12.5%,
6.25%, 3.125% Sn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
6.3 Brillouin Zone of an FCC lattice. . . . . . . . . . . . . . . . . . . . . 141
6.4 Indirect gap energy at the L point varying linearly with hybrid func-
tional (PBE0) mixing fraction . . . . . . . . . . . . . . . . . . . . . . 142
6.5 Gap Difference (E?-EL) varying linearly with lattice constant. Liter-
ature: 5.646A? Our Work: 5.573A? ( 1.3% difference) . . . . . . . . . . 143
6.6 DFT calculated Ge band structure for 2 atom cell using the PBE0
hybrid functional. Dashed box shows the same section of the band
structure given in the literature [8]. The key direct and indirect gaps
are shown with red arrows. . . . . . . . . . . . . . . . . . . . . . . . . 144
6.7 Ge band structure from [8]. The key direct and indirect gaps are
shown with red arrows. . . . . . . . . . . . . . . . . . . . . . . . . . . 145
6.8 E? and EL calculated for different fractions of Sn . . . . . . . . . . . 149
xii
List of Abbreviations
? Mobility
? Velocity
? Mass Density
?(E) Density of States (DoS)
n Mobile Electron Density
p Mobile Hole Density
ND Donor Concentration
N+D Ionized Donor Concentration
N0D Unionized Donor Concentration
NA Acceptor Concentration
N?A Ionized Acceptor Concentration
N0A Unionized Acceptor Concentration
EF Fermi Energy
? Single-Electron Wavefunction
? All-Electron Wavefunction
?k Single-Electron Time-Independent Wavefunction
?k Single-Electron Temporally and Spatially Dependent Wavefunction
Ai[z] Airy Function
CMOS Complimentary Metal-Oxide-Semiconductor
DD Drift-Diffusion
DFT Density Functional Theory
DMOS Doubly-diffused Metal-Oxide-Semiconductor (Transistor)
FET Field Effect Transistor
GaN Gallium Nitride
Ge Germanium
MC Monte Carlo
MOSFET Metal-Oxide-Semiconductor Field Effect Transistor
SiC Silicon Carbide
SiO2 Silicon Dioxide
Sn Tin
xiii
Chapter 1: Introduction
1.1 Post Silicon Materials
As silicon nears the end of its road map, the semiconductor industry is looking
to other materials to continue technological scaling. Despite silicon?s many good
properties, other less mature materials exist which exhibit better performance in
certain regimes. In particular, for high speed integrated circuits electron mobility
is the key figure-of-merit which dictates how fast electrons can move within the
material. Though the mobility of Si is fairly high, it is surpassed by materials like
GaAs and Ge. High mobility materials also facilitate device scaling by allowing
physically smaller devices to maintain reasonable ON currents in transistors.
Additionally, Si has a relatively small bandgap among semiconductors regu-
larly researched (Figure 1.1 [3]). Small bandgap materials perform poorly in high
temperature and high power applications due to their large intrinsic carrier concen-
trations which can overwhelm doped carrier distributions and small critical fields
which cause avalanche breakdown respectively. One popular wide bandgap material
choice is SiC, which boasts not only a large bandgap value, but a higher thermal
conductivity than Si, while still being able to grow a native oxide. Research of SiC
power devices has been taking place since the 1990s and commercial devices have
1
been available to consumers since 2003 [13]. Further research is currently being
done to improve reliability and performance. Current SiC devices suffer from the
presence of interface traps which cause threshold voltage instability [14, 15] as well
as surface roughness which decreases channel mobility [1, 14].
In this thesis, I will first discuss my contribution to the study of SiC surface
mobility and will follow with a discussion of my atomistic calculations of Ge-Sn
alloys. In short, my work involves employing density function theory (DFT) to
solve the many-body Schrdinger equation which tells me the electron density and
potential everywhere within my constructed crystal systems as well as the band
structure. I have created a novel method to extract a perturbation potential due to
defects and/or the presence of an interface from this DFT solution. I then calculate
what we refer to as an atomic-roughness scattering rate using the calculated atomic-
level perturbation potential which is then fed into a Monte Carlo simulation. The
Monte Carlo simulation builds up electron velocity and energy distributions due
to different scattering rates within a semiconductor and provides as output a field-
based mobility. By using this atomic-roughness scattering we avoid approximations
typically used when calculating theoretical mobilities in SiC.
The DFT solution to the many-body Schrdinger equation also provides a band
structure for the given system. I use this band structure to predict the transition of
Ge into a direct-gap material when alloyed with Sn.
2
Figure 1.1: Plot of popular semiconductor elements and compounds relating atomic
bond length to energy bandgap. The scatter data shows a general trend of smaller
bond distances correlate to larger bandgaps. (Figure from [3])
3
1.2 Overview of Group IV(a) Semiconductors
1.2.1 Properties
Both silicon and germanium are good candidates for integrated electronics
due to their ability to thermally grow a native oxide. This allows for the creation of
good interfaces with few defects to be used for MOSFET gate oxides, beneath which
the sensitive channel resides. The compound semiconductor silicon carbide can also
grow a native oxide but it is known to be much more defected due in part to the
presence of carbon which must leave the material during the oxidation process.
Elements and compounds higher in a given column of the periodic table (Fig-
ure 1.2) tend to have larger bandgaps and decrease as you descend in the column.
Large bandgap materials are excellent candidates for high power and high temper-
ature devices because of their large breakdown electric fields and their low intrinsic
carrier concentrations respectively. On the other end of the spectrum, small bandgap
materials perform poorly in these situations, they tend to have increased mobility
which lends them to use in low power and high speed integrated circuits.
Table 1.1 lists various electrical and material properties of the important col-
umn IV semiconductors. The materials in the table are arranged from left to right
as they appear descending in the periodic table. Different properties tend to in-
crease or decrease across this table, almost without exception, which is a result
commonly seen in the periodic table arrangement of the elements. In essence, wide
bandgap materials have larger thermal conductivity and breakdown fields, while
4
13 IIIA 14 IVA 15 VA
5 10.811 6 12.011 7 14.007
2 B C N
Boron Carbon Nitrogen
[He]2s22p1 [He]2s22p2 [He]2s22p3
13 26.982 14 28.086 15 30.974
3 Al Si P Metal
Aluminium Silicon Phosphorus Metalloid
[Ne]3s23p1 [Ne]3s23p2 [Ne]3s23p3 Non-metal
31 69.723 32 72.64 33 74.922 Z mass
4 Ga Ge As Symbol
Gallium Germanium Arsenic Name
[Ar]3d104s24p1 [Ar]3d104s24p2 [Ar]3d104s24p3 GS Config.
49 114.82 50 118.71 51 121.76
5 In Sn Sb
Indium Tin Antimony
[Kr]4d105s25p1 [Kr]4d105s25p2 [Kr]4d105s25p3
Figure 1.2: Metalloid section of the periodic table where most of the important
semiconductor elements reside. Elements in each column have valence 3, 4, and 5
from left to right, and increase in radius from top to bottom.
small bandgap materials have higher mobility and permittivities.
The FCC forms of the first few group IV semiconductors (C, Si, Ge, ?-Sn) and
all polytypes of SiC are indirect gap materials. In contrast to this, a large number of
III-V compound semiconductors are known to be direct gap materials (GaAs, GaSb,
InP, InAs, InSb, AlN, GaN, InN). Though generally true, this fact is not readily
apparent without performing detailed band structure calculations.
5
Table 1.1: Material property comparison of Diamond, 4H-SiC, Si, and Ge at 300K
(descending in column IV from left to right).
Property Diamond [16] 4H-SiC Si Ge
Lattice Structure FCC Hexagonal FCC FCC
Bandgap Eg (eV) 5.45 3.2 1.1 0.66
Direct Gap E? (eV) 7.3 5-6 3.4 0.8
Breakdown Field (MV/cm) 10 3.0 0.6 0.1
Intrinsic Carrier Concentration (cm?3) 10?27 10?8 1010 2?1013
Bulk Electron Mobility (cm2/Vs) 2200 800 1200 3900
Bulk Hole Mobility (cm2/Vs) 1600 115 420 1900
Thermal Conductivity (W/cm K) 20 3-5 1.5 0.58
Relative Permittivity 5.7 9.7 11.9 16
1.2.2 Bonding
Also known as the Carbon Group, elements in this column of the periodic
table (Figure 1.2) have 4 valence electrons and readily form covalent bonds with 4
neighboring atoms in a tetrahedral arrangement due to sp3 hybridization. Electrons
in these hybridized orbitals exhibit 25% s and 75% p character, which result in
the 4-lobed tetrahedral shape that maximizes the separation of the electrons in the
surrounding bonds from each other. A spatial visualization of these orbitals is shown
in Figures 1.3 and 1.4 below.
In addition to the group IV elements forming tetrahedral bonding structures,
it is not unreasonable to assume that compounds of the form AB with one atom
of valence 3 and one of valence 5 will do the same but with higher ionicity bonds
instead. These tetrahedrons can be arranged differently with respect to each other
to give the crystal lattices in which these materials are commonly found. Elemen-
6
Figure 1.4: Actual form of sp3
tetrahedron (level set of the wave-
function) showing overlaid or-
bitals. Yellow and blue colors de-
Figure 1.3: Cartoon visualization note opposite signs the wavefunc-
of the four sp3 orbital lobes. tion.
tal column IV crystals tend to adhere to the face-centered cubic (FCC) diamond
lattice structure, while III-IV materials generally form in the analogous zinc-blende
structure (FCC) or more rarely the wurtzite (hexagonal) structure [6].
In general, as the atomic radius of an element increases (descending within a
column), bonds between nearest neighbor atoms become weaker due to spreading
of higher energy valence wavefunctions. Consequentially the atoms become more
separated and in turn exhibit smaller bandgap. This effect can be understood by
thinking about bringing N identical atoms close enough together for their atomic
orbitals to have significant overlap. When this happens, the energy levels of the
individual atoms (normally degenerate for like-atoms) must split into bands of N
closely-spaced states so as to not violate the Pauli exclusion principle. The amount
of band splitting is directly related to the degree of wavefunction overlap. This
7
causes lower energy orbitals which are confined to the inner core region to split less
than the higher energy valence orbitals with a larger spatial extent.
The manner in which this occurs leads to the formation of forbidden energy
gaps and determines the value of the bandgap based on the equilibrium atomic
spacing and the band occupancies. In addition to the atomic energies broadening
into bands, the bands begin to bend and in certain situations merge as orbitals
hybridize before ultimately separating again once inter-atomic distance becomes
very small.
Energy levels corresponding to wavefunctions of favorable (bonding) charge
distributions will tend to decrease with decreasing interatomic distance. Conversely,
antibonding energies where electrons are localized away from bonding sites will be
increased due to the energetically unfavorable configuration. These antibonding
wavefunctions are analogous to high energy molecular orbitals which occur in chains
of atoms. By taking linear combinations of the N atomic orbitals, A chain of N atoms
gives rise to N possible molecular orbitals which contain from 0 to n-1 nodes between
atoms where the density goes to zero. The 0 node solution corresponds to the fully
bonding wavefunction and the n-1 node to the fully antibonding wavefunction. The
band formation for a chain of 6 atoms is shown in Figure 1.5;
Similarly, in a crystal, the lowest energy level at the bottom of each band
corresponds to the in-phase bonding state combination of atomic orbitals and the
top level to the fully out-of phase antibonding configuration. Levels in the middle
of the band consist of states which are various other combinations of atomic orbitals
8
Figure 1.5: Molecular wavefunctions for a linear chain of 6 atoms. Degen-
erate energies split as interatomic distance decreases due to the interaction of
neighboring atomic orbitals. Light and dark circles represent positive and neg-
ative orbital lobes respectively. By Tem5psu (Own work) [CC BY-SA 4.0
(http://creativecommons.org/licenses/by-sa/4.0)], via Wikimedia Commons
which smoothly transition between these two extremes.
Once atoms are forced closer to and beyond their equilibrium distance, the hy-
bridized bands ultimately begin to split from each other due to the large Coulombic
repulsion of the ions [17]. The widths of the allowed energy bands are limited by the
small number of near-neighbor atoms which are close enough for their wavefunctions
to significantly overlap. Adding additional atoms to the system only results in fur-
ther subdividing states within the band such that it closely resembles a continuum
of states [18]. The band gap of a material is determined by the degree of separa-
tion between the last occupied band refered to as the valence band, and the next
higher unoccupied band known as the conduction band. In the case of FCC column
IV materials, the valence and conduction bands consist of the bonding sigma band
9
and the antibonding sigma? band respectively, the formation of which is depicted in
Figure 1.6.
1.3 Silicon Carbide
1.3.1 Background and Applications
The world of power electronics is ever growing and at present particularly due
to the recent concerns of alternative energy technologies and the commercial emer-
gence of fully-electric vehicles. As Si-based power devices approach their material
limits, alternative materials are being researched to take advantage of their differ-
ent material properties [19]. Typically, wide bandgap semiconductors are the main
go-to for use in high power, high temperature devices for their beneficial electrical
and thermal properties. SiC is of particular interest due to its ability to form a
stable native oxide and the commercial emergence of good quality substrates [20].
Since 2014, the SiC wafer market has grown each year by approximately 21% and
is predicted to continue this trend until at least 2020 to reach $110 million [21]. As
a whole, the SiC power semiconductor market is projected to exceed $1.6 billion in
2022 [22].
SiC crystals exist in numerous polytypes, meaning each crystal shares the
same close-packed layers but differ in the periodic stacking order of these layers. Of
the various polytypes, the 6H and 4H variants have seen the most use in electrical
device applications. Currently, 4H devices are the dominant among commercial
devices due to its higher carrier mobility and low dopant ionization energy [23].
10
Figure 1.6: Splitting and occupancy of 2s (3s,4s) and 2p (3p,4p) states of carbon
(silicon, germanium) as a function of inter-atomic spacing. The distance between
the top of the filled band and the bottom of the unfilled band shows formation of
various bandgap values at experimental inter-atomic distances. Inset images depict
spatial visualizations of the kinds of orbitals found in each band. Orbitals shown: 2s
orbital (bottom right), 2p orbital (top right), overlaid hybrid sp3 orbitals (center),
sp3 sigma bonding orbital (bottom left), and sigma? antibonding orbital (top left).
Color represents sign of wavefunction. [4]
11
In addition to SiC?s high thermal conductivity, saturation velocity, and breakdown
field, it performs better than Si on-chip due to both its lower on-state resistance,
smaller chip area, and higher switching frequency [19,23]. The price paid in increased
device costs is outweighed by the improved characteristics and performance gains
achieved using SiC.
One key application of SiC is for use in the construction of more energy efficient
and power dense DC/DC, DC/AC, AC/DC, and AC/AC converter circuits [23]. In
particular, DC/AC inverter circuits are needed in electric vehicles and photovoltaic
systems to convert the stored battery power into usable power. DC/DC converters
enjoy an increased operating voltage range and a drastic size reduction due to the
need for a smaller inductor. For a given peak-to-peak ripple in a DC-DC converter,
the inductance value needed is inversely proportional to the switching frequency of
the transistor [19]. The properties of SiC allow for high switching frequencies so the
size requirement of the inductor is be reduced, thus shrinking the module?s overall
size [24].
1.3.2 Crystal Structure and Properties
SiC crystals are arranged such that each carbon atom is covalently bonded
to four silicons and vise versa in a tetrahedral manner. Because these tetrahedra
may be stacked differently, SiC forms what are known as polytype variants. The
tetrahedrons create close-packed planes of atoms which vary in their stacking ar-
rangement along the normal direction. The order and periodicity of the stacking
12
defines the polytype of the crystal which is denoted <number><letter>-SiC, where
<number>is the period of the stacking arrangement, and <letter>abbreviates the
Bravais lattice of the crystal (C for cubic, H for hexagonal, R for rhombohedral).
Figure 1.7 illustrates the stacking arrangement for the most common polytypes,
although there are more than 200 known /citeLiu.
Figure 1.7: SiC polytypes depicting only key atoms forming the bonding chain in
the stacking direction. Primitive cells are shown for the hexagonal polytypes and
a non-primitive hexagonal cell is used for 3C. Polytypes differ only in the stacking
order and periodicity of the close-packed atomic planes. For the cubic variant, this is
the {111} plane and for the hexagonal variants, the {0001} plane. (Figure from [5]).
Because the different polytypes of SiC have different atomic arrangements,
they differer in their electrical characteristics and electronic structure. Table 1.2
compares the properties of the two most commonly used polytypes of SiC to the
industry standard Si. As shown in the table, the 4H variant has a slightly larger
bandgap and a greatly increased bulk mobility compared to 6H [22]. The breakdown
field of SiC is approximately 5 times larger than in Si, allowing it to hold-off the
13
much larger voltages seen in power applications. Additionally, SiC exhibits a 2-3 fold
larger thermal conductivity so heat can be dissipated from the device more easily,
preventing thermal runaway and reducing the need for large costly cooling systems
[19]. The increased thermal conductivity in addition to a significantly higher melting
temperature also increases the maximum operating temperature of the device and
consequentially, a larger maximum current density.
Table 1.2: Electrical property comparison of SiC to Si at 300K. Data from [25]
Property 6H-SiC 4H-SiC Si
Bandgap Eg (eV) 3.0 3.2 1.1
Breakdown Field (MV/cm) ?c-axis: 3.2 ?c-axis: 3.0 0.6
N = 1017D cm
?3 ?c-axis: >1 ?c-axis: 2.5
Intrinsic Carrier Concentration (cm?3) 10?5 10?8 1010
Bulk Electron Mobility (cm2/Vs) ?c-axis: 60 ?c-axis: 800 1200
N = 1016 cm?3D ?c-axis: 400 ?c-axis: 800
Bulk Hole Mobility (cm2/Vs) 90 115 420
ND = 10
16 cm?3
Thermal Conductivity (W/cm K) 3-5 3-5 1.5
Relative Permittivity 9.7 9.7 11.9
Because the most commonly used polytype of SiC in power devices is the
4H variant, the remainder of this thesis will discuss 4H-SiC exclusively. Figure 1.8
shows the hexagonal lattice of 4H-SiC along with the Miller indices of associated
directions/planes.
The (0001) Si-face of 4H-SiC is by far the most commonly available wafer type.
On this face, Si atoms terminate 100% of the face but other faces have different Si-to-
C ratios and packing density. The various atomic structures exhibited by each face
causes different electrical transport and oxidation characteristics. Opposite to the Si-
14
[0001]
[11?00]
[112?0]
Figure 1.8: Wigner-Seitz cell of hexagonal real-space lattice structure and relevant
crystal planes
face is the C-face (0001?) which is terminated completely by C atoms. Perpendicular
to these are the a-face (112?0) and the m-face (11?00) which are investigated for their
use in vertical devices and have an equal number of Si and C atoms. Figure 1.9
shows the stacking arrangement and key faces of the 4H-SiC structure.
In terms of process variation for different faces, the Si-face shows the slowest
oxide growth rates which are approximately 12 times slower than the C-face. The
a-face growth rate is slightly slower than the C-face but still significantly faster
than the Si-face which takes 12 hours to grow 50nm of oxide at 1150?C [22]. As
a consequence of different surface morphologies, some faces show greater electron
mobility than others. Specifically, the Si-face exhibits the lowest channel mobility
(< 10 cm2/Vs) compared to the a-face mobility (> 25 cm2/Vs) depending on the
15
Figure 1.9: 4H-SiC lattice depicting the stacking arrangement and crystal faces
16
degree of surface passivation [22]. With NO surface passivation these mobilities can
be improved to 80 and 125 respectively. A mobility of 90 cm2/Vs has been achieved
on passivated C-face [22]. Despite having the lowest mobility, the Si-face is still
commonly used for its superior threshold voltage stability compared to its other
faces [26].
1.3.3 Issues in SiC Power Devices
Traditionally, SiC MOSFETs have suffered from low channel mobility due to
poor quality SiC-SiO2 interfaces [19]. In particular, large surface roughness and
a high density of interface states increases carrier scattering in the channel, caus-
ing the surface mobility to drop approximately 1-2 orders of magnitude from the
bulk value [3]. The low mobility of SiC MOSFETs causes devices to have a small
transconductance and thus poor gate control of current. Various defects have been
detected and characterized by many researchers including [14, 27, 28]. In addition
to causing scattering, some defects can cause carrier trapping and threshold volt-
age instability, leading to reliability issues which is a heavily researched topic [15].
These defects are not passivated by post-metalization H2 anneals as they are in Si,
but NO has been shown to improve interface defect density [14,22].
1.4 Germanium Photonics: Benefits and Challenges
In addition to SiC I also investigate the optical properties of Ge and its po-
tential use in beyond-silicon integrated electronics at the end of the Si road map.
17
Germanium is a natural replacement for CMOS due to it?s ability to grow a native
thermal oxide, a property crucial to form the gate oxide in MOS transistors. Cur-
rently, research has been done seeking to integrate Ge into Si CMOS devices. Ge
is a decent candidate for integration due because it has the same lattice structure
as Si with a different lattice constant. This integration does not come without dif-
ficulty however. The lattice mismatch necessitates a high concentration of misfit
dislocations at the interface and further challenges exist regarding compatibility of
processing growth temperatures [29].
Germanium is known to have an indirect band gap of 0.66eV, allowing it
to function as a relatively inefficient detector of short-wave infrared (SWIR) light.
Instead of trying to integrate a Ge detector into Si CMOS circuitry, the photonics
can be developed on a purely germanium substrate so that both the optics and
supporting electronics are both made of Ge. In doing this, we can exploit the
beneficial electrical properties of Ge such as its superior (3x) mobility to Si.
In this work, our main focus will be on the study of the Ge band structure and
how the conduction band minima change with respect to varying concentrations of
Sn in Ge-Sn alloy. By introducing a certain fraction of Sn to the Ge, the bandgap
can be transformed from a indirect in nature to direct, allowing for more efficient
photon detection. The key to this process is finding the minimum amount of Sn
required to transition so that the gap energy is still as close to its original value as
possible.
18
1.5 Research Accomplishments
SiC:
? Performed a self-consistent DFT calculation of a 125 atom ideal 4H-SiC/SiO2
interface supercell using Quantum ESPRESSO [30].
? Developed a novel technique to extract an interface scattering perturbation
potential from DFT calculations. This scattering is due to the aperiodicity
caused by the presence of the interface and the atomic roughness therein.
? Calculated a 2D quantized inversion layer scattering rate based on the ex-
tracted perturbation potential for use in surface Monte Carlo scattering sim-
ulation.
? Reproduced bulk 4H-SiC field-dependent velocity characteristics from Monte
Carlo simulations.
? Calculated and compared mechanism-dependent component mobilities impor-
tant to MOS device mobility including my own proposed atomic-roughness
mobility which will be present even in an ideal SiO2/SiC interface.
? Combined and adapted physics-based models for doped semiconductor systems
to calculate theoretical mobile hole to acceptor concentration ratio (p/NA)
? Gathered, analyzed, and created empirical parameterizations for concentration-
dependent data of nearly all existing published Hall mobility, resistivity, and
ionization energy data on Al-doped 4H-SiC. I used the ionization energy data
as input into my physics-based model. The Hall mobility and resistivity data
were used verify my physics-based p/NA result by using a method which was
19
first applied to Si.
? Created a more readily usable closed-form expression of my temperature and
concentration-dependent p/NA calculation. The original physics-based model
involves iterative numerical techniques to solve a nonlinear system of integrals
which is computationally expensive so I developed a genetic algorithm to find
the optimal parameters for an expression which reproduces the same two-
dimensional p/NA function but is much more easily implemented.
? Used a variation of my genetic algorithm to parameterize the fully temperature
and concentration dependent Al-doped 4H-SiC Hall mobility function based
on data from nearly all existing published values. This result should be valid
over a larger temperature and doping range than those currently published.
Ge:
? Reproduced key gap energies in the band structure of Ge using a hybrid DFT
calculation.
? Calculated band structure for various GexSn(1?x) alloys with various ratios.
? Determined the fraction at which the band structure of the material becomes
direct.
1.6 Future Work
Currently, the scale at which my DFT calculations can be performed is sig-
nificantly limited by the practical memory and processing power constraints of the
computer on which I am performing the calculations. This limits the size of the
20
supercell to the order of about 100 atoms for practical computations. With a more
powerful computer, it would be possible to simulate significantly larger supercells
which could accommodate interesting defect structures. The issue with trying to
simply add a defect to a smaller supercell is that this cell is periodic and thus any
defect you add will also be periodic. Without a sufficiently large cell, this defect will
have a higher concentration than the real-world analogue which can cause incorrect
results due to the defect interacting with itself. Additionally, to simulate aperiodic
structures like amorphous oxide, surface roughness, defects, etc., the correlation
length of the structure is bounded by the periodicity of the supercell. Therefor,
with a larger supercell, a ?less periodic? structure can be simulated.
To extend my research of scattering at the SiC/SiO2 interface, two main ape-
riodic structures are prime candidates for a DFT-based study:
First, it would be interesting to simulate a non-pristine interface which con-
tains a real atomic configuration of the miscut roughness and step bunching which
recreates an AFM-measured SiC surface. Performing statistical calculations of the
RMS roughness height and the autocovariance of the supercell?s atomic configura-
tion should be matched to the real measurements. A reverse Monte Carlo algorithm
could be used to continually generate random atomic step configurations for the su-
percell until the toughness statistics match, then an amorphous oxide can be added
and structurally relaxed to create the final interface. Even better (but significantly
more computationally intensive) would be to simulate the growth of the oxide from
the rough miscut surface using a molecular dynamics simulation. Finally, with this
supercell structure, the surface-roughness scattering rate of the miscut surface could
21
be calculated using the methods described in my research and this would give in-
sight as to the impact the miscut roughness has on the SiC MOS channel mobility.
Additionally, it would be interesting to directly measure experimentally the chan-
nel mobility of MOS structures using identically- produced SiC wafers with varying
miscuts.
Second, basal plane dislocations (BPDs) are known to be device-impairing
defects in MOS structures. However, device fabricators generally try to transform
these dislocations into supposedly less harmful threading edge dislocations (TEDs).
The actual effect of the TEDs has yet to be fully studied and characterized, so a
DFT study to extract the scattering rate of TEDs would be useful to inform device
designers and fabricators how their presence is affecting the MOS channel mobility.
In terms of my work on incomplete ionization, I have extracted and fit with my
genetic algorithm the Hall mobility literature data as a function of both temperature
and NA, however the resistivity data is a much more complicated function and
I only fit it as a function of NA at T = 300K. To fully check agreement with
my temperature-extended theoretical model results of incomplete ionization, a fit
of the temperature and NA dependent resistivity literature data would need to be
performed.
22
Chapter 2: DFT and Electron Transport Modeling and Theoretical
Framework
2.1 Many-Body Schrodinger Equation
Determining the properties of a multi-atomic system traditionally involves
solving the many-body Schodinger equation. This, however, quickly proves infeasible
when the number of atoms increases as the wavefunction quickly becomes a function
of far too many coordinates (3? (n + N) for n electrons and N nuclei). The time-
independent Schrodinger equation itself is an eigenvalue problem which has allowed
energy states as eigenvalues and their corresponding wavefunctions as eigenvectors.
( ?N ~2 ?n 2 2 ?n ?N n n? ?2 ? ~ ?2 ? q ??? Zj ?
2 ??
2M Rjj 2m
ri 4? ?
j=1 ?i=1 )
0
? i=1 j=1 ~r(i ?R
~
j?
1 q 1
+ +
2 4?0 ? |~r ? ~r ? |i=1 i =6 i ii
2 N N1 q ??? Z Z
)
j j? ? ? = E? r~1, . . . , r~ , R~n 1, . . . , R~N (2.1)
2 4?0 ~j=1 j? 6=j Rj ?
?
R~ j??
Because the wavefunction of a system contains within it all of the information
knowable about the system, the left-hand side of the equation is a large operator
which acts on the wavefunction to extract different energy components. This opera-
tor is also known as the Hamiltonian of the system. The Hamiltonian can be broken
down into a combination of smaller quantum-mechanical operators which are, in
23
order: nuclear kinetic energy, electron kinetic energy, electron-nuclei coulombic po-
tential energy (attraction), electron-electron coulombic potential energy (repulsion),
nuclei-nuclei coulombic potential energy (repulsion). In this equation Mj, Zj, and
Rj are the mass, atomic number, and position of the j
th nuclei and m, q, ri are the
mass, charge, and position of the rth electron.
2.1.1 Born-Oppenheimer Approximation
To begin simplifying the problem, the equation is usually solved in the context
of the Born-Oppenheimer approximation, which assumes that because the nuclei are
far more massive than the electrons, they move on different time scales and as such,
the total wavefunction for the system can be separated into an electronic component
and a purely nuclear (vibrational,rotational) component.
?tot = ?electronic ??nuclear (2.2)
The electronic solution is found by solving the electronic Schrodinger equation with
the positions of the nuclei fixed (usually in their equilibrium configuration). By
slowly varying the location of the nuclei and re-solving the electronic equation, the
electron energy as a function of nuclei positions can be extracted. This term acts
as the nuclear potential term and in combination with the nuclear kinetic energy
term, is used to form the nuclear Schrodinger equation who?s solution results in the
nuclear portion of the wavefunction.
24
2.2 Density Functional Theory
2.2.1 Origins of DFT
The motivation to find a theory such as DFT comes from the immense difficulty
of solving a large quantum mechanical system exactly. The majority of the difficulty
in solving such systems stems from the high dimensionality of the problem. For a
system of N electrons and M nuclei, the total wavefunction is a function of 3(N+M)
spatial variables where N and M are on the order of Avogadro?s number for a real
crystal. A system of this size is effectively impossible to solve given the finite memory
and calculation speed of todays computers; even smaller systems of a few hundreds
of atoms are impractical to solve within any reasonable amount of time [31]. These
restrictions sparked a need for a practical way to solve these systems via a reduction
of dimensionality.
Thomas and Fermi in 1927 were the first to propose that the electron density
was the fundamental variable defining the quantum many-body system. In making
this assumption, calculation of the complicated high-dimensional wavefunction can
be ignored and instead the simple 3-dimensional density is used in the calculation.
Though their formulation in general proved too crude for accurate calculations, their
use of the density as a fundamental variable laid the foundation for the later greatly
successful density functional theory (DFT).
Around the same time the Thomas-Fermi model was being developed, another
important technique emerged which aimed to approximately solve the many-body
25
Schrdinger equation. First introduced in 1928 by Hartree [32], this technique became
known as the self-consistent field method due to need to iterate solutions to achieve
convergence. In this framework, the wavefunction for each electron ?i is calculated
assuming each electron experiences a total potential created by the atomic nucleus
and all other electrons in the system. By starting with an initial guess for the
electron density, an approximate potential and Hamiltonian is determined from
which the wavefunctions can be solved. These wavefunctions in turn determine
an electron density and iteration is performed until the system is solved in a self-
consistent manner. The total wavefunction of the complete system in Hartree?s case
was approximated to be the product of all the one-electron wavefunctions [31].
?(r1, r2, . . . , rN) = ?1(r1) ? ? ??1(rN) (2.3)
The original Hartree equation is shown in Equation 2.4 where the two po-
tential terms are the external potential Vext(r) created by the atomic nuclei and
the potential due to the mean field of all the other electrons. The equation must
be solved iteratively because the latter term is function of the other wavefunction
solutions. ( ?? )N1 |? (r?)|2? ?2 j+ Vext(r) + dr? ?i(r) = i?i(r) (2.4)
2 |r? r?|
j=6 i
The total energy of the system, shown in Equation 2.5, is a simple sum of the
individual electron energies minus the double counting from the electron-electron
Coulombic potential energy term.
?N 1 ?N ? ? |? (r)|2H ? i |?j(r?)|2E = n drdr? (2.5)
2 |r? r?|
n=1 i=6 j
26
Later, Fock and Slater in 1930 showed that by replacing the simple product
with a determinate of the single-electron wavefunctions, the resulting total wave-
function automatically satisfies the antisymmetric requirements of the Pauli exclu-
sion principle and treats particle exchange exactly [31, 33, 34]. The determinant
construction of the total wavefunction is shown below in Equation 2.6.
???? ?? ???1(r1) ? ? ? ?N(r1) ???1
?(r) = ? ??? .. . . . ?? . . .. ???? (2.6)N ! ???1(r ) ? ? ? ? (r )?N N N ?
In 1964, the density functional formalism set forth by Hohenberg and Kohn gave
the work done by Thomas and Fermi mathematical grounds on which to stand. The
Hohenberg-Kohn Theorem showed that the ground state of an atomic system is
fully determined by the density distribution n(r) of the electron [31]. Furthermore,
the energy can be written as a functional of the density, meaning it takes in the
density function as input and returns the total energy (a single value) as output.
More specifically, they showed that there exists a universal energy functional F [n(r)]
which applies to any and all systems. This powerful proof, along with the proof
of a variational principle, allows us to find the ground state by hunting for the
density which minimizes the energy. In principle we can even ignore calculating the
wavefunction entirely, though some techniques still involve solving for it. By simply
working with the density, the dimensionality of the system is reduced from the order
of Avagadro?s number to only 3 spatial dimensions. Though this theory lays the
27
groundwork for greatly simplifying the problem, it unfortunately does not give any
insight into form of the universal functional. The details of each of these theories
are discussed further in Appendix A.
2.2.2 Kohn-Sham Method
Developed in 1965, this formulation of DFT acts as a practical implementation
of the Hohenberg-Kohn Theorem (Appendix A.3) by defining a set of component
energies, each with a clear physical origin, which sum to approximate the universal
energy functional F [n(r)]. The central feature of this method involves defining a
fictitious system of noninteracting electrons moving in an effective external potential,
which gives rise to the same density as the true interacting system. All of the
single-particle wavefunctions satisfy a Schrdinger-like equation known as the Kohn-
Sham equation (Eqn. 2.7 in atomic units), where the potential term is the Kohn-
Sham effective potential - a functional of the density derived from the terms in
the energy functional. The lowest N eigenstates of this equation give the N single-
particle wavefunctions ?i known as the Kohn-Sham orbitals. The ground state
wavefunction of this system can then be exactly written as a Slater determinant of
the single-particle wavefunctions to account for anti-symmetry because there is no
electron-electron interaction (Eqn. 2.8). The total density can be calculated from
the single-electron orbitals using Equation 2.9.
?1?2?i(r) + VKS(r)?i(r) = i?i(r) (2.7)
2
28
???? ?? ??? (r ) ? ? ? ? (r )1 ???
1 1 N 1
. . . ???
?
?(r) = ? ? .. . . .. ???? (2.8)N ! ???1(rN) ? ? ? ?N(r ?N)?
?N
n(r) = ??i (r)?i(r) (2.9)
i
Taking the expectation of the universal functional with the full Kohn-Sham
Slater determinate wavefunction gives us insight to the forms of terms in the total
Kohn-Sham energy functional EKS[n] with which we aim to approximate the true
functional F [n].
EKS[n?(r)] = TS[n] + Eext[n?] + EH [n] + EX [n] +?EC?[n] (2.10)N
KS 1
? 1 n(r)n(r?)
E [n(r)] = ? ?? 2i (r)? ?i(r)dr + Vext(r)n(r)dr + drdr?2
i
1 ?? ?
2 |r? r?|
?? ?i (r)?i(r?)?j(r?)?j(r)? dr?dr + EC [n(r)]
2 |r? r?|
i,j
Here, TS[n] represents the ?single particle? or non-interacting kinetic energy
which is why it can be calculated using a simple sum over all of particle?s wavefunc-
tions. Eext[n] is the external energy which comes from any external potential, such
as the Coulombic attraction to ionic cores of atoms. EH [n] is the Hartree energy
caused by the Coulombic repulsion of all other electrons. EX [n] is the exchange
energy due to Pauli?s principle and the antisymmetry of the wavefunction. These
are the terms which have tractable expressions. The final term EC [n] represents the
correlation energy, which can be thought of as an error term. The correlation energy
29
contains all energy differences between our constructed non-interacting system and
the true system energy functional F [n] and accounts for about 10% of the total
energy of a system. This correction energy accounts for the self-interaction error
within EH as well as the difference in the kinetic energy between the fully interacting
and non-interacting system. The error in EH for the interacting system comes from
the fact that when the integral is solved for a 1-electron system, it does not identi-
cally equal 0, meaning the electron is interacting with itself so a cancellation term
is needed. Exchange energy can also been seen as the effective repulsion of electrons
with parallel spins due to the Pauli exclusion principle, and the Correlation energy
due to interaction of electrons with anti-parallel spins.
The total energy functional of the Kohn-Sham method depends on the approx-
imations made in formulating EC . Additionally, given the computational expense
associated with evaluating EX , the two terms are often combined and approximated
together as one EXC term. Although DFT is exact in theory, we do not know the
form of the universal functional - specifically the EC term, so approximations to this
introduce error and cause the method to be an approximation in practice.
Working backwards from the energy functional, we can generate the effective
Kohn-Sham potential needed in the Kohn-Sham equation by taking the functional
derivative of the external, Hartree, exchange, and correlation energies with respect
to the density. These potentials represent the Coulombic ion potential, the Coulom-
bic electron-electron interaction potential, and quantum-correction potential respec-
30
tively.
?EH ?EXC
VKS(r) = Vext(r) + + ? = Vext(r) + VH(r) + VXC(r) (2.11)?n(r) ?n(r)
n(r?)
VKS(r) = Vext(r) + dr? + V| ? | XC
[n(r)]
r r?
Because the exchange-correlation potential VXC depends on the density n(r),
itself depending on the orbitals ?i(r), which in turn depend on the potential VKS,
solutions to the Kohn-Sham equation must be found iteratively to achieve self-
consistency.
2.3 DFT Calculation Results
Figure 2.1 shows the calculated band structure from my DFT simulation us-
ing the open source software Quantum Espresso [30] with the PBE functional [35]
and pbe-hgh psudopotentials from http://www.quantum-espresso.org. The result-
ing constant-energy ellipses of these minima calculated from this DFT computation
are shown in Figure 2.2.
Notoriously, DFT tends to underestimate the band gap of semiconductors due
to the approximations associated with the form of the energy functional (which
actually contains discontinuities), and possibly due to an inevitability of the Kohn-
Sham theory itself [36]. However, band structures resulting from DFT calculations
still give useful insight into the structure of the actual bands and are often used in
theory-based calculations.
For the purposes of creating a scattering theory, the important conduction
band valleys are the two which lay nearly degenerate at the M point. The resulting
31
20
15
10
5
0
-5
-10
-15
? M K ? A L H A
Figure 2.1: Theoretical band structure of 4H-SiC calculated with DFT.
gap energies and extracted effective masses of the DFT calculation are provided in
Tables 2.1 and 2.2 respectively, where they are compared to the experimental values.
In Chapter 3, this band structure calculation will be developed into a band model
to be used in scattering rate calculations and Monte Carlo scattering simulations.
32
E [eV]
k  [2 ?/a]
z
0.1
0
-0.1
0.5
k  [2 ?/a]
y
0
0.6
0.4
0.2
-0.5 0
-0.2
-0.4
-0.6 k  [2 ?/a]
x
Figure 2.2: Constant energy ellipses for the lowest conduction band minimum at
the M point, calculated using DFT.
Table 2.1: Comparison of calculated energy gap values of 4H-SiC to those found in
the literature.
Gap [eV] Calc. Lit. [7, 8]
Eg (EM) 2.23 3.23
EsM 0.13 0.12
EL 2.65 4
E? 5.12 5-6
Ecr 0.07 0.08
Eso ?0 0.007
Table 2.2: Comparison of extracted effective masses in the M valley from my DFT
calculation to the experimentally measured values.
Mass [m0] Calc. Exp. [37]
EM?? 0.56 0.58
EM?K 0.39 0.31
EM?L 0.31 0.33
33
Chapter 3: Atomic-Level Electron Transport
3.1 Perturbation Scattering Theory
Now that we have set up a way to solve the quantum many-body problem
through DFT, we need to lay the framework required to understand electron trans-
port within a crystal. By understanding the electron transport of a material, electri-
cal characteristics can be predicted and possibly improved by changing the atomic
structure. In this chapter, I will first discuss existing perturbation scattering the-
ory including Fermi?s Golden Rule from which my results are derived and known
scattering rate results. Then I will explain the technique for how I calculated my
own atomic-roughness scattering rates for the best-case SiC-SiO2 interface based on
density functional theory calculations.
The primary difference between traditional Si MOS and emerging SiC MOS
technology is the quality of the semiconductor-oxide interface. The poorer quality of
the SiC/SiO2 interface leads to significantly lower mobility due to surface roughness
scattering and trapped interface charge scattering [38,39]. The treatment of scatter-
ing in SiC proceeds similarly to that of Si, but with a different number of equivalent
valleys and with focus on different dominant scattering mechanisms. In this study, I
apply the standard method of perturbation scattering theory to determine the scat-
34
tering rates for the various important mechanisms in SiC. Additionally, I introduce
a new mechanism deemed intrinsic atomic-roughness scattering not previously in-
vestigated. This new mechanism is calculated in the form of a deformation potential
scattering using data extracted from DFT simulations to model the interruption of
atomic potential periodicity due to the presence of the interface.
When an electron is moving in a crystal, any number of deviations from a
perfect crystal at 0K can cause a perturbation and thus a scattering event. Fig. 3.1
shown below illustrates the family tree of scattering mechanisms which are consid-
ered in various semiconductors.
Opical
Intervalley
Acoustic
Phonon Polar
Opical
Nonpolar
Intravalley
Piezoelectric
Acoustic
Deformation Potential
Scattering Ionized
Impurity
Neutral
Defect Space-Charge
Alloy
Carrier-Carrier
Figure 3.1: Classification of scattering mechanisms in semiconductors. Taken from
[6]
These deviations can be treated using the method of perturbation theory,
35
which treats phonon, impurity, and other interactions as small deviations to the
local potential an electron experiences. The perturbation potential changes the
Hamiltonian of the system slightly and perturbation theory allows us to determine
the resulting eigen energies and wavefunctions from the original solutions. The
theory can then be extended to give a scattering rate for an electron in a given state
based on the perturbation.
The first order time-dependent quantum perturbation theory is derived as
follows. First we start with a general time-dependent Schrdinger equation which all
single-electron wavefunctions must obey.
??k(r, t)
i~ = H0?k(r, t) (3.1)
?t
We also start with an initial unperturbed Hamiltonian H0 with associated
(0)
energies Ek and wavefunctions ?k, for which the solution?s are known, as they
appear in the time-independent Schrdinger equation.
(0) (0) (0)
H0?k = Ek ?k (3.2)
We can write the corresponding time-dependent solution of the initial wave-
function as:
(0) (0) (0)
? (r, t) = ? (r)e?iEk t/~k k (3.3)
We now define a perturbed Hamiltonian H as it differs from the initial Hamil-
tonian H0 by an arbitrary scaling constant ? times a perturbation Hamiltonian H
?.
H = H + ?H ?0 (3.4)
36
Solutions to any Hamiltonian form a basis, so we can write the solution of
the perturbed system as a linear combination of the orthonormal set of solutions
of the unperturbed Hamiltonian. The coefficients needed in the sum are in general
time-dependent and represent the time variation in the wavefunction due to the
perturbation.
? ?
(0) (0) (0)
?(r, t) = ck(t)?k (r, t) = ck(t)?k (r)e
?iEk t/~ (3.5)
k k
Substituting this form of the wavefunction into the perturbed Schrdinger equa-
tion and canceling terms yields:
? ?ck(t) (0) ? (0) ?i~ ? (r)e iEk t/~ ? (0) (0)k = ? H ck(t)?k (r)e?iEk t/~ (3.6)?t
k k
Multiplying both sides of this equality by the complex conjugate of the wave-
(0)?
function for another state ? iEk?t/~k? (r)e , we then integrate both sides over all r
resulting in Equation 3.7. Here I?ve replaced the unperturbed wavefunction in state
k with |k? and the wavefunction for state k? with ?k?| using bra-ket notation.
?ck?(t) ? (0)
i~ = ? ck(t) ?k?|H ? |k? ei[Ek??Ek ]t/~ (3.7)
?t
k
The expression ?k?|H ? |k?, also known as the matrix element of the pertur-
bation potential represents the transition amplitude between the two states. It?s
explicit form is shown in Equation 3.8 and resembles an expectation energy calcu-
lation of the perturbation but uses two different states in the expression.
?
?k?|H ? | ? (0)? ? (0)k = ?k? (r)H ?k (r)dr (3.8)
37
Assuming the coefficients ck(t) vary slowly with time for a weak perturbation,
we can express them as a power series in ?.
(0) (1) (2)
ck(t) = ck + ?ck (t) + ?
2ck (t) + ? ? ? (3.9)
By substituting Equation 3.9 into Equation 3.7 and equating the powers of ?
on both sides we develop a system of equations proportional to increasing powers of ?
which represent different order approximations which can be successively evaluated.
(0)
?c
i~ k?? = 0 (3.10)?t(1)?c ? (0) ? (0)i~ k? = ?k?|H |k? c ei[Ek? Ek ]t/~
?t k
k
(2)
?c ?
~ k? ? | ? | (1)
(0)
i = k? H k? c ei[Ek??Ek ]t/~
?t k
k
... (3.11)
Assuming the electron exists in a single unperturbed state such that at time
(0)
t = 0, ck (0) = 1 and all other coefficients are 0. With this assumption, to first-i
order, the coefficient is determined by Equation 3.12, where the perturbation is
assumed to be harmonic in nature such that H ? = V ?e?i?t.
? t (0)
(1) 1 ?
ck? (t) = ?k?|
? | ? i[EV k e k??Ek ?i~?]t /~i dt? (3.12)
i~ i0
? (0)i[Ek? Ek ?~?]t/~
(1) 1 e i ? 1
c ?k? (t) = ?k?|V |ki?i~ (0)i(Ek? ? Ek ? ~?)/~i
(1) 1 ? | ? | ? i?t sin(?t)ck? (t) = k? V ki e ti~ ?t
? (0)? = (Ek? Ek ? ~?)/2~ (3.13)i
The probability of the electron being in state k? at time t is given by |ck?(t)|2
38
and the transition rate from state ki to state k
? is given by
|ck?(t)|2
S(ki,k?) = lim (3.14)
t?? t [ ]
|?k?|V ? |k ?|2 2i sin(?t)
S(ki,k?) = lim t
t?? [ ~2] ?t2
The squared sinc function term sin(?t) in the limit as t goes to infinity acts
?t
like a Dirac ?-function because the width of the peak becomes very thin. Addition-
ally, instead of integrating to 1, it integrates to ?/t when ? it taken from ?? to
?. This observation allows us to replace the term in Equation 3.14 with ?(?)?/t,
which acts as an energy conservation enforcement term in Equation 3.15. The re-
sulting equation illustrates Fermi?s golden rule, showing the transition probability
from one state to another due to a perturbation is proportional to the matrix ele-
ment squared. The upper sign in the ?-function corresponds to the emission of a
photon with energy ~?, and the lower to the absorption.
2? ( )
S(ki,k?) = |?k?|V ? | ?|2 (0)ki ? Ek? ? Ek ? ~? (3.15)~ i
This result holds true for simple harmonic perturbations which are otherwise
constant in time. The calculation also only considers elastic interactions of the
electron and phonon. From this transition probability the scattering rate can be
determined by multiplying by the number of states per unit volute (?/(2?)3) and
integrating over all final states. The delta function will pick out all eligible final
states which satisfy energy co?nserv?atio?n.
? 2? ? ?
W (k?i) =? ? S(k ,k?)(k
? 2
(i ) sin(?)dk
?
) d?d? (3.16)(2?)3 0 0 0? 2? ? ? |? | ? | 2 (0)W (ki) = k? V ki?| ? Ek? ? Ek ? ~? (k?)2sin(?)dk?d?d?~(2?)2 i0 0 0
39
Figure 3.2 shows an example of the scattering system in consideration, where
initially the electron is in state k. The electron scatters and absorbs a phonon with
wavevector q and ends in state k? which is separated by an angle ? from k.
k?z
k
q
k?
?
?
?k
?k
k?y
k?x
Figure 3.2: Illustration of initial vector k scattering to a different state k? via the
absorption of a phonon with wavevector q. Note: ? is in general not in the same
plane as ?k, which is simply used to define the position of k.
3.2 Scattering and Transport Model
Due to the charge of the electron, when an electric field is applied to a semicon-
ductor, electrons in the conduction band experience a force and thus an acceleration
in the direction opposite to that of the field. This motion is interrupted by scat-
tering events which we treat as instantaneous events that change the wavevector of
the electron. The scattering at each event is caused by one of the many mechanisms
shown previously in Figure 3.1 which consist of various phonons representing lattice
distortions as well as other lattice perturbations. To study the transport of electrons
40
in this system, we rely on the semi-classical Monte Carlo method to stochastically
simulate an electron as it is accelerated and scattered through the band structure of
the crystal. By following the motion of a single electron as it scatters and changes
energy, we can build up statistics of velocity and energy which ideally represent the
ensemble of electrons simultaneously moving within the crystal. The Monte Carlo
method is considered semi-classical because though the scattering mechanisms are
treated quantum-mechanically using Fermi?s Golden Rule, the free-flight transport
between events is treated using classical particle dynamics. As the electron drifts in
the electric field, the wavevector changes in time, obeying Equation 3.17 and where
the velocity is given by Equation 3.18.
dk ?qF
= (3.17)
dt ~
1 dE(k)
V (k) = (3.18)
~ dk
The band structure on which the electrons move form energy surfaces within
the Brillouin Zone (BZ) of the crystal. The BZ is defined to be the Wigner-Seitz cell
constructed for the crystal?s reciprocal lattice. By construction, 1st BZ represents all
points in reciprocal space closest to an arbitrary reciprocal lattice site. The shape
of this zone partitions all of reciprocal space in a periodic fashion and contains
every unique wavevector (k-vector) which correspond to the allowed states within
the crystal. Associated with these states are quantized energies forming the band
structure. By reducing the zone by all symmetries of the crystal, we obtain what is
41
known as the irreducible wedge which contains the non-redundant information of the
energy bands. The edges and vertices of this wedge correspond to highly symmetric
points and directions of the reciprocal lattice. To plot the full band structure, a four
dimensional plot would be required, so instead a representative sample of the zone
is taken by following the edges of the irreducible wedge to form a path between the
high symmetry points. The allowed energies at each k-vector are plotted against
the distance traveled along the path. Fortunately, the band extrema almost always
lay along these symmetry directions so the band gap as well as important valleys of
the conduction band and other features of the energy landscape are revealed in this
manner.
In contrast to the familiar Si which has a diamond crystal structure and fcc
lattice, the most electrically significant polytypes of SiC have a wurtzite crystal
structure which correspond to a hexagonal lattice. The BZ of the hexagonal lattice
is simply a hexagonal prism. A depiction of the hexagonal BZ is shown in Figure
3.3 along with it?s irreducible wedge and high-symmetry points which are visited in
band structure diagrams.
The band structure for 4H-SiC is known to have its conduction band minimum
at the M point, making it an indirect gap material with 6 equivalent minima in the
BZ. In general, it is common to approximate the energy dispersion relation near the
conduction band minima to be parabolic such that energy E is proportional to the
squared magnitude of the wavevector k. However, in wide bandgap semiconductors
like SiC where electric fields allow carriers to gain energy such that they reside away
42
A
L
H
?
M
K
Figure 3.3: Brillouin Zone of a hexagonal lattice.
from the band minimum, a non-parabolicity factor ? is introduced as in Equations
3.19 and 3.20 to better account for the true shape of the bands.
~2k2
?(E) = E(1 + ?E)(= ? )(spherical) (3.19)2m
~2 k2 2k2
?(E) = E(1 + ?E) = l? +
t
? (ellipsoidal) (3.20)2 ml mt
For many semiconductors including SiC, the constant energy surfaces are el-
lipsoidal and have a different masses depending on the direction of the field relative
to the valley. These masses are denoted as the longitudinal effective mass ml aligned
radially out from the center of the BZ to the location of the valley, and transverse
effective mass mt for transport in the plane perpendicular to the longitudinal direc-
tion.
From my full band structure calculation showed previously in Figure 2.1, we
can extract a simplified band model containing only the important energy minima
of the conduction band and the maximum of the valence band. These important
43
bands are shown in Figure 3.4, taken from [8].
Figure 3.4: Simplified band illustration for 4H-SiC with associated energies from
literature [7] (Figure from [8]).
Due to the small difference in energy of the two lowest bands in the M valley,
we will initially treat this as a doubly degenerate valley as a first approximation,
so scattering rates are multiplied by a factor of 2 because of the doubling of the
density of states. Additionally, the calculated overlap integrals of the wavefunctions
from the minimum at any one M valley to the other 5 equivalent M valleys are small
(0.005, 0.162, 0.016, 0.162, 0.005) and expected to be similarly small for points near
the minimum [40]. This indicates that there is little optical intervalley scattering
compared to intravalley scattering so we will approximate the value of the overlap
to be 1 when scattering within a valley and 0 when scattering to other valleys.
44
The Monte Carlo simulation I performed in this work uses deformation poten-
tial values fitted by Hjelm et al. [41] as well as the corresponding phonon tempera-
tures and other physical constants which are reproduced in Table 3.1. The phonon
temperatures are related to the phonon energies and frequencies via Equation 3.21.
Eph = ~?ph = kBTph (3.21)
Table 3.1: Phonon and material properties taken from the model used by Hjelm [40]
Parameter Value
Acoustic Deformation Potential DA (eV) 19
Polar-optical Phonon Temperature (K) 1393
Nonpolar Optical Deformation Potential DtK (eV/cm) 12e8
Nonpolar Optical Phonon Temperature (K) 989
Density ? (g/cm3) 3.2
Velocity of Sound ?s (cm/s) 13.73e5
M Valley Longitudinal Effective Mass ml 0.29
M Valley Transverse Effective Mass mt 0.42
Nonparabolicity of M valley ? (eV-1) 0.4
3.3 Scattering Rate Calculations
Electrons moving under the influence of an applied field accelerate until they
are scattered into a different state from one of numerous mechanisms. Each scat-
tering event corresponds to the emission or absorption of a phonon which allows
for the change in the electron?s momentum. Because phonons are bosons, their oc-
cupancy is determined by Bose-Einstein statistics and the associated distribution
function in Equation 3.22 often appears in the scattering rate formulae for different
45
mechanisms.
1
Nq(x) = (3.22)
ex ? 1
As the electron travels through the crystal scattering, its energy is constantly
changing which affects the number of states available for it to scatter into and
thus the probability of scattering in the first place, making each scattering rate a
function of electron energy. Additionally, the various scattering mechanisms consid-
ered all have different perturbation potentials which correspond to different energy-
dependent scattering rate formulae. For simple formulae, the integrals may have a
closed form which can be trivially implemented into the Monte Carlo code. How-
ever, some of the rate formulae involve complex integrals which must be evaluated
numerically. The calculations can then be stored into a lookup table and referenced
at each iteration to speed up the computation. In the following sections I only
include scattering rate formulas for the dominant mechanisms in SiC.
3.3.1 Acoustic
Acoustic phonons represent acoustic vibration modes of atoms in the crystal.
These phonons are generally lower frequency than other phonons and arise due to
thermal excitations.
The scattering probability for acoustic phonons is calculated to account for
inelastic scattering with ellipsoidal equienergy surfaces and nonparabolic bands. The
derivation of which is given in the landmark paper by Jacoboni and Reggiani [42],
46
and the result shown in Equations 3.23 and 3.24.
1/2 3 2
P ab
m (kBT ) D
A (E) =
d [ A??(E)
?1/2?
25/2?~4v4s? ? ]x2,a x2,a
(1 + 2?E) N (x?)x?2dx? ?q + 2?kbT Nq(x )x
?3dx? (3.23)
x1,a x1,a
1/2
em md (k
3
BT ) D
2
PA (E[) = A?(E)?1/2?25/2??~4v4s? ? ]x2,e x2,e
(1 + 2?E) [N (x? ?2 ? ? ?3q ) + 1]x dx ? 2?kbT [Nq(x ) + 1]x dx? (3.24)
x1,e x1,e
Here, the upper equation is the probability for phonon absorption and the
lower for phonon emission. The limits for integration are derived from the energy and
momentum conserving delta function and are shown in Table 3.2. The associated
integrals are evaluated numerically to take into account the exact form of Nq.
Table 3.2: Limits of integration used in acoustic phonon scattering probability cal-
culations. Here, C(?) = 4Es/(k
2
BT (1? 4?Es)) and Es = mdvs/2
A[[b?sorption ? ] Emission Conditionx1,a = C(?) ?Es(1 + 2?E)? ?] x = N/A ? < Es? 1,e 1?4?Es
x2,a = C(?) [ Es(1 + 2?E) + ?] x[2,e = N/Ax = 0 x = 0 ] ? > Es?1,a ? ?1,e ? 1?4?Es
x2,a = C(?) ? + Es(1 + 2?E) x2,e = C(?) ? ? Es(1 + 2?E)
3.3.2 Non-Polar Optical Scattering
Non-polar optical scattering corresponds to excitations of a polar optical vi-
bration mode of the crystal. These modes are generally excited by interaction with
light.
47
Analytically, the scattering rate for optical phonons is given by Equation 3.25,
a result also derived in the Jacoboni and Reggiani paper [42]. Here, the resulting
equation is formulated to account for ellipsoidal, nonparabolic bands. This scatter-
ing rate is used to account only for intravalley scattering and intervalley scattering
is regarded as small enough to be neglected.
[ ]
(D K)2
(3/2) ?
Pop(E) = ?t
md Nop ?(E ? ~?op) [1 + 2?(E ? ~?op)] (3.25)
2?~3??op Nop + 1
3.3.3 Polar Optical Scattering
Polar optical phonons are excitations of the polar optical modes of a lattice
and occur when different atoms in the primitive unit cell of an ionic semiconductors
vibrate out of phase with each other. This oscillation causes an alternating electric
field and changing potential that can scatter carriers. The analytical expression
shown in Equation 3.26 is taken from [43].
? ( ) [ ]
e2 md?pop 1 1 1 + 2?E ? Nq
Ppop(E) = ? ? ? F0(E,E ?) (3.26)
2~ ?? ?0 ?(E) Nq + 1
Here, E ? = E ? ~?pop corresponding to energy after emission or absorption of
the phonon. The other terms are calculated as:
{ ? ?1/2 1/2 ? }
F (E,E ? ?
?
1 ? (E) + ? (E ) ?
0 ) = C A ln ?? ?+B?1/2(E)? ?1/2(E ?) ?
A = {2(1 + 2?E)(1 + ?E ?) + ?[?(E) + ?(E ?)]}2
B = ?2??1/2(E)?1/2(E ?){4(1 + ?E)(1 + ?E ?) + ?[?(E) + ?(E ?)]}
C = 4(1 + ?E)(1 + ?E ?)(1 + 2?E)(1 + 2?E ?)
48
3.4 Novel Interface Atomic Roughness Scattering
3.4.1 Existing Surface Roughness Model
In this work, I investigate the inclusion of a novel scattering mechanism to
model the most fundamental scattering due to the break in periodicity at an ideal
semiconductor-oxide interface which will be referred to as atomic roughness scat-
tering. Presently, the rather large scattering contribution in SiC at the interface is
caused by surface roughness scattering [44]. The surface roughness is modeled as a
series of bunched steps with size range around 1 to 5 nanometers [27]. Steps form
due to the presence of a crystallographic miscut of around 4-8? which enables the
growth of homoepitaxial layers in the c-direction without the need for the presence of
screw dislocations [45]. The steps in conjunction with the gate-induced surface field
perpendicular to the interface cause a perturbation to the potential that electrons
in the channel see as defined by Equation 3.27.
d?(z)
??(x, y, z) = ?z(x, y) (3.27)
dz
Here, the model assumes the simple form of the surface perpendicular field
times the height of the steps. From the form of the perturbation, we can see that at
higher surface fields the level of interaction and thus the rate of scattering increases.
At high surface fields this mechanism is predicted to dominate in 4H-SiC MOSFETs
[27]. In the typical implementation of this mechanism, a roughness power spectrum
49
is introduced, and assumes the form given by [44]:
[ ??2L2S(q) = ] (3.28)
1 + q
2L2
2
In this standard but relatively simplistic model, ? represents the average step
amplitude and L is the correlation length or the standard deviation of the step
separation. This power spectrum corresponds to an exponential autocovariance
function given in Equation 3.29 [44, 46] and the transition rate in Equation 3.30.
The values of ? and L are fitted to experimental data by [27] to be 3.5nm and 7nm
respectively and account for the distribution of nano and macro-scale steps.
<? = ??(r), ?(r? r?)? = ?2e?r/L (3.29)
2?2e(2?2L2? F 2?)(z)S(k, k , z) = ?(E ? E ?) (3.30)
~ q2L21 +
2
3.4.2 My Perturbation and DFT-Based Model: Calculation and Scat-
tering Rate Results
In our investigation of surface roughness scattering, we have decided to study
the fundamental scattering for an ideal 4H-SiC/SiO2 interface that exists despite
the presence of a miscut. Scattering arising due to the miscut can be integrated into
the theory later as an additive potential. In our model, inversion layer electrons
have their energies quantized into subbands of a quantum well near the interface so
that motion is restricted to the plane of the interface. The scattering mechanism is
treated using the well-known perturbation potential form given earlier in Equation
3.16. In this work, we developed a novel technique for extracting the actual form of
50
this potential from a series of DFT calculations, rather than an approximation, and
should be applicable to any conceivable interface. To extract the perturbation, we
first perform a relaxation calculation on a 4H-SiC/SiO2 interface supercell (Figure
3.5) to determine the energetically favorable atomic configuration. Then we perform
another relaxation of a pure 4H-SiC bulk structure. The converged DFT calculations
result in a value of the potential everywhere inside each of the supercells. We
take the potential field from each calculation and align the bulk potential with
the potential far away from the interface, then by subtraction we determine the
potential difference the electron experiences due to the presence of the interface.
The perturbation potential function varies in the plane of the interface based on the
location of the bulk atoms and the oxide atoms. The perturbation decays to zero
with depth into the bulk of the semiconductor as the effect of the presence of the
interface on the periodic bulk potential diminishes.
Figure 3.6 shows the extracted potential at the interface which we take to be
the average height of the first oxygen layer. The potential here contains all con-
stituent potential components found in the Kohn-Sham Hamiltonian: ionic, Hartree
(electron-electron), and exchange-correlation.
Subtracting this interface potential from the periodic bulk potential, we obtain
a perturbation potential, shown in Figure 3.7. The deep wells in the figure, shown in
blue, correspond to the oxygen atom locations. The regularly spaced yellow potential
barriers are the inverted carbon atom wells, showing where the carbon atoms of the
51
Figure 3.5: 4H-SiC/SiO2 interface supercell created by [9], used to extract calculated
potential.
bulk would be if the bulk continued past the interface. Since the scattering rate
depends on the square of the perturbation potential, the choice of sign is arbitrary.
To evaluate the matrix element given in Equation 3.8, we need to assume a
form for the initial and final states of the electron. Because this scattering is occur-
ring in the channel at the surface of a MOSFET device, there exists a mesoscopic
scale potential well due to band bending near the semiconductor surface and the
presence of a band offset between the semiconductor and insulator. The shape of
the well is regularly approximated in literature to be triangular and the channel
52
electrons which reside there are confined only in the z dimension towards the bulk
crystal [44,47]. The electrons in the well act similarly to a particle in a box, residing
in discrete energy states known as subbands which occur on top of the crystal band
structure, shown in Figure 3.8.
From these assumptions, we use separation of variables to break the Schrdinger
equation into an equation with only x and y dependence, and an equation with
purely z dependence. The resulting x and y dependent equation takes the form of
a free electron equation so the electronic wavefunction in the x,y (interface) plane
is taken to be a plane wave. In the z dimension, the standard analytical solution
for a particle in an infinite triangular well is used [47]. The solutions are known to
consist of Airy functions (Ai[z]), shifted to contain the appropriate number of zeros
for each subband energy level [48]. The full form of the wavefunction is shown in
Equation 3.31 and associated subband energies (referenced to the conduction band
energy EC) in Equation 3.32.
1
[ ? = ? ei(kxx+kyy)(n ) ( ?n(z) )] (3.31)? (1/3)
2m?eF En
?n(z) = An ? Ai ? z ? z ? ?(z ? z)~2 i ieF
An = ?( ) [1 ]~2 (1/3) ? ( ( ))? Ai ?( 3?? n?) 1
(2/3)
2m eF 2 4
1 ? t3
A(i[?] = ) co(s (+ ?t )dt? 0 3
2~2 2 (1/3)
)(2/3)
e F 3? 1
En = ? n? (3.32)2m 2 4
Here, ? is the interface area, kx and ky are the 2d wavevector components of
53
the electron in the plane of the interface, A is for Airy function normalization, F is
the surface electric field value, zi is the interface location, n is the subband index,
and ? is the Heaviside function. The wavefunctions here are the solutions obtained
for an infinite triangular well of slope F with an infinite conduction band offset. The
actual offset for a 4H-SiC/SiO2 system is experimentally determined to be 2.7eV,
which is quite large compared to the subband energies for typical surface field values
(E1 = 0.18eV@F = 1MV/cm) so the approximation is quite reasonable.
Substituting our assumed form of the wavefunctions into Equation 3.15 results
in the transition rate S(k,k?) from state k to state k? shown in Equation 3.33. The
delta function in this equation contains the total energy conservation condition and
assumes that the energy of the phonon in the exchange is negligible, as is typically
done for deformation potential interactions. To obtain the total scattering rate
W2D(kq, n) out of state k, we sum over all eligible final states k? as shown in Equation
3.34.
?? ?
2? ? ?? ?? ? ?2
S(k , k?, n, n?) = ?e?iq q ((kxx+kyy)??n?(z)?V (x, y, z)ei(kxx+)kyy)?n(z)~ ?? ? zi
~2 |k?|2 ~2q ? |k
2
q|
? ? + En? ? ? En dzdydx (3.33)2m 2m
???
W (k) = S(k ?q, kq, n, n
?)dk?q (3.34)
(2?)2
n?
The delta function in Equation 3.33 which enforces energy conservation in turn
imposes allowed magnitudes of k?q, shown in Equation 3.35, depending on the initial
54
and final subband energies. ?
2m?
k?q = (kq)
2 + (E
~2 n
? En?) (3.35)
After computing the transition rate for a given initial state kq, n, we see that
the allowed final states lay in concentric rings in the kx, ky plane. This happens
because when the electron scatters to a higher subband, it gains potential energy
and must lose kinetic energy to conserve total energy, thus its |kq| must shrink.
Conversely, when the electron scatters to a lower energy subband, it loses potential
energy and gains kinetic energy, creating the larger rings. It is also possible for
electrons to have kinetic energy too small, such that the energy required to scatter
to the next higher band is larger than the amount of kinetic energy available to
be lost, making this a forbidden transition. Figure 3.9 shows the transition rates
which are calculated for all allowed k?q, given the initial state in the second subband
(n = 2) with initial kinetic energy of 1eV.
To determine the total scattering rate out of our given initial state, we simply
integrate the transition rates over the rings of allowed states, and then sum over
all rings. The calculation has been performed for various surface fields and shows
the 2D atomic roughness scattering rate as a function of electron energy in Figure
3.10. The discontinuities in the scattering rate arise when the electron gains enough
energy to access each higher subband which increases the scattering probability by
increasing the number of available states. It is worth noting that the scattering rate
calculated with the traditional combined empirical surface model is proportional
to the square of the gate field at the interface (F 2?) [39, 44, 49?52], and the same
55
holds for my first-principles calculated scattering rates. This result is confirmed by
plotting the log of the scattering rate against the log of the field, shown in Figure
3.11, where the exponent is given by the slope of resulting straight lines. From the
results of my scattering rate calculations, we see that the intrinsic atomic-roughness
scattering rates come out to be between 1?1011 and 5?1011 per second for a range
of typical gate-induced surface fields from 0.1 MV/cm to 1.0 MV/cm.
3.4.3 Comparison of Intrinsic and Extrinsic Mobilities
Component mobilities due to different mechanisms add according to Matthiessen?s
rule. In Equation (3.36) total mobility ?Tot is separated into its intrinsic (?int) and
extrinsic (?ext) components. The intrinsic mobility is itself comprised of material-
dependent mobilities such as the bulk (?B), surface phonon (?SP ), and atomic-scale
surface roughness (?AR). Extrinsic mobility components are associated with process-
dependent issues such as the Coulombic scattering from interface traps (?C), and
the surface roughness (?SR) mainly due to nanometer-scale interface steps resulting
from a crystallographic miscut included to promote homoepitaxial growth. The bulk
mobility term is included with the intrinsic terms because we are focused on inter-
face issues, independent of the substrate doping. Note that in previous works the
CESRM scattering approach inherently lumps together all non-Coulombic interface
effects, and we calculate this value for comparison. The extrinsic miscut steps and
the bunching thereof are the largest roughness features present in calculations fit to
experimentally measured interfaces, so to a good approximation, we consider this
56
Table 3.3: Important component mobilities at the SiC/SiO2 interface for T=300K
? cm2/V s E? MV/cm
0.1 0.5 1
?B 1070
?SP 1470 560 400
? 4.4? 107AR 1.8? 106 4.8? 105
?int 620 370 290
?SR 350 14 4
?C 429 447 591
?ext 193 14 4
?Tot 147 13 4
term to be purely extrinsic [51].
1 1 1
= + = (3.36)
?Tot ?int ?ext
1 1 1 1 1
( + + ) + ( + + ? ? ? )
?B ?SP ?AR ?C ?SR
To show their relative contributions to the total mobility, the other predomi-
nant scattering mobilities present at the SiC/SiO2 interface are quantified and then
compared. These mobilities are presented in Table 3.3 and consist of values cal-
culated based on models and data from Potbhare et al. [1]. The details of these
calculations are shown in Appendix C. Note that the total mobility measured ex-
perimentally is highly dependent on oxide quality, the interface state density and the
surface roughness profile. The values calculated here are representative of a typical
device and are also consistent with the component mobilities found by Uhnevionak
et al. [52]. In table 3.3, intrinsic mobilities are separated from extrinsic mobility
components.
57
For the highest gate fields we used, the total intrinsic mobility ?int becomes
approximately 290 cm2/Vs, significantly higher than the corresponding extrinsic
mobility ? 2ext of 4 cm /Vs. At lower fields, the intrinsic mobility is approximately
620 cm2/Vs compared to the extrinsic mobility of 193 cm2/Vs. The important ob-
servation here is that by reducing the mobility limiting effects of the extrinsic terms,
mobility improvements on the order of 4 to 70 times can be expected, depending
on the field. These benefits, however, hinge on future improvements in fabrication
techniques with special care taken to minimize the extrinsic scattering effects. One
key method of reducing extrinsic surface roughness is likely to come from eliminat-
ing the 4? to 8? Si-face miscut which gives rise to surface steps and step bunching.
Though the practice of including a miscut is ubiquitous for growing homoepitaxial
layers in industry, it is expected to be a major source of mobility limiting surface
scattering, especially at high gate fields.
On-axis homoepitaxial growth of SiC has been a topic of research with some
successes, in particular the C-face (0001) of 4H-SiC seems to be the easiest to control
the polytype [53], though success has also been seen on the Si-Face (0001) where
spiral growth is mediated by screw-dislocations [54,55]. These on-axis growth tech-
niques are also shown to reduce the basal plane dislocation density and improve
epitaxial layer quality [53,54] but island growth using these techniques may prevent
extrinsic surface roughness from being eliminated completely. Other works have
investigated the effect of transport versus miscut angle and also found that mobility
generally improves when miscuts are reduced. While the more conventional [56] Si-
face is considered in this work, experiments have been done with electron transport
58
along on the C-face [56,57] and a-face (1120) [58, 59].
In the C-face experiment by Fukuda et al. [57], higher surface roughness was
confirmed to exist by TEM images for larger miscut angles. MOSFETs fabricated on
substrates with a vicinal 0.8?off-angle cut showed peak field-effect channel mobility
of 92.3 cm2/Vs compared to 84.5 cm2/Vs for 8?off-angle devices [57]. Hijikata et
al. [56] also concluded from C-face experiments that smaller off-angle devices should
lead to better device performance.
Harada et al. [58] performed mobility measurements on trench MOSFETs
fabricated on 0.7? C-face vicinal wafers where the channel plane was aligned at
various angles relative to the a-face. A peak mobility of 160 cm2/Vs was achieved
when the channel was aligned to the the a-face, compared to the lowest mobility of
110 cm2/Vs angled half way between a-faces.
Experiments were also performed on trench devices fabricated on wafers with
a 8?surface miscut-angle. The deviation from vertical trench sidewalls compounds
with the wafer miscut angle resulting in the channel plane deviating further from
the exact a-face plane. Under these conditions Harada found the peak mobility to
reduce from 160 cm2/Vs to 70 cm2/Vs [58]. Yano et al. [59] found similar results,
reducing their peak mobility from 66 cm2/Vs to a low of 6 cm2/Vs for the greatest
channel misalignment.
The aforementioned experimental work argues that channel misalignment to
important crystallographic planes reduces channel mobility. These results are con-
gruent with our calculations which predict much higher mobility if miscut-induced
surface roughness is reduced.
59
3.5 Conclusion
In this work I classified the traditional mobility terms into intrinsic and extrin-
sic components, where the intrinsic are purely due to inherent material properties,
and the extrinsic are due to electron scattering arising from device fabrication. In
particular, the intrinsic components are separated from the extrinsic components
of surface roughness mobility, which have previously been computed by lumping
these components together. The intrinsic surface roughness is due to atomic-scale
potential-variation and the presence of the interface, while the extrinsic is largely
due to the wafer miscut generating nanometer-scale steps and step bunches.
I found that the intrinsic mobility of the 4H-SiC/SiO2 interface is much larger
than the extrinsic mobility. This is especially apparent with respect to surface
roughness scattering, where extrinsic mobility due to miscut roughness is orders of
magnitude lower than intrinsic surface mobility due to atomic-scale roughness. This
result is found to be largely insensitive to the magnitude of the perturbation, where
even if the perturbation is doubled, the extrinsic surface roughness dominates the
degradation of mobility. This result argues for further developing on-axis (or close
to on-axis) MOS technology for SiC. Experimental results support this conclusion
as well [56?59].From these experiments, evidence is clear that deviations from ideal
crystallographic planes create longer-range disorder and increased surface roughness.
In summary, we conclude that the interface of SiC/SiO2 does not intrinsically lead
to low mobility and that process-induced imperfections are the cause of the interface
mobility degradation experimentally observed.
60
V [eV]
0 -50
-100 -100
-200 -150
-300
-200
10
5
5
0 -250
y [nm] 0 -5 x [nm]
Figure 3.6: Interface potential taken at the semiconductor-insulator interface. The
blue (well) areas depicted are within the oxygen atoms.
61
?V [eV] 200
100
500
0
0 -100
-200
-500 -300
-400
-1000 -500
10
5 -600
5
0 -700
y [nm] 0 -5 x [nm]
Figure 3.7: Interface perturbation potential taken at the semiconductor-insulator
interface. The blue (well) areas are the oxygen atoms, and the yellow (peak) areas
are where the carbon sites are located if the bulk continued beyond the interface.
62
Figure 3.8: Real-space band extrema cartoon depicting the well structure which
forms at the SiC-SiO2 interface when a channel is formed. The subbands are enu-
merated 1, 2, 3... and exist as states found only within the well near the surface.
On the right, the triangular approximation to the well is made and the correspond-
ing approximate wavefunctions are shown which increase in number of nodes with
energy.
63
Figure 3.9: Calculated transition rates for electrons due to atomic-roughness scat-
tering with DFT-extracted interface perturbation potential. Initial state has 1eV of
kinetic energy and is in the second subband of a triangular well with surface field of
1MV/cm. Scattering within the second subband is shown in red, scattering to the
third subband in orange, and to the first subband in blue.
64
Figure 3.10: 2D atomic roughness scattering rates for electrons in the 1st (blue), 2nd
(red), and 3rd (orange) subbands. The x axis shows total energy of the electron, so
each curve begins at the corresponding subband energy for the given field. Insets
show the relative shape of the triangular well that each plot represents. Bottom
right plot contains both energy and field variation.
65
Figure 3.11: Field-dependent scattering within the 1st subband for electrons with
various kinetic energies.
66
Chapter 4: Incomplete Ionization in p-Type SiC
4.1 Background and Introduction
This work is motivated in part to further the understanding of p-type SiC
devices, and ultimately to aid in the development of SiC CMOS devices. Presently,
more attention is being paid to making advances in n-type SiC devices. In this work,
I sought to investigate issues related to p-type SiC. In particular, I analyze the de-
gree of incomplete ionization which becomes an important effect at higher doping
concentrations that exist in p-type and CMOS devices. My major contributions in
this area include aggregating and analyzing almost all existing published resistivity,
Hall mobility, and ionization energy data on Al-doped 4H-SiC, and using this infor-
mation to verify my physics-based formulation of the incomplete ionization. I also
extend this incomplete ionization fraction calculation to higher temperatures and
developed an easy-to-use parameterization of the temperature and concentration-
dependent ionization fraction by using a genetic algorithm. As a small aside, I have
also used the genetic algorithm to optimally parameterize the experimental Hall
mobility data over a larger range of temperatures and doping concentrations and
using more data than those that are currently published.
Dopant atoms are used in semiconductor devices to add mobile electrons (using
donors) and mobile holes (using acceptors) which act as charge carriers. To simplify
67
device simulations/calculations, the assumption is sometimes made that all dopant
atoms are ionized. As a consequence, this assumption implies that all acceptors
have captured an electron from the valence band (NA = N
?
A ) and all donors have
donated an electron to the conduction band (N = N+D D ) and as such, each dopant
atom contributes one charge carrier (barring electron-hole recombination effects and
assuming this is an uncompensated semiconductor). In a real physical system, the
energy needed to ionize these dopant atoms mainly comes from thermal energy
which follows a Maxwell-Boltzmann distribution under no applied field. Based off
this distribution, some portion of the dopant atoms will inevitably remain unionized
due to the non-zero probability of electrons having energy less than the ionization
energy. The degree of incomplete ionization depends on not only the temperature,
but also on the ionization energy of the state, and the conduction/valence band
density of states for donors/acceptors respectively.
For the case of p-type 4H-SiC, unfortunately, all acceptor dopant species have
large ionization energies. As a result of this fact, a greater degree of incomplete ion-
ization is expected in p-type SiC than in its n-type counterpart or in other materials
like Si and Ge [60]. Consequently, this effect must be treated more carefully in SiC
in order to accurately model mobile hole densities for given doping and temperature
conditions.
Of the acceptor dopants, aluminum (Al) is arguably the most favorable as it
has the highest solid solubility [10, 61?66] and lowest ionization energy at around
190 ? 230meV [10, 65?72]. The other two main acceptor dopants are Ga which
has an ionization energy of approximately 250 ? 270meV and B with two levels -
68
one around 300meV from Si lattice site occupation and one around 650meV due to
defect-related complex known as the D-center [10, 65?70, 73?79]. SiC doping levels
are illustrated in Figure 4.1.
Figure 4.1: Doping levels in 4H-SiC.
For comparison, the most common n-type dopants in SiC also introduce two
levels each: 40? 70meV and 90? 130meV for nitrogen (N) as well as 40? 50meV
and 80? 100meV for phosphorous (P) [68,71,72,79?88]. P-type dopants in Ge have
an ionization energy of merely 10meV and for Si 45?73meV [10,68,89?91]. Various
dopant levels in Ge, Si, and GaAs are shown in Figure 4.2.
69
70
Figure 4.2: Doping levels in Ge, Si, and GaAs.
4.2 Standard Incomplete Ionization Formulation
In materials like Si, where incomplete ionization is often neglected, if it is
included, it is generally done so in an approximate manner to ease computational
complexity. The standard approximations made usually include using a single dis-
crete impurity energy level (generally also treated to be independent of the doping
concentration) and applying Maxwell-Boltzmann quantum-state occupation statis-
tics in lieu of the true Fermi-Dirac statistics for fermions. These approximations
become worse for higher doping concentrations when the dopant wavefunctions start
to overlap and the semiconductor approaches degeneracy.
Under these assumptions we can calculate the ratio of mobile holes to ac-
ceptor atoms p/NA by invoking the charge neutrality condition and manipulating
the resulting equation to remove the dependence on the Fermi level (an unknown
quantity) and obtain an algebraic expression.
To begin, we will first make assumptions about the physical system which
result in a single impurity level. This means we treat the dopant atoms as far apart
and non-interacting i.e. no wavefunction overlap. This assumption is better for
low-doped samples because the average distance between nearest-neighbor dopant
(?1/3)
atoms goes as NA . The discrete dopant energy levels form an acceptor density
of states in the bandgap, which, at sufficiently low doping concentrations, takes the
approximate form of a delta-function at the dopant energy level EA.
?i(E) = NA?(E ? EA) (4.1)
71
This density of states must integrate to the total acceptor density NA because there
are NA atoms and thus NA total discrete energy states. We will see later that this
acceptor density of states function spreads out and becomes more complex when
more physical and quantum-mechanical effects are accounted for.
Now we can compute the total number of ionized acceptors by multiplying the
acceptor density of states by the occupancy probability function, then integrating
over all energies. For the occupancy of acceptor states, we use a modified Fermi-
Dirac function which includes the degeneracy gA = 4 of acceptor state occupancy.
The four-fold degeneracy comes from the choice of spin-up or spin-down electron
occupying the acceptor state, as well as the choice to create a heavy or light hole in
the valence band due to the degeneracy at the valence band maximum (?-point) in
4H-SiC. This is described in further detail in section B.2.
? ?
N?A = ?i(E)fA(E)dE (4.2)
??
1
fA(E) = ( ) (4.3)
1 + g exp E?EFA kBT
Carrying out the integration:
? ?
? NA?(E ?(EA)N = )A dE (4.4)
?? 1 + gA exp
E?EF
kBT
N?
NA
A = ( ) (4.5)
1 + gA exp
EA?EF
kBT
The analogous calculation is also performed to calculate the number of holes
in the valence band, except we use the valence band density of states ?v(E) defined
for ?? ? E ? EV . The 1 ? f(E) term comes from the fact that holes stem from
72
unoccupied valence band states rather than filled states.
4?(2m?)3/2?p
?v(E) = E3 V ? E (4.6)h
1
f(E) = ( ) (4.7)
? 1 + exp E?EFkBTEV
p = ??v(E)(1? f(E))dE (4.8)?? ?
4?(2m?)3/2 EVp EV(? Ep = )dE (4.9)
h3 ?? 1 + exp EF?E
kBT
In this standard formulation, we approximate the Fermi-Dirac occupancy dis-
tribution f(E) with the Maxwell-Boltzmann which is a valid when EF?EV & 3kBT
i.e. non-degenerate doping scenarios when the Fermi level is at least a few kBT away
from the valence band edge. This will allow us to obtain a closed-form solution to
an otherwise non-analytic expression.
( )
? 1( ) ? E ? EF1 f(E) = exp (4.10)
1 + exp EF?E kBT
kBT
4?(2m?)3/2
? EV ? ( )
? p ? E ? EFp EV E exp dE (4.11)
h3 ?? kBT
Equation 4.11 has an exact solution which can be written using NV , the effec-
tive density of states in the valence band which depends on the hole effective mass
in the valence band m?p: ( )
p ? EV ? EFNV(exp ) (4.12)kBT
2?m? (3/2)pkBT
NV = 2 (4.13)
h2
Now we must apply the charge-neutrality condition. For our system, we as-
sume no donor counter doping, and since intrinsic carrier concentration ni is so small
73
in 4H-SiC, we can also neglect the equilibrium electron concentration.
p+N+D = n+N
?
A (4.14)
n2
p = i +N?A (4.15)p
p ? N?A (4.16)
The resulting charge neutrality equation states that the hole concentration
must equal the ionized acceptor concentration. This result stands to reason as
the majority of holes are created by promoting electrons from the valence band
to the acceptor states rather than to the conduction band due to the much larger
energy difference. In general, both of these terms are functions of the unknown
Fermi-level EF , as exhibited in Equations 4.5 and 4.12. By making the Maxwell-
Boltzmann occupation statistics approximation (low doping limit), however, the
governing equations can be manipulated to explicitly remove the EF dependence and
provide a closed-form solution of only known quantities. We can rewrite Equation
4.5 in terms of the hole concentration given in Equation 4.12, and apply the charge
neutrality condition in Equation 4.16 leaving us with a quadratic equation for the
solution of the hole concentration p.
? ( NANA = ) ( ) (4.17)
1 + g exp EV ?EFA exp
EA?EV
kBT kBT
N?
NA
A = p = ( ) (4.18)
1 + gAp exp EA?EV
NV kBT
NA
p = (4.19)
1 + p
2?
74
Finally solving the quadratic and introducing the auxiliary variable of known
quantities ?:
?
p = ?? + ?(2 + 2?N)A (4.20)
NV ?EA
? = exp ? (4.21)
2gA kBT
?EA = EA ? EV (4.22)
As a reminder, NA is the acceptor doping density, NV is the valence band
effective density of states, ?EA is the acceptor ionization energy, and gA is the
acceptor state degeneracy equal to 4 in 4H-SiC.
For other materials where ni is not small (or for very high temperatures), one
must solve a cubic equation to calculate the hole concentration:
p3 + 2?p2 ? (2?N 2 2A + ni )p? 2?ni = 0 (4.23)
For those curious, this has the closed form solution:
?
1 [? ](1/3) 3? 2 x 2?p = 4x3 + y2 + y ? [
3 3 2 3 ? ] ? (4.24)(1/3)4x3 + y2 3+ y
x = ?6NA? ? 4?2 ? 3n2i (4.25)
y = ?36N 2A? ? 16?3 + 36?n2i (4.26)
A comparison of the quadratic solutions (dashed) to the full cubic solutions (solid)
are shown in Figure 4.3. The more accurate cubic solutions have three characteristic
regions. For the lowest temperatures, we see the freeze-out region where incomplete
ionization occurs due to the lack of thermal energy needed to ionize all of the ac-
ceptor atoms and create mobile holes. In the middle temperature region, the mobile
75
hole concentration equals the total acceptor doping and remains fairly constant.
This region of stability extends over a larger temperature range for higher acceptor
doping. This is the region where devices should normally aim to operate within.
Finally, at extremely high temperatures, the intrinsic carrier concentration starts
to dominate and the material is said to be intrinsic because the hole concentration
equals the intrinsic carrier concentration independent of the doping.
Plotted in a different manner, Figure 4.4 shows the plot of p/NA which shows
incomplete ionization as a fraction. From this plot we can also determine as general
trend that incomplete ionization occurs more severely for higher doping concentra-
tions. This can be a problem for MOS-type devices which often have body/channel
doping in the range 1016cm?3 to 1018cm?3 where the conductivity of that region
could be significantly less than predicted assuming complete ionization.
For degenerate doping, some issues arise that this simplified model can not
handle. Most importantly, the Maxwell-Boltzmann statistics will overpredict occu-
pancy of valence band states because it allows more than one electron to occupy a
single spin-state which is forbidden for Fermions. To solve this problem we need to
use the full Fermi-Dirac distribution in the hole concentration calculation. Addi-
tionally, as doping increases, so too does the number of mobile holes in the valence
band. These holes, however, are not uniformly distributed throughout the crys-
tal. Instead, they will arrange themselves preferentially around acceptor ions via
Coulombic attraction, leading to an increasing amount of screening of the acceptor
atom?s potential [92]. This, in turn, leads to a lower ionization energy because the
76
10 18
10 16
10 14
n
i
14
10
14
10
10 12
16
10
16
10
18
10
18
10
10 10
0 500 1000 1500 2000 2500
T [K]
Figure 4.3: Hole concentration calculated using the standard incomplete ionization
model for variousNA doping. Dashed lines show the quadratic solution which ignores
ni and the solid lines show the cubic solution which includes it. Calculated for 4H-
SiC with ?EA = 200meV.
average depth of the potential well decreases, as shown in Figure 4.5. This doping
concentration dependent ionization energy can also be included into the simplified
incomplete ionization model.
Increasing doping concentration also leads to interactions among dopant atoms
?1/3
because the average distance between atoms shrinks as NA . Once these atoms
become close enough such that their valence state wavefunctions overlap, the degen-
eracy of their individual ionization energies must be lifted - leading to a quantum-
mechanical broadening of the acceptor density of states. This effect combined with
the lowering of the average ionization energy due to screening can lead to overlap
77
-3
p [cm ]
1
0.8
0.6
14
10
0.4 1510
16
10
17
10
0.2 18
10
0
0 200 400 600 800 1000
T[K]
Figure 4.4: p/NA calculated for different doping concentrations as a function of
temperature for 4H-SiC with ?EA = 200meV.
of the acceptor states with the valence band. Lastly, the as the doping concentra-
tion increases the material as a whole becomes more disordered. These randomly
distributed dopants will randomly form clusters leaving local regions of greater-
than-average density of acceptors and other regions of less-than-average density.
The Coulombic charge created by these ionized atoms causes a shift in the local
potential which will change with position about its equilibrium mean value. This
locally fluctuating potential which changes on a mesoscopic scale causes the local
band structure to shift above and below its equilibrium state. Because the density
of states is a value averaged over the entire volume of the crystal, this effect leads to
a smearing of the states. Finally, at high enough doping concentration, the acceptor
78
p/N
A
Hydrogenic
Screened
r
Figure 4.5: Coulombic potential of impurity atom core. the solid line is the bare
hydrogenic 1/r potential and the dashed line is the screened potential under the
effects of mobile carriers redistributing themselves.
79
V(r)
density of states can evolve into a quasi-continuous band which supports a conduc-
tion mechanism parallel to valence band transport. This type of conduction is often
referred to as impurity conduction and includes various mechanisms: Mott conduc-
tion, Anderson conduction, nearest neighbor hopping, and variable range hopping.
These effects will be included in our advanced theoretical model of incomplete ion-
ization in p-type 4H-SiC.
It is also possible for dopants to substitute inequivalent lattice sites which have
different ionization energies due to the different forms of the local lattice potential.
In 4H-SiC, Al preferentially substitutes Si atoms in the lattice but there are two
inequivalent Si atom sites which are denoted as the hexagonal h and the cubic k site.
For visualization of this effect, refer to Figure 4.6. We know that all SiC lattices are
formed by combining arrangements of tetrahedrally bonded Si and C atoms. Each Si
is bonded to 4 C atoms and vice-versa. For 4H-SiC when two tetrahedra are joined
in an interpenetrating manner (vertex to center), there are two arrangements which
occur. The hexagonal site is formed when the configuration makes a more mirror-
symmetric structure where the external atoms lie on the vertices of triangular prism.
The cubic site is formed when one tetrahedra is rotated 60?about the connection
axis leaving the external atoms to lie on the vertices of a triangular antiprism. The
number and type of inequivalent lattice sites present in a SiC crystal depends on its
polytype. For 3C-SiC, every site is a cubic site and in 6H-SiC, there is one hexagonal
site and two different cubic sites as shown in Figure 4.7.
To include the effects of different dopant site ionization energies, the total
80
Figure 4.6: Local arrangement of atoms leading to inequivalent lattice sites in 4H-
SiC.
Figure 4.7: Inequivalent sites (h=hexagonal, k=cubic) labeled in SiC polytypes.
The black circles represent carbon atoms and the white circles represent Si atoms.
ionized dopant concentration is split into a sum. The fraction of atoms in each
doping site is denoted with fh and fk and are usually taken to be
1 . Equation 4.27
2
81
also leads to a cubic solution for the hole concentration.
fhNA fkNA
N?A = ( ) + ( ) (4.27)
1 + g exp EAh?EFA 1 + g exp
EAk?EF
kBT A kBT
fhNA fkNA
p = +
1 + p ( 1 + p ) (4.28)2?h 2?k
NV ?EAx
?x = exp ? (4.29)
2gA kBT
[? ] fk = 1? fh (4.30)?1 (1/3) 3 2 x 2(?h + ?k)
p = ? 4x3 + y2 + y ?
3 [? ] ? (4.31)3 2 3 (1/3) 34x3 + y2 + y
x = 6(fhNA(?k ? ?h)?NA?k + 2?h?k)? 4(?h + ?k)2 (4.32)
y = 36(fhNA(?
2 ? ?2k h)?NA?2k) + 72N 2 2 3 3A?h?k + 24(?h?k + ?h?k)? 16(?h + ?k)
(4.33)
For materials where there are more than two inequivalent lattice sites, equation 4.27
can be easily extended by adding additional terms for each inequivalent site leading
to increasing polynomial degree solutions. For polynomial degree greater than 3,
a numerical solution must be obtained to solve for the acceptor concentration and
care must be taken to ensure that the ?positive? physical solution between 0 and NA
is the one obtained. This is applicable in n-type 6H-SiC where nitrogen incorporates
into three different sites: one hexagonal site and two inequivalent cubic sites where
fh = f
1
k1 = fk2 = leading to a quartic equation for p.3
This modification can also be used to model systems where multiple acceptor
dopant species are present such as p-type 4H-SiC doped with a deep B well and a
shallow Al adjustment at the surface. The ionization energy difference between Al
and B is quite large (> 100meV ) so the solution to the cubic equation may differ
82
significantly from the solution of either dopant species individually depending on
the ratio of their concentrations. For the situation of different doping species, the
fraction f must be solved for by using the individual doping concentrations NAl and
NB: fAl = NAl/(NAl +NB) and fB = 1? fAl.
4.3 Doping-Dependent Ionization Energy
For a more accurate calculation of incomplete ionization, one must include the
effects of doping-dependent ionization energy which lowers the energy required to
ionize the dopant atoms as the doping concentration increases. This effect turns
out to significantly affect the resulting p/NA ratio for high doping, and in particular
semiconductors like SiC where the low-doped ionization energy is quite high.
The doping-dependent ionization energy can be experimentally measured through
various techniques including CV, thermal admittance spectroscopy, Hall and 4-point
resistivity fits, donor-acceptor pair luminescence, free-carrier concentration spec-
troscopy, Raman spectroscopy, deep-level transient spectroscopy, IR absorption, etc.
The experimental data from the literature confirms that increasing doping concen-
tration lowers the ionization energy, by moving the acceptor density of states towards
the valence band. This effect is caused by the increased screening from mobile holes
in the valence band which distribute themselves such that the Coulombic energy of
the system is decreased, resulting in a lower acceptor ionization energy [92]. The
same effect is explained for the case of mobile electrons screening donor states in Si
in other works [92?95].
83
To develop an empirical model of the doping-dependent ionization energy in
Al-doped 4H-SiC, a comprehensive review of experimentally measured ionization
energies was first performed. The experimentally determined energies from the lit-
erature are presented in Figure 4.8. The spread in the ionization energies at any
given doping concentration may be partially explained by the varied techniques used
to extract the energies. Some techniques may be capturing different parts of the ac-
ceptor density of states which may not necessarily be a single discrete energy level.
Other contributing factors to these discrepancies might include different fabrica-
tion techniques and differing degrees of counter-doping among the various research
groups. For this work, the ionization energy is considered to be the mean (center)
of the acceptor density of states. The empirical model of the ionization energy is
extracted from the experimental data by minimizing the least-squares error.
A well-behaved logistical equation is used to parameterize the doping-dependent
ionization energy (Equation 4.34) [96?98] and the experimental data is used to de-
termine the optimal parameters. This model contains three parameters, ?EA0 as
the isolated dopant ionization energy, NE as the reference doping where the ion-
ization energy is half of its isolated value, and a characterizing the degree of the
doping dependence. The ionization energy based on this model is shown in black in
Figure 4.8. This is the parameterized form of ?EA(NA) that is used in the rest of
this work.
?(EA0?EA(NA) = )a (4.34)
1 + NA
NE
84
Figure 4.8: Doping density-dependent fit to experimental acceptor ionization energy
for 4H-SiC doped with Al. The blue line is the Pearson-Bardeen model and the black
line is the logistical parameterization.
85
In earlier literature works, ionization energies may also be fit using the Pearson-
Bardeen model which takes the form ?EA(NA) = ?EA0 ? ?(NA)1/3 instead of the
form in Equation 4.34. The Pearson-Bardeen model is based on the average impu-
rity separation distance (N1/3). The Pearson-Bardeen formula is not used in my
incomplete ionization model because at high enough doping, it predicts that the
ionization energy eventually becomes negative and continues past the valence band
edge. Instead, the logistical expression in Equation 4.34 is used so that the center of
the impurity band asymptotically approaches the valence band edge. This property
gives better control of how the high-doped regions behave so that they can be tuned
with other parameters of the incomplete ionization model in accordance to the ex-
perimental data. The logistical formula may also be potentially more physically
sound because the change in ionization energy is mainly due to screening and thus
should tend towards zero [92]. Figure 4.8 shows the logistical parameterization in
black and the fit using the Pearson-Bardeen parameterization is in blue. Both mod-
els agree well for doping concentrations less than 3?1020cm?3. Later discussion will
address doping concentrations at and above this level, where a different mechanism
dominates the p/NA fraction such that the exact ionization energy in this region is
not as important.
4.4 Acceptor Density of States
As mentioned earlier, the acceptor density of states is not always best de-
scribed by a single energy level, but instead by a density of states function which
86
is spread in energy space. This is particularly true in the case of higher doping
when the interactions between dopant atoms become significant. When the doping
concentration is increased, the number of dopants in close proximity to each other
also increases, causing greater quantum-mechanical interaction. The overlap of the
dopant-state wavefunctions along with the Pauli-Exclusion principle lifts the degen-
eracy of the single dopant energy levels and broadens them from series of localized
?-function states (all of equal energy) into a quasi-continuous doping band [93, 99].
The width of this band is calculated by adapting work from Lee et al. which was
performed on donor states [100].
In this work, the acceptor density of states spread by wavefunction overlap
?i0(E) is taken as a single rectangle function centered at EA(NA) which integrates
to the total acceptor density NA (with respect to the energy variable E). A single
impurity band is used in this model because the vast majority of experimental evi-
dence reports no significant ionization energy variation (to within analysis resolution
of a few meV) amongst Al substitution of inequivalent silicon lattice sites (hexagonal
and cubic) in 4H-SiC [65, 69, 70, 85, 101?103]. The energy bandwidth (width of the
rectangle function) of the donor states is characterized by the parameter B which is
calculated using the tight binding model, hydrogenic s-orbitals, and assuming that
the dopants are randomly distributed. Based on the work of Lee et al. [100], we
implement the following set of equations for ?i0(E) and B for acceptors in SiC.
87
???
??NA/B, ?
B ? E ? EA(N ) ? BA
? 2 2?i0(E) = ? (4.35)? 0, otherwis? (e )
2 ?4B = 2 J(R)4?NAR exp ?N R3A dR (4.36)
0 3
q2?
J(R) = (1 + ?R) exp(??R) (4.37)
4?r0
Under the tight binding approximation, the dopant band width B is related to
the energy transfer integral J(R) between nearest neighbor dopant states a distance
R apart. J(R) is calculated by using hydrogenic s-orbitals to evaluate the energy
transfer integral. In this calculation, ? = (1/a )(E /E )1/20 A 0 is the scaled inverse
radius for the acceptor state, a0 is the Bohr radius, and E0 is the ground state
energy of the hydrogen atom [91]. By assuming the dopants are uniformly randomly
distributed, the nearest neighbor distance R will follow a Poisson distribution [104].
The bandwidth B is then calculated by finding twice the average value of J(R) and
the band is centered symmetrically about EA [99].
4.5 Disorder Effects and Band Tailing
In addition to the impurity density of states spreading due to wavefunction
overlap, the band is further spread due to the effects of lattice disorder and subse-
quent potential fluctuations arising from the randomly distributed ionized dopants.
At any given point in a doped crystal, the potential will differ from that of a pure
crystal due to the Coulomb potential produced the dopant ions. Random fluctu-
ations within the local concentration of the dopant atoms causes a potential fluc-
88
tuation about its mean value on a mesoscopic scale. A less significant potential
deviation is also created by the local strain induced by the dopant ions and clus-
ters changing the periodic crystal structure, but we will be neglecting this smaller
effect in this work. Depending on the local atomic configuration of the dopant ions,
the potential throughout the crystal and thus the energy reference for the quantum
states will fluctuate about its mean value. [93,94,100,104,105].
These potential fluctuations lead to a distortion or smearing of the states near
the conduction and valence band extrema as the energy reference of these states
varies locally but the macroscopic density of states formulae represent spatial aver-
ages over the entire lattice. If the potential varies slowly with respect to the dopant
wavefunction, the local impurity states are also shifted with respect to this energy
- leading to impurity band tailing [92, 94, 100, 105, 106]. This effect is incorporated
by separately averaging the standard valence band density of states formula and
the rectangular acceptor density of states over the local potential energy distribu-
tion inside the crystal (denoted as P (V )) [92?94, 100, 104, 105]. The averaging is
accomplished by means of convolution integrals, given in Equations 4.43 through
4.46.
89
( )
1 2
P (V ) = ( exp )?
V
(4.38)
(2?)1/2? 2?2
1/2
N?
? = A q4? (4.39)
8?22 2r0
??2(= ??2 ?2h + ?)i (4.40)1/2
r0kBT
( ?h = ) ( )[ ] (4.41)q2p1/2 ?1/3
r0kBT 4 4
?i = + ? ?N2 A (4.42)q p(1 + p/NA) 3 3
The standard deviation of the Gaussian potential distribution is characterized
by the screening length ?. The screening length is itself comprised of hole and ionic
components. The hole screening length comes from the standard Debye length in a
semiconductor. The ionized impurity screening length also comes from the Debye
length but is increased by the average distance between acceptor atoms.
With this modification, the density of impurity and valence band density of
states become ?i(E) and ?v(E) respectively. The resulting density of states looks a
slightly smeared version of the input - causing tail states to appear at the band edges,
as shown in Figure 4.9. The variation of the local potential corresponds to the local
occupancy of the acceptor states causing regions of higher-than-average acceptor
occupancy (local negative charge) and lower-than-average acceptor occupancy (local
positive charge) due to the spatially-independent Fermi energy. Because the local
potential energy of the system affects the local shift in the conduction and valence
bands in the same manner, the local band gap Eg remains unchanged [94].
90
? +?
?i(E) =? ?i0(E ? V )P (V )dV (4.43)??E?EV
?v(E)?= ?v0(E ? V )P (V )dV (4.44)?? ????NA/B, ?
B ? E ? EA(NA) ? B
? 2 2?i0(E) = ? (4.45)0, otherwise
4?(2m?)3/2?p
?v0(E) = EV ? E (4.46)
h3
With these resulting dopant density of states ?i(E) and valence band density
of states ?v(E), the only remaining unknown in my model is the Fermi level EF
which can be solved for numerically and then used to compute the ratio p/NA.
The full visualization of the density of states under various circumstances or
levels of approximation is depicted in Figure 4.10. For the acceptor density of states,
in the case of low doped samples there will be a single energy level at energy EA
which will spread out due to wavefunction overlap as the doping density increases.
Both the acceptor and valence band density of state functions then smear out when
disorder effects begin to appear due to randomly arranged dopant clusters and lattice
disorder at even higher doping concentrations.
4.6 Theoretical Incomplete Ionization Model
The theoretical incomplete ionization model described in this work is derived
from physical principles and includes the effects discussed in the previous sections:
doping-dependent ionization energy, dopant band spreading due to wavefunction
91
Figure 4.9: P-doped semiconductor density of states diagram. Left: Shows the
position-dependent variation of the local band structure (negative of potential).
Right: Dashed blue lines show the standard valence band density of states and the
rectangular acceptor band density of states. After the inclusion of disorder effects,
the resulting density of states with band tailing is shown with solid blue lines. Note:
Image is not to scale - depending on the doping concentration, the doping band will
be significantly closer to and may even overlap the valence band density of states
tail. The doping band width B as well as its height will also vary considerably
depending on the doping density.
92
Figure 4.10: Evolution of the density of states under various conditions: Row a)
Acceptor density of states starts as a single energy level at low doping, spreads due
to wavefunction overlap at higher concentrations, then smears due to disorder Row
b) Valence band density of states smearing due to disorder effects row c) Combined
picture of the valence band and acceptor density of states
93
overlap, and density of states smearing due to disorder. To solve for the ionization
ratio, charge neutrality is invoked in order to relate the hole concentration p to
the ionized acceptor concentration N?A which are in turn both described by integral
equations involving an occupation probability function and a density of states. In
this work uncompensated p-type doping is considered, so N+D is set equal to zero
in the charge neutrality equation. Additionally, because this model is used to treat
4H-SiC, the intrinsic carrier concentration is extremely low (from 300? 800K ni ?
10?8 ? 1010cm?3) so n = n2i /p is assumed to be negligible.
Forms of impurity band conduction have been observed in SiC at high doping
concentrations and in compensated samples [99,107], suggesting the need for a mod-
ification to the standard charge neutrality and incomplete ionization equations. To
include this effect, the total number of mobile holes p is rewritten as a sum of ?nor-
mal? valence band conducting holes pv and holes conducting via an impurity band
mechanism pi. Adopting this new form of p, the charge neutrality equation must also
be rewritten because only the valence band holes contribute positive charge. The
fraction of unionized acceptor states which contribute impurity conducting holes is
quantified by the fraction (1 ? ?). With these modifications, the ratio p/NA given
by Equation 4.51 is no longer technically the ?ionization fraction? because there are
more mobile holes than ionized acceptors due to the additional impurity conduction
holes (although this distinction is not always made).
94
p = pv + pi (4.47)
pv +N+ = nc +N?D A (4.48)
pi = (1? ?)N0A (4.49)
NA = N
0
A +N
?
A (4.50)
p pv N?
= ? + (1? ?) =6 A (4.51)
NA NA NA
Impurity conduction effects have also been observed in Si [97,98,108], which is
where this modification to the incomplete ionization formula was originally applied
[97,98]. In these papers, the parameter ? is obtained using the empirical expression
in Equation 4.52 which is also applied in this work.
1
?(NA) = ( ) (4.52)d
1 + NA
Nb
Here, Nb characterizes the doping level at which the impurity conduction be-
gins to dominate, and d characterizes the rate at which the conduction increases in
the regime of the impurity conduction mechanism.
To ensure that pv is correct in the case of high doping when the Fermi Level
approaches the valence band edge, the full Fermi-Dirac integral (Eq. 4.53) is used to
calculate the valence band hole carrier concentration in stead of using the Maxwell-
Boltzmann approximation. The unknown Fermi level EF is determined by solving
the charge neutrality equation (Eq. 4.55). As mentioned earlier, there is no counter
doping present, and the free electron concentration is neglected. The integral used to
compute the ionized acceptor states (Eq. 4.54) uses a modified Fermi-Dirac function
95
due to the degeneracy gA associated with filling the state in different ways. In the
case of 4H-SiC, gA = 4 due to the valence band maximum exhibiting a heavy-hole
light-hole degeneracy, which allows for the creation of either a heavy-hole or a light-
hole in the valence band, in addition to the choice of spin ?up? or spin ?down? for an
electron filling the acceptor state [65, 66,72,98,102,103,108?118].
? EV
v ?v((E)p = )dE (4.53)
? ?? 1 + exp EF?EkBT?
? ?i(EN = () )A dE (4.54)
?? 1 + gA exp
E?EF
kBT
pv = N?A (4.55)
From the charge neutrality condition and integral equations, we have three
equations and three unknowns pv, N?A , and EF . Plugging in equations 4.53 and 4.54
into equation 4.16, allows for the Fermi level EF to be determined numerically using
trapezoidal integration combined with the method of bisection for a given doping
NA. Once the correct Fermi level is found, we evaluate the integrals in equations
4.53 and 4.54 to calculate pv and N?A . These values are then used in Equation 4.51
to find the ratio of mobile holes to acceptor doping p/NA. For solving a compensated
system (doped with both NA and ND), the slightly modified equations are discussed
in the appendix.
I would like to emphasize for the sake of contrast with the model described
in the next section that the only empirical inputs into this theoretical model are
the logistical parameterization of the doping-dependent ionization energy ?EA(NA)
96
(Equation 4.34) and the impurity conduction parameter ? (Equation 4.52) which
are both extracted to match experimental data.
4.7 Experimental Incomplete Ionization Verification
In the previous section p/NA was calculated from a largely theoretical stand-
point. This section discusses the separate experimental technique used to verify the
theoretically calculated p/NA. This experimental technique described by Altermatt
et al uses measured values of Hall mobility and resistivity over a range of doping
densities from the literature [96?98]. With this data, a second p/NA ratio is ulti-
mately extracted and compared to the results from the theoretical technique. The
core of this experimental method relies upon equating p/NA to a ratio of mobili-
ties ?Cond/?Hall where the Hall mobility is measured directly and the conductivity
mobility is derived from resistivity measurements.
Resistivity data (taken with 4-point measurements) is converted into a cor-
responding conductivity mobility value ?Cond using Equation 4.56. The derived
conductivity mobility is then divided by the measured Hall mobility ?Hall to obtain
experimentally obtained values of p/NA. The hall correction factor r assumed to be
1 in this work [103,119?122].
97
? 1?Cond (4.56)
q?NA
?RH 1
?Hall = = (4.57)
r rq?p
?Cond rp p
= ? (4.58)
?Hall NA NA
Because the fraction ?Cond/?Hall is used to obtain p/NA, each of these mobil-
ities must be representable as continuous functions of the doping NA. To achieve
this, each set of data is separately fit to an empirical function using the least-squares
method. The experimentally-derived p/NA is then compared to the p/NA calculated
using the theoretical model.
For clarity, the theoretical model does rely on a small degree of experimen-
tal input, namely, the experimental ?EA(NA) and ?(NA) which were presented in
the previous section. The experimental data used in the theoretical model is, how-
ever, largely unrelated to the data which forms the experimental model because
the different data sets are measuring completely different physical properties. This
should allow the experimental model to act as a valid check on the result from the
theoretical model.
98
?Hall0
?Hall(NA) = ( ) (4.59)b
( 1 +? NANH )
?Cond0 ?? + ?2 + 2?NA
?Cond(NA) = ( ( )c) + ?Imp (4.60)
1 + (NA NN AC )
NV ?EA(NA)
?(NA) = exp ? (4.61)
2gA ( kBT )
?(NA ?NImp)2
?Imp(NA) = ?Imp0 exp (4.62)
2?2Imp
The Hall mobility is fit using the standard Caughey-Thomas [123] form often
used in literature [71, 72, 112, 114, 115, 117, 122, 124?126]. The fit for conductivity
mobility is derived partially from physical considerations, the details of which are
described in the Appendix (Section B.1.1). With the analytical forms for ?Hall and
?Cond given by Equations 4.59-4.62, their ratio provides an analytical form for p/NA.
The result of this calculation is presented in Figure 4.14 later in this work.
The parameters used in the theoretical and experimental models are presented
in Table 4.1. With the parameters in this table, the Hall mobility, ionization energy,
and the hole concentration can be determined for a given acceptor doping density.
Later in this work I will present a more usable approximate expression for calculating
the p/NA as a function of both acceptor concentration and temperature.
The plots of the physically-motivated empirical fits of ?Hall and ?Cond are
shown in Figures 4.11 and 4.12 along with the respective experimental data from
literature. The experimentally measured resistivity values are also plotted in Figure
4.13, which highlights the change in the conduction mechanism in the high doping
regime. In this figure, at an Al doping concentration of approximately 1020cm?3 a
99
Table 4.1: Parameters characterizing Al-doped 4H-SiC
Parameter Value
?EA0 [meV] 214.86
?EA(NA) N
?3
E [cm ] 8.12? 1019
a 0.632
? 2Hall0 [cm /Vs] 109.6
?Hall(NA) N [cm
?3
H ] 2.92? 1018
b 0.6335
?Cond0 [cm
2/Vs] 109.83
N [cm?3C ] 2.92? 1018
c 0.5891
? [cm2Imp0 /Vs] 1.355
?Cond(NA)
NImp [cm
?3] 2.65? 1020
?Imp [cm
?3] 9.06? 1019
NV (300K) [cm
?3] 2.49? 1019
gA 4
N [cm?3] 4.5? 1020b
?(NA)
d 2.9
different conduction mechanism (likely a form of impurity or hopping conduction)
seems to dominate as the slope of the trend changes. This concentration at which
this slope change occurs agrees well with an estimated value of the critical concen-
tration for impurity conduction [127] Ncrit = (2.2a0)
?3 exp(1 ? r) ? 1 ? 1020cm?3
for 4H-SiC, where r = 9.76 [71]. The black dashed line in Figure 4.13 is generated
using the resistivity equation provided by Heera et al. [10] with our own data as
input.
From the data presented in Figure 4.11, the spread in experimental mobility
100
Figure 4.11: Hall mobility fit to experimental measurements of 4H-SiC doped with
Al
101
Figure 4.12: Conductivity mobility derived from experimental measurements of re-
sistivity
102
Figure 4.13: Resistivity fit to experimental measurements of 4H-SiC doped with Al.
Dashed black line comes from using the predicted resistivity from Heera [10] using
our values for mobility and ionization energy.
103
values is immediately obvious - likely due to the large differences in fabrication
methods and parameters. This spread is also readily apparent in the resistivity
data plotted in Figure 4.13. Interestingly, it turns out that all of the scattered data
that lies above the dashed trend line in this plot for N 18A > 10 cm
?3 comes from
devices doped using ion implantation as opposed to in situ doping via epitaxial
growth. This suggests that it is difficult to repair the lattice damage caused by
implantation and/or fully activate the implanted dopants.
To ensure a robust fit to the experimental data using Equations 4.59-4.62,
outlier points are excluded on a paper-by-paper basis. The remaining data points,
which are included in the least-squares error calculation, are outlined by black circles
in Figures 4.11-4.13. The standard deviation ??err of a 10 data point moving
window is also plotted, shown as blue dashed lines.
Conductivity mobility (resistivity) data for extremely low doping concentra-
tions is mostly unavailable due to difficulties obtaining completely pure samples,
but is required because of how it will affect the quality of the fitting function in
this region. For this reason, we add resistivity data points to the low doping region
that are calculated using the simple charge neutrality condition valid for low dop-
ing [10] described in the Appendix 3 (Section B.1.1). These resistivity data points
are labeled as Analytical in Figures 4.12 and 4.13.
The experimental incomplete ionization fraction is calculated by dividing the
conductivity mobility fit curve (Equation 4.60) by the Hall mobility fit curve (Equa-
tion 4.59) as described in Equation 4.58. The resulting p/NA fraction is presented
as the solid blue line in Figure 4.14. The standard deviation error bars ??err calcu-
104
Figure 4.14: Incomplete ionization ratio p/NA for 4H-SiC doped with Al at T =
300K. Experimental data fit (?Cond(NA)/?Hall(NA)) is shown in solid blue. Our
theoretical model, calculated with Equation 4.51 is shown in dashed red.
105
lated using a sliding window of 10 data points are also plotted as blue dashed lines
in this figure. The upper and lower bounds are calculated by dividing the ??err of
the Conductivity curve by the ??err of the Hall curve. The discrete points shown
are the ?extracted? discrete p/NA values calculated using: {? ? Cond ?Cond} /?Hall({NA } )
(red circles) and ? ({NHallCond A }?)/{? }?Hall (blue circles) where {}? indicates the set
of data points from the discrete measurements, and the functions indicate data cal-
culated from the empirically fitted curves (Equations 4.59 and 4.60). In addition
to the standard deviation lines, these data points help give some indication of the
spread in the data about the direct division of the two fitted functions.
Result Comparison: Now that p/NA has been extracted using experimen-
tal measurements, this result is used to confirm the values calculated using the
theoretical model presented in Sections 4.3-4.6. The theoretical model results are
plotted as the red dashed line in Figure 4.14 and for comparison, the experimental
model predicts the solid blue line. Upon comparing the theoretical results with the
empirically fitted, good agreement is found to within several percent.
4.8 Temperature Dependence
Because the theoretical model from Sections 4.3-4.6 is fully temperature de-
pendent, it can be used to predict the p/NA ratio at elevated temperatures, which
will often be the more relevant operating regime for SiC devices. Figure 4.15 shows
our predicted p/NA for elevated temperatures, although these values are not yet
confirmed with experimental data.
Due to the non-analytical form of p/NA derived using the theoretical model
106
described in Sections 4.3-4.6, and the significant number of parameters needed to
express the experimental data, a more usable form of p/NA for device designers is
desirable. The results of the temperature-dependent theoretical model are closely
approximated by the analytical expression in Equation 4.63. Optimal parameters
used in this expression are obtained by utilizing a genetic algorithm. The form
of this expression is derived from an expression presented by Kuzmicz [108]. This
parameterization was obtained for the following doping and temperature ranges,
respectively: 1014cm?3 ? NA ? 3?1020cm?3 and 300K ? T ? 800K. The benefit of
using this expression enables scientists and engineers can simply evaluate Equation
4.63 for their given temperature and doping conditions directly instead of having to
iteratively solve the nonlinear system of integrals described in Sections 4.3-4.6 for
each doping and temperature point.
p [ ]? 1? A exp ?(B |ln(NA/N C0)|) (4.63)
NA
Starting with the expression used by Kuzmicz [108] for Si, a temperature
dependence to the parameter C is added to help better match the theoretical model
results of p/NA for p-type 4H-SiC. The values of A, B, C, and N0 were optimized
for Al-doped 4H-SiC using a genetic algorithm with a squared error fitness function.
This parameterization is plotted in Figure 4.15 with a dashed line and agrees well
with the theoretical model results with a significant reduction in computational
effort. The full temperature dependent parameters A, B, C, and N0 were found to
107
be: ( )?0.2343
T
A(T ) = 0?.9611 ( ) (4.64)300???? 0.7748T0.0978 ( ) NA < N0300B(T,NA) = ???? (4.65)T0.3702? 0.0(347 ) otherwise3000.968
N (T ) = 2.7248? 1019 T ?30 ( ) cm (4.66)300
T
C(T ) = 5.4738? 1.2263 (4.67)
300
Plotting the theoretical model data in a different manner, Figure 4.16 shows
the total mobile hole concentration p calculated by multiplying the theoretical model
results in Figure 4.15 by NA. From these high temperature calculations, the tem-
perature dependence of hole concentration seems to diminish as the impurity band
conduction mechanism begins to dominate at around 2? 1020cm?3 for all tempera-
tures considered. Further experimental measurements at elevated temperature are
needed to confirm this prediction.
4.9 Results and Conclusion
From the gathered resistivity data, a consistent relationship between acceptor
doping concentration NA and resistivity is observed for epitaxially grown layers
doped with Al. In contrast, samples which are doped by implantation consistently
exhibit on the order of 10? higher resistivity compared to epitaxially grown layers
of the same doping concentration but even among this data there is a large spread
in values. This result is likely attributable to residual damage sustained during the
implantation process as well as incomplete high temperature activation of dopant
108
800
700
600
500
400
300
Figure 4.15: Incomplete ionization ratio p/NA calculated using our Model at el-
evated temperatures T = 300K ? 800K (solid line) compared to our empirical
parameterization (dashed).
109
800
700
600
500
400
300
Figure 4.16: Hole concentration in 4H-SiC doped with Al calculated using our Model
at elevated temperatures T = 300K ? 800K.
110
atoms into lattice sites. Both of these issues are apparently unable to be fully
rectified by the annealing process which generally follows implantation.
At an Al doping concentration greater than 1020cm?3 the start of a parallel
form of conduction is observed which is likely a type of impurity band conduction
known as variable range hopping [128, 129]. The parallel conduction mechanism
causes the resistivity to rapidly decrease with increasing doping concentration. This
increase in conduction is quickly counteracted by the solid solubility limit of Al
in 4H-SiC [130] which is around 2 ? 1020cm?3. Beyond this limit the resistivity is
observed to increase with further doping instead of decrease. The specific mechanism
is not investigated in this work but may potentially be attributed to multiple effects
including the formation of precipitated clusters of Al [64, 130, 131] which would
increase scattering due to the non-crystalline structure as well there is a potential
that holes conducting via impurity states may exhibit a higher effective mass.
Both the theoretical model and experimental data indicate that Al doping
concentrations between 1017cm?3 and 1020cm?3 typically result in mobile hole con-
centrations equal to 10% or less of the doping concentration. This poses a challenge
for semiconductor device designers trying to achieve low resistivity p-type regions
in devices. The minimum p/NA ratio predicted by the theoretical model is below
30% even at temperatures as high as 800K.
111
Chapter 5: Transport Model Validation and Genetic Algorithm Ap-
plication
5.1 The Monte Carlo Method
In principle, to simulation carrier motion within a semiconductor we would
need to solve the Boltzmann transport equation (BTE). The BTE, shown in Equa-
tion 5.1 is a nonlinear differential equation which provides us with a probability
distribution function fk in 3 dimensions of space, 3 of momentum, and one of
time [132?134].
( )
?fk
( ) = ?
1?kE? ? ?
eF
rfk + ? ?
?fk
~ ~ k
fk + (5.1)
?t ? ?t scatterM?fk V
= [fk?(1? fk)Si(k?, k)? fk(1? f ? 3 ?3 k?)Si(k, k )] d k?t scatter (2?) i
In this equation, the distribution function is evolved in time due to the effects
of drifting and scattering in an applied field. Due to the complexity associated
with solving for this function in a 6-dimensional phase-space directly, we instead
employ the use of the Monte Carlo Method to obtain a solution instead. The aim
of the this technique is to model electron transport by generating numerical results
statistically through a repeated random sampling process. This method allows us
to take calculated scattering rates and follow an electron on it?s trajectory as it
travels through a crystal, scattering via random mechanisms with appropriately
112
weighted probabilities. During the flight of the electron through the crystal, we also
collect statistics about its energy, velocity, and valley occupation. Generally, it is
assumed that each scattering particle acts independently from all other carriers, so
the statistics generated from following one carrier undergoing many events should
be representative of the ensemble of many carriers all mutually scattering.
Figure 5.1 lays out the flowchart the Monte Carlo algorithm follows. To sim-
plify, first a total scattering rate is calculated for an electron given its initial energy.
Then, a flight time is randomly determined from an exponential distribution. This
flight time represents the amount of time the electron is drifting under the effect
of the electric field and accelerating uninhibited. During the flight the electron?s
k vector changes according to 3.17. At the end of the flight the electron scatters
according to a randomly chosen mechanism with probability proportional to each
mechanism?s individual scattering rate. After deciding the mechanism, the elec-
tron?s energy changes by an amount determined by the phonon energy associated
with the mechanism. Finally, the electron?s k vector is updated to the new k? di-
rection after scattering interaction with a phonon of wavevector q. Since there are
generally many vectors which all have the same energy, the new vector is again
chosen randomly from a distribution associated with the decided mechanism, and
is based off the analytical form of the scattering rate. Certain mechanisms tend to
prefer directional, anisotropic scattering and others are more isotropic. Finally, the
entire process is then started anew and repeated many times to accrue sufficient
data to produce meaningful averages and smooth distributions.
113
Start
Define Params.
- Electric Field (F)
- Material Props.
Sample Data
Initialize - E, Vx
i=0, E=3k
2 B
T
Free Flight
- tf = tf ? dti=i+1
- kx = kx ? qFdt/~
- Calc. New ~k and E
Calculate Scattering Rates T
Random Time of Flight
t > dt
tf = ??T ln(1 ? r f1)
Choose Scatter
F
Mech. using r2
Determine Final State
Based on Mech., and
?k, ?k from r3, and r4
Update State
- k = kfinal
- Calc. New Energy
Sample Data F
i>MaxFlights
- E, Vx
T
Output Histograms
- Velocity Distribution
- Energy Distribution
- k Distribution
Stop
Figure 5.1: Flowchart of Monte Carlo Algorithm for a single field configuration using
4 random numbers r1, r2, r3, and r4
114
Calculations are performed using the algorithm for different strength electric
fields. Averages are then taken for each calculation and plotted against the strength
of the field. Of particular interest are the plots of average energy and average drift
velocity (in the direction of the field). The derivative of the velocity with respect to
field strength produces a plot of the differential field-effect mobility as a function of
field. From the velocity plot, we can see material properties such as possible velocity
overshoot and eventual saturation as the field gets too large.
Figure 5.2 shows the scattering rates of all mechanisms which we consider
in the Monte Carlo simulation. The work done currently only includes the 3D
mechanisms and is in good agreement with work previously done in the literature
[40]. To maintain realism within our approximation, three subbands are chosen to
account for a well configuration at the onset of inversion where there is 1.4 eV of
band-bending in the SiC ( 1016cm?3). The lowest three subbands for a well with
surface field of 2MV/cm are 0.37, 0.66, and 0.89 eV above the conduction band
minimum, allowing the electrons to have enough room to gain kinetic energy while
remaining inside the well for relatively small applied lateral electric field. More
subbands can be included depending on the actual well configuration being studied.
The 3D Monte Carlo simulation was performed for 4H-SiC with transport on
the (0001) interface. Figure 5.3 shows the results of the time-averaged electron
velocity determined as a function of applied lateral field. The velocity is taken to
be the component only in the direction of the field so that the mobility extracted
115
15
10
Pop em
Pop ab
Opem
Opab14
10 Ac em
Ac ab
2D 1st
2D 2nd13
10 2D 3rd
12
10
11
10
10
10
0 1 2 3
Energy [eV]
Figure 5.2: Scattering rates for electrons with different total energy. Rates asso-
ciated with phonon emission are shown as dashed lines and solid lines are used
for absorption. The 2D Atomic Roughness scattering rates out of the first three
subbands are also shown for a surface field of 2MV/cm.
directly relates to the field as Equation 5.2.
?drift = ??F (5.2)
For low electric fields, the drift velocity is known to increase linearly with
applied field. From the results obtained, we see a velocity overshoot peak of ap-
proximately 1.7e7 cm/s at around 700kV/cm. The resulting value of this velocity
peak agrees with its theoretically expected value which can be derived by equat-
ing the kinetic energy of the electron to the energy of the most probable phonon,
116
Scattering Rate [1/s]
8
10
7
10
6
10
5
10
4
10
0 1 2 3 4
10 10 10 10 10
F [kV/cm]
Figure 5.3: Average electron velocity (in the direction of the applied lateral field)
plotted against field strength. Monte Carlo simulation results are solid lines and
data-fitted curve from [11] is shown with the dashed line. Curves shown are for
different temperatures: Blue 300K, Red 600K
and solving for the velocity the electron achieved before phonon interaction. To get
the average velocity between scattering events simply divide the emission velocity
by two. For SiC, polar-optical phonon emission accounts for the highest scattering
rate and thus is the most probable interaction type. Using the values from the MC
simulation results in a velocity saturation prediction of 1.7e7 cm/s.
m???
2
emit ? ~?pop (5.3)
2
1 2~?pop
vsat = ? = 1.7e7cm/s2 m
Experimental data exists only up to 400kV/cm, where the velocity is fitted to
saturate at 2.2e7 cm/s [11] (red curve in Figure 5.3), which agrees fairly well with
our simulation results. Beyond this point, the results of the Monte Carlo show that
the velocity overshoots and begins to decrease as it saturates due to the increasing
117
V [cm/s]
scattering rate of inelastic acoustic phonons. In an inelastic scattering event, kinetic
energy of the electrons goes into the energy of the phonon, reducing the electron?s
speed after the event.
Because this calculation only includes M valley, all electrons exist in the lowest
conduction band and thus have potential energy defined to be zero. The energy then
shown in the Figure 5.4 represents both average total and average kinetic energy of
the electrons in the material. From the plot we see that at the breakdown field of
SiC (3MV/cm), the electrons have 1.7eV of kinetic energy which is reasonably small
compared to the SiC/SiO2 conduction band offset barrier of 2.7eV [135]. Addition-
ally, from the distributions presented in Figure 5.5, at 3MV/cm the vast majority
of electrons are still below 2.7eV but at this point the validity of the simulation
begins to break down. At 700kV/cm where the average velocity is shown to peak,
the energy distribution shows that almost all electrons have less than 1eV of energy.
From the three energy distributions shown, we see that there is quite a drastic field
dependence on the shape and spread of electron energies. At low fields the distribu-
tion is similar to an exponential-decay and as field increases, the distribution peaks
at higher and higher energies. The tail also follows the same trend and extends
to higher energy as the distribution broadens. At the breakdown field, we see the
energy distribution has transformed to a nearly Gaussian shape because the electron
is scattering so frequently and constantly changing energy.
Absolute field effect mobility is calculated by dividing the electron velocity at
a given field by the value of the field. To extract differential mobility the derivative
118
1
10
0
10
-1
10
-2
10
0 1 2 3 4
10 10 10 10 10
F [kV/cm]
Figure 5.4: Plot of average electron energy for given lateral field strength. Zero
energy corresponds to electrons at the conduction band minimum. Only M-Valley
scattering is considered in this simulation.
of the average velocity vs field plot is taken. The results of both calculations are
shown in Figures 5.6 and 5.7 respectively.
5.2 Conclusions and Future Work
Results of this Monte Carlo simulation closely match the experimental velocity
data available for SiC from [11], and the MC simulation results of [40]. To extend this
work, I plan to add the effects of 2D interface scattering within the quantum well to
account for atomic and surface roughness scattering. Additionally, these simulations
were performed assuming transport along the (0001) face (Si-face). Because my
method for extracting an atomic roughness perturbation potential is general, the
next step planned is to create and relax an a-face aligned interface supercell using
119
Energy [eV]
4 4
?10 F=100 kV/cm ?10 F=700 kV/cm 4?10 F=3000 kV/cm
6 8 5
4
6
4
3
4
2
2
2
1
0 0 0
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
-1 10 -1 10 10k  [m ] ?10 k  [m ] ?10 k  [m -1 ] ?10
x x x
4 4 4
?10 F=100 kV/cm ?10 F=700 kV/cm ?10 F=3000 kV/cm
8 8 5
4
6 6
3
4 4
2
2 2
1
0 0 0
0 1 2 3 4 0 1 2 3 4 0 1 2 3 4
Energy [eV] Energy [eV] Energy [eV]
Figure 5.5: Distributions of wavevector component in the direction of applied lateral
field (top) and electron energy (bottom) for various lateral field magnitudes.
600
500
400
300
200
100
0
0 1 2 3 4
10 10 10 10 10
F [kV/cm]
Figure 5.6: absolute field mobility for electrons plotted against applied lateral field,
taken as ?(F )/F .
120
Mobility [cm 2 /Vs]
600
500
400
300
200
100
0
-100
0 1 2 3 4
10 10 10 10 10
F [kV/cm]
Figure 5.7: Differential electron field mobility as a function of field, extracted from
Monte Carlo simulation results.
DFT. This work will allow for a comparison to be made between various interfaces
to predict the mobility achievable in an ideal-case scenario where there is no surface
roughness scattering due to the presence of a miscut.
5.3 Genetic Algorithm For Parameter Extraction
5.3.1 Method Overview
To develop a more usable form of the p/NA results from Section 4.8 without
the need to solve the complete iterative system of integrals, a genetic algorithm
was designed an applied to determine optimal parameters for an empirical formula.
In general, a genetic algorithm is a stochastic technique which uses concepts from
biological evolution in order to optimize a function. For my purposes, the genetic
algorithm is used to optimize a set of parameters for an empirical function to best
121
Differential Mobility [cm 2 /Vs]
fit experimental data.
Each parameter set is denoted as an ?individual? and each parameter?s value
as a ?gene? in the context of the genetic algorithm. An individual?s ?chromosome?
is the combination of all of its genes. In the algorithm, many individuals make
up a ?generation? and each are tested by calculating a fitness function for a given
individual Fi (Eq. 5.4). The fitness function creates an objective ranking of the
individuals by determining the ?goodness-of-fit? for each parameter set. In this case
the fitness function for each individual is taken to be the least-squares error between
the empirical function using that individual?s parameter set fi,empirical(x) and the
experimental data set fdata(xj). A weighting function w(x) may also be added to
ensure that fitting preferentially favors certain more sensitive or important regions
or minimizes the impact of regions where experimental data is unreliable. The sum
in Equation 5.4 is taken over the j experimental data points.
??[ ]
Fi = w(xj) (f 2i,empirical(xj)? fdata(xj)) (5.4)
j
In the genetic algorithm, to limit memory usage, the population size is fixed
such that each generation has the same number of individuals. Each generation
of individuals is primarily made from ?breeding? the individuals from the previous
generation to create new individuals. Mimicking the biological process of sexual
reproduction, different genes are randomly taken from each parent and combined
to create offspring. After each generation is created, old individuals die out and
are replaced by their children and the cycle is repeated until a satisfactory fitness
122
level is achieved. The basic flow diagram of the genetic algorithm is shown below
in Figure 5.8.
Start
Initial Population
Calculate Fitness
min(F ) < err Createi
i F New Generation
T
Stop
Figure 5.8: Flowchart of the basic Genetic Algorithm process.
123
Create New Generation
Breeding Elite
Select Promote
Crossover Stock
Mutate Random
Figure 5.9: Operations performed to create a new generation of individuals in a
genetic algorithm.
In the algorithm flowchart Figure 5.8, the initial population is created by
assigning each individual?s genes (parameter values) by selecting randomly from
within the bounds of reasonable values set for each parameter. The fitness of each
individual in the generation is then evaluated to assign a qualitative ranking of
how well each individual fits the experimental data. The best fitness among the
generation is then used to determine the if the algorithm is to stop by comparing
it to a desired threshold value. If sufficient fitness is not achieved among any of
the individuals of this generation, a new generation is created and the process is
repeated. Figure 5.9 diagrams the process of creating a new generation showing its
component groupings.
Every generation is based off of the previous generation except for the initial
124
population. The main bulk of the generation is created with a breeding process which
selects parent individuals from the previous generation and ?breeds? them to create
individuals in the new generation. Parents are selected randomly from the previous
generation but preference is given to individuals with better fitness. This effect can
be easily achieved by implementing a hard cutoff where the best X individuals act as
the ?elite parent pool? from which the parents of a set number of children in the next
generation are chosen. The remaining children are created by choosing completely
random parents without any bias based on fitness. This ensures that the best
X parents will have more representation (through their offspring) in the following
generation. For simplicity of implementation, every pair of parents has one child,
and individuals are not monogamous. Additionally, there is a chance that duplicate
pairs may appear in which case that pair of parents ends up with multiple children.
After selection of the parents, their child is generated with a two-step process:
crossover and mutation, the details of which will be described later. Additionally,
during the selection process, the best X individuals can also be promoted directly
from the current generation into the new generation in a process known as ?Elitist
Selection?. By adding this feature to the generation creation process, the current
best solution is never forgotten, ensuring that the fitness function monotonically
improves with each generation.
Earlier, it was mentioned that individuals with better fitness were selected
more frequently to be parents. Doing this ensures that better individuals produce
more children and, in effect, biases the algorithm to stay near solutions which are
known to be good via the breeding process. However, this can eventually lead to
125
convergence issues. It is possible and indeed likely that many of the best individuals
will start to become near-clones of each other i.e. the best individuals will have
very similar genes. In effect, this causes the gene pool to stagnate and keeps the
algorithm stuck in a local minimum of the fitness function rather than moving
toward the global minimum. To help combat a stagnated gene pool, new ?stock?
individuals are introduced to each generation as a supplement to all of the offspring
made by breeding the previous generation?s individuals. These ?stock? individuals
have no parents and are created instead with purely random genes (like the initial
population individuals) to promote genetic diversity further down the generations.
The likelihood that these individuals themselves give highly fit solutions is low but
they may introduce genes that aren?t currently available elsewhere in the population
that can lead to highly fit individuals after breeding.
5.3.2 Breeding Process Details
The breeding process selects two individuals as parents and mixes their genes
to create an offspring. The selection of parent individuals is performed randomly
but bias is given towards selecting individuals with better fitness. To generate a
child from a given a pair of parents, the crossover operator is applied. Crossover
is performed by generating a random bit array, with length equal to the number of
genes in an individual. These bits then determine from which parent each gene of
the child is to be taken i.e. for a ?1? take parent 1?s gene and for a ?0? take parent
2?s gene. Written mathematically: Gi,Child = riGi,P1 + (1? ri)Gi,P2 where ri is the
126
ith bit in the crossover array and the Gi terms are the i
th genes of the child and
parent individuals. Due to the large number of possible crossover arrays, even if
two individuals share the same parents, they will likely differ significantly from each
other.
Parent 1: 0.53 22.2 3.87 ? ? ? 3.5? 1018 1059 11.91
Parent 2: 1.96 45.1 1.46 ? ? ? 7.2? 1017 1360 25.34
Crossover:
1.96 22.2 1.46 ? ? ? 3.5? 1018 1059 25.34
Bit Array: 0 1 0 ? ? ? 1 1 0
Figure 5.10: Random binary bit array used to generate an offspring from two parent
individuals
After crossover, each new offspring is then subject to a mutation stage which adds
some random spreading or broadening to the genes (parameter values). This can
help to fine-tune the parameters in later generations once most of the major con-
vergence has been achieved mainly due to crossover in the early generations. A
simple way to add mutation is to multiply each of the child?s genes by a random
Gaussian variable centered at 1 using set standard deviations ?i for each parameter
N (? = 1, ?2i ). As a note, for parameters which may span many orders of magnitude
(such as doping concentrations, reference doping values, electron concentrations, ox-
ide charge, etc.), it may be better to apply mutation to the log of the parameter
127
and then re-exponentiate to recover the final mutated value.
Crossover:
1.96 22.2 1.46 ? ? ? 3.5? 1018 1059 25.34
Mutation:
Child: 2.05 21.6 1.15 ? ? ? 3.8? 1018 1057 25.82
Figure 5.11: Random mutation applied after the crossover to finally create the
offspring.
An alternative breeding strategy which could also be used is to form each child?s
genes via a random linear combination of it?s parents genes instead of choosing the
exact value of either parent. The crossover formula Gi,Child = riGi,P1 + (1? ri)Gi,P2
still applies, but now each ri is a random variable which can vary continuously from
0 to 1 and is chosen from a distribution. In this strategy, ri acts as more of a mixing
fraction than a selector. For simplicity the distribution may be uniformly random
to create more hybridized individuals, or it may have a more bimodal shape that is
biased towards 0 and 1 to give more of a crossover effect (Fig. 5.12). Once again,
for parameters that span many orders of magnitude this linear mixing should be
performed on the log of the gene and then re-exponentiated to obtain the mixed
value.
128
PDF
5
4
3
2
1
0
0 0.2 0.4 0.6 0.8 1
r
Figure 5.12: Example bimodal probability density function which could be used to
pick random gene mixing fractions.
Parent 1: 0.53 22.2 3.87 ? ? ? 3.5? 1018 1059 11.91
Parent 2: 1.96 45.1 1.46 ? ? ? 7.2? 1017 1360 25.34
Mixing:
Child: 1.79 23.1 2.30 ? ? ? 2.6? 1018 1137 18.89
Mixing Array: 0.12 0.96 0.35 ? ? ? 0.82 0.74 0.48
Figure 5.13: Child created using linear combinations of its parents? genes.
The values and rates of convergence achievable will depend on the problem
being solved as well as the chosen form of the fitting function and the exact im-
129
plementation of the algorithm (size and proportions of individual groups in each
generation, crossover and mutation functions, etc.). A good solution should be ob-
tainable for most reasonable implementations given that the parameterized function
has enough parameters to reasonably capture the shape of the data and the evolution
is given enough generations to converge.
5.3.3 Genetic Algorithm Applied to Doping and Temperature-Dependent
Hall Mobility
In researching incomplete ionization in Al-doped 4H-SiC, I gathered a large
number of references and extracted thousands of measurements of Hall mobility
taken under a wide range of doping and temperature conditions. As these con-
ditions change, the hole concentration and the dominant scattering mechanisms
vary. The combination of the way these two varying quantities interact give rise to
dramatic variations in Hall mobility as a function of temperature and doping. Con-
sequently, the shape of the Hall mobility function is complex and, given the large
amount of experimental data, a genetic algorithm is the ideal choice for finding the
optimal parameterization of an empirical Hall mobility fit. In this section I present
my method and results for the temperature and doping-dependent Hall mobility
function of Al-doped 4H-SiC.
Multiple runs of the algorithm have been performed starting from a new ran-
dom population each time in order to demonstrate the convergence of the algorithm.
The fitness of the current best individual is plotted against the generation number in
130
Figure 5.14 for each of these separate runs. Initially, the completely random initial
populations (generation 1) give rise to a range of random fitness values from about
10 to 13 and within 10 generations, all lie between 8 and 9. Rapid improvements are
made within the first 100 or less generations in all runs. Occasionally, large jumps
in fitness can be seen which is when a ?breakthrough? is made. This generally occurs
when a new individual is bred which overtakes the best elite candidate rather than
the best candidate (or near-clones of it) making small improvements through slight
mutation.
12
11
10
9
8
0 1 2 3 4
10 10 10 10 10
Generation
Figure 5.14: Fitness across generations for different runs of the algorithm.
It is interesting to note that for this particular fitting function and data set,
about half of the runs got temporarily stuck in the same local minimum, where
fitness stagnates at ? 8.15 for a varying number of generations. It turns out that
131
Fitness
this is a local minimum which contains nearly-optimized parameters for the majority
of the fitting function, i.e. the mobility matches the data for most doping and
temperature conditions, with the exception of the simultaneous high doping and
low temperature region. As will be described later, the shape of the function in this
specific region is mainly dominated via multiplication with a hyperbolic tangent
function containing 4 parameters. Because the hyperbolic tangent goes to unity
in every other region, these parameters are difficult for the algorithm to optimize.
Changing the parameters in the hyperbolic tangent does not noticeably change the
value of fitness as most of the data points exist outside this region where the fitting
function is mostly insensitive to changes in these parameters. Only after the majority
of the function has settled to near-optimized values does the fitness function have the
sensitivity to determine how good or bad various values of the hyperbolic tangent
parameters are. After this point, the algorithm can start to select better individuals
which fit the data in this specific region and continue forward to the final fitting
function which works across the entire space.
The Hall mobility function with the best fitness achieved is shown below in
Figure 5.15, plotted against the scattered experimental data points taken from lit-
erature which is used to evaluate the fitness function. Vertical projection lines show
the error between the fitting function and the experimental data. Not visible in the
figure are the approximately equal number of data points which are underneath the
fit surface which have an approximately equal amount of error to those above it.
132
10 3
3
10
10 2
2
10
10 1
1
10
0 0
10 10
15
10 1610 17 20010 18 30010 19 500
400
10 20 60010 21 70010 800
-3
N  [cm ] Temperature [K]
A
Figure 5.15: Genetic Algorithm fit of the doping and temperature-dependent Hall
mobility of Al-doped 4H-SiC.
The empirical function used is a modified form of the Caughey-Thomas mo-
bility model, with parameters given in Table 5.1:
( )??
? T0 T
?Hall(NA, T ) = B ? ( ref(NA, T ) ) ? (5.5)(? ?T )
( [ 1 + NA( ( )N0 exp )TT0 ])A
B 1 NA 1(NA, T ) = tanh ? T ?B + (5.6)
2 Nref 2
133
2
 [cm /Vs]
Hall
Table 5.1: Optimized parameters for empirical Hall mobility in Al-doped 4H-SiC
obtained using a genetic algorithm
Parameter Value
?0 [cm
2/Vs] 320.1
Tref [K] 190.15
? 2.926
N [cm?30 ] 3.8? 1017
T0 [K] 108.6
? 0.959
? [K?1] 9.5? 10?4
? [K?1] 0.0139
B [K] 318.86
Nref [cm
?3] 2.15? 1023
A 0.691
134
Chapter 6: Germanium Modeling
6.1 Introduction
Germanium is a group IV semiconductor commonly used in Short Wave In-
frared (SWIR) optical devices due to its relatively small band gap of 0.66eV. Like
silicon in the period above it, the conduction band minimum of germanium does
not lie at the same point in k space as the valence band maximum, making it an
indirect gap material and thus reducing its absorption efficiency. Unlike Si however,
the direct gap of Ge is only slightly larger than its indirect gap energy. With clever
bandgap engineering Ge is able to transition to a direct gap material. One such
method showing promise is alloying Ge with Sn in various ratios. Using DFT we
can calculate the effects the alloy has on the band structure for different percentages
of tin and thus predict the percentage needed to transition germanium into a direct
gap material.
Because the bandgap of Ge is indirect, it is an inefficient absorber of light
meaning that incoming photons with energy close to the bandgap energy penetrate
deeper into the material before being absorbed. If the absorbing material is made
thick compared to the minority carrier diffusion length, generated carriers will be
more likely to recombine before crossing the junction and thus will not be collected
135
- decreasing quantum efficiency. Additionally, thicker materials increase the chance
of defects which can act as further sights of recombination. By using direct-gap
semiconductors as the absorbing material, the active layer can be made thin and
quantum efficiency can be kept high.
Absorption of light into a semiconductor is described by the Lambert-Beer Law
(Equation 6.1), which uses a parameter known as the absorption coefficient which
describes how the intensity of light with a given wavelength changes with depth into
the material due to absorption. In Equation 6.1, I0 is the monochromatic intensity
at the surface, ? is the absorption coefficient (a function of wavelength), and x is
the depth into the semiconductor. The absorption coefficient ? is also related to a
material property known as the extinction coefficient ? by Equation 6.2.
I(x) = I e??x0 (6.1)
4?f?
? = (6.2)
c
The absorption coefficient amalgamates details of the band structure and turns
conduction band minima into kinks in the coefficient plot. As photon energies
increase and, more states are available for valence band electrons to be excited to
once the energy passes higher conduction band minima, making an excitation and
thus an absorption more probable. Indirect minima show a quadratic dependence
of the absorption on the energy and direct minima show an approximately square
root dependence [136]. Figure 6.1 shows the absorption coefficient as a function
of photon energy for select semiconductors which absorb in the visible and near-
infrared. Here, we see that absorption starts in Ge and Si at 0.66eV and 1.15eV for
136
room temperature measurements, corresponding to their bandgap energies. For Ge
absorption towards 0.8eV we see a transition into a rapidly increasing absorption due
to the direct valley. The final bend corresponds to the second indirect minimum
at 0.85eV. Further confirmation that the 0.8eV minimum is the direct gap value
can be seen by looking at the low temperature curve. At low temperatures, it is
less likely for an electron to be excited into an indirect minimum because there
are less phonons available to assist the momentum transfer needed. Because of
this, absorption starts closer toward the energy of the direct gap as the indirect
transitions become increasingly unlikely.
Figure 6.1: Absorption coefficient for various semiconductor materials. (Plot from
[12])
137
As with most semiconductors and their alloys, the valence band maximum of
germanium is located at the gamma point in the Brillouin Zone. The conduction
band minimum for germanium is located at the L point for the conventional FCC
lattice structure and is 140meV less than the direct gap energy [137]. This relatively
small difference has been the motivation for band structure research and engineering
to achieve a direct-gap form of the material for use in optical devices. The most
promising techniques applied to achieve such a structure include applying strain
and alloying with various other elements [138]. The obvious choice of alloy material
has been tin due to its location in the period immediately below germanium. The
alpha allotrope variant of tin has the same diamond crystal structure as germanium
but acts like a semimetal with a negative band gap at the gamma point [138]. An
elementary application of the simple linear form of Vegards law gives the indication
that the transition from direct to indirect gap should occur at approximately a 21%
uniform tin concentration [139]. Experimental and calculated results both show the
presence of a bowing parameter which is needed to fit the non-linear experimental
data for how the band gaps change with varying tin concentration [140]. In general,
it is known that with increasing Sn concentration, both the direct and indirect
gaps of Ge1?xSnx shrink but the direct gap does so at a faster rate. If the exact
concentration of Sn needed to cause this transition can be determined, direct gap
devices can be fabricated while still maintaining as much of the gap as possible.
138
6.2 DFT-Based Analysis of Ge
Band structure calculations in the past have predominately been performed
using single-electron empirical pseudopotential methods with extrapolation to fit
data to GeSn alloys [141]. We have performed calculations for specific compositional
percentages of GeSn using Density-Functional Theory (DFT). Figure 6.2 shows the
cells on which we have used to perform calculations, containing 12.5%, 6.25%, and
3.125% Sn. These cells were constructed using repeated 8 atom cubic unit cells.
Figure 6.2: Atomic supercells of Ge1?xSnx used in DFT calculations for 12.5%,
6.25%, 3.125% Sn
Before we calculate the predicted band structure of the alloys however, we
first needed to reproduce the experimentally known results of a pure Ge crystal.
To do this we must chose not only an appropriate functional and pseudopotential,
but also ensure the calculation results are converged with respect to the numerical
discretization built in to the code. The accuracy of DFT results are highly dependent
on the combination of the system being studied, the pseudopotential used, and the
functional applied to the calculation [142]. Pseudopotentials act as an effective ionic
139
potential that each valence electron interacts with, and acts to stabilize calculations
by treating sharp core potential as a smoother approximation within some defined
radius. In doing this the high frequency components of core wavefunctions are
supressed while leaving the electrically and chemically-important valence portion
unchanged. By making this approximation, larger systems of atoms become solvable
by reducing the number of plane waves needed in the wavefunction expansion within
calculation. Different methods exist for generating pseudopotentials and accuracy
is generally dependent on the configuration of the system being studied. After
testing numerous pseudopotentials generated using different functionals and valence
occupancies, we were able to use a pseudopotential generated using the Perdew-
Burke-Ernzerhof (PBE) functional to obtain an accurate band structure of pure
germanium crystal. The self-consistent calculation was performed using the hybrid
PBE0 functional which allowed us to correct the indirect and direct bandgap energies
to their experimental values by applying the appropriate mixing fraction of exact
Hartree-Fock exchange energy. Calculations in DFT transform the set of Kohn-
Sham equations for non-interacting particles by expanding the potentials and wave
functions onto a finite basis of plane waves, the highest energy of which is set by a
cutoff energy threshold. By adjusting this cutoff energy the number of plane waves
in the expansion can be adjusted. There basis set should be sufficiently large to
ensure good representation of the wave functions but this increase comes at the
cost of increased computation time and memory usage. In addition to requiring a
sufficient cutoff energy, calculations also should be performed at enough k points
to adequately sample the BZ. During the calculation, integrals over the entire BZ
140
are needed to calculate the charge density at each iteration. This integration is
approximated by taking a weighted sum over a finite set of special points in the
irreducible wedge of the BZ with weights corresponding to the number of equivalent
points in the entire BZ [143]. The BZ of Ge is the same truncated octahedron as
that of Si, shown in Figure 6.3.
L
U
?
X
K W
Figure 6.3: Brillouin Zone of an FCC lattice.
In addition to increasing the cutoff energy, increasing the number of k points
increases the computation time. Ideally, we would use the Monkhorst-Pack k point
grid as a way to select points in an unbiased manner and perform a relatively
computation-heavy self-consistent field (SCF) calculation to obtain an accurate rep-
resentation of the system potential. This potential can then be used as an input to
an inexpensive non-self-consistent field (NSCF) calculation with k points selected
along the path of the high symmetry points of the desired band structure plot.
Unfortunately, Quantum Espresso does not support the ability to perform NSCF
calculations using hybrid functionals, so our calculations had to be performed with
141
the desired k point path directly in the SCF calculations. To minimize integration
errors, we used a long path length as well as a high k point density to attempt to
cover the majority of the irreducible wedge.
While trying to tune the parameters used to obtain the experimental band
gaps for the pure Ge crystal cell, we noticed that both the direct and indirect gap
energies varied linearly with the mixing fraction as shown in Figure 6.4.
Figure 6.4: Indirect gap energy at the L point varying linearly with hybrid functional
(PBE0) mixing fraction
The lattice constants of materials are not reliably reproducible with DFT in
general, and the degree of error changes with both material and functional used.
More specifically, the PBE functional used in this work is known to regularly over-
estimate the lattice constant in solids by up to 2% [144]. According to the results
by [142], the lattice constant of Ge is overestimated by about 1.9% for the PBE
functional. We started with calculating the band structure from the experimen-
tally reported lattice constant of 5.646A? and then performed calculations for other
142
slightly deviated lattice parameters. By adjusting the lattice constant of the cell we
effectively apply a hydraulic strain to the crystal. Figure 6.5 shows the resulting plot
of the energy gap difference with respect to lattice constant which exhibits a linear
relationship. We adjust the lattice constant to give the experimental gap difference
of 140meV resulting in a lattice constant 1.3% strained from the experimental value
of 5.646A? [138].
Figure 6.5: Gap Difference (E?-EL) varying linearly with lattice constant. Litera-
ture: 5.646A? Our Work: 5.573A? ( 1.3% difference)
The final band structure obtained of the pure germanium crystal is shown in
Figure 6.6, with an indirect gap of 0.66eV and a direct gap of 0.80eV achieved using
a lattice constant of 5.573A?. Our calculation was able to produce a band structure
which matched the gap energies reported in literature [8] for the lowest direct and
indirect valleys, with the rest of the valleys in moderately good agreement. The band
structure diagram from literature is given in Figure 6.7 and a comparison of the gap
143
Figure 6.6: DFT calculated Ge band structure for 2 atom cell using the PBE0 hybrid
functional. Dashed box shows the same section of the band structure given in the
literature [8]. The key direct and indirect gaps are shown with red arrows.
energies in Table 6.1. This calculation was performed on a simple 2 atom cell with
an FCC Bravais lattice and the results give confidence in using this pseudopotential
Table 6.1: Comparison of calculated energy gap values of Ge to those found in the
literature.
Gap [eV] Calc. Lit. [8]
Eg (EL) 0.659 0.66
E?1 0.799 0.8
?E 0.788 0.86
EX 0.977 1.2
E?2 2.926 3.22
Eso ?0 0.29
144
Figure 6.7: Ge band structure from [8]. The key direct and indirect gaps are shown
with red arrows.
and functional for further, larger calculations on GeSn alloys.
6.3 DFT-Based Analysis of GeSn Alloys
To generate supercells with variable percentage of tin, we change the cell from
an FCC Bravais lattice with a 2 atom basis to a simple cubic lattice with an 8 atom
basis or a tetragonal lattice with a 16 or 32 atom basis, each with a single tin atom.
As with the 2 atom cell, we tune the lattice constant and mixing fraction of an 8
atom cell which can be repeated to generate the other cells. The mixing fraction
used to tune the gap values of the pure germanium 8 atom cell will be kept constant
when tin is added and when the larger cells with tin are used.
Performing computations using larger supercells comes with the cost of signifi-
145
cantly increased computation time. The Virtual Crystal Approximation (VCA) is a
work-around to this problem, wherein the characteristics from the pseudopotentials
of both Ge and Sn are merged into a hybrid pseudopotential representing a non-
existent atom which approximates a percent concentration of each atom, allowing
for small atomic cells but introducing non-physical atoms [139, 141]. Because of
the potentially limiting transferability inherent to certain kinds of pseudopotentials,
care must be taken to ensure that pseudopotentials maintain accuracy. The effects of
these approximations are evident in published results on the subject, where the var-
ious transition percentage predictions range greatly from 6-21% or more depending
on the method used [138]. With access to the High-Performance Computing Clus-
ter Deepthought2 at the University of Maryland, the computationally expensive,
large supercell DFT calculations can be run massively parallelized with Quantum
ESPRESSO [30] to obtain accurate band structures for different cell sizes and frac-
tions of Sn. We believe that performing these ab initio calculations provides us more
accurate information about the band structure of the system than is obtained using
the VCA.
Calculation time is system dependent and increases greatly with an increase
in the number of k points used, the number of bands in the calculation, and the
number of plane waves. As an approximation, the time to complete a calculation
increases proportional to the number of atoms per unit cell and is given by Equation
6.3 [143].
Tcpu ? NiterNk ? (O(N 2bNpw) +O(NbNpw logN 2pw) +O(NbNpw)) (6.3)
146
Where Niter is the number of iterations required to achieve self-consistency,
Nk is the number of k points specified, Nb is the number of bands, and Npw is
the number of plane waves in the expansion. The number of k points required to
achieve convergence generally decreases with an increase in the supercell size because
of the inverse size relationship of real-space and reciprocal-space. In contrast, the
number of bands required will increase as it is calculated by taking the number
of atoms in the supercell and multiplying by the number of valence electrons per
atom. The number of plane waves required to achieve convergence will also generally
increase with an increasing supercell size. The number of iterations to achieve
self-consistency is more difficult to predict but can be assumed to fall within 5 to
20. To complete calculations in a reasonable amount of time due to this highly-
nonlinear time scaling, we use the multiple parallelization levels available in the
Quantum Espresso PWscf program. Quantum Espresso is set up with different
levels of parallelism forming a hierarchical structure. The top level divides the
processors into pools, each of which takes care of the calculation at a group of k
points. The next level, known as plane wave parallelization, distributes the wave
function coefficients across the processors in each pool of plane wave processors,
offering one of the biggest calculation speedups. Once the speedup for this level
saturates, the final level of parallelization can extend the processor scaling. This level
divides each group of plane wave processors in to task groups, each which perform
the calculation on a group of electronic states. The plane wave groups can also be
partitioned into linear algebra groups which parallelize diagonalization and matrix
multiplication by distributing across groups of a square number of processors [143].
147
Through various trials we have been able to find optimal distributions of pro-
cessors for each parallelization level and bring calculation times down significantly.
For the two atom cell using hybrid functional DFT with 1500 k points and a cutoff
energy of 100Ry, the calculation time was able to be reduced from multiple hours
in a serial hybrid functional calculation to under 5 minutes in parallel using 144
processors. Similarly to the two atom cell, the 8, 16, and 32 atom cell calculations
were parallelized to reduce their computation time while maintaining accuracy with
respect to cutoff energy and number of k points. The band structures were obtained
from these calculations and the direct and indirect gap energies were extracted for
the various compositional percentages of the germanium-tin alloy. Plotting these
energies, we extract the Sn percentage needed to transition from an indirect to a
direct gap material shown in Figure 6.8. The transition at 8.5% tin is in agreement
within the range predicted by various other methods [138] but we believe this value
to be close to the true value.
In performing these calculations, we have ensured their convergence by in-
creasing the number of k points and the cutoff energy used until the total energy of
each system studied varied less than 0.01 Ry. The systems studied were represen-
tative of uniformly distributed tin in a GeSn crystal which we take to approximate
a real world uniform distribution. From these calculations, we obtained the band
structure for different compositional percentages of Sn and have extracted their var-
ious direct and indirect gap energies. By plotting these, we found the crossing point
which indicates the transition from an indirect gap to occur at 8.5% Sn. At this
148
Figure 6.8: E? and EL calculated for different fractions of Sn
Sn concentration, and presumably all percentages higher than this, the material has
undergone the transition into a direct gap material. By manufacturing the alloy
directly at this transition percent, the band structure will be direct but will also
have the added benefit of maintaining as much of the initial gap value as possible.
This improvement will allow for optical Ge devices to be manufactured with greater
quantum efficiency.
149
Appendix A: DFT Evolution and History
A.1 Hartree-Fock Method
Fock?s extension of Hartree?s theory, the Hartree-Fock method, attempts to
solve the quantum many-body problem by first assuming that the total wavefunc-
tion can be approximated by a single Slater determinant. Electrons are treated
as moving in an external potential from the ions plus a mean field formed by the
other electrons in the system. Via application of the variational principle, a set of
N-coupled equations of the form Eqn. A.1 can be derived and solved iteratively to
produce the approximate orbitals. The effects of electron correlation (interaction of
opposite spin electrons) are neglected in this formulation which results in a system-
atic error. However, the exchange effects of electrons with like-spin are accounted
(for in this theory. ?? )N N ?
?1?2 |?j(r?)|
2 ? ? ?
? i
(r?)?j (r?)?j(r)
+ Vext(r) + dr? ?i(r) dr? = i?i(r)
2 |r? r?| |r? r?|
j 6=i j=6 i
(A.1)
In this formulation, a trial wavefunction is taken as an initial guess. Due to
the variational theorem, the energy expectation of the system taken using the trial
wavefunction will be greater than the energy of the true ground state of the system.
150
The energy is obtained by taking the expectation of the Hartree-Fock Hamiltonion
with the Slater d?eterminant total wavefunction?, resulting in Equation A.2.
1 ? 1 ? ?? ?j(r)??i (r?)?j(r)?HF ? ?2 i(r?)E = ?j(r) ?j(r)dr + drdr? (A.2)2? 2 |r? r?|j=1 ? j=1,i=1? ?1 ?j(r)??i (r?)?i(r)?j(r?) ?? drdr? + ??
2 |r? r?| j
(r)?j(r)Vext(r)dr
j=1,i=1 j=1
Though this technique proved to be an improvement over solving the original
many-body Schrdinger equation, there was still need for a practical application of the
SCF technique due to the complexity still associated with the number of electrons.
A.2 Thomas-Fermi Model
In 1927 Thomas and Fermi created the first model for a quantum electronic
system based purely on the electron density. With this model, instead of solving
for the N wavefunctions of an N-electron system self-consistently, the problem was
reduced to a function of only the 3 spatial dimensions of the density. Though
initially formulated to apply to single atoms, the theory was later also applied to
molecules, crystals and metals [145]. The simplified problem under this formulation
allows solutions to be directly obtained from its variational Euler equation.
In their assumptions they consider an electron gas which satisfies Fermi statis-
tics where electron-electron and electron-nucleus interactions are treated classically
as purely Coulombic in nature. The kinetic energy contribution at each point in the
system is determined by using the exact kinetic energy of a homogeneous electron
gas with density equal to the value at that point, making this a local density ap-
proximation [31]. This kinetic energy per volume is then integrated over the entire
151
system to yield total kinetic energy, shown in equation A.3.
2 ( ) ?3 ~ 2
T [n(r)] = 3?2 3 [n(r)]5/3d3r (A.3)
5 2me
The potential energy terms used in the Thomas-Fermi (TF) model are shown
below in Equation (A.4). From the symmetric form of the electron-electron repulsion
term we can see that despite double counting being accounted for by the 1 term,
2
there exists a self-interaction energy which is erroneously included.
?
Ue?N = n(r)VN(r)d
3r (A.4)
?M Zje2
VN(r) = ?
? |r?Rj=1 j|
1 n(r)n(r?)
Ue?e = e
2 d3rd3r?
2 |r? r?|
Furthermore, the TF model originally did not account for electron exchange
energy which is required to obey the Pauli exclusion principle nor does it account
for the effects of electron correlation. A correction to approximate exchange energy
derived from a homogeneous electron gas model was eventually added by Dirac in
1930 [31].
Due mainly to the simplistic kinetic energy approximation and in part to the
other approximations, the total energy functional given by the TF model (Eqn.
(A.5)) ends up being quite inaccurate and predicts the total energy of a bonded
system to be larger than its constituent parts, making molecules unstable [145].
This is obviously a serious issue and a poor result in most situations so other more
152
sophisticated theories were later developed.
( E)TF [?n(r)] = T [n] + E?e?N [n] + Ee?e[n] ? ? (A.5)22 3
ETF
3h 3 1 n(r)n(r?)
[n(r)] = [n(r)]5/3d3r + n(r)VN(r)d
3r + e2 d3rd3r?
10me 8? 2 |r? r?|
Thomas-Fermi existed as a predecessor to DFT by creating approximate en-
ergy functionals of the electron density in an attempt to simplify the quantum
many-body problem without the mathematical rigor and justification which came
later from the work of Hohenberg and Kohn. This later work showed that a func-
tional approach could be made to reproduce the exact energies of the system instead
of mere approximations [146].
A.3 Hohenberg-Kohn Theory
In their revolutionary 1964 paper [147], Hohenberg and Kohn developed a
theorem forging the groundwork for DFT. Their proofs gave deep insight to the
quantum many-body problem and to potential solution techniques. First, they
proved that there is a one-to-one mapping between the external ionic potential
Vext(r) and the ground state density, as well as a one-to-one between the ground state
density and the total wavefunction (Eqn. A.6). This property allows us to determine
the potential from the density, from which we can construct a Hamiltonian. Since
the Hamiltonian tells us everything about the system including the wavefunction
solutions, the density determines all states of the system [31]. Their second proof
shows that there exists a universal kinetic+interaction functional of the density
F [n(r)] (applying systems with any Vext(r)) which obeys a variational principle.
153
This functional differs from the total ground state energy functional by the external
ionic potential energy as shown in Equation A.7.
?
n(r) = N ??(r, r2, . . . , rN)?(r, r2, . . . , rN)dr2 ? ? ? drN (A.6)
?
EHK [n(r), Vext(r)] = Vext(r)n(r)dr + F [n(r)] (A.7)
EHKGS [n(r)] = minE
HK [n, Vext]
n(r)
?
N = n(r)dr (A.8)
Because this functional obeys a variational principle, it enables us to hunt for
densities which minimize this functional under the constraint that the density inte-
grates total number of electrons in the system (Eqn. A.8). The configuration which
minimizes the energy functional must then be the ground state density from which
we can determine everything about the system using the first proof. By following
this process, the problem is drastically simplified due to the reduction in dimension-
ality. Instead of searching for an unfathomably high dimensional wavefunction, we
can simply work with the density, a function of only 3 spatial dimensions.
Though this theory lays the groundwork for greatly simplifying the problem, it
unfortunately does not give any insight into form of the universal functional. Further
theories go on to approximate the functional in various forms and generally use
equations with clear physical origins which are relatively simple to evaluate. Other
methods solve the problem self-consistently by recasting the problem from an energy
154
minimization into a Schro?dinger-like equation which must be solved consistently such
as in Kohn-Sham DFT.
155
Appendix B: Semiconductor Physics
B.1 Low Doping Limit of Hole Concentration
The following derivation reiterates the assumptions made to compute the hole
concentration in the low-doping limit of non-compensated semiconductor which were
presented in Section 4.2. This result is used in the derivation of the conductivity
mobility parameterization in Section B.1.1.
To begin, charge neutrality is used to compute the hole concentration. For the
system considered, it is assumed that there is no significant donor counter doping.
Additionally, because intrinsic carrier concentration is so small in 4H-SiC, ni is ne-
glected. The resulting charge neutrality equation states that the hole concentration
must equal the ionized acceptor concentration p = N?A . In general, both of these
terms are functions of the unknown Fermi-level EF , as is shown in Section 4.6. In the
low doping limit, however, the governing equations can be manipulated using some
approximations which allow us to remove EF as an independent unknown variable
and provide a closed-form analytical solution.
In the case of low doping, the Fermi level is within the bandgap far from
the valence band edge so the tail of the true Fermi-Dirac occupancy function for
holes is well approximated by the Maxwell-Boltzmann distribution. We apply this
156
approximation to the full integral in Equation B.1 which results in Equation B.2.
4?(2m?)3/2
? ?EV
p EV(? Ep = )dE (B.1)
h3? ?? 1 + exp (EF?EkBT )
4?(2m?)3/2 EV ?
? p ? E ? EFp EV E exp dE (B.2)
h3 ?? kBT
Equation B.2 has an exact solution which can be written as:
( )
EV ? EF
p ? NV(exp ) (B.3)kBT
2?m? (3/2)pkBT
NV = 2 (B.4)
h2
Where NV is the effective density of states in the valence band which depends
on the hole effective mass m?p and the temperature T . Next, we create and evaluate
the integral for the ionized acceptor density N?A using the modified Fermi-Dirac
statistics by including the acceptor state degeneracy gA. For low doped samples,
dopant atoms are far apart and thus all impurity states introduced have virtually
the same energy (assuming small or no inequivalent site energy-dependence). This
creates a density of impurity states in the form of a delta-function at the ionization
energy EA i.e. ?i(E) = NA?(E ? EA).
? ?
? NA?(E ? EA)N = ( )A dE (B.5)
?? 1 + gA exp
E?EF
kBT
N?
NA
A = ( ) (B.6)
1 + gA exp
EA?EF
kBT
157
Finally, we can rewrite Equation B.6 in terms of the hole concentration given
in Equation B.3, leaving us with a quadratic equation for the solution of the hole
concentration p.
N? = ( NA) ( )A (B.7)
1 + g exp EV ?EF exp EA?EVA kBT kBT
NA
N? = p = ( )A (B.8)
1 + gAp exp EA?EV
NV kBT
NA
p = p (B.9)
?1 + 2?
p = ?? + ?(2 + 2?N)A (B.10)
NV ??EA? = exp (B.11)
2gA kBT
?EA = EA ? EV (B.12)
Here, ? is an auxiliary variable of known physical parameters, NA is the ac-
ceptor doping density, NV is the valence band effective density of states, ?EA is
the acceptor ionization energy, and gA is the acceptor state degeneracy equal to 4
in 4H-SiC.
B.1.1 Conductivity Mobility Parameterization
To develop a semi-empirical expression which parameterizes this work?s defi-
nition of conductivity mobility (Equation B.13) we need to create an expression to
158
represent the resistivity term, as it is the only unknown parameter.
? 1?Cond (B.13)
q?NA
Using the standard expression for resistivity (Equation B.14), the two un-
knowns are now the hole concentration p and the hole mobility ?p.
1
? = (B.14)
qp?p
The hole mobility ?p is taken to have the same form as the Hall mobility
?p = ?
c
Cond0/(1+(NA/NC) ) and p is calculated u?sing the solution to the low doping
charge neutrality equation of the form p = ?? + ?2 + 2?NA calculated in Section
4.2. This expression is then substituted for ? into the definition of conductivity
mobility (Eq. B.13). This provides the first part of our empirical conductivity
mobility in Equation 4.60. An additional impurity mobility term is then added
which has the form of a Gaussian because of the approximate shape of the observed
resistivity data for large NA.
B.2 Modified Fermi-Dirac for Dopant States
To calculate the number of ionized acceptor states, we must multiply the ac-
ceptor impurity density of states ?i(E) by the probability of electron occupancy of
said states. The statistical occupation function is essentially a Fermi-Dirac distribu-
tion, except there is a degeneracy factor which appears because impurity states can
only be singly charged as the acceptor may only accept one electron and a donor can
159
only loose one electron due to the large Coulombic energy associated with double
ionization.
For an acceptor, the gained electron can either be spin up or spin down without
changing the energy of the system (to a good approximation), which creates two
degenerate configurations. Additionally, since we are dealing with acceptor states,
the ?accepted? electron is coming from the valence band maximum located at the
?-point, where there is a degeneracy of the light-hole and heavy-hole bands in SiC.
This additional degeneracy multiples the previous degeneracy, leading to a 4-fold
degeneracy which should be included into the Fermi-Dirac function for acceptor
state occupancy.
As a starting point, the mathematical derivation using the micro-canonical
ensemble for the standard Fermi-Dirac proceeds as follows:
Thermodynamics states that the most likely configuration for a system will be
the macrostate with the largest number of microstates. First we will enumerate the
total number of microstates of the system, then using Lagrange multipliers, we will
apply the constraints of total particle number and total energy to find the occupancy
function which defines this macrostate.
At each energy Ei, we want to fill the ni degenerate states with mi electrons
where 0 ? mi ? ni. We can represent the number of electrons in each state using
an occupancy function fi (a fraction < 1) which multiplies the ni states mi = fini.
Here, fi represents the probability of occupying fini states at energy Ei.
To determine fi, we first enumerate the number of ways these states and
electrons can be arranged - knowing that each state can only hold one electron. The
160
number of distinct configurations Ci for each energy Ei is calculated by dividing the
number of ways to rearrange all of the states, by the number of ways to rearrange
the filled states, and by the number of ways to rearrange the empty states. This is
equivalent to enumerating the number of ways to fill ni states with mi electrons:
ni! ni!
Ci = = (B.15)
mi!(ni ?mi)! (fini)!(ni ? fini)!
Because the electrons and the states at a given energy are indistinguishable
amongst themselves, the two factorial terms appear in the denominator. To find
the most likely distribution of f across all energy states, we must multiply the
number of configurations at each energy to get the total number of configurations
or microstates for the entire system. Then, to obtain a more tractable solution, we
apply Sterling?s approximation ln(n!) ? n ln(n) ? n which is valid due to the large
number of electrons and states.
? ? ni!
C = Ci = (B.16)
?(fi i ini)!(ni ? fini)!
ln(C) = lnCi = (B.17)
? i
(ni ln(ni)? (fini) ln(fini)? (ni ? fini) ln(ni ? fini))
i
The system is also constrained by the total number of electrons N and the
total energy U defined by:
161
?
N =? fini (B.18)i
U = Eifini (B.19)
i
Finally we can determine the distribution f that maximizes the number of con-
figurations, which thermodynamics ensures will be the most probable distribution
in thermal equilibrium. To do this, we maximize the configuration function C (or
ln(C)) subject to our constraint equations using the method of Lagrange multipliers.
[
? ? ? ]
ln(C)? a fjnj ? b E( ) j
fjnj = 0 (B.20)
?fi j j
1? fi
ln ? a? bEi = 0 (B.21)
fi
1
fi = (B.22)
1 + exp(a+ bEi)
Further analysis of the units in the variational equation and relating the vari-
ational energy term to the change in energy due to the change in entropy allows us
to notice that b must be 1/kBT to satisfy dU =
1d(ln(C)) = Td(kB ln(C)) = TdSb
from thermodynamics. Analogously it can be shown that a must be ?EF/kT from
(?a/b)dN = ?dN , where EF is called the Fermi energy and ? is the electro-chemical
potential of the system. Finally, we rewrite our distribution function, known as the
Fermi-Dirac distribution:
1
F (E) = ( ) (B.23)
1 + exp E?EF
kBT
162
This equation is valid for the occupancy of Fermions (namely electrons) which
have only one way to occupy a state - such as in the conduction band of SiC and
other semiconductors. Because the density of states used in conjunction with this
formula already accounts for (i.e. ?labels?) states with ?up? and ?down? spin, con-
duction band (and valence band) states can only be unoccupied or singly occupied
and the spin of the electron is tied to the state. When determining the distribution
function for acceptor or donor state occupancy, there are different ways in which the
state may be occupied. These extra ways of occupying the state lead to a degeneracy
(depending on how many extra ways) which must be included in the distribution
formulation. The density of states associated with a doped impurity does not have
spin association built in to it?s states. For example, in the case of donors, a single
doping level at ED may have a Gaussian density of states gi(E) centered around
ED which integrates to the donor concentration ND. With this density of states,
each state is either singly occupied with either a spin-up or spin-down electron (in-
dependent of how every other state is occupied) or singly ionized with no electron
present. Each state is not ?labeled? as a ?spin-up? or ?spin-down? state beforehand
as is the case with the conduction band and valence band density of states via the
factor of 2 included for spin in their values. Therefore each electron occupying a
donor state may do so in two ways: with either spin ?up? and with spin ?down?. This
choice manifests itself in increasing the total number of possible ways to fill the ni
donor states with the mi donor electrons. The total number of ways to assign the
states is thus doubled for each of the mi electrons which occupy the donor states
due to the choice of each being filled with a spin ?up? or spin ?down? electron. Any
163
electron chosen to not occupy a donor state will instead occupy a conduction band
state. The total number of electrons in the system is limited-by and equal-to the
total donor concentration.
So, for electrons occupying donor states:
2min ! 2finii ni!
Ci = = (B.24)
mi!(ni ?mi)! (fini)!(n? fini)!
Carrying out the same derivation leads to the equation:
( )
ln 2 ? 1? fi ? a? bEi = 0 (B.25)
fi
which results in a modified Fermi-Dirac function to represent the probability
of occupancy of donor states:
1
FD(E) = ( ) (B.26)
1 + 1 exp E?EF
2 kBT
This occupancy function can be used directly with the donor density of states
to calculate the unionized donor concentration N0D (i.e. number of donor states
which are still are occupied by electrons). To calculate the often more useful ionized
donor concentration N+D we must use 1? FD(E):
1
1? FD(E) = ( ) (B.27)
1 + 2 exp EF?E
kBT
Similarly, for acceptors, there are two choices to make when filling each accep-
tor state with a hole. The electrons still in the valence band may not only be spin
?up? or spin ?down? but also may be residing in the heavy-hole or a light-hole band.
164
This effect occurs due to the light-hole/heavy-hole valence band degeneracy at the
gamma point which is common among many semiconductors (including SiC). These
choices multiply the number of ways to choose and fill mi of ni acceptor states with
electrons by 4 for each of the ni ?mi holes:
4ni?min ! 4ni?finii ni!
Ci = = (B.28)
mi!(ni ?mi)! (fini)!(n? fini)!
resulting in the modified Fermi-Dirac function for the electron occupancy of
acceptor states for calculating N?A :
1
FA(E) = ( ) (B.29)
1 + 4 exp E?EF
kBT
B.3 Valence Band Density of States
To derive the density of states in the valence band, we assume that the holes
are free particles confined to a crystal with side lengths L. From Bloch?s theorem,
we know that the wavefunction may be expressed as:
?k(r) = e
ik?ruk(r) (B.30)
kxL = 2?nx (B.31)
kyL = 2?ny (B.32)
kzL = 2?nz (B.33)
n = . . . ,?2,?1, 0, 1, 2, . . . (B.34)
165
with the stipulation that k must obey periodic boundary conditions with the
crystal boundary. From this, we know that in reciprocal space, there is one state
per (2?/L)3. Using the dispersion relation for a free particle:
~2|k|2
Ek = ? (B.35)2m
From this we can calculate the number of states N within a spherical volume
of radius |k|. The number must be divided by 8 because we must account for the
equivalence of ?kx,y,z values, which represent a phase shift in the wavefunction but
are actually the same state. To account for the spin degeneracy of the states we
multiply by two.
(
| |3 ?
)3/2
2 (4/3)? k V V 2m E
N = = |k|3 = (B.36)
8 (2?/L)3 3?2 3?2 ~2
Here, V is the volume of the crystal. To get the density of states DoS(E) (per
volume) we take the derivative of N with respect to E and divide by V :
(
dN 1 2m?
)3/2?
?v0(E) = = E (B.37)
dE 2?2 ~2
4?(2m?)3/2?
= E
h3
B.3.1 Dopant DoS Spreading
To approximate the band width of the acceptor impurity states, we use the
tight binding model and assume hydrogenic s-like wavefunctions [91] (?0(r)) to eval-
uate the energy transfer integral for states associated with atoms at locations Ri
166
and Rj.
?
?3
?0(r ?Ri) = ( exp)(?? |r ?Ri|) (B.38)? 1/2
1 EA
? ? = (B.39)a0 E0
q2
J (|Ri ?Rj|) = ? (r ?R )? (r ?R )d3r (B.40)
4?r0 |r ?Ri|
0 i 0 j
q2?
J(R) = (1 + ?R) exp(??R) (B.41)
4?r0
Here, R = |Ri ?Rj| is the distance between nearest neighbor dopant atoms,
a0 is the Bohr radius, ? is the scaled inverse radius for the acceptor state, and E0 is
the ground state energy of the hydrogen atom. Assuming the dopants are uniformly
randomly distributed in the crystal (reasonable assumption for a low to moderately
doped box-like profile), the probability that the nearest neighbor dopant lies at a
distance between R and R + dR from a given dopant is given by an exponential
distribution.
( )
4
4?N 3A exp ? ?NAR R2dR (B.42)
3
This probability is derived from the Poisson probability distribution of finding
m atoms in a volume w with uniform density n (taken from Kane [104]).
(nw)m
P (m,w) = e?nw (B.43)
m!
Averaging the energy transfer integral weighted by the nearest neighbor dis-
tance probability distribution leads to the total bandwidth of the impurity levels
167
B.
? ( )
?J(R)? = J(R)4?N R2 ?4exp ?N 3A AR dR (B.44)
3
B = 2 |?J(R)?| (B.45)
B.4 Compensated Systems
For the case of compensated systems (containing both NA and ND concentra-
tions), we must go back to the full charge-neutrality equation.
p+N+D = n+N
?
A (B.46)
p+N+ 2 ?D = ni /p+NA (B.47)
p2 + (N+D ?N
?
A )p? n
2
i (B.48)
Here, p, N+ ?D , and NA are all functions of the unknown Fermi level EF :? EV ?V ((E)p = )dE (B.49)
??? 1 + exp EF?EkBT?
N?
?A(E)
= ( )A dE (B.50)
??? 1 + gA exp E?EFkBT?
+ ?D(E()N )D = dE (B.51)
?? 1 + gD exp
EF?E
kBT
Here, ?V (E), ?A(E), and ?D(E) are the valence band, acceptor, and donor density
of states. As before, this equation can be solved using the quadratic formula for p.
When these integrals are substituted into Equation B.48, it can be rewritten such
168
that only two integrals need to be solved for each iteration of the numerical solution.
?
?D?(EF ) + (D?(EF ))2 + 4n2
p(E ?) = iF (B.52)2?
? ? + ? ? ?D(E() ) ? ?A(E()D (EF ) ND NA = )dE (B.53)
?? 1 + g EF?ED exp 1 + gA exp
E?EF
kBT kBT
This equation can be solved numerically for EF by using trapezoidal integra-
tion combined with the method of bisection for a given doping NA and ND.
169
Appendix C: Standard Component Mobility Formulations
C.1 Bulk Mobility
Bulk mobility is modeled using the dopant and temperature dependent Caughey-
Thomas model described by Equation (C.1) [123].
? (300)?max
?B =
T (C.1)
1 + (D(T ))?
Nref
Where the values for bulk mobility parameters are: ? 2max is 1071 cm /Vs; ? is
2.4; N is 1.9e17 cm?3ref ; and ? is 0.4.
C.2 Surface Phonon Mobility
Surface phonons, which partially account for the reduction in MOSFET surface
mobility, are calculated using the acoustic phonon deformation potential and surface
field. The analytical form is described by Equation (C.2) [148], where ?bulk is the
bulk material density of SiC; vs is the velocity of sound in SiC; m
?, mc, m? are the
(2D) density of states, conductivity, and perpendicular effective masses, respectively;
Dac is the acoustic phonon deformation potential; e is the elementary charge; ~ is
170
the reduced Planck?s constant; and kB is Boltzmann?s constant.
A B
?SP = + (C.2)
F? 1/3TF?
3 ~3?bulkv2
A = s
2m?m( 2cDac )
e~3? v2 (1/3)bulk s 9~2B =
m?mcD2ackB 4em?
C.3 Combined Empirical Surface Roughness Mobility
The CESRM is described using the analytical fit given in Equation (C.3) [1].
Step roughness in epitaxially grown 4H-SiC miscut 8?from the (0001) plane is on
average comprised of bunched steps resulting in an rms roughness height of between
0.38 nm [149] and 3.5 nm [51]. Compared to the single bilayer step height of 0.24
nm it is clear that the surface roughness typically accounted for in literature on this
matter is predominantly affected by the large-scale step bunching effects rather than
the smaller atomic-scale roughness limit we distinguish in this paper. Experimental
measurements performed by Kimoto et al. for 3.5?off-angle wafers confirm that the
most probable step consists of 4 Si-C bilayers for the Si-face [150]. For relevant 4H-
SiC MOSFETs we use ?SR = 3.5? 1012 V/s extracted by Potbhare et al. [1, 51] In
the calculation, ? is the RMS surface height variation, L is the roughness correlation
171
length, qsc is the screening wavevector, and ? is the scattering angle.
?SR
(?SR = ) (C.3)F 2?
~3 1
?SR =
2m ?cm e?2L2 ?SR
? ?SR =?/2 ( sin)3((?) )d?
0 ?sin(?) + ? qsc 1 + sin2(?)L2 4m kBT
? ~2 ~28m kBT/
C.4 Coulomb Mobility
Trapped and fixed charges caused by atomic defects at the interface create
charged scattering sites which induce a Coulombic scattering effect for the channel
electrons. The Coulomb mobility is treated using Equation (C.4) [1] where z is the
depth into the device. The parameters for fixed interface sheet charge density Nf ,
interface trapped sheet charge density Nit, surface inversion sheet charge density
Ninv, and inversion layer depth Zav are given in Table C.1.
Inversion layer charge as a function of surface field was derived from the work
of Arnold [2] and matched to the corresponding interface state density calculated
by Potbhare et al. [1]. The inversion layer depth Zav was taken to be the expected
depth of the electron density using the lowest energy (1st) Airy function solution
?1(z). The resulting Coulomb mobilities are a function of z so the values reported
in Table 3.3 are evaluated at their corresponding Zav.
172
16?2~kBT
?C =( ? (C.4)? m e3(Nf +Nit)f(z)pi/2 )
q2
f(z) = 1? sc
( ?? 8m kBT~2 sin2(?) +0 )q2sc
?
exp ? 8m kBT2z ? sin
2(?) + q2sc d?~2
e2Ninv
qsc =
SiCZavkBT
Table C.1: Extrinsic Coulomb Scattering Mobility Parameters [1, 2]
E? [MV/cm] 0.1 0.5 1
N [cm?2f ] 1.3? 1012
N [cm?2] 1? 1012 3.35? 1012 3.55? 1012it
N [cm?2inv ] 4? 1010 1.5? 1012 3.65? 1012
Zav [nm] 3.77 2.21 1.75
173
Bibliography
[1] Siddharth Potbhare, Neil Goldsman, Gary Pennington, Aivars Lelis, and
James M McGarrity. Numerical and experimental characterization of 4 h-
silicon carbide lateral metal-oxide-semiconductor field-effect transistor. Jour-
nal of Applied Physics, 100(4):044515, 2006.
[2] Emil Arnold. Charge-sheet model for silicon carbide inversion layers. IEEE
Transactions on Electron Devices, 46(3):497?503, 1999.
[3] Shizuo Fujita. Wide-bandgap semiconductor materials: For their full bloom.
Japanese journal of applied physics, 54(3):030101, 2015.
[4] P.A. Tipler and R.A. Llewellyn. Modern Physics. W.H. Freeman and Com-
pany, 5 edition, 2008.
[5] Y Umeno, Y Kinoshita, and T Kitamura. Ab initio dft simulation of ideal
shear deformation of sic polytypes. Modelling and Simulation in Materials
Science and Engineering, 15(2):27, 2006.
[6] S. Adachi. Properties of Group-IV, III-V and II-VI Semiconductors. John
Wiley & Sons, Ltd, 2005.
[7] C Persson and Ulf Lindefelt. Relativistic band structure calculation of cubic
and hexagonal sic polytypes. Journal of Applied Physics, 82(11):5496?5508,
1997.
[8] Michael E Levinshtein, Sergey L Rumyantsev, and Michael S Shur. Properties
of Advanced Semiconductor Materials: GaN, AIN, InN, BN, SiC, SiGe. John
Wiley & Sons, 2001.
[9] Shahrzad Salemi. Electronic Structure of SiC/SiO2 by Density Functional
Theory. PhD thesis, University of Maryland, College Park, 2012.
[10] V. Heera, D. Panknin, and W. Skorupa. p-Type doping of SiC by high dose Al
implantation - problems and progress. Applied Surface Science, 184(1-4):307?
316, 2001.
174
[11] Imran A Khan and James A Cooper. Measurement of high-field electron trans-
port in silicon carbide. IEEE Transactions on Electron Devices, 47(2):269?273,
2000.
[12] GE Stillman, VM Robbins, and N Tabatabaie. Iii-v compound semicon-
ductor devices: Optical detectors. IEEE Transactions on Electron Devices,
31(11):1643?1655, 1984.
[13] B Jayant Baliga. Sic power devices: From conception to social impact. In 2016
46th European Solid-State Device Research Conference (ESSDERC), pages
192?197. IEEE, 2016.
[14] Devanarayanan Perinthatta Ettisserry. INTEGRATED MODELING OF RE-
LIABILITY AND PERFORMANCE OF 4H-SILICON CARBIDE POWER
MOSFETS USING ATOMISTIC AND DEVICE SIMULATIONS. PhD the-
sis, University of Maryland, College Park, 2015.
[15] Anant Agarwal and Sarah Haney. Some critical materials and processing issues
in sic power devices. Journal of Electronic Materials, 37(5):646?654, 2008.
[16] Karl E Spear and John P Dismukes. Synthetic diamond: emerging CVD
science and technology, volume 25. John Wiley & Sons, 1994.
[17] Robert Eisberg and Robert Resnick. Quantum physics of atoms, molecules,
solids, nuclei, and particles. Quantum Physics of Atoms, Molecules, Solids,
Nuclei, and Particles, 2nd Edition, by Robert Eisberg, Robert Resnick, pp.
864. ISBN 0-471-87373-X. Wiley-VCH, January 1985., page 864, 1985.
[18] R. Serway, C. Moses, and C. Moyer. Modern Physics. Thomson Brooks, 5
edition, 2005.
[19] Hangseok Choi et al. Overview of silicon carbide power devices. Fairchild
semiconductor, 2016.
[20] Robert K Willardson and Eicke R Weber. SiC materials and devices, vol-
ume 52. Academic Press, 1998.
[21] Yole Developpement. Sic, gan and other wider-bandgap materials present new
choices for power electronics. Market Focus: Power Electronics 91, 9(8), 2015.
[22] Gang Liu, Blair R Tuttle, and Sarit Dhar. Silicon carbide: A unique platform
for metal-oxide-semiconductor physics. Applied Physics Reviews, 2(2):021307,
2015.
[23] Lijuan Li, Canbing Li, Yijia Cao, and Feng Wang. Recent progress of sic
power devices and applications. IEEJ Transactions on Electrical and Elec-
tronic Engineering, 8(5):515?521, 2013.
175
[24] Tristan Evans, Toshio Hanada, Yuki Nakano, and Takashi Nakamura. De-
velopment of sic power devices and modules for automotive motor drive use.
In 2013 IEEE International Meeting for Future of Electron Devices, Kansai,
pages 116?117. IEEE, 2013.
[25] Wai-Kai Chen. VLSI Handbook Second Edition. Taylor & Francis Group,
2007.
[26] Takahide Umeda, Mitsuo Okamoto, Ryouji Kosugi, Shinsuke Harada, Ryo
Arai, Yoshihiro Sato, Takahiro Makino, and Takeshi Ohshima. Sic mos inter-
face states: Difference between si face and c face. ECS Transactions, 58(7):55?
60, 2013.
[27] Siddharth Potbhare. Modeling and characterization of 4H-SiC MOSFETS:
High field, high temperature, and transient effects. PhD thesis, University of
Maryland, College Park, 2008.
[28] J. Fan and P.K. Chu. Silicon Carbide Nanostructures. Engineering Materials
and Process, 2014.
[29] Subal Sahni. Highly integrated Germanium Photo-detectors and III-V Hybrid
Lasers for Silicon Photonic Applications. PhD thesis, University of California
Los Angeles, 2007.
[30] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto
Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococ-
cioni, Ismaila Dabo, Andrea Dal Corso, Stefano de Gironcoli, Stefano Fabris,
Guido Fratesi, Ralph Gebauer, Uwe Gerstmann, Christos Gougoussis, An-
ton Kokalj, Michele Lazzeri, Layla Martin-Samos, Nicola Marzari, Francesco
Mauri, Riccardo Mazzarello, Stefano Paolini, Alfredo Pasquarello, Lorenzo
Paulatto, Carlo Sbraccia, Sandro Scandolo, Gabriele Sclauzero, Ari P Seitso-
nen, Alexander Smogunov, Paolo Umari, and Renata M Wentzcovitch. Quan-
tum espresso: a modular and open-source software project for quantum sim-
ulations of materials. Journal of Physics: Condensed Matter, 21(39):395502
(19pp), 2009.
[31] Robert O Jones. Density functional theory: Its origins, rise to prominence,
and future. Reviews of modern physics, 87(3):897, 2015.
[32] Douglas R Hartree. The wave mechanics of an atom with a non-coulomb
central field. part i. theory and methods. In Mathematical Proceedings of
the Cambridge Philosophical Society, volume 24, pages 89?110. Cambridge
University Press, 1928.
[33] Vladimir Fock. Na?herungsmethode zur lo?sung des quantenmechanischen
mehrko?rperproblems. Zeitschrift fu?r Physik, 61(1-2):126?148, 1930.
[34] John C Slater. Note on hartree?s method. Physical Review, 35(2):210, 1930.
176
[35] John P Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient
approximation made simple. Physical review letters, 77(18):3865, 1996.
[36] John P Perdew. Density functional theory and the band gap problem. Inter-
national Journal of Quantum Chemistry, 28(S19):497?523, 1985.
[37] D Volm, BK Meyer, DM Hofmann, WM Chen, NT Son, C Persson, Ulf Lin-
defelt, O Kordina, E So?rman, AO Konstantinov, et al. Determination of the
electron effective-mass tensor in 4h sic. Physical Review B, 53(23):15409, 1996.
[38] P. Friedrichs, T. Kimoto, L. Ley, and G. Pensl. Silicon carbide volume 2:
Power devices and sensors, 2009.
[39] G Pennington and N Goldsman. Simulation of electron transport in (0001)
and (1120) 4H-SiC inversion layers. Journal of Applied Physics, 106(6):063701,
2009.
[40] Mats Hjelm. Monte Carlo simulations of homogeneous and inhomogeneous
transport in silicon carbide. PhD thesis, Royal Institute of Technology, Sweden,
2004.
[41] Mats Hjelm, Hans-Erik Nilsson, A Martinez, KF Brennan, and E Bellotti.
Monte carlo study of high-field carrier transport in 4h-sic including band-to-
band tunneling. Journal of applied physics, 93(2):1099?1107, 2003.
[42] Carlo Jacoboni and Lino Reggiani. The monte carlo method for the solution
of charge transport in semiconductors with applications to covalent materials.
Reviews of modern Physics, 55(3):645, 1983.
[43] Carlo Jacoboni and Paolo Lugli. The Monte Carlo method for semiconductor
device simulation. SSprinter-Verlag, Vienna, Austria, 1989.
[44] Shinya Yamakawa, Hiroaki Ueno, Kenji Taniguchi, Chihiro Hamaguchi, Kazuo
Miyatsuji, Kazuo Masaki, and Umberto Ravaioli. Study of interface roughness
dependence of electron mobility in si inversion layers using the monte carlo
method. Journal of applied physics, 79(2):911?916, 1996.
[45] Wolfgang J Choyke, Hiroyuki Matsunami, and Gerhard Pensl. Silicon carbide:
recent major advances, volume 1. Springer Science & Business Media, 2013.
[46] SM Goodnick, DK Ferry, CW Wilmsen, Z Liliental, D Fathy, and OL Kri-
vanek. Surface roughness at the si (100)-sio 2 interface. Physical Review B,
32(12):8171, 1985.
[47] Tsuneya Ando, Alan B Fowler, and Frank Stern. Electronic properties of
two-dimensional systems. Reviews of Modern Physics, 54(2):437, 1982.
[48] L Farhang Matin, H Hasan Bouzari, and F Ahmadi. Solving schrodinger
equation specializing to the stark effect in linear potential by the canonical
function method. Journal of Theoretical and Applied Physics, 8(3):140, 2014.
177
[49] SM Goodnick, DK Ferry, CW Wilmsen, Z Liliental, D Fathy, and OL Kri-
vanek. Surface roughness at the Si (100)-SiO2 interface. Physical Review B,
32(12):8171, 1985.
[50] Han Fu, KV Reich, and BI Shklovskii. Surface roughness scattering in multi-
subband accumulation layers. Physical Review B, 93(23):235312, 2016.
[51] Siddharth Potbhare, Neil Goldsman, Aivars Lelis, James M McGarrity,
F Barry McLean, and Daniel Habersat. A physical model of high temperature
4H-SiC MOSFETs. IEEE Transactions on Electron devices, 55(8):2029?2040,
2008.
[52] Viktoryia Uhnevionak, Alexander Burenkov, Christian Strenger, Guillermo
Ortiz, Elena Bedel-Pereira, Vincent Mortet, Fuccio Cristiano, Anton J Bauer,
and Peter Pichler. Comprehensive study of the electron scattering mechanisms
in 4H-SiC MOSFETs. IEEE Transactions on Electron Devices, 62(8):2562?
2570, 2015.
[53] K Kojima, H Okumura, S Kuroda, and K Arai. Homoepitaxial growth of 4H-
SiC on on-axis (0001) C-face substrates by chemical vapor depositon. Journal
of crystal growth, 269(2-4):367?376, 2004.
[54] J Hassan, JP Bergman, A Henry, and E Janze?n. On-axis homoepitaxial growth
on Si-face 4H?SiC substrates. Journal of Crystal Growth, 310(20):4424?4429,
2008.
[55] Stefano Leone, Franziska C Beyer, Henrik Pedersen, Olof Kordina, Anne
Henry, and Erik Janze?n. High growth rate of 4H-SiC epilayers on on-axis
substrates with different chlorinated precursors. Crystal Growth & Design,
10(12):5334?5340, 2010.
[56] Yasuto Hijikata, Hiroyuki Yaguchi, Sadafumi Yoshida, Y Takata,
K Kobayashi, H Nohira, and T Hattori. Off-Angle dependence of charac-
teristics of 4H-SiC-Oxide interfaces. In Materials science forum, volume 527,
pages 1003?1006. Trans Tech Publ, 2006.
[57] Kenji Fukuda, Makoto Kato, Shinsuke Harada, and Kazutoshi Kojima. High
inversion channel mobility of 4H-SiC MOSFETs fabricated on c (000-1) epi-
taxial substrate with vicinal (below 1?) off-angle. In Materials science forum,
volume 527, pages 1043?1046. Trans Tech Publ, 2006.
[58] Shinsuke Harada, Sachiko Ito, Makoto Kato, Akio Takatsuka, Kazutoshi Ko-
jima, Kenji Fukuda, and Hajime Okumura. Isotropic channel mobility in
UMOSFETs on 4H-SiC C-face with vicinal off-angle. In Materials Science
Forum, volume 645, pages 999?1004. Trans Tech Publ, 2010.
178
[59] Hiroshi Yano, H Nakao, Tomoaki Hatayama, Yukiharu Uraoka, and Takashi
Fuyuki. Increased channel mobility in 4h-sic UMOSFETs using on-axis sub-
strates. In Materials science forum, volume 556, pages 807?810. Trans Tech
Publ, 2007.
[60] T Kimoto and J.A. Cooper. Appendix A: Incomplete Dopant Ionization in
4H-SiC, pages 511?515. John Wiley & Sons, Ltd, 2014.
[61] Yu.A. Vodakov, E.N. Mokhov, M.G. Ramm, and A.D. Roenkov. Doping pe-
culiarities of SiC epitaxial layers grown by sublimation sandwich-method. In
Springer Proc. Phys, volume 56, pages 329?334, 1992.
[62] Yu.A. Vodakov, E.N. Mokhov, R.C. Marshall, J.W. Faust, and C.E. Ryan.
Silicon Carbide, 1973. University of South Carolina Press, Columbia, SC,
pages 508?19, 1974.
[63] Y.M. Tairov and Y.A. Vodakov. Group IV materials (mainly SiC). In J.I.
Pankove, editor, Electroluminescence, pages 31?61. Springer Berlin Heidelberg,
1977.
[64] M.K. Linnarsson, U. Zimmermann, J. Wong-Leung, A. Scho?ner, M.S. Jan-
son, C. Jagadish, and B.G. Svensson. Solubility limits of dopants in 4H?SiC.
Applied Surface Science, 203:427?432, 2003.
[65] M.V. Rao, J.B. Tucker, M.C. Ridgway, O.W. Holland, N. Papanicolaou, and
J. Mittereder. Ion-implantation in bulk semi-insulating 4H?SiC. Journal of
Applied Physics, 86(2):752?758, 1999.
[66] T. Troffer, M. Schadt, T. Frank, H. Itoh, G. Pensl, J. Heindl, H.P. Strunk, and
M. Maier. Doping of SiC by Implantation of Boron and Aluminum. Physica
Status Solidi (a), 162(1):277?298, 1997.
[67] I.G. Atabaev, Kh.N. Juraev, and M.U. Hajiev. Spectral Dependence of Op-
tical Absorption of 4H-SiC Doped with Boron and Aluminum. Journal of
Spectroscopy, 2018, 2018.
[68] J. Lutz, H. Schlangenotto, U. Scheuermann, and R. De Doncker. Semiconduc-
tor Power Devices Physics, Characteristics, Reliability. Springer International
Publishing, 2 edition, 2018.
[69] S. Greulich-Weber. EPR and ENDOR Investigations of Shallow Impurities in
SiC Polytypes. Physica Status Solidi (a), 162(1):95?151, 1997.
[70] M. Ikeda, H. Matsunami, and T. Tanaka. Site effect on the impurity levels in
4H, 6H, and 15R SiC. Physical Review B, 22:2842, 09 1980.
[71] T Kimoto and J.A. Cooper. Appendix C: Major Physical Properties of Com-
mon SiC Polytypes, pages 521?524. John Wiley & Sons, Ltd, 2014.
179
[72] M. Bakowski, U. Gustafsson, and U. Lindefelt. Simulation of SiC high power
devices. Physica Status Solidi (a), 162(1):421?440, 1997.
[73] S.G. Sridhara, L.L. Clemen, R.P. Devaty, W.J. Choyke, D.J. Larkin, H.S.
Kong, T. Troffer, and G. Pensl. Photoluminescence and transport studies of
boron in 4H SiC. Journal of Applied Physics, 83(12):7909?7919, 1998.
[74] T. Troffer, C. Ha??ler, G. Pensl, K. Ho?lzlein, H. Mitlehner, and J. Vo?lkl. Boron-
related defect centers in 4H silicon carbide. In International Conference on
Silicon Carbide and Related Materials 1995, number 142, pages 281?284. IOP
Publishing Ltd., 1996.
[75] A.A. Lebedev. Deep level centers in silicon carbide: A review. Semiconductors,
33(2):107?130, 1999.
[76] S.G. Sridhara, L.L. Clemen, R.P. Devaty, W.J. Choyke, D.J. Larkin, H.S.
Kong, T. Troffer, and G. Pensl. Photoluminescence and transport studies of
boron in 4H SiC. Journal of Applied Physics, 83(12):7909?7919, 1998.
[77] M.M. Anikin, A. Lebedev, A.L. Syrkin, and A.V. Suvorov. Investigation of
deep levels in SiC by capacitance spectroscopy methods. Soviet Physics of
Semiconductors, 19(1):69?71, 1985.
[78] T. Troffer, G. Pensl, A. Scho?ner, A. Henry, C. Hallin, O. Kordina, and
E. Janze?n. Electrical characterization of the gallium acceptor in 4H-and 6H-
SiC. Materials Science Forum, 264:557?560, 1998.
[79] G.L. Harris and INSPEC (Information service). Properties of Silicon Carbide.
EMIS Datareviews Series. INSPEC, Institution of Electrical Engineers, 1995.
[80] A.-B. Chen and P. Srichaikul. Shallow Donor Levels and the Conduction Band
Edge Structures in Polytypes of SiC. Physica Status Solidi (b), 202(1):81?106,
1997.
[81] S. Kagamihara, H. Matsuura, T. Hatakeyama, T. Watanabe, M. Kushibe,
T. Shinohe, and K. Arai. Parameters required to simulate electric characteris-
tics of SiC devices for n-type 4HSiC. Journal of Applied Physics, 96(10):5601?
5606, 2004.
[82] A.O. Evwaraye, S.R. Smith, and W.C. Mitchel. Shallow and deep levels in
ntype 4HSiC. Journal of Applied Physics, 79(10):7726?7730, 1996.
[83] T. Kimoto, A. Itoh, H. Matsunami, S. Sridhara, L. L. Clemen, R. P. Devaty,
W. J. Choyke, T. Dalibor, C. Peppermller, and G. Pensl. Nitrogen donors and
deep levels in highquality 4HSiC epilayers grown by chemical vapor deposition.
Applied Physics Letters, 67(19):2833?2835, 1995.
180
[84] W. Go?tz, A. Scho?ner, G. Pensl, W. Suttrop, W. J. Choyke, R. Stein, and
S. Leibenzeder. Nitrogen donors in 4Hsilicon carbide. Journal of Applied
Physics, 73(7):3332?3338, 1993.
[85] N. Donato and F. Udrea. Static and Dynamic Effects of the Incomplete Ion-
ization in Superjunction Devices. IEEE Transactions on Electron Devices,
65(10):4469?4475, 2018.
[86] M. Laube, F. Schmid, G. Pensl, G. Wagner, M. Linnarsson, and M. Maier.
Electrical activation of high concentrations of N+ and P+ ions implanted into
4HSiC. Journal of Applied Physics, 92(1):549?554, 2002.
[87] R. Wang, I.B. Bhat, and T.P. Chow. Epitaxial growth of n-type SiC us-
ing phosphine and nitrogen as the precursors. Journal of Applied Physics,
92(12):7587?7592, 2002.
[88] S. Rao, T.P. Chow, and I. Bhat. Dependence of the Ionization Energy of
Phosphorous Donor in 4H-SiC on Doping Concentration. Materials Science
Forum, 527-529:597?600, 01 2006.
[89] S.M. Sze and K.K. Ng. Physics of Semiconductor Devices. John Wiley &
Sons, 2006.
[90] G.W. Neudeck and R.F. Pierret. Advanced Semiconductor Fundamentals.
Pearson Education, New York, 2003.
[91] H. Brooks. Theory of the Electrical Properties of Germanium and Silicon.
volume 7 of Advances in Electronics and Electron Physics, pages 85 ? 182.
Academic Press, 1955.
[92] E.F. Schubert. Doping in III-V Semiconductors. Cambridge University Press,
1993.
[93] H. Ikeda and F. Salleh. Influence of heavy doping on Seebeck coefficient in
silicon-on-insulator. Applied Physics Letters, 96(1):012106, 2010.
[94] D.S. Lee and J.G. Fossum. Energy-Band Distortion in Highly Doped Silicon.
IEEE Transactions on Electron Devices, 30(6):626?634, 1983.
[95] J.R. Lowney, A.H. Kahn, J.L. Blue, and C.L. Wilson. Disappearance of im-
purity levels in silicon and germanium due to screening. Journal of Applied
Physics, 52(6):4075?4080, 1981.
[96] A. Schenk, P.P. Altermatt, and B. Schmithusen. Physical Model of Incomplete
Ionization for Silicon Device Simulation. In 2006 International Conference on
Simulation of Semiconductor Processes and Devices, pages 51?54. IEEE, Sep.
2006.
181
[97] P.P. Altermatt, A. Schenk, and G. Heiser. A simulation model for the density
of states and for incomplete ionization in crystalline silicon. I. Establishing
the model in Si: P. Journal of Applied Physics, 100(11):113714, 2006.
[98] P.P. Altermatt, A. Schenk, and G. Heiser. A simulation model for the density
of states and for incomplete ionization in crystalline silicon. II. Investigation
of Si:As and Si:B and usage in device simulation. Journal of Applied Physics,
100(11):113715, 2006.
[99] V.I. Fistul. Heavily Doped Semiconductors. Plenum Press, New York, 1969.
[100] T.F. Lee and T.C. McGill. Variation of impurity-to-band activation energies
with impurity density. Journal of Applied Physics, 46(1):373?380, 1975.
[101] I.G. Ivanov, A. Henry, and E. Janze?n. Ionization energies of phosphorus and
nitrogen donors and aluminum acceptors in 4H silicon carbide from the donor-
acceptor pair emission. Physical Review B, 71:241201, Jun 2005.
[102] A. Parisini and R. Nipoti. Analysis of the hole transport through valence
band states in heavy Al doped 4H-SiC by ion implantation. Journal of Applied
Physics, 114(24):243703, 2013.
[103] G. Pensl and W.J. Choyke. Electrical and optical characterization of SiC.
Physica B: Condensed Matter, 185(1-4):264?283, 1993.
[104] E.O. Kane. Thomas-Fermi Approach to Impure Semiconductor Band Struc-
ture. Physical Review, 131(1):79, 1963.
[105] T.N. Morgan. Broadening of Impurity Bands in Heavily Doped Semiconduc-
tors. Physical Review, 139(1A):A343, 1965.
[106] A.L. Efros. Electron localization in disordered systems (the anderson transi-
tion). Soviet Physics Uspekhi, 21(9):746, 1978.
[107] G. Busch and H. Labhart. U?ber den Mechanismus der elektrischen
Leitfa?higkeit des Siliciumcarbids. Helvetica Physica Acta, 19:463, 1946.
[108] W. Kuz?micz. Ionization of impurities in silicon. Solid-State Electronics,
29(12):1223?1227, 1986.
[109] S. Asada, T. Okuda, T. Kimoto, and J. Suda. Hall scattering factors in p-
type 4H-SiC with various doping concentrations. Applied Physics Express,
9(4):041301, mar 2016.
[110] H. Matsuura, K. Sugiyama, K. Nishikawa, T. Nagata, and N. Fukunaga. Oc-
cupation probability for acceptor in Al-implanted p-type 4HSiC. Journal of
Applied Physics, 94(4):2234?2241, 2003.
182
[111] H. Matsuura. The influence of excited states of deep dopants on majority-
carrier concentration in a wide-bandgap semiconductor. New Journal of
Physics, 4:12?12, mar 2002.
[112] A. Koizumi, J. Suda, and T. Kimoto. Temperature and doping dependencies
of electrical properties in Al-doped 4H-SiC epitaxial layers. Journal of Applied
Physics, 106(1):013716, 2009.
[113] A. Koizumi, N. Iwamoto, S. Onoda, T. Ohshima, T. Kimoto, K. Uchida, and
S. Nozaki. Compensation-dependent carrier transport of Al-doped p-type 4H-
SiC. In Materials Science Forum, volume 679, pages 201?204. Trans Tech
Publ, 2011.
[114] H. Tanaka, S. Asada, T. Kimoto, and J. Suda. Theoretical analysis of Hall
factor and hole mobility in p-type 4H-SiC considering anisotropic valence band
structure. Journal of Applied Physics, 123(24):245704, 2018.
[115] M. Ruff, H. Mitlehner, and R. Helbig. SiC Devices: Physics and Numerical
Simulation. IEEE Transactions on Electron Devices, 41(6):1040?1054, 1994.
[116] L. Kasamakova-Kolaklieva, L. Storasta, I. Ivanov, B. Magnusson, S. Contreras,
C. Consejo, J. Pernot, M. Zielinski, and E. Janze?n. Temperature-Dependent
Hall Effect Measurements in Low Compensated p-Type 4H-SiC. Materials
Science Forum, 457-460:677?680, 2004.
[117] H. Matsuura, M. Komeda, S. Kagamihara, H. Iwata, R. Ishihara,
T. Hatakeyama, T. Watanabe, K. Kojima, T. Shinohe, and K. Arai. De-
pendence of acceptor levels and hole mobility on acceptor density and tem-
perature in Al-doped p-type 4H-SiC epilayers. Journal of Applied Physics,
96(5):2708?2715, 2004.
[118] G. Wagner, W. Leitenberger, K. Irmscher, F. Schmid, M. Laube, and G. Pensl.
Aluminum Incorporation into 4H-SiC Layers during Epitaxial Growth in a
Hot-Wall CVD System. Materials Science Forum, 389-393:207?210, 01 2002.
[119] T. Troffer, M. Schadt, T. Frank, H. Itoh, G. Pensl, J. Heindl, H. P. Strunk, and
M. Maier. Doping of SiC by Implantation of Boron and Aluminum. Physica
Status Solidi (a), 162(1):277?298, 1997.
[120] R. Nipoti, R. Scaburri, A. Halle?n, and A. Parisini. Conventional thermal
annealing for a more efficient p-type doping of Al implanted 4H-SiC. Journal
of Materials Research, 28(1):1722, 2013.
[121] H. Fujihara, J. Suda, and T. Kimoto. Electrical properties of n- and p-type 4H-
SiC formed by ion implantation into high-purity semi-insulating substrates.
Japanese Journal of Applied Physics, 56(7):070306, jun 2017.
183
[122] Y. Negoro, T. Kimoto, H. Matsunami, F. Schmid, and G. Pensl. Electrical
activation of high-concentration aluminum implanted in 4H-SiC. Journal of
Applied Physics, 96(9):4916?4922, 2004.
[123] D.M. Caughey and R.E. Thomas. Carrier Mobilities in Silicon Empirically
Related to Doping and Field. Proceedings of the IEEE, 55(12):2192?2193,
1967.
[124] D. Stefanakis and K. Zekentes. TCAD models of the temperature and doping
dependence of the bandgap and low field carrier mobility in 4H-SiC. Micro-
electronic Engineering, 116:65?71, 2014.
[125] M. Roschke and F. Schwierz. Electron mobility models for 4H, 6H, and 3C
SiC [MESFETs]. IEEE Transactions on Electron Devices, 48(7):1442?1447,
2001.
[126] W.J. Schaffer, G.H. Negley, K.G. Irvine, and J.W. Palmour. Conductivity
Anisotropy in Epitaxial 6H and 4H SiC. MRS Proceedings, 339:595, 1994.
[127] E. McIrvine. Phenomenology of impurity conduction in semiconductors. Jour-
nal of Physics and Chemistry of Solids, 15:356?358, 1960.
[128] H. Matsuura, A. Takeshita, T. Imamura, K. Takano, K. Okuda, A. Hidaka,
S. Ji, K. Eto, K. Kojima, T. Kato, S. Yoshida, and H. Okumura. Depen-
dence of conduction mechanisms in heavily Al-doped 4H-SiC epilayers on Al
concentration. Applied Physics Express, 11(10):101302, 2018.
[129] R.M. Hill. Variable-Range Hopping. Physica Status Solidi (a), 34(2):601?613,
1976.
[130] M.K. Linnarsson, M.S. Janson, U. Zimmermann, B.G. Svensson, P.O.A?. Pers-
son, L. Hultman, J. Wong-Leung, S. Karlsson, A. Scho?ner, H. Bleichner, and
E. Olsson. Solubility limit and precipitate formation in Al-doped 4H-SiC epi-
taxial material. Applied Physics Letters, 79(13):2016?2018, 2001.
[131] M.K. Linnarsson, P.O.A?. Persson, H. Bleichner, M.S. Janson, U. Zimmer-
mann, H. Andersson, S. Karlsson, R. Yakimova, L. Hultman, and B.G. Svens-
son. Precipitate formation in heavily Al-doped 4H-SiC layers. In Materials
Science Forum, volume 353, pages 583?586. Trans Tech Publications Ltd.,
Zurich-Uetikon, Switzerland, 2001.
[132] Kazutaka Tomizawa. Numerical Simulation of Submicron Semiconductor De-
vices. Artech House, 1993.
[133] Hongchin Lin and Neil Goldsman. An efficient solution of the boltzmann
transport equation which includes the pauli exclusion principle. Solid-state
electronics, 34(10):1035?1048, 1991.
184
[134] Hongchin Lin, Neil Goldsman, and ID Mayergoyz. An efficient determinis-
tic solution of the space-dependent boltzmann transport equation for silicon.
Solid-state electronics, 35(1):33?42, 1992.
[135] V.V. Afanasev, M. Bassler, G. Pensl, M.J. Schulz, and E. Stein von Kamienski.
Band offsets and electronic structure of sic/sio2 interfaces. Journal of Applied
Physics, 79(6):3108?3114, 1996.
[136] Jenny Nelson. The Physics of Solar Cells. World Scientific Publishing Com-
pany, 2003.
[137] David S Sukhdeo, Hai Lin, Donguk Nam, Ze Yuan, Boris M Vulovic, Suyog
Gupta, James S Harris, Birendra Dutt, and Krishna C Saraswat. Approaches
for a viable germanium laser: tensile strain, gesn alloys, and n-type doping.
In 2013 Optical Interconnects Conference, pages 112?113. IEEE, 2013.
[138] S Gupta. Germanium-Tin (GeSn) Technology. PhD thesis, Stanford Univer-
sity, California, 2013.
[139] Suyog Gupta, Blanka Magyari-Ko?pe, Yoshio Nishi, and Krishna C Saraswat.
Achieving direct band gap in germanium through integration of sn alloying
and external strain. Journal of Applied Physics, 113(7):073707, 2013.
[140] Pairot Moontragoon, Zoran Ikonic?, and Paul Harrison. Band structure calcu-
lations of si?ge?sn alloys: achieving direct band gap materials. Semiconductor
science and technology, 22(7):742, 2007.
[141] Suyog Gupta, Blanka Magyari-Ko?pe, Yoshio Nishi, and Krishna C Saraswat.
Band structure and ballistic electron transport simulations in gesn alloys. SIS-
PAD, 2012.
[142] Philipp Haas, Fabien Tran, and Peter Blaha. Calculation of the lattice con-
stant of solids with semilocal functionals. Physical Review B, 79(8):085104,
2009.
[143] Andrea Dal Corso. A pseudopotential plane waves program (pwscf) and some
case studies. In Quantum-Mechanical Ab-initio Calculation of the Properties
of Crystalline Materials, pages 155?178. Springer, 1996.
[144] Lianhua He, Fang Liu, Geoffroy Hautier, Micael JT Oliveira, Miguel AL Mar-
ques, Fernando D Vila, JJ Rehr, G-M Rignanese, and Aihui Zhou. Accuracy
of generalized gradient approximation functionals for density-functional per-
turbation theory calculations. Physical Review B, 89(6):064305, 2014.
[145] Nandor L Bala?zs. Formation of stable molecules within the statistical theory
of atoms. Physical Review, 156(1):42, 1967.
[146] Nikita S Shah. Ab Initio Molecular Dynamics Calculation Of Vibrational
Properties Of Refractory Carbides. PhD thesis, Bhavnagar, 2013.
185
[147] Pierre Hohenberg and Walter Kohn. Inhomogeneous electron gas. Physical
review, 136(3B):B864, 1964.
[148] Claudio Lombardi, Stefano Manzini, Antonio Saporito, and Massimo Vanzi.
A physically based mobility model for numerical simulation of nonplanar de-
vices. IEEE Transactions on Computer-Aided Design of Integrated Circuits
and Systems, 7(11):1164?1171, 1988.
[149] Gary Pennington, Siddharth Potbhare, Neil Goldsman, James M McGarrity,
and Aivars Lelis. Impact of surface steps on the roughness mobility in 4h-
sic. In Semiconductor Device Research Symposium, 2005 International, pages
143?144. IEEE, 2005.
[150] Tsunenobu Kimoto, Akira Itoh, Hiroyuki Matsunami, and Tetsuyuki Okano.
Step bunching mechanism in chemical vapor deposition of 6H?and 4H?SiC
{0001}. Journal of Applied Physics, 81(8):3494?3500, 1997.
[151] N.S. Saks, A.K. Agarwal, S.H. Ryu, and J.W. Palmour. Low-dose aluminum
and boron implants in 4H and 6H silicon carbide. Journal of Applied Physics,
90(6):2796?2805, 2001.
[152] E.M. Handy, M.V. Rao, O.W. Holland, P.H. Chi, K.A. Jones, M.A. Derenge,
R.D. Vispute, and T. Venkatesan. Al, B, and Ga ion-implantation doping of
SiC. Journal of Electronic Materials, 29(11):1340?1345, 2000.
[153] M. Ikeda, H. Matsunami, and T. Tanaka. Site effect on the impurity levels in
4H, 6H, and 15R SiC. Physical Review B, 22(6):2842, 1980.
[154] S.S. Li and W.R. Thurber. The dopant density and temperature dependence
of electron mobility and resistivity in n-type silicon. Solid-State Electronics,
20(7):609?616, 1977.
[155] S.R. Smith, A.O. Evwaraye, W.C. Mitchel, and M.A. Capano. Shallow accep-
tor levels in 4h- and 6h-sic. Journal of Electronic Materials, 28(3):190?195,
Mar 1999.
[156] I.G. Ivanov, B. Magnusson, and E. Janze?n. Analysis of the sharp donor-
acceptor pair luminescence in 4H-SiC doped with nitrogen and aluminum.
Physical Review B, 67:165211, Apr 2003.
[157] H. Matsuura, K. Aso, S. Kagamihara, H. Iwata, T. Ishida, and K. Nishikawa.
Decrease in Al acceptor density in Al-doped 4H-SiC by irradiation with 4.6
MeV electrons. Applied Physics Letters, 83(24):4981?4983, 2003.
[158] P. Achatz, J. Pernot, C. Marcenat, J. Kacmarcik, G. Ferro, and E. Bustar-
ret. Doping-induced metal-insulator transition in aluminum-doped 4H silicon
carbide. Applied Physics Letters, 92(7):072103, 2008.
186
[159] J. Pernot, S. Contreras, and J. Camassel. Electrical transport properties of
aluminum-implanted 4HSiC. Journal of Applied Physics, 98(2):023706, 2005.
[160] V. Heera, K.N. Madhusoodanan, W. Skorupa, C. Dubois, and H. Romanus. A
comparative study of the electrical properties of heavily Al implanted, single
crystalline and nanocrystalline SiC. Journal of Applied Physics, 99(12):123716,
2006.
[161] N.I. Kuznetsov and A.S. Zubrilov. Deep centers and electroluminescence in
4H?SiC diodes with a p-type base region. Materials Science and Engineering:
B, 29(1-3):181?184, 1995.
[162] J.-M. Bluet, J. Pernot, J. Camassel, S. Contreras, J.L. Robert, J.F. Michaud,
and T. Billon. Activation of aluminum implanted at high doses in 4H?SiC.
Journal of Applied Physics, 88(4):1971?1977, 2000.
[163] W. Kaindl, M. Lades, N. Kaminski, E. Niemann, and G. Wachutka. Experi-
mental characterization and numerical simulation of the electrical properties
of nitrogen, aluminum, and boron in 4H/6H-SiC. Journal of Electronic Ma-
terials, 28(3):154?160, Mar 1999.
[164] H. Itoh, T. Troffer, and G. Pensl. Coimplantation Effects on the Electrical
Properties of Boron and Aluminium Acceptors in 4H-SiC. In Silicon Carbide,
III-Nitrides and Related Materials, volume 264 of Materials Science Forum,
pages 685?688. Trans Tech Publications Ltd, 12 1997.
[165] J. Weisse, M. Hauck, T. Sledziewski, M. Tschiesche, M. Krieger, A.J. Bauer,
H. Mitlehner, L. Frey, and T. Erlbacher. Analysis of compensation effects
in aluminum-implanted 4H-SiC devices. In Materials Science Forum, volume
924, pages 184?187. Trans Tech Publ, 2018.
[166] W.J. Schaffer, H.S. Kong, G.H. Negley, and J.W. Palmour. Hall effect and CV
measurements on epitaxial 6H- and 4H-SiC. Institute of Physics Conference
Series, 137:155?159, 1994.
[167] N.S. Saks, A.V. Suvorov, and D.C. Capell. High temperature high-dose im-
plantation of aluminum in 4H-SiC. Applied Physics Letters, 84(25):5195?5197,
2004.
[168] G.A. Lomakina, Yu.A. Vodakov, E.N. Mokhov, V.G. Oding, and Kholuyan.
G.F. Comparative Study of the Electrical Properties of Three Polytypes of
Silicon Carbide. Soviet Physics Solid State, 12(10):2356?2359, 1971.
[169] I.R. Arvinte. Investigation of dopant incorporation in silicon carbide epilayers
grown by chemical vapor deposition. PhD thesis, Universite Co?te d?Azur, 11
2016.
187
[170] F. Schmid, M. Krieger, M. Laube, G. Pensl, and G. Wagner. Hall Scatter-
ing Factor for Electrons and Holes in SiC, pages 517?536. Springer Berlin
Heidelberg, Berlin, Heidelberg, 2004.
[171] S. Ji, K. Kojima, Y. Ishida, H. Tsuchida, S. Yoshida, and H. Okumura. Low
Resistivity, Thick Heavily Al-Doped 4H-SiC Epilayers Grown by Hot-Wall
Chemical Vapor Deposition. In Silicon Carbide and Related Materials 2012,
volume 740 of Materials Science Forum, pages 181?184. Trans Tech Publica-
tions Ltd, 3 2013.
[172] S. Contreras, M. Zielinski, C. Konczewicz, L.and Blanc, S. Juillaguet,
R. Mu?ller, U. Ku?necke, P. Wellmann, and J. Camassel. Results of SIMS,
LTPL and temperature-dependent Hall effect measurements performed on Al-
doped ?-SiC substrates grown by the M-PVT method. Materials Science
Forum, 527-529:633?636, 2006.
[173] Y. Tanaka, N. Kobayashi, H. Okumura, R. Suzuki, T. Ohdaira, M. Hasegawa,
M. Ogura, S. Yoshida, and H. Tanoue. Electrical and Structural Properties
of Al and B Implanted 4H-SiC. Materials Science Forum, 338-342:909?912,
2000.
[174] M. Obernhofer, M. Krieger, F. Schmid, H.B. Weber, G. Pensl, and A. Scho?ner.
Electrical and Structural Properties of Al-Implanted and Annealed 4H-SiC. In
Silicon Carbide and Related Materials 2006, volume 556 of Materials Science
Forum, pages 343?346. Trans Tech Publications Ltd, 8 2007.
[175] P. Kwasnicki. Evaluation of doping in 4H-SiC by optical spectroscopies. PhD
thesis, Universite? Montpellier II - Sciences et Techniques du Languedoc, 12
2014.
[176] J. Pernot, J. Camassel, S. Contreras, J.-L. Robert, J.-M. Bluet, J.F. Michaud,
and T. Billon. Control of Al-implantation doping in 4H?SiC. Materials Science
and Engineering: B, 80(1):362 ? 365, 2001.
[177] J.-M. Bluet, J. Pernot, T. Billon, S. Contreras, J.F. Michaud, J.-L. Robert,
and J. Camassel. Al and Al/C High Dose Implantation in 4H-SiC. In Silicon
Carbide and Related Materials - 1999, volume 338 of Materials Science Forum,
pages 885?888. Trans Tech Publications Ltd, 5 2000.
[178] G. Liaugaudas, D. Dargis, P. Kwasnicki, H. Peyre, R. Arvinte, S. Juillaguet,
M. Zielinski, and K. Jarasiunas. Optical Characterization of p-Type 4H-SiC
Epilayers. Materials Science Forum, 821-823:249?252, 06 2015.
[179] S. Asada, J. Suda, and T. Kimoto. Analytical formula for temperature depen-
dence of resistivity in p-type 4H-SiC with wide-range doping concentrations.
Japanese Journal of Applied Physics, 57(8):088002, jun 2018.
188
[180] S. Ji, K. Kojima, Y. Ishida, S. Saito, T. Kato, H. Tsuchida, S. Yoshida, and
H. Okumura. The growth of low resistivity, heavily Al-doped 4H?SiC thick
epilayers by hot-wall chemical vapor deposition. Journal of Crystal Growth,
380:85?92, 2013.
[181] M. Rambach, A.J. Bauer, and H. Ryssel. Electrical and topographical charac-
terization of aluminum implanted layers in 4H silicon carbide. Physica Status
Solidi (b), 245(7):1315?1326, 2008.
[182] T. Shirai, A. Danno, K.and Seki, H. Sakamoto, and T. Bessho. Solution
Growth of p-Type 4H-SiC Bulk Crystals with Low Resistivity. In Silicon
Carbide and Related Materials 2013, volume 778 of Materials Science Forum,
pages 75?78. Trans Tech Publications Ltd, 5 2014.
[183] T. Kimoto, H. Yano, S. Tamura, N. Miyamoto, K. Fujihira, Y. Negoro, and
H. Matsunami. Recent Progress in SiC Epitaxial Growth and Device Process-
ing Technology. Materials Science Forum, 353-356:543?548, 2001.
[184] T. Miyazawa, K. Nakayama, A. Tanaka, K. Asano, S. Ji, K. Kojima, Y. Ishida,
and H. Tsuchida. Epitaxial growth and characterization of thick multi-layer
4H-SiC for very high-voltage insulated gate bipolar transistors. Journal of
Applied Physics, 118(8):085702, 2015.
[185] K. Tone and J.H. Zhao. A comparative study of C plus Al coimplantation and
Al implantation in 4H and 6H-SiC. IEEE Transactions on Electron Devices,
46:612 ? 619, 04 1999.
[186] C. Persson, A.F. da Silva, and B. Johansson. Metal-nonmetal transition in
p-type SiC polytypes. Physical Review B, 63(20):205119, 2001.
189