Electronic Structure of Heterogeneous Materials [623783]

Electronic Structure of Heterogeneous Materials
Application to Optical Properties
Nat´alia Leit ˜ao Marques Morais
(Licenciada)
Dissertac ¸ ˜ao para obter o grau de Mestre em
Engenharia Fisica Tecnol ´ogica
J ´ uri
Presidente:
Orientador: Jos´e Lu´ıs Rodrigues J ´ulio Martins
Vogal:
September 2014

Imagination is more important than knowledge. For knowledge is limited to all we now know and understand, while
imagination embraces the entire world, and all there ever will be to know and understand.
Albert Einstein

Acknowledgments
I would like to thank the Instituto de Engenharia de Sistemas e Computadores, Microsistemas & Nanotec-
nologias (INESC-MN) of Lisbon, that very kindly hosted me in as I predict to be the beginning of my career as a
researcher in Physics, subject that I’ve always been passionate about since at a very young age.
I would like to thank my coordinator Jos ´e Luis Martins for guiding me on this project that concluded this
Mestrado em Engenharia F ´ısica Tecnol ´ogica (MEFT) 5 years Master program and providing me an opportunity
to perform a high-quality research in the domain of Condensed Matter Physics and also to his research fellow
Carlos Reis to provide some other help I needed.
iii

Abstract
The objective of this work is to find a description of the group IV elements Silicon, Carbon and Germanium, in
bulk, calculating its electronic structure and optical properties. To calculate the band structure we will use pseu-
dopotentials. We will use an Empirical Pseudopotential Method to find the better fitting to the pseudopotentials to
the experimental known data about the band structure to each of this elements. We start by fitting the pseudopo-
tentials to an ab initio pseudopotential generator [1], to find the first acceptable set of pseudopotentials. From that
we further adjust the potentials to the experiment, and find the better fit. After that the optical properties of the bulk
Si, C and Ge are calculated. The purpose is to generate a pseudopotential to each of this elements that simulates
correctly the properties and can be transferable to supercells of Si-Ge-C.
Keywords
Nanotechnologies, Simulation of Materials, Condensed Matter Physics, Solid State Applications, Pseudopo-
tentials
v

Resumo
O objectivo teste projecto ´e encontrar uma descric ¸ ˜ao para um s ´olido cristalino de elementos de grupo IV Sil ´ıcio,
Carbono e Germ ˆanio, e calcular a estrutura electr ´onica e propriedades ´opticas. Para calcular a estrutura de ban-
das v ˜ao se usar pseudopotencials. Vai ser usado o M ´etodo de Pseudopotencial Emp ´ırico (EPM) para encontrar
o melhor ajuste dos pseudopotenciais aos dados experimentais conhecidos da estrutura de bandas de cada um
destes elementos. Antes disso, comec ¸a-se por ajustar os pseudopotenciais a um gerador de pseudopotenci-
aisab initio , para que se encontre uma regi ˜ao aceit ´avel de pseudopotencias. A partir da ´ı, melhora-se o ajuste,
ajustando os potenciais `a experi ˆencia. Depois as propriedades ´opticas dos cristais de Si, C e Ge podem ser cal-
culadas. O objectivo ´e de gerar um pseudopotencial para cada um destes elementos que simule correctamente
as propriedades, que depois seja transfer ´ıvel para super-redes de Si-Ge-C.
Palavras Chave
Nanotecnologias, Simulac ¸ ˜ao de Materiais, F ´ısica da Mat ´eria Condensada, Aplicac ¸ ˜oes de F ´ısica do Estado
S´olido, Pseudopotenciais
vii

Contents
1 Introduction xvii
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii
1.2 State of The Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii
1.2.1 Si, Ge and C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix
1.2.2 Group IV Semiconductor Compounds and Alloys . . . . . . . . . . . . . . . . . . . . . . . . . xx
1.2.2.A Si and Ge heterostructures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx
1.2.2.B Si and Ge heterostructures with C impurities . . . . . . . . . . . . . . . . . . . . . . xxiv
1.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxvi
2 Theoretical Introduction xxvii
2.1 The many-body problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxviii
2.2 The Adiabatic Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxix
2.3 Separable Schr ¨odinger equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxix
2.4 The Hartree-Fock Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxx
2.5 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxi
2.6 Hohenberg-Kohn Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxii
2.7 Kohn-Sham Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxiii
2.8 Pseudopotentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxv
2.9 Ab initio pseudo-potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxv
2.10 Empirical Pseudopotential methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxvii
2.11 Non-local and Spin-Orbit Pseudopotentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxviii
2.12 Screening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xlii
2.12.1 Definitions of the dielectric function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xlii
2.12.2 Screening in a metal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xliii
2.13 Optical properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xlv
2.14 Imaginary part of the dielectric function 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xlvi
3 Results xlviii
3.0.1 Simple test – free electron . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xlix
3.1 Silicon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . l
3.1.1 Description only with a local pseudopotential . . . . . . . . . . . . . . . . . . . . . . . . . . . l
3.1.2 Description with non-local pseudopotential . . . . . . . . . . . . . . . . . . . . . . . . . . . . liv
3.2 Carbon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lix
3.3 Germanium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxiv
ix

4 Conclusions and Future Work lxxi
Bibliography lxxv
5 Appendix 1
.1 Spin-Orbit projectors for the pseudopotential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
.2 Clebsch-Gordon coefficients for mixing states lands=1
2. . . . . . . . . . . . . . . . . . . . . . . . 2
.3 Matrix elements of the momentum matrix operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
.3.1 Orthogonality of the basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
.3.2 Local Pseudopotencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
.3.3 Nonlocal Potencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
.3.4 Spin-Orbit contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
x

List of Figures
1.1 The figure shows a diamond crystal structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix
1.2 The energy-band structure of a) Si and b) Ge are calculated with a tight-binding model. The top of
the valence bands is set at zero energy. [2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix
1.3 The energy-band structure of a) Si and b) Ge, calculated with a tight-binding model. The top of the
valence bands is set at zero energy. [3] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx
1.4 This models are consistent with the proposed structure of SiGe. Top: one-eighth of the ordered unit
cell Bottom: [101] projection through the unit cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi
1.5 The unit cell of the (Ge) 2/(Si) 2superlattice and its band-structure are shown . . . . . . . . . . . . . . xxi
1.6 The figure shows a) electronic band structure of the (Ge) 6/(Si) 4supperlattice on Si [001] substrate
and band structure of the free-standing b) (Ge) 5/(Si) 5c) (Ge) 4/(Si) 6superlattice along lines of high
symmetry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii
1.7 Transition energies of various n+m= 10 superlattices s as a function of lateral strain in the Si layer
are represented . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii
1.8 Band structures E(~k)for the principal symmetry directions of the diamond lattice for (a) Si 0:2Ge0:8
and (b) Si 0:74Ge0:26where calculated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii
1.9 The figures show the predicted single-impurity defect levels of a) T2and b)A1symmetry as a
function of the composition x. Shown are also the conduction and valence band edges . . . . . . . . xxiii
1.10 Band structures of the 12magic sequence grown on Si 0:4Ge0:6where calculated . . . . . . . . . . xxiv
1.11 The figure shows the omparison of between Si 6Ge4superlattice and the magic sequence of the
direct absorption spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiv
1.12 Direct dipole-allowed band-gaps where measured . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiv
1.13 Valence-band offsets for compressive strained Si 1xGex, and Si 1xyGexCy(x=10%, 20%, and
30%) and tensile strained Si 1yCyand Si 1xyGexCy(y=1%, 2%, and 3%) are plotted as a function
of the effective lattice mismatch (expressed in “effective” Ge or C concentrations, respectively) . . . xxv
1.14 Conduction-band offsets for compressive strained Si 1xGex, and Si 1xyGexCy(x=10%, 20%,
and 30%) and tensile strained Si 1yCyand Si 1xyGexCy(y=1%, 2%, and 3%) are plotted as a
function of the effective lattice mismatch (expressed in “effective” Ge or C concentrations, respectively)xxv
1.15 The band-gap narrowing for ternary Si 1xyGexCyalloys, strained on Si(001) is represented as a
function of lattice mismatch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxvi
2.1 The figure shows the Venn’s diagram corresponding a potential in the space of all potentials V(~ r)to
a ground state electron density in the space of the ground state electron densities GS(~ r). . . . . . xxxi
2.2 The figure shows a schematic plot of a pseudopotential in reciprocal space with the G’s that corre-
spond toG2= 3;8;11withG2in units of (2=a)2[4] . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxviii
xi

3.1 The figure shows a) the free electron bands in the fcc lattice and b) the Brillouin zone for the fcc latticexlix
3.2 The graphic shows the functions that compose the local pseudopotential and the pseudopotential
itself, unscreened and screened . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . l
3.3 The figure shows the LDA band structure of bulk Silicon calculated with the program of reference
[5]. Experimental values are indicated by the double arrows. . . . . . . . . . . . . . . . . . . . . . . . lii
3.4 The Figure shows the a) calculated density of states (blue line), the photo emission spectroscopy
and inverse photo emission data obtained from reference [6] (yellow line) and b) the calculated joint
density of states with the parameters from Table 3.1, for Silicon . . . . . . . . . . . . . . . . . . . . . lii
3.5 Real part 1and imaginary part 2of the dielectric function of Silicon are calculated using the pa-
rameters on Table 3.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . liii
3.6 Comparing the calculated dielectric function of Silicon, using the parameters on Table 3.1, with
experimental results in reference [7] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . liv
3.7 The figure shows the calculated reflectance for Silicon using Table 3.1 (blue line) and the experi-
mentally obtainced from reference [7] (yellow line) . . . . . . . . . . . . . . . . . . . . . . . . . . . . liv
3.8 The figure shows the fit of the expression (3.9) (line) to the ppseudopotential of Silicon generated
with the program in reference [1] (dots) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lv
3.9 The ballpark figure to fit the local part of the Silicon pseudopotential using function (3.21) is shown.
The green point is the result of the fit, in Table 3.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lvi
3.10 The ballpark figure to fit the local part of the Silicon pseudopotential using function (3.21) is shown.
The green point is the result of the fit, in Table 3.4, the red points are the pairs of values ( Ra,kTF)
for which the eigenvalue of energy is the one of the Ep, the green line is adjusted to these points,
the orange points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy is the one of
theEp+ 0:5eVand the yellow points are the pairs of values ( Ra,kTF) for which the eigenvalue of
energy is the one of the Ep0:5eV.Ep=4:16eVis obtained from reference [1] for Silicon . . . . lvi
3.11 A fit that was made (line), of (3.20) to the s“projector” of Silicon, generated with [1] (points) . . . . . lvii
3.12 The ballpark figure to fit the non-local part of the pseudopotential of Silicon using function (3.23) is
shown. The green point is the result of the fit, in Table 3.6, the red points are the pairs of values
(Rb,B) for which the eigenvalue of energy is the one of the Es, the yellow parabola is adjusted to
these points, the orange points are the pairs of values ( Rb,B) for which the eigenvalue of energy
is the one of the Es+ 0:5eVand the yellow points are the pairs of values ( Rb,B) for which the
eigenvalue of energy is the one of the Es0:5eV.Es=10:83eVis obtained from reference [1] . . lviii
3.13 Band structure of Silicon was a) calculated using LDA, from [5] and b) calculated using the parame-
ters on Table 3.9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lviii
3.14 It is represented a) the calculated density if states (blue line),the photo emission spectroscopy and
inverse photo emission data obtained from reference [6] (yellow line), b) the calculated joint density
of states, c) calculated dielectric function, d) calculated (blue line) and experimental (yellow line,
[7])1, e) calculated (blue line) and experimental (yellow line, [7]) 2f),g) calculated (blue line) and
experimental (yellow line, [7]) reflectance for Silicon with the local pseudopotential of equation (3.7)
and non-local projector for the pseudopotential of equation (3.19) with the parameters written in
Table 3.9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lx
3.15 The figure shows the fit of the expression (3.9) to the ppseudopotential of Carbon generated with
the program in reference [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxi
xii

3.16 The ballpark figure to fit the local part of the pseudopotential using function (3.21) is shown. The
green point is the result of the fit, in Table 3.11, the red points are the pairs of values ( Ra,kTF) for
which the eigenvalue of energy is the one of the Ep, the green line is adjusted to these points, the
orange points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy is the one of the
Ep+0:5eVand the yellow points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy
is the one of the Ep0:5eV.Ep=5:41eVis obtained from reference [1] for Carbon . . . . . . . . lxi
3.17 The fit of the expression (3.20) to the s“projector” of Carbon, generated with the program in refer-
ence [1] is represented . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxi
3.18 The ballpark figure to fit the non-local part of the pseudopotential of Carbon using function (3.23) is
shown. The green point is the result of the fit, in Table 3.13, the red points are the pairs of values
(Rb,B) for which the eigenvalue of energy is the one of the Es, the red parabola is adjusted to these
points, the orange points are the pairs of values ( Rb,B) for which the eigenvalue of energy is the one
of theEs+ 0:5eVand the yellow points are the pairs of values ( Rb,B) for which the eigenvalue of
energy is the one of the Es0:5eV.Es=13:63eVis obtained from reference [1] . . . . . . . . . . lxii
3.19 Band structure of diamond a) form reference [5] and b) calculated with the parameters in Table 3.16 lxiii
3.20 It is represented the a) DOS of Carbon, calculated here (blue line), the photo emission spectroscopy
data from reference [8] (yellow line) divided by a factor of 20, b) the calculated joint density of states,
the c) real part (blue line) and imaginary part (purple line) of the dielectric function, d) the comparison
between the calculated (blue) and experimentally obtained (yellow, [9]) 1, e) comparison between
the calculated (blue) and experimental (yellow, [9] 2and f) the calculated (blue) and experimentally
obtained (yellow, [9]), divided by a factor of 100, reflectance. . . . . . . . . . . . . . . . . . . . . . . . lxiv
3.21 The figure shows the fit of the expression (3.9) to the ppseudopotential of Germanium, generated
with the program in reference [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxv
3.22 The ballpark figure to fit the local part of the pseudopotential using function (3.21) is shown. The
green point is the result of the fit, in Table 3.18, the red points are the pairs of values ( Ra,kTF) for
which the eigenvalue of energy is the one of the Ep, the yellow line is adjusted to these points, the
orange points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy is the one of the
Ep+0:5eVand the yellow points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy
is the one of the Ep0:5eV.Ep=4:05eVis obtained from reference [1] for Germanium . . . . . . lxvi
3.23 The fit of the expression (3.20) to the s“projector” of Germanium generated with the program in
reference [1] is shown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxvi
3.24 The ballpark figure to fit the non-local part of the pseudopotential of Germanium using function
(3.23) is shown. The green point is the result of the fit, in Table 3.13, the red points are the pairs
of values (Rb,B) for which the eigenvalue of energy is the one of the Es, the yellow parabola is
adjusted to these points, the orange points are the pairs of values ( Rb,B) for which the eigenvalue
of energy is the one of the Es+ 0:5eVand the yellow points are the pairs of values ( Rb,B) for which
the eigenvalue of energy is the one of the Es0:5eV.Es=11:92eVis obtained from reference [1] lxvii
3.25 Calculated band structure of Germanium, without the spin-orbit splitting, with the parameters on
Table 3.23 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxvii
3.26 The pseudopotentials a) Vlocal;screen (kx), b)
Vnonlocal (kx)and c)
Vspinorbit (kx)of Germanium are
represented graphically, using the parameters of Table 3.27 . . . . . . . . . . . . . . . . . . . . . . . lxviii
xiii

3.27 The pseudopotentials Vlocal;screen (kx),Vlocal;screen (kx)+
Vnonlocal (kx)andVlocal;screen (kx)+
Vnonlocal (kx)+

Vspinorbit (kx)of Germanium are calculated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxix
3.28 Calculated band structure of Germanium, a) using reference [5] and with the values of some impor-
tant transitions in eV, b) with the program, using the parameters of Table 3.27 . . . . . . . . . . . . . lxix
3.29 It is represented a) the calculated density if states (blue line),the photo emission spectroscopy and
inverse photo emission data obtained from reference [6] (yellow line) divided by a factor of 3, b)
the calculated joint density of states, c) calculated dielectric function, real part (blue) and imaginary
(purple), d) calculated (blue line) and experimental (yellow line, [7]) 1, e) calculated (blue line)
and experimental (yellow line, [7]) 2f),g) calculated (blue line) and experimental (yellow line, [7])
reflectance for Germanium with the local pseudopotential of equation (3.7), the non-local projector
for the pseudopotential of equation (3.19) with the parameters written in Table 3.9 and the spin-orbit
projectors using equations (2.74) and (3.25-3.27) with lmax= 2 . . . . . . . . . . . . . . . . . . . . . lxx
4.1 The figure shows the a) band structure and the b) density of states, calculated for Silicon using the
pseudopotentials obtained here and using a DFT -MGGA calculation . . . . . . . . . . . . . . . . . . lxxii
4.2 The figure shows the a) band structure and the b) density of states, calculated for Carbon using the
pseudopotentials obtained here and using a DFT -MGGA calculation . . . . . . . . . . . . . . . . . . lxxiii
4.3 The figure shows the a) band structure and the b) density of states, calculated for Germanium using
the pseudopotentials obtained here and using a DFT -MGGA calculation . . . . . . . . . . . . . . . . lxxiii
4.4 The figure shows the a) band structure and the b) density of states, calculated for Germanium using
the pseudopotentials obtained here and using a DFT -MGGA calculation . . . . . . . . . . . . . . . . lxxiii
xiv

List of Tables
1.1 Calculated values for the gaps (in eV) are shown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix
3.1 The table shows previously existent parameters for Silicon, that we first check in this research. . . . li
3.2 Resulting band structure subjected to the potential (3.9) of Silicon, using the parameters from Table
3.1 multiplied by a factor fwith0f1. The the opening of the gap is clearly shown . . . . . . . . li
3.3 Experimental and calculated in the current work transitions of Silicon in eVare calculated with the
parameters from Table 3.1 in equation (3.9) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . li
3.4 Results to the fit of equation (3.9) to the ppseudopotential of Silicon, generated with the program in
reference [1], using the default weight function chosen by M ATHEMATICA . . . . . . . . . . . . . . . . lv
3.5 The results to the fitting, using (3.22), to the ppseudopotential of Silicon generated with the program
in reference [1] are shown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lvii
3.6 The obtained parameters of the fitting of (3.20) to the s“projector” generated with [1] are shown . . . lvii
3.7 The obtained final parameters, used to calculate the important energetic transitions of Silicon are
shown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lvii
3.8 Important energy transitions of Silicon where calculated with the parameters in Table 3.7 . . . . . . . lviii
3.9 The final parameters for the pseudopotential of Silicon where obtained after adjusting to the experiment lix
3.10 Experimental and calculated in the current work transitions of Silicon in eVare calculated with the
parameters on Table 3.9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lix
3.11 The table shows the results to the fit of equation (3.9) to the ppseudopotential of Carbon generated
with the program in reference [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lix
3.12 The results to the fitting, using (3.22), to the ppseudopotential of Carbon, generated with the pro-
gram in reference [1] are shown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lix
3.13 Results to the fit of equation (3.20) to the s“projector” of Carbon generated with the program in
reference [1] are shown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxii
3.14 These are the obtained final parameters, used to calculate the important energetic transitions of
Carbon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxii
3.15 This important energy transitions of Carbon were calculated with the parameters in Table 3.14 . . . lxiii
3.16 The discovered final parameters for Carbon, after the adjustment to the experiment are shown here lxiii
3.17 Experimental and calculated in the current work transitions of Diamond in eV, calculated with the
parameters on Table 3.16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxiii
3.18 The table shows the results to the fit of equation (3.9) to the ppseudopotential of Germanium
generated with the program in reference [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxv
3.19 The table shows the results to the fitting to the ppseudopotential of Germanium generated with the
program in reference [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxv
xv

3.20 Results to the fit of equation (3.20) to the s“projector” of Germanium generated with the program in
reference [1] are shown here . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxvi
3.21 This obtained parameters are used to calculate the some transitions of Germanium without the
spin-orbit part . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxvi
3.22 Important energetic transitions of Germanium, calculated without using the spin-orbit splitting, with
the parameters on Table 3.21 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxvii
3.23 Discovered parameters, after the better “adjustment” to the experiment, used to calculate the band
structure of Germanium, without the spin-orbit splitting. The “adjustment” is as close as possible
because without the spin-orbit contribution we cannot fully describe Germanium. . . . . . . . . . . . lxvii
3.24 Important transitions calculated for a bulk Germanium without the spin-orbit splitting, with the pa-
rameters on Table 3.23 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxviii
3.25 This set of parameters are used to calculate important transitions of Germanium, with the spin-orbit
splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxviii
3.26 Important transitions of Germanium, with the spin-orbit splitting, are calculated with the parameters
in Table 3.25 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxviii
3.27 The set of parameters used to calculate the band structure of Germanium, with the spin-orbit splitting
are shown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxix
3.28 Important optical transitions of Germanium where calculated using the parameters of Table 3.27 . . lxix
xvi

1
Introduction
Contents
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii
1.2 State of The Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii
1.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxvi
xvii

1.1 Motivation
In micro and nanotechnologies of today, there is an high interest in semiconductor materials that have a direct
gap that can be grown in Silicon, since it is the material that is widely used in integrated circuits. We are interested
in study materials based in super-lattices of Si-Ge with C impurities, and to simulate its optical properties.
We are interested in simulating cells with many atoms. There are already quite a few methods to do so, each
one with its positive and negative points. The vast majority of electronic structure calculations are done using
the Kohn-Sham equations [10] with the local-density approximation (LDA) [10, 11] for the exchange-correlation
energy and potential (as we will describe later on). However this calculations can lead to results with very bad
agreement with experiment, namely the band gap of semiconductors and insulators. For example, in Si, with LDA,
the band-gap is predicted to be one half of its value, while in Ge the band gap is very small or even disappears.
There are new ways of calculating the exchange and correlations functionals, like the Tran-Blaha model [12].
This gives very much more improved band gaps for a variety of insulators and semiconductors.
The GW approximation (GWA) is used to calculate the self-energy of a many-body system of electrons [13]. The
approximation to be made is that the expansion of the self-enerfy "in therms of the single particle Green’s function
Gand the screened Coulomb interaction Wcan be approximated to the lowest order term. The GW method can
yield very accurate band gaps, but leads to very heavy calculations, particularly if done self-consistently.
The Empirical Pseudopotential Method (EPM) relies on the experimental results to fit a set of parameters used
to describe the potentials that act on the electrons. This has the advantage that, if the programming is efficient,
the calculations can be made very quickly, and therefore a very large number of atoms can be included. There are
although “dangers” in this method, since it is very tempting to use a model with a big number of parameters that fit
very well to the experience but don’t have any physical meaning, because with a large number of parameters we
can fit everything up to 1%or less.
The method we are going to use is the an EPM, with functions of the potentials with only a few parameters, but
with physical meaning, to fit to the experiment. It is possible for the band gap to be adjusted very precisely.
There is no black box! This means that everything is done from first principles, and everything is rederived from
the beginning, since this is a research project with pedagogical purposes. For this reason, the software we are
going to use is M ATHEMATICA , since we can see more clearly the programming we are doing while doing it.
1.2 State of The Art
Silicon and germanium are in the group of materials with indirect band-gap. For this reason they cannot emit
light effectively. It is desirable by the industry to be possible to silicon to be integrated with optical interconnects.
Different ideas are under discussion to meet this goal. One of the options are the strained superlattices of Silicon
and Germanium. By using strained layers it is possible to overcome indirect band behavior by the folding back
mechanism, which might allow the use of Si/Ge SL as a light emitting structure.
There is a limit on the number of strained-layers that can be accommodated on a given substrate – it is called
the critical thickness. For the case of germanium grown epitaxially on silicon, the maximum number of germanium
monolayers which can be deposited is six. Recently, it was shown that the addition of carbon into Silicon and
Germanium layers can be helpful to eliminate this problems. Adding an element with a much smaller radius than
that of silicon to a layer containing Ge atoms ( =rSi= 1:17˚A,rGe= 1:22˚A, butrC= 0:77˚A) also gives the possibility
to manipulate the strain, helping as well with thermal stability.
xviii

1.2.1 Si, Ge and C
Silicon, Germanium and Carbon are Group IV elements in the periodic table. They all crystalize in the diamond
crystal structure (Figure 1.1) which is an face centered cubic Bravais lattice with a two-point basis, with lattice
constantsa= 5:43˚A,a= 5:66˚Aanda= 3:57˚Afor Silicon, Germanium and Carbon, respectively. Carbon also
appears as graphite, which is more stable at ambient conditions.
Figure 1.1: The figure shows a diamond crystal structure
The band structure of Silicon and Germanium was calculated in reference [2] (Figure 1.2) using a tight-binding
model.
Figure 1.2: The energy-band structure of a) Si and b) Ge are calculated with a tight-binding model. The top of the valence
bands is set at zero energy. [2]
For Silicon, the distance between the conduction-band minimum and the point is equal to 0:89(2=a), where
ais the lattice constant. The values of the fundamental and direct gap for both Si and Ge are in Table 1.1 [2].
Silicon Germanium
Direct 3.41 0.90
Fundamental 1.05 0.76
Table 1.1: Calculated values for the gaps (in eV) are shown
As we can see, we have for both Si and Ge indirect band gaps in their electronic structure. Si, being an indirect-
gap semiconductor, silicon is not used in photonics and optoelectronics. The development of a direct-gap group
IV compound would have major advantages in the integration in microelectronics and optoelectronics. Because
recent advances in integrated microelectronics continue to rely on the unique properties of Si, the synthesis of
band-gap engineered, epitaxial, and preferably lattice-matched heterostructures on Si are considered essential for
future generations of high-speed devices. A lot of effort is been made to fabricate materials that could combine
the unique properties that the group IV elements are known to possess. These properties include a wealth of
technological applications (Si), the highest thermal conductivity (diamond), high hole mobility (Ge).
xix

1.2.2 Group IV Semiconductor Compounds and Alloys
When combining elements with different lattice constants, we have to take into account the effects of lattice
mismatch and strain. Lattice mismatch happens when a compound is grown on a different substrate. Both Si
and Ge crystallize in the diamond structure, but their lattice constants differ by about 4.2%. As a result, the Si/Ge
superlattices are under internal stress. This stress produces a distortion of the lattices and creates dislocations.
For thin layers, the growth can be well behaved, in which case the lateral lattice constant is the same in the Si and
Ge layers and equal to that of the substrate.
When two semiconductors like Silicon and Germanium are put together, discontinuities can occur in the band
structure. For a lattice matched interface (no strain), we need just to determine how the band structure of the two
materials line up at the interface. When the materials are strained, we have two problems in the calculation of
the band structure: Hydrostatic strain will produce additional shifts, and uniaxial or biaxial strain splits degenerate
bands. Figure 1.3 shows the different contributions. Thus, the total change in a band can be expressed as [3]

E=
Ea+
Eh+
Es; (1.1)
where
Eaa stands for the material differences (for the unstrained case),
Ehis the shift due to hydrostatic strain,
and
Ess reflects the splitting due to biaxial strain. It should be noted that each one of the contributions can have
different signs, compensating one another.
Figure 1.3: The energy-band structure of a) Si and b) Ge, calculated with a tight-binding model. The top of the valence bands
is set at zero energy. [3]
For growth with no strain, the lattice constant along the growth axis is reasonably given by Poisson’s ratio
=d"trans
d"axial, in which"trans is the transverse strain and "axial, the axial strain. The strain in each layer will then be
given by [2]
"jj=ajj
ai1"?=2
1"jj; (1.2)
in which"jjand"?are the lateral and perpendicular strain, respectively, ai,ajj, anda?are the equilibrium (bulk)
lattice constants of the strained material, of the substrate, and the lattice spacing perpendicular to the interface,
respectively. The lattice constant parallel to the interface ajjis the same along the structure, being forced to be the
same of the substrate.
1.2.2.A Si and Ge heterostructures
Diamond-Cubic SiGe [14] : Bulk Si and Ge form solid solution with diamond structure that has been said to
be as an almost ideal solution. The constructed with a (SiGe)-(GeSi) sequence of atomic layers ain the [111]
xx

direction, as shown in Figure 1.4 [14] and it is a Si 0:5Ge0:5composition. At the left, the ordering is described as
alternating Ge and Si 111 lattice planes, where 3/4 of the bonds are homopolar (Si-Si, and Ge-Ge) and 1/4 are
heteropolar (Si-Ge). This structure has been described as ”microscopically strained” because of the difference in
bond lengths. At the right 3/4 of the bonds are heteropolar (Si-Ge).
Figure 1.4: This models are consistent with the proposed structure of SiGe. Top: one-eighth of the ordered unit cell Bottom:
[101] projection through the unit cell.
Ultrathin (Ge) m/(Si)n[15] strained-layer superlattices grown on Si substrates have attracted considerable in-
terest. Thus, the electronic and optical properties of these superlattices can be changed to specific needs. It is the
possible to obtain a direct or a quasidirect band gap based on two indirect semiconductors.
It is interesting the difference between tetragonal and orthorhombic symmetry. The orthorhombic symmetry
occurs if the indices n and m are even. In Figure 1.5 [15] we have the crystal structure of Si 2Ge2and its electronic
structure on a Si substrate. The calculated lowest transition is indirect ( Eg= 0:90eVat0:95M), while the lowest
direct transitions in being 1.36 and 1.55 eV.
Figure 1.5: The unit cell of the (Ge) 2/(Si) 2superlattice and its band-structure are shown
Electroreflectance and photoluminescence experiments Ge/Si superlattices with n+m= 10 have stimulated
interest in the possibility of the existance of quasidirect transitions in this material. The minimum of the conduction
band occurs at 0:83Xin Si. It is expected that it is folded back to k= 0for a total periodicity of ten.
Besides the proper periodicity, another requirement has to be met in order to obtain a quasidirect transition; the
proper strain distribution and thus ajj. It is needed tensile strain in the Si layers to lower the minimum of the twofold
states (which are folded back to ) below those of the four other states. This is represented in Figure 1.6 [15].
In Figure 1.7 the direct ( E) and the indirect transitions are plotted as a function of composition and strain.
The energy of the quasidirect Etransition remains constant for Si substrates. The biggest variation of indirect
gaps with composition is found for the ENgap. It shows an approximately linear description with increasing Si
xxi

Figure 1.6: The figure shows a) electronic band structure of the (Ge) 6/(Si) 4supperlattice on Si [001] substrate and band
structure of the free-standing b) (Ge) 5/(Si) 5c) (Ge) 4/(Si) 6superlattice along lines of high symmetry.
content [0:14(mn)eV].
Figure 1.7: Transition energies of various n+m= 10 superlattices s as a function of lateral strain in the Si layer are represented
SixGe1xalloys [16] : Alloys of Si xGe1xhave continuously variable lattice parameters and band gap (although
is indirect), and they have potential for practical applications. For example, they have been successfully used to
create heterojunction bipolar transistors with cutoff frequencies bigger than 100 GHz, which is much higher than
for the traditional Si ones. A problem facing Si-Ge technology is the mismatch that causes compressive strain in
SixGe1xlayers grown on Si. According to Vegard’s law, in a Si xGe1xalloy,ajj=xaSi+ (1x)aGe, in whichaSi
andaGeare the cubic lattice constants of Si and Ge structures, respectively. The strain increases with increasing
Ge concentration and as increasing Si-Ge film thickness as well.
The band structure for Si xGe1xalloys was calculated in reference [16] (Figure 1.8)
The band gap, Eg, is indirect, with the valence-band maximum in the point and the conduction-band minimum
changes from L[L= (2=a)(1
2;1
2;1
2)]to near the point X[= (2=a)(1;0;0)]. The change occurs at approximately
x= 0:25(for a temperature of 4K). It is this feature that makes the defect levels of this alloy interesting to study,
since alloys with compositions near the x= 0:25can possibly produce deep levels in the gap for impurities such
as As and P .
Near the band gap, every sp 3bonded impurity with a valence greater than that of tetrahedrally bonded host
by unity is expected to have an s-like level, a triply degenerate p-like deep level. In Figure 1.9 [16] is shown
the predicted single-impurity defect levels of p-like and s-like symmetries. Shown are also in these figures the
conduction-band minima as functions of composition xwhere the zero of energy is taken to be the to of the
valence band maxima for all x.
Magic sequence SiGe 2Si2Ge2Si and the genetic algorithm [17] : It was identified the sequence of Si and Ge
xxii

Figure 1.8: Band structures E(~k)for the principal symmetry directions of the diamond lattice for (a) Si 0:2Ge0:8and (b)
Si0:74Ge0:26where calculated
Figure 1.9: The figures show the predicted single-impurity defect levels of a) T2and b)A1symmetry as a function of the
composition x. Shown are also the conduction and valence band edges
layers with strong transition across the electronic band-gap from amongst all possible superlattices [Si n0Gep0/Sin1Gep1/…/SinNGepN]1
including substrate orientation and strain.
It is used a efficient research method: It is constructed a population of superlattices according to chance
and their relative success, namely, their ability for light-emission at the band-edges. New superlattice candidates
(offspring) are created from the previous population by interchanging random sets of layers in the superlattice
between two parents (crossover), and by flipping random Ge layers into Si layers and vice-versa in a single parent
(mutation). At each generation, the worst individuals in the previous population are replaced by the offspring, thus
guiding the population as a whole towards the global optimum through survival of the fittest. To judge fitness, i.e.,
the strength of the optical transition, it is computed the dipole matrix element between the valence band minimum
and conduction band maximum at of each superlattice candidate, which is directly proportional to the strength of
the optical transition.
The set of results are a variation of a magic sequence composed of =SiGe 2Si2Ge2Si followed by a Ger-
manium buffer layer of n= 1232monolayers. The magic sequence satisfies: wave vector directness and the
dipole matrix element between the valence band maximum and the conduction band minimum is nonzero. The
first condition is satisfied when the structure is grown on substrates Si 1xGexwithx0:4(Figure 1.10).
The second condition is shown by the spectrum of absorption in Figure 1.11 top. It also contains the spectrum
of absorption absorption for the current state-of-the-art superlattice Si 6Ge4. This puts the magic sequence much
xxiii

Figure 1.10: Band structures of the 12magic sequence grown on Si 0:4Ge0:6where calculated
more adjusted for practical applications.
Figure 1.11: The figure shows the omparison of between Si 6Ge4superlattice and the magic sequence of the direct absorption
spectra
We can evaluate the effect of deviation of the best sequences in the “directness” of the band gap. By changing
the substrate and the mutations in the magic sequence, it was constructed Figure 1.12. nis the magic se-
quence with a Ge buffer of nmonolayers, while is the sequence SiGe 2Si2Ge2SiGe 2SiGe 9and
the sequence
SiGe 2SiGe 2Si2Ge2SiGe 2SiGe 6.
Figure 1.12: Direct dipole-allowed band-gaps where measured
1.2.2.B Si and Ge heterostructures with C impurities
A possible solution to the mismatch problem is the incorporation of C which has a lattice constant of 3:57˚Ain a
Diamond crystal structure, which is much smaller than those of Si and Ge. Incorporation of C into SiGe material
should reduce the lattice mismatch because of the smaller size of C, compensating for the larger size of Ge. The
xxiv

linear approximation for the lattice constant between Si, Ge, and diamond is [3]
a(x;y) = (1xy)aSi+xaGe+yaC; (1.3)
resulting in a Ge:C ratio of 8.2 for strain compensation. The incorporation of a third component also adds additional
flexibility in band-gap engineering. This could pose a challenge to the GaAs technologies. We define an “effective
lattice-mismatch” mfeffas [3]
mfeff=a(x;y)aSi
aSi; (1.4)
for ternary Si 1xyGexCyon Si(001) substrates. A positive mfeffthat the material is compressively, mfeff<0is
tensile strain, and mfeff= 0means strain-compensated Si 1xyGexCyalloy. The hydrostatic contribution is [3]

Eh=av;c("?+ 2"jj); (1.5)
whereav;cis the appropriate hydrostatic deformation potential for the valence or conduction band, respectively.
For the material dependent term
Ea,

Ea(x;y) =
Ea(x) +
Ea(y); (1.6)
Figures 1.13, 1.14 and 1.15 summarize the results for the offsets of strained Si 1xyGexCyon Si(001) [3]. The
effective concentration corresponds to the concentration needed for identically strained binary layers.
Figure 1.13: Valence-band offsets for compressive strained Si 1xGex, and Si 1xyGexCy(x=10%, 20%, and 30%) and
tensile strained Si 1yCyand Si 1xyGexCy(y=1%, 2%, and 3%) are plotted as a function of the effective lattice mismatch
(expressed in “effective” Ge or C concentrations, respectively)
Figure 1.14: Conduction-band offsets for compressive strained Si 1xGex, and Si 1xyGexCy(x=10%, 20%, and 30%) and
tensile strained Si 1yCyand Si 1xyGexCy(y=1%, 2%, and 3%) are plotted as a function of the effective lattice mismatch
(expressed in “effective” Ge or C concentrations, respectively)
xxv

Figure 1.15: The band-gap narrowing for ternary Si 1xyGexCyalloys, strained on Si(001) is represented as a function of
lattice mismatch
The band-gap narrowing is obtained by adding the valence and the conduction-band offsets. The band gap
for the alloys is always smaller than that of silicon. The addition of C (Ge) into compressive strained Si 1xGex
(tensile strained Si 1yCy) leads to a smaller change in band-gap narrowing than an equivalent strain reduction in
the binary alloy (lower Ge or C content, respectively).
1.3 Thesis Outline
This thesis is organized as follows:
In chapter 2 we will make a theoretical introduction with the physics needed to this project. Starting by the
by the many body problem, going through the Density Functional Theory (DFT) we get to the conclusion we
can use an independent electron approximation to the problem and we do this with the help of pseudopotentials.
We explain what is a pseudopotential, using not only local pseudopotentials but also the non-local (important in
Carbon) and spin-orbit contributions (important in the heavier Germanium) that are needed to fully describe bulk
group IV elements in question.
In chapter 3 we describe the research for the best fitted pseudopotentials to the experiment, that describe each
one of this elements. We start by checking a previous description used for Silicon with only a local pseudopotential.
After we improve the description by adding a non-local contribution to the pseudopotential. We search for the best
fitted local and non-local parts of the pseudopotential to the experiment and calculate the optical properties. The
same search is done for Carbon. For Germanium we add a spin-orbit contribution and find as well the best fitted
potentials (local, non-local and spin-orbit) to the experimental results. This research work will be done in the
software M ATHEMATICA .
The software used in good for the research and learning process, but has limitations such as the computation
time increases highly with the precision requested for the calculation. In chapter 4 introduced the discovered
pseudopotentials in a F ORTRAN program, where were obtained better results namely in the band structure and the
density of states for each of the bulk elements. An indication of future research was also made.
xxvi

2
Theoretical Introduction
Contents
2.1 The many-body problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxviii
2.2 The Adiabatic Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxix
2.3 Separable Schr ¨odinger equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxix
2.4 The Hartree-Fock Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxx
2.5 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxi
2.6 Hohenberg-Kohn Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxii
2.7 Kohn-Sham Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxiii
2.8 Pseudopotentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxv
2.9 Ab initio pseudo-potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxv
2.10 Empirical Pseudopotential methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxvii
2.11 Non-local and Spin-Orbit Pseudopotentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxviii
2.12 Screening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xlii
2.13 Optical properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xlv
2.14 Imaginary part of the dielectric function 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xlvi
xxvii

We want to calculate the properties of a solid – a condensed matter system. This task can be quite complicated,
if not impossible, to solve. Most of the time (almost all the time) we have to use approximations to calculate the
properties of the system, which means we have always to keep in mind those approximations and what they
change in the final result. If start from the beginning of the problem: the many body problem, because that is what
we have, a solid with a big number of electrons and nucleus. After we are going to prove that an independent
electron approximations with the use of pseudopotentials is enough to describe our solid. What is here written,
from section 2.1 to 2.7, is based on the notes of reference [18]
2.1 The many-body problem
Consider a system of melectrons and nnuclei, each one with spatial coordinates ~ riand~Rjand spiniand
j. In the non-relativistic limit, the time-independent Schr ¨odinger equation can be used to calculate the properties
of the system is
H(~ r1;1;:::;~ rm;m;~R1;1;:::;~Rn;n) =E(~ r1;1;:::;~ rm;m;~R1;1;:::;~Rn;n); (2.1)
where the many particle Hamiltonian is
H=nX
i=11
2Mir2
~Ri+mX
j=11
2r2
~ rj+X
i<jZiZj
j~Ri~Rjj+X
i<j1
j~ ri~ rjjX
i;jZi
j~Ri~ rjj: (2.2)
The subscripts in the Laplacian r2indicate on which coordinate they operate. MiandZiare the mass and the
charge, respectively, of the ith nucleus. The terms in the Hamiltonian are, respectively, the Kinetic Energy of the
nuclei and the electrons, the Potential Energy related to the potential caused by the nuclei and felt by the nuclei,
the potential caused by the electrons and felt by them and the potential caused by the ions and felt by the electrons
(or vice versa). Notice that it is written in an adimentional form, adequate to computational purposes. This units
are called atomic units, a system in which the numerical values of following fundamental physical constants are all
unity by definition: the electron mass, me, the elementary charge, e, the reduced Planck’s constant ~=h
2and the
Coulomb’s constant1
40, where0is the dielectric permittivity of vacuum. In this work, all the equations will be
written in the atomic system of units unless it is said otherwise.
Since electrons are fermions, the wavefunction is anti-symmetric with respect to the electron coordinates,
position and spin (Pauli’s exclusion principle),
(:::;~ ri;i;:::;~ rj;j;:::) =(:::;~ rj;j;:::;~ ri;i;:::): (2.3)
An analytical solution for the equation (2.1) with the Hamiltonian (2.2) is only known for very simple systems
such as the Hydrogen atom and the H+
2ion, and numerical solutions can be obtained for atoms and molecules
with a small number of electrons such as the He atom or the H 2molecule. These are very simple systems. Not
such miracle is to be expected with solids of higher dimensions. If we could solve the exact Hamiltonian for these
more complicated systems we could in principle predict all of it’s properties. In a macroscopic solid, there are
about 1023nuclei and a similar number of electrons [19]. The equation (2.1) to be solved would have something in
the order of 1023variables, which is not possible to be solved with the current computational technology.
Therefore we are compelled to make approximations motivated by physical considerations, such as the ones
described in the following sections.
xxviii

2.2 The Adiabatic Approximation
The first approximation to be made takes into account that the mass of the nuclei is much larger (generally
104105times [19]) than the mass of the electrons, that is Mime. Therefore we can say that the nuclear motion
is much slower (or say the nuclei are fixed) than the motion of the electrons. In this way we can neglect their Kinetic
Energy, obtaining the Hamiltonian Hefor the electronic problem, which depends on the nuclear coordinates,
He=mX
j=11
2r2
~ rj+X
i<jZiZj
j~Ri~Rjj+X
i<j1
j~ ri~ rjjX
i;jZi
j~Ri~ rjj
He (k)(~ r1;:::;~ rm;~R1;:::;~Rn) =U(k)(~R1;:::;~Rn) (k)(~ r1;:::;~ rm;~R1;:::;~Rn); (2.4)
and whose eigenvalues also depend on the position and spin of the nuclei and form a family identified by the
quantum number (k). This is called the Born-Oppenheimer or adiabatic approximation. The energy U(k)can be
considered as the potential energy in the Hamiltonian used Shcr ¨odinger equation that describes its motion
H(k)
N=nX
i=11
2Mir2
~Ri+U(k)(~R1;:::;~Rn)
H(k)
N(k;q)(~R1;:::;~Rn) =E(k;q)(k;q)(~R1;:::;~Rn); (2.5)
where the quantum number qconcerns the vibrational, rotational and translational states. We say that the total
wave function is the product of the solutions of the equations (2.4) and (2.5),
(k;q)(~ r1;:::;~ rm;~R1;:::;~Rn)' (k)(~ r1;:::;~ rm;~R1;:::;~Rn)(k;q)(~R1;:::;~Rn); (2.6)
which is a complete set expansion of the wavefunction  =P
k;qck;q (k)(k;q). The decoupling of the electronic
and nuclear motions can also be obtained using perturbation theory. This is the Born-Oppenheimer aproximation
[20]
2.3 Separable Schr ¨odinger equation
But the problems in the equation (2.2) are almost the same as the ones in equation (2.4) since we still have a
number of electrons in the order of 1023which interact with one another through the Coulomb interaction and thus
we cannot separate the Schr ¨odinger equation in its different variables. We could neglect the Coulomb interaction
but unfortunately we know this is a bad approximation. Another solution is to consider that a given electron is
subjected by a potential depending on the average distribution of the electrons. This results in the Hamiltonian [19]
Hs=mX
j=11
2r2
~ rj+mX
j=1nX
i=1Vat(~ rj~Ri)
=mX
j=1"
1
2r2
~ rj+nX
i=1Vat(~ rj~Ri)#
=mX
j=1Hj; (2.7)
whereVatincludes the average Coulomb electron-electron and the electron-nucleus interaction. The nucleus-
nucleus term is dropped, since it is a constant for the same position of the nuclei. Like this it is possible to separate
the Hamiltonian into a sum of mindependent terms.
xxix

2.4 The Hartree-Fock Method
Considering we have a separable Hamiltonian like the one in equation (2.7), we can construct and anti-
symmetric wave-function using a Slater determinant of orthonormal one electron spin orbitals, i(~ r); i= 1;:::;m ,
D(~ r1;:::;~ rm) =1p
m!1(~ r1)1(~ r2):::  1(~ rm)
2(~ r1)2(~ r2):::  2(~ rm)
…………
m(~ r1)m(~ r2)::: m(~ rm); (2.8)
in whichhijji=ij. In this case we can show that hDjDi= 1.
We use the variational principle in quantum mechanics to calculate de ground state energy E0and wavefunction
0of the many-body system,
E0= minh jHj i
h j i; E 0=h 0jHj 0i
h 0j 0i: (2.9)
In the Hartree-Fock method we search the wavefuntions in the form of slater determinants,
EHF= minhDjHjDi; EHF=hDHFjHjDHFi; (2.10)
where we always have that EHFE0. To minimize the energy in equation (2.10) the Euler-Lagrange method can
be used,


jh
hDjHjDiX
ijij
hijjiij
i
= 0; (2.11)
whereijare the Lagrange multipliers. The expectation value of the energy for the Slater determinant is
hDjHjDi=X
iZ

i(~ r)
1
2r2

i(~ r)d3r
+X
iZ

i(~ r)X
jZj
j~ r~Rjj

i(~ r)d3r
+X
i<jZiZj
j~Ri~Rjj
+1
2X
i;jZZ

i(~ r)
j(~ r0)1
j~ r~ r0ji(~ r)j(~ r0)d3rd3r0
1
2X
i;jZZ

i(~ r)
j(~ r0)1
j~ r~ r0ji(~ r0)j(~ r)d3rd3r0; (2.12)
that is a sum of five contributions: kinetic, external potential, ion-ion, Hartree and exchange. The ion-ion contribu-
tion is a constant. We can define the external ionic potential,
v(~ r) =X
jZj
j~ r~Rjj; (2.13)
and the total electronic charge density,
(~ r) =hDjX
i(~ ri~ r)jDi=X
j
j(~ r)j(~ r): (2.14)
With this definitions we can rewrite the external potential and the Hartree contributions,
Z
v(~ r)(~ r)d3r1
2ZZ(~ r)(~ r0)
j~ r~ r0jd3rd3r0: (2.15)
xxx

Introducing (2.15) into (2.12) and the later in (2.11) we obtain the Hartree-Fock equations,
h
1
2r2+v(~ r) +Z(~ r0)
j~ r~ r0jd3r0i
i(~ r)X
jj(~ r)Z

j(~ r0)1
j~ r~ r0ji(~ r0)d3r0=X
jijj(~ r): (2.16)
Knowing that the Slater determinants are invariant with respect to unitarian transformations i! 0
j=
P
iUij i, we can use this do diagonalize (2.16), obtaining
h
1
2r2+v(~ r) +Z(~ r0)
j~ r~ r0jd3r0i
i(~ r)X
jj(~ r)Z

j(~ r0)1
j~ r~ r0ji(~ r0)d3r0=ii(~ r): (2.17)
The correlation energy Ec=E0EHFis the difference between the exact energy of a system and the Hartree-
Fock energy. This equation is a non-linear eigenvalue differential integral equation in 3 dimensional space. The
Hartree-Fock method provides a connection between the many body wavefuntion to mone body wavefunctions.
2.5 Density Functional Theory
Density functional theory is a method to investigate the electronic structure in the ground state of atoms,
molecules our condensed matter systems. It says that the properties of a many-electron system can be uniquely
determined by the ground state electron density of the system, that depends on 3 spatial coordinates.
But before we have to make sure that the properties of the system can be indeed be uniquely associated to a
determined electron density (ilustration in the Venn Diagram in Figure 2.1).
Figure 2.1: The figure shows the Venn’s diagram corresponding a potential in the space of all potentials V(~ r)to a ground state
electron density in the space of the ground state electron densities GS(~ r)
If it is univocal, in principle, by knowing the electron density, we can obtain the potential that acts on this
electrons and from that we can calculate the ground state eigenfunctions and all the other properties of the system.
First let us assume that two different potentials can lead to the same ground state electron density,
V1(~ r)- 1(~ r)XXXXX z
V2(~ r)- 2(~ r) :GS(~ r),
xxxi

and consider the Hamiltonian
H=T+Vee+Vext; (2.18)
where,
T =P
i1
2r2
~ ri
Vee=P
i<j1
j~ ri~ rjj
Vext =P
iv(~ ri); (2.19)
are the Kinetic, Coulomb Potential and external Potential energies, respectively and v(~ ri)is calculated using equa-
tion (2.13). We calculate de ground state energy using V1,
EGS1=h 1jT+Vee+V1j 1i
=h 1jT+Veej 1i+Z
V1(~ r)GS(~ r)d3r: (2.20)
Using the variational principle,
h 2jT+Vee+V1j 2i=h 2jT+Veej 2i+Z
V1(~ r)GS(~ r)d3r
=h 2jT+Veej 2i+Z
V2(~ r)GS(~ r)d3r+Z
(V1V2)(~ r)GS(~ r)d3r
=EGS2+Z
(V1V2)(~ r)GS(~ r)d3r; (2.21)
so,
EGS1<EGS2+Z
(V1V2)(~ r)GS(~ r)d3r: (2.22)
But starting with the calculation of EGS2and following the same line of thinking we get
EGS2<EGS1+Z
(V2V1)(~ r)GS(~ r)d3r; (2.23)
which results with,
EGS1> EGS2+Z
(V1V2)(~ r)GS(~ r)d3r
EGS1< EGS2+Z
(V1V2)(~ r)GS(~ r)d3r; (2.24)
that is a contradiction. We reach then the conclusion that in the absence of degeneracies, two different potentials
cannot lead to the same electron density in the ground state. This means that any property of the many-body
system is a functional of (~ r). From that comes the name Density Functional Theory.
2.6 Hohenberg-Kohn Theorem
Considering a system of melectrons, the electronic Hamiltonian can be described by Equation (2.18), we
define the set of all normalized anti-symmetric wavefunctions as
A=f j (:::;~ ri;:::;~ rj;:::) = (:::;~ rj;:::;~ ri;:::)andh j i= 1g: (2.25)
xxxii

We can define the ground state energy of the system as a functional of the external potential Vext(~ r)as
E[Vext] = min
2Ah jHj i: (2.26)
The subset AofAis the set of all wavefunctions that correspond to the charge density (~ r),
A=f j 2AandmZ
d3r2:::Z
d3rmj (~ r;~ r2;:::;~ rm)j2=(~ r)g: (2.27)
An universal functional of the charge density is
F[] = min
2Ah jT+Veej i; (2.28)
which is independent of Vext(~ r). The ground state energy can be calculated as a functional of the electron density
for each potential Vext,
E0=EGS[] = min
2Ah jT+Vee+Vextj i
= min
min
2Ah jT+Vee+Vextj i
= min

min
2A(h jT+Veej i+h jVextj i)
; (2.29)
in whichh jVextj i=R
Vext(~ r)(~ r)d3r, so
EGS[] = min
Z
Vext(~ r)(~ r)d3r+ min
2Ah jT+Veej i
= min

F[] +Z
Vext(~ r)(~ r)d3r
= min
EVext[]: (2.30)
Therefore Hohenberg-Kohn theorem states that the minimum of the energy functional EVextis the ground state
energyE0of the system. Note that contrary to the above proof, we did not require the ground state of the system
to be non-degenerate.
2.7 Kohn-Sham Equations
Consider the set of one electron wavefunctions iand occupation numbers fithat have a charge density equal
to,
B=f(f1;:::;fk;1;:::;k)jhijji=ij;kX
i=1fiji(~ r)j2=(~ r)g; (2.31)
where some authors define 0fi1, and we define an energy functional of the charge density that it is called
kinetic energy functional but it is not the true kinetic energy of the system
T0[] = min
(f1;:::;fk;1;:::;k)2BkX
i=1fihij1
2r2jii: (2.32)
The exchange and correlation functional energy is defined as
Exc[] =F[]1
2ZZ(~ r)(~ r0)
j~ r~ r0jd3rd3r0T0[]; (2.33)
whereF[]was defined in equation (2.28). It contains the rest of the many-body contributions to the energy. Also
sinceT0[]is not the exact kinetic energy of the interacting energy but instead the kinetic energy of the ground
xxxiii

state of a system of non interacting electrons with density (~ r), it contains also kinetic energy terms related to the
electron-electron interaction. Rewriting it in a different way, we have
F[] =T0[] +1
2ZZ(~ r)(~ r0)
j~ r~ r0jd3rd3r0+Exc[]: (2.34)
We can calculate the energy functional,
EVext[] =F[] +Z
Vext(~ r)(~ r)d3r
=T0[] +Z
Vext(~ r)(~ r)d3r+1
2ZZ(~ r)(~ r0)
j~ r~ r0jd3rd3r0+Exc[]; (2.35)
and the ground state energy,
EGS[] = min
EVext[]
= min

T0[] +Z
(~ r)Vext(~ r)d3r+1
2ZZ(~ r)(~ r0)
j~ r~ r0jd3rd3r0+Exc[]
= min

min
(f1;:::;fk;1;:::;k)2BkX
i=1fihij1
2r2jii+Z
(~ r)Vext(~ r)d3r+
+1
2ZZ(~ r)(~ r0)
j~ r~ r0jd3rd3r0+Exc[]
= min
min
(f1;:::;fk;1;:::;k)2BkX
i=1fihij1
2r2jii+Z
(~ r)Vext(~ r)d3r+
+1
2ZZ(~ r)(~ r0)
j~ r~ r0jd3rd3r0+Exc[]
; (2.36)
but this is the same as minimizing over all wavefunctions and therefore,
EGS[] = min
(f1;:::;fk;1;:::;k)kX
i=1fihij1
2r2jii+Z
(~ r)Vext(~ r)d3r+1
2ZZ(~ r)(~ r0)
j~ r~ r0jd3rd3r0+Exc[]
:(2.37)
The Euler equation that minimizes the expression above with respect to the one electron wavefunctions iis


ikX
i=1fihij1
2r2jii+Z
(~ r)Vext(~ r)d3r+1
2ZZ(~ r)(~ r0)
j~ r~ r0jd3rd3r0
+Exc[]X
i;ji;j(hijji)ij
= 0; (2.38)
withVextgiven by Equation (2.19), is called the Kohn-Sham equation. The result of the minimization is

1
2r2+v(~ r) +vH(~ r;] +vxc(~ r;]
i(~ r) =ii(~ r);
vH(~ r;] =Z(~ r0)
j~ r~ r0jd3r0;
vxc(~ r;] =Exc[]
(~ r);
(~ r) =kX
j=1fj
j(~ r)j(~ r); (2.39)
wherevHis the Hartree potential and vxcis the exchange and correlation potential. The curve parameter in (~ r;]
means that the function is dependent of a variable, and the rect parameter means that it is a function. Although
xxxiv

the equation resembles the Schr ¨odinger equation for non-interacting particles, the dependence of vHandvxcin
the charge density makes it a non-linear system of equations. A common approximation is the Local Density
Approximation (LDA) that says that it depends only on the charge density in the point of interest vxc(~ r) =vxc[(~ r)].
The way to solve this kind of equations is through iterative, self-consistent methods. The kind of logic is the
following:
1. Guess initial in
2. Calculate vH[in]andvxc[in]
3. Solve the Kohn-Sham equation
4. Calculate out=P
iji(~ r)j2
5. Ifinoutstop. If not in=F(in;out)and start from 2.
Although it was a highly used method ab initio method, the DFT with LDA gives the wrong value for the band
gap, namely half of the value for the band gap of Silicon.
2.8 Pseudopotentials
The pseudopotential model describes a solid as a sea of valence electrons moving in a periodic background
of cores. The space can be divided into two regions: the region near the nuclei, the “core” composed primarily of
tightly bound core electrons which are not very affected by the neighbour atoms; and the valence electron region
which is involved in bonding the atoms together. This results that the atoms in the same group – such as Carbon,
Silicon and Germanium (group IV, for ex.) are treated in mostly the same way – apart from a few “details”. The focus
of the calculation is only on the accuracy of the valence electron wavefunction away from the core. The potential
in the ion core is strongly attractive for the valence electrons, but the requirement for the valence wavefuntions to
be orthogonal to those of the core contributes to an effective repulsive potential for valence states. This results in
a net weakly attractive potential that affects the valence electrons.
2.9 Ab initio pseudo-potentials
We first consider an atom of atomic number Z. An one-electron Hamiltonian can be written as [4]
H=1
2r2+Vion+Vscr; (2.40)
whereVion=Z
ris the ion core potential, that can be taken as a linear superposition of spherical potentials, and
Vscris the screening potential, a potential very important in many body physics. Usually it is divided into two parts
(as it was seen before), the Hartree potential, VH(~ r;], that comes from the Poisson’s equation
r2VH=4e2(~ r); (2.41)
where(~ r)is the valence electron charge density. The other part is the exchange and correlation potential Vxc,
that was also mentioned before. If we use the Local Density Approximation (LDA), then Vxc(~ r) =Vxc[(~ r)]. The
total potential is, thus
VT(~ r) =Vion(~ r) +VH+Vxc(~ r): (2.42)
xxxv

If there is one state for which we know the the wavefuntion and the value of the energy we can invert the
Schr ¨odinger equation to obtain the total potential [4]
VT=1
2r2
+E: (2.43)
This equation is well behaved if is nodeless, since it is highly preferable for the pseudopotential to be smooth
and the wiggles associated with the nodes are undesirable. The quantityr2
is extremely sensitive to numerical
errors when !0. If there are no numerical errors, what normally happens is that if !0, thenr2 !0as
well. If there is an error and r2 doesn’t go to zero when !0, this quantity will diverge. But equation 2.43
is the simplest possible case. For example, we can extract the energy levels of interest by performing an atomic
structure calculation starting from all electron atomic calculations. Within the density-functional theory this is done
by assuming a spherical screening approximation and self-consistently solving the radial Kohn-Sham equation [21]

1
2d2
dr2+l(l+ 1)
2r2+VT(~ r;]
rRnl(r) ="nlrRnl(n); (2.44)
that results in the “all electron” wavefuntions and energies. We have to take into account that equation (2.43) is only
well behaved if the wavefuntions used have no nodes. This can be achieved by the construction of pseudo-wave
functions with no nodes (for this reason, the quantum number nwill be omitted in the further calculations) based
on the wavefunctions of the equation (2.44) as it was done successfully in reference [21]. Other characteristics of
this pseudo-wave funtion are [21]: the normalized atomic radial pseudo-wave-function with angular momentum lis
equal to the normalized radial all-electron wave function after a cutoff radius rcl,
RPP
l(r) =RAE
l(r)forr>rcl; (2.45)
the charge enclosed for the two wavefuntions within rclmust be equal,
Zrcl
0jRPP
l(r)j2r2dr=Zrcl
0jRAE
l(r)j2r2dr; (2.46)
so that the norm of the wavefunction is conserved after normalization; and the valence all-electron and pseudopo-
tential eigenvalues must be equal,
"PP
l="AE
l: (2.47)
A pseudopotential under this conditions is called a “norm-conserving pseudopotential”. Once we have the
pseudo-wave function we can calculate the pseudopotential by inversion of the Schr ¨odinger equation [21],
VPP
l="ll(l+ 1)
2r2+1
2rRPP
l(r)d2
dr2[rRPP
l(r)]: (2.48)
By inverting the Shcr ¨odinger equation for each of the wavefunctions separately, for the same energy level n,
we get with different potentials for each quantum number l,Vl. This is called non-locality of the pseudopotential.
The pseudopotential is decomposed into a sum over angular momentum components [4],
VT=V0P0+V1P1+V2P2+:::; (2.49)
where theP`projects out the `thangular momentum component,
Pl=lX
m=lj`mih`mj; (2.50)
xxxvi

whereh~ rj`mi=Y
`m(
)andY
`m(
)is centered on the atom. Another complication is to take into account the spin
orbit effects in heavier elements (like Ge). Non-locality and spin-orbit considerations will be further developed in
later sections.
2.10 Empirical Pseudopotential methods
The empirical pseudopotential method relies on experimental results for the construction of the pseudopotential
and the predictions made with the pseudopotentials should converge as best as possible with experience. The
history of the empirical pseudopotential methos
Lets assume first that the pseudopotential is local, i.e., independent of `. The Schr ¨odinger equation for a
periodic system is [22]

1
2r2+V(~ r)
~k(~ r) =E(~k) ~k(~ r): (2.51)
In a crystal, the potential V(~ r)is periodic in the lattice. We can use a plane wave expansion that will only have
plane waves with the periodicity of the lattice. With the local approximation, the pseudopotential can be written as
V(~ r) =X
~GV(~G)S(~G)ei~G
~ r=X
~GU(~G)ei~G
~ r; (2.52)
where~Gis a reciprocal lattice vector, V(~G)are the form factors and S(~G)is the structure factor,
S(~G) =1
NaNaX
i=1ei~G
~ i: (2.53)
Once the form factors are decided, we solve (2.51). We can assume that the wave functions ~k(~ r)can be
expanded in plane waves, with no loss of generality and solve the secular equation, which is the Schr ¨odinger
equation (2.51) in the reciprocal space [22] [23],
detjH(~k;~G~G0)E(~k)Ij= 0; (2.54)
where
H(~k;~G~G0) =1
2(~k~G)2~G;~G0+V(~G~G0)S(~G~G0): (2.55)
The form factors depend only on the magnitude of j~G~G0jif the pseudopotential can be taken as spherically
symmetric, which is generally the case for tetrahedral semiconductors [22]. In the Chelikowsky pseudopotential,
for diamond or zinc-blende semiconductors, generally only three form factors are enough to determine the pseu-
dopotential, those for G2= 32
a
2;82
a
2;112
a
2. The factor G2= 0is not important since it only gives the
level zero of energy and S(~G) = 0 forG2= 42
a
2;122
a
2. So of the six smallest reciprocal lattice vectors, only
the form factors G2= 32
a
2;82
a
2;112
a
2are required to specify the crystalline potential [4][22] (Figure 2.2).
These three values are fitted to optical transition energies and the whole band structure follows from them.
Before solving the equation (2.55), the values of the form factors are needed. With the EPM the form factors
are obtained by fitting them to the experiment. Usually they are adjusted to the optical data. The method is similar
with the one discussed to DFT:
1. Estimate initial V(~G)
2. Solve secular equation
xxxvii

Figure 2.2: The figure shows a schematic plot of a pseudopotential in reciprocal space with the G’s that correspond to
G2= 3;8;11withG2in units of (2=a)2[4]
3. Calculate band structure and optical properties
4. Compare with experiment
5. If it agrees with experiment stop. If not change V(~G)and start from 2.
In reference [24], a semi-empirical pseudopotential is used. First an ab initio method is used, in which spherical
atomic potentials (with only the local part) v(r), such that the solutions of
(
1
2r2+Vnonlocal (~ r) +X
X
Rv()(j~ r~Rj))
~ i= ~i~ i; (2.56)
whereis the chemical atomic type and ~Rstands for all possible atomic positions of , including those related
by lattice translations, will have large overlaps with the LDA solutions iandi, so that they reproduce the LDA
results for bulk systems with a good approximation. This means that it will suffer from poor reproduction of the
observed excitation energies. Therefore, the potential described in the reciprocal space by
v()(q) =X
iC(n)e(qan)2=b2
n; (2.57)
whereC(n),anandbnare free parameters (the author used 20 of each), that will be further adjusted, this time
empirically, to reproduce the experimentally observed excitation energies. Notice the big number of parameters
used in the fit.
This work is based in the Empirical Pseudopotential Method with some differences. The parameters that are
adjusted to the experiment are not the form factors but the parameters of a function we define as the pseudopoten-
tial. Also, as this work is predicted to be used in superlattices, we calculate the form factors for the pseudopotential
in a lattice of equally distant points in the reciprocal space. If we want to describe a superlattice, we need to fit the
whole curve of the potential, because in a supperlattice, the ~Gvectors may not be constant. Like in the work of
reference [24], we will be adjusting an empirical expression to experimental results, but the number of parameters
used will be much less and each one will have a physical meaning.
2.11 Non-local and Spin-Orbit Pseudopotentials
If we take into account the spin-orbit effects, we can obtain pseudo-wave-functions RPP
`j(r), with energies "`j
and normalizationR1
0r2jRPP
`j(r)j2dr= 1that are constructed from the respective all-electron wave-functions. From
xxxviii

the inversion of the radial Schr ¨odinger equation we obtain the correspondent ionic pseudopotentials VPP
`j(r). The
indexjtakes the values `1
2, except for `= 0, where the only allowed value is j=1
2. It is in this distinct
pseudopotentials for j=`1
2andj=`+1
2that the effect of the spin-orbit is included in the calculations,
as the major spin-orbit effect is in the core region, since the dominant contribution comes from the motion of
electrons in the Coloumb potential in the innermost region of the atomic cores. To restrict the non-local part of
the pseudopotential to the core region we define a local potential VL(r)that is arbitrary in the core region, and is
identical to the pseudopotentials VPP
`j(r)outside the core region ( VL=VPP
`jforr >rc). We define the non local
pseudopotential as

VNL
`j(~ r;~ r0) =VPP
`j(~ r;~ r0)VL(~ r): (2.58)
The non-local part of the pseudopotential for ` > ` maxcan be neglected as long as the local part and `max
reasonably chosen.
It is convenient also to separate the pure spin-orbit part from the average non-local pseudopotential, because
the spin-orbit can often be treated as a small perturbation, which is often not the case of the non-local component.
We therefore define for `>0the degeneracy weighted average

VNL
`(~ r;~ r0) =`
2`+ 1
VNL
„1
2(~ r;~ r0) +`+ 1
2`+ 1
VNL
„+1
2(~ r;~ r0); (2.59)
and the spin-orbit part

VSO
`j(~ r;~ r0) =
VNL
`j(~ r;~ r0)
VNL
`(~ r;~ r0): (2.60)
In its semi-local form the action of the pseudopotential operator on a spinor is, in a compact notation
VPP=VL+X
`X
mj`miVNL
`h`mj+X
`X
jX
mjj`jmjiVSO
`jh`jmjj; (2.61)
wherej`miandj`jmjiare angular momentum states, the first with just orbital components the second with the
composition of the orbital and spin angular momenta. Semi-local means that it is non-local in the angular but not
radial coordinates.
For pratical implementations it is better to have the full explicit expression with all the integrals in spherical
coordinates, the spherical harmonics Y`m(;)and the explicit form of the Clebsch-Gordon coefficients for the
composition of an angular momentum `with a spin1
2(explicit derivation in Appendix .2),
xxxix

VPP +1
2(r;; )
1
2(r;; )
=VL(~ r;~ r0) +1
2(r;; )
1
2(r;; )
+`maxX
`=0`X
m=`
VNL
`(~ r;~ r0)
Y`m(;)R2
0d0R
0sin(0)d0Y`m(0;0) +1
2(r;0;0)
Y`m(;)R2
0d0R
0sin(0)d0Y`m(0;0) 1
2(r;0;0)!
+`maxX
`=1
VSO
„+1
2(~ r;~ r0)0
Y``(;)Z2
0d0Z
0sin(0)d0Y``(0;0) 1
2(r;0;0)
+`maxX
`=1
VSO
„+1
2(~ r;~ r0)Y„(;)
0Z2
0d0Z
0sin(0)d0Y„(0;0) +1
2(r;0;0)
+`maxX
`=1`1
2X
mj=(`1
2)
VSO
„+1
2(~ r;~ r0)0
@q
2`+1+2mj
4`+2Y`mj1
2(;)q
2`+12mj
4`+2Y`mj+1
2(;)1
A
r
2`+ 1 + 2mj
4`+ 2Z2
0d0Z
0sin(0)d0Y`mj1
2(0;0) +1
2(r;0;0)
+r
2`+ 12mj
4`+ 2Z2
0d0Z
0sin(0)d0Y`mj+1
2(0;0) 1
2(r;0;0)!
+`maxX
`=1`1
2X
mj=(`1
2)
VSO
„1
2(~ r;~ r0)0
@q
2`+12mj
4`+2Y`mj1
2(;)q
2`+1+2mj
4`+2Y`mj+1
2(;)1
A
r
2`+ 12mj
4`+ 2Z2
0d0Z
0sin(0)d0Y`mj1
2(0;0) +1
2(r;0;0)
+r
2`+ 1 + 2mj
4`+ 2Z2
0d0Z
0sin(0)d0Y`mj+1
2(0;0) 1
2(r;0;0)!
: (2.62)
From the computational point of view, the semi-local form of the pseudopotential is less efficient than the full
non local form. The procedure of Kleinman and Bylander allows the construction of a fully non-local potential
[21, 25, 26],
VKB=X
`;m
VNL
`PP
`;mED
PP
`;m
VNL
`
D
PP
`;m
VNL
`PP
`;mE; (2.63)
where PP
`;j;mj(r;; ) =RPP
`j(r)Y`m(;).
In the absence of spin terms, we first define a function, referred to as “projector” in the literature (more details
in Appendix .1). We write it between “” because it is now exactly a projector. One of the reasons for that is that
ha`mja`mi6= 1. The function associated with the “projector” is
a`m(r;; ) =1p
jb`j
VNL
`(r)R`(r)Y`m(;); (2.64)
where
b`=Z1
0r2drR`(r)
VNL
`(r)R`(r); (2.65)
and the non-local KB pseudopotential operator is
VKB (r;; ) =`maxX
`=0`X
m=`a`m(r;; )sgn(b`)Z1
0r02dr0Z2
0d0Z
0sin(0)d0a
`m(r0;0;0) (r0;0;0); (2.66)
or in shorter notation
xl

VKB=`maxX
`=0`X
m=`ja`misgn(b`)ha`mj; (2.67)
where sgn (b`)determines if the potential is attractive or repulsive. The sgn b`in between the ketand the brais
another difference from a real projector, as the one in equation (2.50). Considering spin, we have the “projectors”
a`jmj=a`jmj+1
2(r;; )
a`jmj1
2(r;; )
; (2.68)
whose definition again include the relevant Clebsch-Gordon coefficients. Notice that the a`jmjare the spinors and
a`jmjmsthe components of the spin. For the case mj=(`+1
2)we have
a„+1
2`+1
2(r;; ) =1q
jb„+1
2j
VNL
„+1
2(r)R„+1
2(r)
Y„(;)
0
a„+1
2(`+1
2)(r;; ) =1q
jb„+1
2j
VNL
„+1
2(r)R„+1
2(r)0
Y``(;)
;
for the other case we have
a„+1
2mj(r;; ) =1q
jb„+1
2j
VNL
„+1
2(r)R„+1
2(r)0
@q
2`+1+2mj
4`+2Y`mj1
2(;)q
2`+12mj
4`+2Y`mj+1
2(;)1
A
a„1
2mj(r;; ) =1q
jb„1
2j
VNL
„1
2(r)R„1
2(r)0
@q
2`+12mj
4`+2Y`mj1
2(;)q
2`+1+2mj
4`+2Y`mj+1
2(;)1
A: (2.69)
where in the Spherical Harmonics Y`mwe wrotem=mjms, and
b`j=Z1
0r2drR`j(r)
VNL
`j(r)R`j(r); (2.70)
and the non-local KB pseudopotential operator is
VKB +1
2(r;; )
1
2(r;; )
=`maxX
`=0`+1
2X
j=max(`1
2;1
2)jX
mj=ja`jmj+1
2(r;; )
a`jmj1
2(r;; )
sgn(b`)
Z1
0r2drZ2
0d0Z
0sin(0)d0(a`jmj+1
2(r;; ) +1
2(r;; ) +a`jmj1
2(r;; ) 1
2(r;; ));(2.71)
or in shorter notation
VKB=`maxX
`=0`+1
2X
j=max(`1
2;1
2)jX
mj=jja`jmjisgn(b`)ha`jmjj: (2.72)
Notice that in the KB formalism it is not possible to separate exactly the spin-orbit part from the non-spin
non-local part because R`(r)6=R`j(r), although since they are similar, such a separation should be a good
approximation. In the empirical pseudopotential method, we do not have to worry about that detail. For the non
local part we have
a`m(r;; ) =f`(r)Y`m(;); (2.73)
and for the spin-orbit part we have for the case mj=(`+1
2)
xli

a„+1
2`+1
2(r;; ) =f`;`+1
2(r)Y„(;)
0
a„+1
2(`+1
2)(r;; ) =f`;`+1
2(r)
0
Y``(;)
;
and for the other case
a„+1
2mj(r;; ) =f`;`+1
2(r)0
@q
2`+1+2mj
4`+2Y`mj1
2(;)q
2`+12mj
4`+2Y`mj+1
2(;)1
A
a„1
2mj(r;; ) =f`;`1
2(r)0
@q
2`+12mj
4`+2Y`mj1
2(;)q
2`+1+2mj
4`+2Y`mj+1
2(;):1
A (2.74)
wheref`j(r)are empirical functions. From the expression of the KB pseudopotential those functions must
behave asr`for smallr. It must also decay very rapidly for r>rc. In the recipe in reference [21], the f`jare zero
forr>rc.
The representation in equation (2.72) is different than the one that is normally used. What it is more usual is to
work with the representation ~L
~Sto the spin-orbit coupling [27].
2.12 Screening
2.12.1 Definitions of the dielectric function
The dielectric funtion (!;~k)has significant consequences for the physical properties of solids. The dipole
moment density is defined.
~D=0~E+~P; (2.75)
in S.I. units. In this section S.I. units will always be used. In this equation ~Eis the electric field, ~Pthe polarization
and0the dielectric permittivity in space. In the case of linear optics, the induced polarization depends linearly on
the electric field strength in a manner that can often be described by the relationship.
~P=0~E; (2.76)
and thus
~D=0~E+0~E=0(1 +)~E=0~E; (2.77)
with= 1 +being the dielectric constant, also known as relative permittivity. The vector ~Dis related to the
external applied charge density extin the same way as ~Eis related to the total charge density =ext+ind,
whereindis the charge density induced in the system by ext.
The divergence relation of the electric field [28] is
r
~D=r
0~E=ext (2.78)
r
~E=
0=ext+ind
0: (2.79)
The relation between the fourier components of ~Dand~Eis
xlii

~D(~k) =(~k)~E(~k); (2.80)
then we can rewrite (2.78) and (2.79).
r
~E=r
X~E(~k)ei~k
~ r=X(~k)
0ei~k
~ r(2.81)
r
~D=r
X
(~k)~E(~k)ei~k
~ r=X
ext(~k)ei~k
~ r: (2.82)
We divide (2.82) by (2.81). Each equations are satisfied term by term, so
(~k) =ext(~k)
(~k)= 1ind(~k)
(~k): (2.83)
The electrostatic potentials 'and'extdefined by ~E=r'and~D=r'extsatisfy
r2'=
0(2.84)
r2'ext =ext
0; (2.85)
so [28]
'ext(~k)
'(~k)=ext(~k)
(~k)=(~k); (2.86)
where'(~k)is the total or screened potential from the external potential 'ext(~k)
'(~k) ='ext(~k)
(~k); (2.87)
and we need now to find the expression for (~k). In the electron gas, in the limit ~k!0,(!;0)describes the
collective excitations of the Fermi sea. In the other limit !!0,"(0;~k)describes the electrostatic screening of the
electron-electron and electron-lattice interactions in crystals [28]
2.12.2 Screening in a metal
In a metal we have an uniform gas of electrons of charge concentration n0ein a background of positive
chargen0e. The dielectric function for real systems is very complicated. Here we will only derive and work with
an approximation. A derivation from reference [28] will be followed. We consider a sinusoidal variation of positive
charge density in the xdirection
+(x) =n0e+ext(k) sin(kx); (2.88)
whereext(k) sin(kx)gives rise to an electrostatic field, the external field applied to the electron gas. We know the
relations between the electric field, the potential and the charge density ~E=r'andr
~E=
0, which gives
rise to the Poisson equation r2'=
0. For the positive charge we have
'='ext(k) sin(kx);  =extsin(kx); (2.89)
and using the Poisson equation
xliii

k2'ext(k) sin(kx) =ext(k)
0sin(kx),k2'ext=ext(k)
0: (2.90)
The electron gas will be subjected to this external potential 'extand to the induced electrostatic potential
'ind(k) sin(kx)of the deformation of the electron gas itself. And the electron charge density is
(x) =n0e+ind(k) sin(kx): (2.91)
The amplitude of the total electrostatic potential '(k) ='ext(k) +'ind(k)is related to the total charge density
(k) =ext(k) +ind(k)by the Poisson equation
k2'ext=ext(k)
0: (2.92)
In the Thomas-Fermi approximation is assumed that a local internal chemical potential can be defined as a
function of the electron concentration at that point. In a reagion where there is no electrostatic contribution:
=0
F=~2
2m(32n0)2
3; (2.93)
at absolute zero, with Fbeing the Fermi energy. In a region where there is electrostatic potential '
=F(x)e'(x)'~2
2m[32n(x)]2
3e'(x)=~2
2m[32n0]2
3: (2.94)
This results in the local value of the Fermi level being
F(x) =e'(x) +0
F; (2.95)
and by a Taylor series expansion of Fas a function of the electron concentration nwe have:
F(n(x)) =0
F+dF
dn
n=n0[n(x)n0]: (2.96)
Using equations (2.95) and (2.96) together:
dF
dn0[n(x)n0] =e'(x): (2.97)
From (2.93) we can calculate the derivatived"F
dn0=2
3"F
n0and thus:
2
3F
n0[n(x)n0] =e'(x) =e'(x),n(x)n0=3
2n0e'(x)
F; (2.98)
beingn(x)n0the induced part of the electron concentration and thus
ind(k) =e[n(x)n0] =3
2n0e2
F'(k) =3
2n0e2
F(k)
k20: (2.99)
From equation (2.83) we discover that
(k) = 1ind(k)
(k)= 1 +k2
TF
k2withk2
TF=3
2n0e2
"F0: (2.100)
xliv

2.13 Optical properties
To calculate the optical properties such as Reflectivity, we first compute the imaginary part of the dielectric
function in the limit ~k!0, from which we can get the real part from a Kramers-Kroning transformation. This
transformation relates the real and imaginary part of the generalised susceptibility, , that is related to the dielectric
function by = 1 +, and therefore
1= 1 +1;
2=2; (2.101)
where the index 1 denotes the real part and 2 the imaginary part. By using the Kramers-Kroning transformation
we get the real part of the susceptibility from the imaginary part [29]
1(!) =2
PZ+1
02()
2!2d; (2.102)
wherePdenotes the principal part of the integral. The real part of the dielectric function is therefore calculated
from the imaginary part 2as
1(!) = 1 +2
PZ+1
02()
2!2d: (2.103)
The optical measurements that give the fullest information on the electronic system are measurements of the
reflectivity of light at normal incidence on single crystals [28]. The reflectivity coefficient is given by:
r(!) =Erefl
Einc(2.104)
whereEincis the incident electric field and Ereflis the reflected electric field, as its subscripts indicate. By
definition, the refractive index n(!)and the extinction coefficient K(!)are related to the dielectric function (!) =
1(!) +i2(!)as [28]
p
(!)n(!) +iK(!) (2.105)
If the incident traveling wave with wavevector kis traveling in the xdirection, then the ycomponent is
Ey
inc=Aei(kx!t); (2.106)
and the reflected wave in vacuum can be written as
Ey
refl=Bz
refl=A0ei(kx!t); (2.107)
For the transmited wave in the dielectric medium we find
Ey
trans =ckBz
trans
!=1=2Bz
trans =A00ei(kx!t): (2.108)
from the Maxwell equation cr~H=@~E
@tand the dispersion relation !2=c2k2for electromagnetic waves.
The boundary conditions at the border ( x= 0) are such that both EyandBzshould be continuous:
Ey
inc+Ey
refl=Ey
trans
By
inc+By
refl=By
trans,AA0=A00
A+A0=1=2A00,A0
A=11=2
1=2+ 1; (2.109)
xlv

and
r(!) =Erefl
Einc=A0
A=1=21
1=2+ 1=n(!) +iK(!)1
n(!) +iK(!) + 1: (2.110)
But the quantity that is actually measured in the experiments in the reflectance R, the ratio of the reflected
intensity to the incident intensity:
R=E
reflErefl
E
incEinc=rr=n(!)iK(!)1
n(!)iK(!) + 1n(!) +iK(!)1
n(!) +iK(!) + 1
=(n(!)1)2+K2(!)
(n(!) + 1)2+K2(!)(2.111)
2.14 Imaginary part of the dielectric function 2
To describe the interaction between an external electromagnetic field and Block electrons in a semiconductor
crystal we will use s semi-classical approach. This means that the field is treated classically but the electrons
are described using quantum mechanical wave functions. In this chapter we will use S.I. units. The unperturbed
Hamoltonian is
H0=p2
2m+V(~ r): (2.112)
Using the Coulomb gauge invariance we say that the scalar potential  = 0 and that the vector potential ~A(~ r;t)
behaves asr
~A= 0. In this Gauge, the electric and magnetic fields are given by
~E=1
c@A
@tand~B=r~A: (2.113)
We obtain the quantum mechanical Hamiltonian describing the motion of an electron in an external electro-
magnetic field replacing ~ p!~ p+ (e~A
c)as
H=1
2m"
~ p+
e~A
c!#2
+V(~ r): (2.114)
The first term can be expanded as
1
2m
~ p+e~A
c!2
=p2
2m+e
2mc~A
~ p+e
2mc~ p
~A+e2A2
2mc2: (2.115)
Knowing that ~ p=~
irwe calculate the action of ~ p
~Ain a function f(r)
(~ p
~A)f(r) =~A
(~
irf) +~
ir
~A
f: (2.116)
From the Coulomb Gauge, r
~A= 0, thereforee
2mc

~ p
~A=e
2mc
~A
~ p. We say that the amplitude of the field
j~Ajis very small, so we neglect the terme2A2
2mc2, which depends quadratically on the field. So the Hamiltonian will
approximately be given by [30]
H=H0+e
mc~A
~ p=H0+HeR; (2.117)
whereHeRis the electron-radiation interaction Hamiltonian. If we consider to the field to be small we can treat this
as a perturbation, so we can apply the time-dependent perturbation theory. We will use the Fermi Golden Rule
to calculate the transition probability per unit volume Rfor an electron in the valence band state~kv;nvE
to the
conduction band~kc;ncE
, [30]
xlvi

R=2
~X
~kc;~kvD
~kc;ncHeR~kv;nvE2
(Ec(~kc)Ev(~kv)~!): (2.118)
We need first to evaluate first the matrix element
D
~kc;ncHeR~kv;nvE2
=e
mc2D
~kc;nc~A
~ p~kv;nvE2
: (2.119)
We will say that ~A=A~ e,in which~ eis an unit vector in the direction of ~A,
D
~kc;ncHeR~kv;nvE2
=e
mc2
jAj2D
~kc;nc~ e
~ p~kv;nvE2
; (2.120)
In Appendix .3 we conclude that for a transition to happen we have that ~kc=~kv=~k,
D
~kc;ncHeR~kv;nvE2
=D
~k;ncHeR~k;nvE2
=e
mc2
jAj2D
~k;nc~ e
~ p~k;nvE2
=e
mc2
jAj2jPcvj2;(2.121)
where~Pcvis the momentum matrix operator and we write Ain terms of the Amplitude of the incident electric field
E=j~E(~ q;!)j, [30]
A=E
2qei(~ q
~ r!t): (2.122)
The eletric dipole transisiton probability Ris then given by the Fermi Golden Rule in equation (2.118), [30]
R=2
~e
m!2E(!)
22X
kjPcvj2(Ec(~k)Ev(~k)~!): (2.123)
The power loss by the field due to absorption in unit volume is given by R~!and is equal to the rate of decrease
in the energy of the incident beam per unit volume, [30]
dI
dt=dI
dxdx
dt
=Ic
n=2!I
n2; (2.124)
whereis the absorption coefficient and nthe refractive index of the medium. Iis related to the field amplitude by
[30]
I=n2
8jE(!)j2: (2.125)
Doing!=cq, we obtain
2(!) =2e
m!2X
~kjPcvj2(Ec(~k)Ev(~k)~!) (2.126)
xlvii

3
Results
Contents
3.1 Silicon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . l
3.2 Carbon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lix
3.3 Germanium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lxiv
xlviii

The work that has been developed was the development of a pseudopotential method to calculate the band
structure of Si, Ge and C, all group IV elements. Since they have the same number of valence electrons, using the
pseudopotential method, the description for this elements is similar, except for Ge, that is an heavier element, we
have to consider spin-orbit effects. All the code is in the web at https://fenix.tecnico.ulisboa.pt/homepage/
ist167936/dissertacao-de-mestrado
3.0.1 Simple test – free electron
The solid has the diamond crystal structure, and thus the primitive lattice vectors are those for the fcc lattice
~ a1=
01
21
2

a ~ a 2=1
201
2

a ~ a 3=1
21
20

a; (3.1)
and the reciprocal primitive vectors are the primitive vectors for an bcc lattice
~b1=1 1 1
2
a~b2=11 1
2
a~b3=1 11
2
a: (3.2)
The whole crystal has the size N1~ a1N2~ a2N3~ a3. The lattice constant used is a= 10:26a:u:which is the
same for Silicon. To solve the secular equation (2.54) we use a plane wave basis
D
~ r~ki1i2i3;~GI1I2I3E
=1pN1N2N3Vcellei(~ki1i2i3+~GI1I2I3)
~ r(3.3)
in which
~ki1i2i3=i1
N1~b1+i2
N2~b2+i3
N3~b3 (3.4)
withij= 0;:::;Nj1and
~G=I1~b1+I2~b2+I3~b3 (3.5)
is a reciprocal lattice vector to calculate the matrix elements of the kinetic energy operator. As it is demonstrated
in Appendix .3.1 this basis functions are orthogonal among them. This means that the kinetic energy operator will
be a diagonal matrix
1
2r2~k;~GE
=1
2j~k+~Gj2~k;~GE
(3.6)
We start with the band structure calculation of a bulk element with coordination number 4 and consider that the
potential that acts on the valence electrons is V= 0. This means that the Hamiltonian that describes this system
has only the kinetic energy. Figure 3.1 a) shows the band structure for a free electron in a fcc Brillouin zone, where
Figure 3.1: The figure shows a) the free electron bands in the fcc lattice and b) the Brillouin zone for the fcc lattice
xlix

we see a folded parabola, as expected for free electrons. In the Figure is represented E(~k)in the Brillouin zone of
Figure 3.1 b).
3.1 Silicon
3.1.1 Description only with a local pseudopotential
We start by choosing a simple analytical form for the unscreened local pseudopotential. The family of pseu-
dopotentials is
Vlocal(r) =4
rerfr
Ra
+16
q2z(pRa)3er2
R2a: (3.7)
This expression is chosen empirically. The 4
r(Figure 3.2, blue line) is the Coulomb potential of the nucleus
with the core electrons for all group IV atoms. The erf
r
Ra
(Figure 3.2, yellow line) term, when multiplied by the
Coulomb term, it smooths the potential in the origin, so it doesn’t go to infinity as r!0(Figure 3.2, green line)
and erf
r
Ra
= 1forr>Ra. The value of4
rerf
r
Ra
at the origin is8pRa. We want that the potential to go be
well behaved when r!+1, so we add and exponential to make this effect (Figure 3.2, red line), multiplied by an
apparently strange prefactor that will be understood later. The result is the purple line of Figure 3.2.
Figure 3.2: The graphic shows the functions that compose the local pseudopotential and the pseudopotential itself, unscreened
and screened
Performing a Fourier transform of this function we get with the local pseudopotential in the reciprocal space,
Vlocal(k) =Z2
0Z
0Z+1
0Vlocal(k)j0(kr)r2sin()drdd =
= 4Z+1
0Vlocal(k)j0(kr)r2dr=
= 161
q2z1
k2
ek2R2
a
4: (3.8)
We understand now that the parameter qzis when we have Vlocal(qz) = 0 . After we divide by equation (2.100)
to obtain the screened potential that we next are going to use in the secular equation,
Vlocal;screen (k) =Vlocal(k)
(k)= 16k2
q2z1ek2R2
a
4
k2+k2
TF; (3.9)
where the parameters Ra,qzandkTFare the parameters to be adjusted to the experiment. We know that the
equation (2.100) is used in metals but we are going to use it here with the semiconductors Si and Ge and the
insulator diamond-C because it has only 1 parameter that we are going to control and adjust. Since we are using
kTFas an empirical parameter, the form is not important. What we achieved is that when k!0,Vlocal;screen (k)is
well behaved. Notice that the potential has spherical symmetry, so it will only depend on j~kj. In practice we actually
l

RakTFqz
0.9 0.35 1.48
Table 3.1: The table shows previously existent parameters for Silicon, that we first check in this research.
just calculate the Fourier component Vlocal;screen(j~Gj), where and ~Ga vector in the reciprocal lattice. In Figure
3.2, the orange line, Vlocal;screen(r)is shown.
So now we add to the program discussed in section 3.0.1 the local potential operator to be included in the
calculations. In Appendix .3.2 we prove that, using a plane wave basis set, the non-zero matrix elements are those
with the same vector ~k, but we also prove the selection rule ~G00=~G0~Gand calculate the analytical expression
for the matrix elements of the local operator.
First we check previously obtained values, written in Table 3.1
First we multiply the pseudopotential function by a factor fwith0f1. We already saw the result to f= 0
in the previous section. The purpose is to observe the evolution of the effects of turning on the potential in the
band structure. The results are in Table 3.2
f= 0 f= 0:1 f= 0:3
f= 0:5 f= 0:75 f= 1
Table 3.2: Resulting band structure subjected to the potential (3.9) of Silicon, using the parameters from Table 3.1 multiplied
by a factorfwith0f1. The the opening of the gap is clearly shown
With the increasing fwe can observe the lifting of the degeneracies in Land, as well as the translation of
the bands to higher energies.
The band structure obtained with f= 1is at first sight quite similar to the one obtained by reference [5] that
we can see in Figure 3.3. This graphic has noted the experimental values to some transitions in bulk Silicon from
Landoldt-Bornstein. In the calculated band gap with [5] does not have the same value as the experimental results,
so the Figure is not at scale. We compare them with the ones obtained here, displayed in Table 3.3. We can see
that some values are similar but others are quite far and we see in the band structure that the second set of p
bands is higher than the second sband, which should be reversed, fact that is confirmed by experience. We can
T0T1T2T3T4T5T6
Experience 12.5 1.2 2.04 3.35 4.15 1.17 2.9
Calculated 12.9 1.3 1.42 3.19 2.39 1.07 3.0
Table 3.3: Experimental and calculated in the current work transitions of Silicon in eVare calculated with the parameters from
Table 3.1 in equation (3.9)
li

Figure 3.3: The figure shows the LDA band structure of bulk Silicon calculated with the program of reference [5]. Experimental
values are indicated by the double arrows.
confirm this is happening by looking at the values of transitions T3andT4.
In the first programs, we first calculated the energy values only for the ~kpoints that would give us the band
structure. But now we want to calculate the density of states and the optical properties, which requires the calcu-
lation of the eigenvalues and vectors on a cubic grid. We calculated a grid of 111111points. First we calculate
the density of states,
D(E) = lim
N1;N2;N3!11
N1N2N3X
nX
i1X
i2X
i3(EEn(~ki1;i2;i3)); (3.10)
where
(E) = lim
!0D(E) = lim
!01
p
2e1
2(E
)2
: (3.11)
The joint density of states is also calculated
J(E) = lim
N1;N2;N3!11
N1N2N34X
vX
cX
i1X
i2X
i3[E(Ec(~ki1;i2;i3)Ev(~ki1;i2;i3))]; (3.12)
where the indices vandcare correspondent to the valence and conduction bands, respectively. This is a measure
of the transitions between energy levels. The results for the density of states and joint density of states are in
Figure 3.4. We observe a noise from using the Gaussian function. To minimize this effect we would have first to
Figure 3.4: The Figure shows the a) calculated density of states (blue line), the photo emission spectroscopy and inverse
photo emission data obtained from reference [6] (yellow line) and b) the calculated joint density of states with the parameters
from Table 3.1, for Silicon
increaseNito decrease after. This is not very practical when we use a software like M ATHEMATICA , because
it increased a lot the number of calculations, increasing as well the computation time (see Chapter 4 for DOS
lii

without noise). In Figure 3.4 a), blue line, we can see clearly the band gap after 0eV. The joint density of states is
linked to the optical transistors, calculating only the direct electronic transitions. We can see that the graphic starts
to be different from zero around T3= 3:19eV, which is the direct transition of lowest energy. T5,T1andT2are
not described here. We can compare with the results obtained by reference [6] in Figure 3.4 a), yellow line, and
notice that there is no noise from the Gaussian use, as before, but the peaks are approximately in the same place,
although they have different intensities. The different intensities for E < 0can be explained with the existence of
background noise to the lower energies It is also calculated the area under the curve from 13eVto0eVthat
from reference it is known that it should be 4, that is what it is obtained.
Next we calculate the imaginary part of the dielectric function 2(!)using the expression (2.126). The mo-
mentum matrix operator jPcvj2=D
~k;c~ p~k;vE2
, with~ p=i~ris calculated using the fact that the normalized
wave-functions are a linear combination of the plane wave basis functions (equation (3.3)), with the coefficients
that come out of the eigenvectors in the diagonalization of the Hamiltonian matrix operator.
i1i2i3(~ r) =1pN1N2N3VcellX
I1;I2;I3cI1I2I3ei(~ki1i2i3+~GI1I2I3)
~ r(3.13)
With this we calculate
h i1i2i3j~rj i1i2i3i=1
N1N2N3VcellZZZ
crystalX
I1;I2;I3c
I1I2I3ei(~ki1i2i3+~GI1I2I3)
~ r~r0
@X
I0
1;I0
2;I0
3cI0
1I0
2I0
3ei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ r1
Ad3r=
=1
N1N2N3VcellX
I1;I2;I3X
I0
1;I0
2;I0
3c
I1I2I3i(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)cI0
1I0
2I0
3
ZZZ
crystalei(~ki1i2i3+~GI1I2I3)
~ rei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ rd3r=
=1
N1N2N3VcellX
I1;I2;I3X
I0
1;I0
2;I0
3c
I1I2I3i(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)cI0
1I0
2I0
3
N1N2N3Vcelli1;i0
1i2;i0
2i3;i0
3I1;I0
1I2;I0
2I3;I0
3=
=X
I1;I2;I3jcI1I2I3j2i(~ki1i2i3+~GI1I2I3): (3.14)
Using equation (2.103) we calculate the real part of the dielectric function 1(!). The results for 1and2are
shown in Figure 3.5
Figure 3.5: Real part1and imaginary part 2of the dielectric function of Silicon are calculated using the parameters on Table
3.1
We can compare this results to the ones obtained in reference [7] in Figure 3.6. We can see that the peaks
are mostly in the same place but the differences arise because excitonic effects in the calculations. We can also
liii

see that the imaginary "2should start around 3eVbut what was calculated starts around 2:5eV, indicating that the
highest valence band and lowest conduction band are less separated in than they should be.
Figure 3.6: Comparing the calculated dielectric function of Silicon, using the parameters on Table 3.1, with experimental
results in reference [7]
From equation 2.105 we can obtain
1(!) =n2(!)K2(!)2(!) = 2n(!)K(!); (3.15)
and calculate
K(!) =1p
2=r
1(!) +q
2
1(!) +2
2(!)n(!) ="2(!)
2K(!): (3.16)
The reflectivity coefficient r(!)and reflectance R(!)can be obtained from equations (2.110) and (2.111),
respectively. In Figure 3.7 we can see the result of the calculation of R(!)and the experimental Reflectance form
reference [7]. The figure agrees with experiment except for the curve between 2 and 3eVin a).
Figure 3.7: The figure shows the calculated reflectance for Silicon using Table 3.1 (blue line) and the experimentally obtainced
from reference [7] (yellow line)
3.1.2 Description with non-local pseudopotential
The results in the above section could be more satisfactory if we had a more “physical” description of Silicon.
That is to say to differentiate the effect of the potential applied in the selectrons from the one applied in the p
electrons. This means we are going to consider the non-local part of the pseudopotential. In our case we are
going to use `max= 0. This means that the non-local part is going to act on the selectrons while the local part
liv

affects both sandpelectrons. The matrix elements are calculated using the result in Appendix .3.3. Since the
Empirical Pseudopotential Method is used, the projectors a`m(r;; )are calculated using expression (2.73).
For the function f`(r), consistent with the local form of the pseudopotential and satisfying the desired conditions,
we are going choose again a Gaussian multiplied by rl,
f`(r) =B`r`er2
R2
b; (3.17)
whereB`andRbare constants to be adjusted. The reason for rlcomes from equation 2.64, when r!0,
R`(r)!rl. Expression (3.17) has an analytical Fourier transform,
F`(k) =Z+1
0r2j`(kr)B`r`er2
R2
bdr=R`
bB`pR3
b
4Rbr
2`
eR2
bk2
4; (3.18)
This time we are going to search not only for the parameters Ra,qzandkTFto the local potential from equation
(3.7) but also the parameters RbandB`to the non-local potential from equation (3.17). We are going to start for
scratch, which means we are lost in a 5-dimensional space of 5 parameters to adjust, without knowing where to
start with. So in the beginning, we are going to be based in the pseudopotentials generated with the program from
reference [1] and fit the functions to the pseudopotentials to obtain the first parameters we will work with. Since we
have`max= 0, we will use only Y00=1p
2in equation (2.73) and the projector a00is
a00(r) =1p
2B0er2
R2
b; (3.19)
and its 3DFourier transform
A00(k) =1
4Ble1
4k2R2
bpR3
b: (3.20)
We start by fitting the local part of the potential at the ppseudopotential of reference [1] to equation (3.9) (Figure
3.8). The purpose is to start to discover the acceptable region for the values of the parameters to be adjusted. The
results are on Table 3.4.
Figure 3.8: The figure shows the fit of the expression (3.9) (line) to the ppseudopotential of Silicon generated with the program
in reference [1] (dots)
RaqzkTF
0.93 2.17 0.65
Table 3.4: Results to the fit of equation (3.9) to the ppseudopotential of Silicon, generated with the program in reference [1],
using the default weight function chosen by M ATHEMATICA
From here we fix the result to qzdefinitely and work with the rest of the parameters. Notice from equation (3.9)
thatqzis where we have Vlocal(k) = 0 . Next we calculate the quantity
lv

minimization (Ra;kTF) =Z
(PPgenp(k)Vlocal;screen (k;Ra;kTF))2k2dk (3.21)
where PPgenp(k)is theppart of the pseudopotential generated with [1] and Vlocal;screen (k;Ra;kTF)is equation
(3.9), where we change the parameters RaandkTFto draw a ballpark figure, as we can see in Figure 3.9. The
green point corresponds to the result of the fit of Table 3.4.
Figure 3.9: The ballpark figure to fit the local part of the Silicon pseudopotential using function (3.21) is shown. The green
point is the result of the fit, in Table 3.4
In the next step we solve the Schrodinger equation and find out the pairs of values ( Ra,kTF) for which the
eigenvalue of energy is the one of the Epfor bulk Silicon. This value is Ep=4:16eVand comes from reference
[1]. It is also calculated the pairs for which the eigenvalue of energy is Ep0:5eV. In Figure 3.10 all these results
are condensated. The green line is a fit of the pairs (Ra;kTF)for whichEp=4:16eVto a linear function.
Figure 3.10: The ballpark figure to fit the local part of the Silicon pseudopotential using function (3.21) is shown. The green
point is the result of the fit, in Table 3.4, the red points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy
is the one of the Ep, the green line is adjusted to these points, the orange points are the pairs of values ( Ra,kTF) for which
the eigenvalue of energy is the one of the Ep+ 0:5eVand the yellow points are the pairs of values ( Ra,kTF) for which the
eigenvalue of energy is the one of the Ep0:5eV.Ep=4:16eVis obtained from reference [1] for Silicon
In another step, we do last adjustments to the local potential by putting it in a self consistent Schrodinger
equation, adjusting the eigenvalue of energy to Es=10:83eV(value from reference [1]),
lvi


1
2r2+Vlocal;screen (r)
i+1+rAKBs(r)Z
AKBs(r) i(r)dr=Es i+1; (3.22)
whereAKBs(r)is thesprojector generated with [1]. This results with the values of Table 3.19. We use this results
to fit the non-local part of the pseudopotential.
RaqzkTF
0.93 2.17 0.54
Table 3.5: The results to the fitting, using (3.22), to the ppseudopotential of Silicon generated with the program in reference
[1] are shown
Next we fit equation (3.19) to the spseudopotential generated with [1] to obtain the first parameters RbandB0
to start with. The graphic is in Figure 3.11 and the results in Table 3.6
Figure 3.11: A fit that was made (line), of (3.20) to the s“projector” of Silicon, generated with [1] (points)
RbB0
0.85 10.2
Table 3.6: The obtained parameters of the fitting of (3.20) to the s“projector” generated with [1] are shown
The ballpark figure is draw with the same principle as before and as we did previously,
minimization (Ra;kTF) =Z
(AKBs(r)A00(k;Rb;B))2k2dk; (3.23)
and we calculate the pairs (Rb;B)for which the eigenvalue of energy of the Schrodinger equation is Es=10:83eV
– and we saw that the value correspondent to Rb= 0:85isB= 10:2- and also Es0:5eV. The results are in
Figure 3.12. The yellow line is a quadratic tendency line.
Finally we set the parameters for the local and non-local pseudopotentials (Table 3.7) and with this values
RaqzkTFRbB0
0.93 2.17 0.54 0.85 10.2
Table 3.7: The obtained final parameters, used to calculate the important energetic transitions of Silicon are shown
we calculate the transitions of Table 3.3. Also we calculate the same values adding each of the parameters and
incrementh1= 0:01forRa,qz,kTFandRbandh2= 0:1forB0to understand the consequences of changing
each of the parameters. The results are in Table 3.8. We can still see that transitions T3andT4are still reversed.
From now we can manually change the values of the parameters to better adjust to the experimental values of
this transitions. We take into account that changing RaorkTFwe seepbands changing (fundamentally lowering
or rising), since the local potential describes the ppotential for Silicon and changing RbandB0will change the s
lvii

Figure 3.12: The ballpark figure to fit the non-local part of the pseudopotential of Silicon using function (3.23) is shown. The
green point is the result of the fit, in Table 3.6, the red points are the pairs of values ( Rb,B) for which the eigenvalue of energy is
the one of the Es, the yellow parabola is adjusted to these points, the orange points are the pairs of values ( Rb,B) for which the
eigenvalue of energy is the one of the Es+ 0:5eVand the yellow points are the pairs of values ( Rb,B) for which the eigenvalue
of energy is the one of the Es0:5eV.Es=10:83eVis obtained from reference [1]
Transition value 0Ra+ 0:01ktf+ 0:01Rb+ 0:01B+ 0:1
T0= 12:5 11.4367 11.4469 11.4533 11.2423 11.3854
T1= 1:2 1.24625 1.25111 1.25243 1.24625 1.24625
T2= 2:04 1.85174 1.85578 1.83495 2.10205 1.92583
T3= 3:35 3.22261 3.20447 3.19349 3.22261 3.22261
T4= 4:15 2.70572 2.74958 2.72366 3.16722 2.84453
T5= 1:17 1.53095 1.49858 1.47963 1.60367 1.55232
T6= 1:07 2.81646 2.83021 2.83327 2.81646 2.81646
Table 3.8: Important energy transitions of Silicon where calculated with the parameters in Table 3.7
bands since the non-local part of the potential describes the spotential for Silicon. We see on Table 3.8 the effects
of changing each of the parameters to better guide us.
After the searching, the final band structure for Silicon is on Figure 3.13, the final values of the parameters on
Table 3.9 the values of the transitions on Table 3.10. The shape of the band structure is very similar to the one
Figure 3.13: Band structure of Silicon was a) calculated using LDA, from [5] and b) calculated using the parameters on Table
3.9
calculated with LDA. We can see that the pandsbands that where exchanged have now taken their places, as
lviii

RaqzkTFRbB0
0.972 2.17 0.62 1.06 6.1
Table 3.9: The final parameters for the pseudopotential of Silicon where obtained after adjusting to the experiment
T0T1T2T3T4T5T6
Experience 12.5 1.2 2.04 3.35 4.15 1.17 2.9
Calculated 10.7 1.3 2.29 3.91 3.83 1.17 3.0
Table 3.10: Experimental and calculated in the current work transitions of Silicon in eVare calculated with the parameters on
Table 3.9
we can see in the values of the transitions T3andT4. The value for T5is well fitted and is the most important
since is the one who gives the value of the indirect gap of Silicon. We now calculate the same properties as we
did before in the first description of Silicon, namely the density of states, D(E), the joint density of states, J(E),
the real and imaginary parts of the dielectric function, 1(!)and2(!)and the reflectance R(!). We calculated a
grid of 999points. All this results are Figure 3.14. In the real part of the dielectric function we see once again
that the peak is in the same place but the difference come for not having considered the excitonic effects. In the
imaginary part, the peaks are dislocated to lower energies and the function starts growing at a value of !to high
than it was supposed to. The reflectance agrees with experiment except for the curve between 2and4eVin the
calculatedR(!)and the first peak is in a 1eVhigher energy than it should.
3.2 Carbon
With Carbon we use both a local and non-local pseudopotential and use again lmax = 0. We fit the local
pseudopotential generated from reference [1] to the expression 3.9 to obtain the parameters RaandkTF. This
time, theqzparameter was directly obtained from solving Vlocal;screen (k) = 0 , from the pseudopotential generated
from the program [1]. The graphic of the fit is in Figure 3.15 and the results are in Table 3.11. We can see in the
Figure that this fit is not very good because the potential described by the fitted function has higher value in almost
all the points than the one gave by reference [1].
RaqzkTF
0.19 5.69 0.76
Table 3.11: The table shows the results to the fit of equation (3.9) to the ppseudopotential of Carbon generated with the
program in reference [1]
In the same way with the previous section, a ballpark figure was obtained with the computing of expression
(3.21) and we solve the Schrodinger equation and find out the pairs of values ( Ra,kTF) for which the eigenvalue
of energy is the one of the Epfor bulk Carbon. This value is Ep=5:41eVand comes from reference [1]. It is
also calculated the pairs for which the eigenvalue of energy is Ep0:5eV. In Figure 3.16 all these results are
condensate. The yellow line is a fit of the pairs (Ra;kTF)for whichEp=5:41eVto a quadratic function. The
green point corresponds to the result of the fit.
The last adjustments to the local potential using (3.22), used to adjust the non-local part, are in Table 3.12 and
the results of the fit of the non-local pseudopotential are in Figure 3.17 and Table 3.13.
RaqzkTF
0.19 5.69 0.70
Table 3.12: The results to the fitting, using (3.22), to the ppseudopotential of Carbon, generated with the program in reference
[1] are shown
lix

Figure 3.14: It is represented a) the calculated density if states (blue line),the photo emission spectroscopy and inverse photo
emission data obtained from reference [6] (yellow line), b) the calculated joint density of states, c) calculated dielectric function,
d) calculated (blue line) and experimental (yellow line, [7]) 1, e) calculated (blue line) and experimental (yellow line, [7]) 2f),g)
calculated (blue line) and experimental (yellow line, [7]) reflectance for Silicon with the local pseudopotential of equation (3.7)
and non-local projector for the pseudopotential of equation (3.19) with the parameters written in Table 3.9
lx

Figure 3.15: The figure shows the fit of the expression (3.9) to the ppseudopotential of Carbon generated with the program in
reference [1]
Figure 3.16: The ballpark figure to fit the local part of the pseudopotential using function (3.21) is shown. The green point is
the result of the fit, in Table 3.11, the red points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy is the one
of theEp, the green line is adjusted to these points, the orange points are the pairs of values ( Ra,kTF) for which the eigenvalue
of energy is the one of the Ep+ 0:5eVand the yellow points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy
is the one of the Ep0:5eV.Ep=5:41eVis obtained from reference [1] for Carbon
Figure 3.17: The fit of the expression (3.20) to the s“projector” of Carbon, generated with the program in reference [1] is
represented
lxi

RbB0
0.54 47.0
Table 3.13: Results to the fit of equation (3.20) to the s“projector” of Carbon generated with the program in reference [1] are
shown
The ballpark figure is draw with the same principle as before using (3.23). We calculate the pairs (Rb;B)
for which the eigenvalue of energy of the Schrodinger equation is Es=13:63- and we saw that the value
correspondent to Rb= 0:54isB= 46:5- and alsoEs0:5eV. The results are condensated in Figure 3.18. The
red line is a quadratic tendency line.
Figure 3.18: The ballpark figure to fit the non-local part of the pseudopotential of Carbon using function (3.23) is shown. The
green point is the result of the fit, in Table 3.13, the red points are the pairs of values ( Rb,B) for which the eigenvalue of energy
is the one of the Es, the red parabola is adjusted to these points, the orange points are the pairs of values ( Rb,B) for which the
eigenvalue of energy is the one of the Es+ 0:5eVand the yellow points are the pairs of values ( Rb,B) for which the eigenvalue
of energy is the one of the Es0:5eV.Es=13:63eVis obtained from reference [1]
Finally we set the parameters for the local and non-local pseudopotentials (Table 3.14).
RaqzkTFRbB0
0.19 5.69 0.70 0.54 46.5
Table 3.14: These are the obtained final parameters, used to calculate the important energetic transitions of Carbon
In Figure 3.19 a) we see the band structure of Carbon, generated by reference [1] with some experimental
measured and theoretical predicted energetic transitions on it. This values where obtained from GW values, from
the LB (Landoldt-Bornstein), from PRB22 and PRB87 (Physical Reviews B), and Himpsel. In Table 3.15 we use
the values from Table 3.14 to calculate the important transisitons of Carbon, as we did in the previous section for
Silicon in Table 3.8. For some transitions we had two informations about its value, one experimental (PRB22) and
other theoretical (GW).
Using the same logic of the previous section, we adjust further the parameters of Table 3.14 until the values of
the transitions and the shape of the band structure is the most satisfactory as possible. The final parameters are
on Table 3.16 and the calculated band structure in Figure 3.19 b). The important transitions are on Table 3.17.
Once again we gave more importance to the transition of the indirect gap T10but we noticed that by better
lxii

Transitions 0Ra+0.01kTF+0.01Rb+0.01B+0.1
T0=23. 21.5054 21.5642 21.5587 20.4208 21.4666
T1=17.3 15.8529 15.8996 15.884 14.6045 15.8073
T2=14.4 13.5541 13.5996 13.6105 13.2196 13.541
T3=2.9 2.87054 2.87965 2.88081 2.87054 2.87054
T4=10.3 9.17417 9.15452 9.12433 9.17417 9.17417
T5=7.3 6.30523 6.2949 6.26954 6.30523 6.30523
T6=14.4 13.6916 13.7382 13.7046 16.729 13.8165
T7=14. 12.9692 13.0095 13.0017 12.0654 12.9347
T8=7.0 6.4844 6.50641 6.51223 6.4844 6.4844
T9=6.3 4.95551 4.90878 4.87322 5.68798 4.98377
T10=5.48 4.40137 4.36195 4.33045 5.09497 4.42815
Table 3.15: This important energy transitions of Carbon were calculated with the parameters in Table 3.14
RaqzkTFRbB0
0.21 5.73 0.70 0.555 47.04
Table 3.16: The discovered final parameters for Carbon, after the adjustment to the experiment are shown here
Figure 3.19: Band structure of diamond a) form reference [5] and b) calculated with the parameters in Table 3.16
T0T1T2T3T4T5T6T7T8T9T10
Experience 23.0 17.3 14.4 2.9 10.3 7.3 14.4 14.0 7.0 6.3 5.48
Calculated 19.8 13.9 13.1 2.9 9.1 6.3 18.8 11.6 6.5 6.1 5.48
Table 3.17: Experimental and calculated in the current work transitions of Diamond in eV, calculated with the parameters on
Table 3.16
adjusting this value, the others started to be worse adjusted than they where previously. The gap in Lwas poorly
opened. Is was also calculated D(E),J(E),1(!),2(!)andR(!). We calculated a grid of 999points. The
results are in Figure 3.20. We see in the calculated DOS that the peak in Figure 3.20 a) is dislocated to higher
energies compared the PE data. The calculations are not in accordance with the experience. In the real part
of the dielectric function, the peaks about in the same place but they have different intensities. The peak of the
calculated imaginary part of the dielectric function in is in higher energies than the experimentally obtained. The
same happens for the Reflectance.
lxiii

Figure 3.20: It is represented the a) DOS of Carbon, calculated here (blue line), the photo emission spectroscopy data from
reference [8] (yellow line) divided by a factor of 20, b) the calculated joint density of states, the c) real part (blue line) and
imaginary part (purple line) of the dielectric function, d) the comparison between the calculated (blue) and experimentally
obtained (yellow, [9]) 1, e) comparison between the calculated (blue) and experimental (yellow, [9] 2and f) the calculated
(blue) and experimentally obtained (yellow, [9]), divided by a factor of 100, reflectance.
3.3 Germanium
For Germanium we use not only the local part of the pseudopotential (3.7), the non-local part with `max= 0
and projector (3.19) but also the spin-orbit contribution, because Germanium is an heavier element (bigger atomic
numberZ), the spin-orbit contribution for the Hamiltonian
HSO=Ze2
2m2ec21
r3~S
~L; (3.24)
is starting to be measurable. The matrix elements for the spin-orbit operator are calculated in Appendix .3.4. Since
the Empirical Pseudopotential Method is used, the projectors a`jmj(r;; )are calculated using equation (2.74).
For the function f`jwe choose one with the same shape as (3.17) for the same reasons.
f`j=C`jr`er2
R2c: (3.25)
lxiv

Since for`= 0the only allowed value of jisj=1
2, if the degeneracy averaged perturbation is zero, we have
C„+1
2=r
`
`+ 1C„1
2; (3.26)
and
sgn(b„+1
2) =sgn(b„1
2) =1; (3.27)
and we use `max= 2because with spin orbit we describe also the contribution of the core electrons, which in
Germanium also take part the d(`= 2) orbitals.
For searching for the local and non-local parts of the pseudopotential we do the same as for Carbon. So we
start by fitting the local part of the pseudopotential (3.9) at reference [1]. The results are on Figure 3.21 and Table
3.18.
Figure 3.21: The figure shows the fit of the expression (3.9) to the ppseudopotential of Germanium, generated with the
program in reference [1]
RaqzkTF
0.94 1.83 0.59
Table 3.18: The table shows the results to the fit of equation (3.9) to the ppseudopotential of Germanium generated with the
program in reference [1]
In the same way with the previous section, a ballpark figure was obtained with the computing of expression
(3.21) and we solve the Schrodinger equation and find out the pairs of values ( Ra,kTF) for which the eigenvalue
of energy is the one of the Epfor bulk Carbon. This value is Ep=4:05eVand comes from reference [1]. It is
also calculated the pairs for which the eigenvalue of energy is Ep0:5eV. In Figure 3.16 all these results are
condensate. The yellow line is a fit of the pairs (Ra;kTF)for whichEp=4:05eVto a linear function. The results
are in Figure 3.22.
Using (3.22), we fix the parameters on Table 3.19, fit the non-local part of the potential (Figure 3.23, Table 3.20),
RaqzkTF
0.94 1.83 0.52
Table 3.19: The table shows the results to the fitting to the ppseudopotential of Germanium generated with the program in
reference [1]
and the ballpark figure is draw with the same principle as before, using equation (3.23), with Es=11:92eVfrom
[1] (Figure 3.24).
Finally we set the parameters for the local and non-local pseudopotentials (Table 3.21), and use it to calculate
lxv

Figure 3.22: The ballpark figure to fit the local part of the pseudopotential using function (3.21) is shown. The green point is
the result of the fit, in Table 3.18, the red points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy is the one of
theEp, the yellow line is adjusted to these points, the orange points are the pairs of values ( Ra,kTF) for which the eigenvalue
of energy is the one of the Ep+ 0:5eVand the yellow points are the pairs of values ( Ra,kTF) for which the eigenvalue of energy
is the one of the Ep0:5eV.Ep=4:05eVis obtained from reference [1] for Germanium
Figure 3.23: The fit of the expression (3.20) to the s“projector” of Germanium generated with the program in reference [1] is
shown
RbB0
0.71 12.1
Table 3.20: Results to the fit of equation (3.20) to the s“projector” of Germanium generated with the program in reference [1]
are shown here
RaqzkTFRbB0
0.94 1.83 0.52 0.71 10.5
Table 3.21: This obtained parameters are used to calculate the some transitions of Germanium without the spin-orbit part
some transitions of Germanium, without the spin-orbit part (Table 3.22), indicated in Figure 3.28. Only some of the
transitionsTnindicated in the Figure are calculated in this case, because the ones that are not calculated have that
value due to the spin-orbit splitting that we are not taking into account. Also we only fit the values in Land, since
they are more precise than the others. After we finely adjust the new set of parameters (Table 3.23) to calculate
the band structure of Figure 3.25. We observe that the shape is quite similar to the one in Figure 3.28 except for
the spin-orbit splitting and the transitions of Table 3.22 are calculated again in Table 3.24 with the parameters of
Table 3.23.
lxvi

Figure 3.24: The ballpark figure to fit the non-local part of the pseudopotential of Germanium using function (3.23) is shown.
The green point is the result of the fit, in Table 3.13, the red points are the pairs of values ( Rb,B) for which the eigenvalue of
energy is the one of the Es, the yellow parabola is adjusted to these points, the orange points are the pairs of values ( Rb,B) for
which the eigenvalue of energy is the one of the Es+ 0:5eVand the yellow points are the pairs of values ( Rb,B) for which the
eigenvalue of energy is the one of the Es0:5eV.Es=11:92eVis obtained from reference [1]
Transitions 0Ra+ 0:01kTF+ 0:01Rb+ 0:01B+ 0:1
T2= 3:0 3.04158 3.02267 3.0137 3.04158 3.04158
T4= 12:912.0827 12.0869 12.0909 11.9919 12.0636
T5= 7:9 8.33772 8.28299 8.25422 8.34293 8.33899
T6= 4:2 4.1031 4.07357 4.06538 4.1031 4.1031
T7= 0:74 0.25232 0.24454 0.23122 0.35451 0.27527
T8= 1:4 1.14544 1.15008 1.1511 1.14544 1.14544
T9= 7:7 7.12262 7.12819 7.13326 7.05045 7.10688
T10= 10:610.2042 10.2014 10.2048 10.1029 10.1827
Table 3.22: Important energetic transitions of Germanium, calculated without using the spin-orbit splitting, with the parameters
on Table 3.21
RaqzkTFRbB0
0.94 1.83 0.52 0.71 12.0
Table 3.23: Discovered parameters, after the better “adjustment” to the experiment, used to calculate the band structure of
Germanium, without the spin-orbit splitting. The “adjustment” is as close as possible because without the spin-orbit contribution
we cannot fully describe Germanium.
Figure 3.25: Calculated band structure of Germanium, without the spin-orbit splitting, with the parameters on Table 3.23
We see that T6is much better adjusted but others, T4,T5,T8,T9are a little bit less. After we introduce the
spin-orbit splitting. We set Rc=Rband adjust C`jin a way that the splitting of the most energetic valence bands
lxvii

T2T4T5T6T7T8T9T10
Experience 3.0 12.9 7.9 4.2 0.74 1.4 7.7 10.6
Calculated 3.0 11.8 8.4 4.1 0.61 1.5 7.0 9.9
Table 3.24: Important transitions calculated for a bulk Germanium without the spin-orbit splitting, with the parameters on Table
3.23
RaqzkTFRbB00RcC`j
0.94 1.83 0.52 0.71 12.0 0.71 0.279
Table 3.25: This set of parameters are used to calculate important transitions of Germanium, with the spin-orbit splitting
isT0= 0:296. This value is C`j= 0:279for all`andj. So we use the parameters of Table 3.25 to calculate the all
the transitions Tnin Table 3.26. As we can see from table 3.26 the introduction of the spin-orbit splitting changed
the band structure and changed the values of Tn. Using Table 3.26 we adjust again the parameters of Table 3.25
to the ones in Table 3.27.
With this parameters we can plot and compare the separate potentials Vlocal(~k),
Vnonlocal (~k)and
Vspinorbit (~k)
in Figure 3.26 with ~kgoing in the direction of X=1 0 0

. Notice the widely different vertical scales. We can see
Figure 3.26: The pseudopotentials a) Vlocal;screen (kx), b)
Vnonlocal (kx)and c)
Vspinorbit (kx)of Germanium are represented
graphically, using the parameters of Table 3.27
that the most important contribution for the electronic structure of the solid is the local pseudopotential, followed
by the non-local part and the spin-orbit pseudopotential, as it is to be expected, since from equations (2.58) and
(2.60), these potentials are treated as perturbations. We can calculate the sum of the potentials above, plotted in
Figure 3.27. We can see that the non-local pseudopotential still has a significant effect. The effect of the spin-orbit
contribution is very small but is still enough to do the splitting of the last two valence pbands and the first two
conduction pbands as we can see in Figure 3.28.
We calculate the optical transitions, including those of 0:85XandXindicated in Figure 3.28, in Table 3.28 and
observe that most of the values are much better fitted than previously.
After, with the program in we calculated D(E),J(E),1(!),2!andR(!)in the same way we did before.
Transitions Ra+ 0:01kTF+ 0:01Rb+ 0:01B00+ 0:1Rc+ 0:01C`j+ 0:1
T0=0.296 0.295937 0.294133 0.295015 0.295937 0.295937 0.325729 0.38345
T1=0.888 0.0409927 0.0784698 0.0588591 0.309062 0.0927434 0.0166953 0.146117
T2=3.0 2.61454 2.59893 2.58871 2.61454 2.61454 2.57129 2.40364
T3=0.21 0.227238 0.225022 0.225656 0.227238 0.227238 0.249606 0.406467
T4=12.9 12.0424 12.0455 12.0514 11.9349 12.023 12.0667 12.0834
T5=7.9 8.1037 8.05216 8.02233 8.10862 8.10473 8.07875 8.05465
T6=4.2 3.80045 3.77378 3.7645 3.80046 3.80045 3.77077 3.69779
T7=0.74 0.356431 0.353207 0.336776 0.484888 0.381147 0.331019 0.302728
T8=1.4 1.27053 1.27433 1.27579 1.27052 1.27053 1.28317 1.22654
T9=7.7 7.14404 7.14962 7.15642 7.06527 7.12927 7.16898 7.19284
T10=10.6 10.1278 10.1236 10.1291 10.0082 10.106 10.1522 10.1703
Table 3.26: Important transitions of Germanium, with the spin-orbit splitting, are calculated with the parameters in Table 3.25
lxviii

RaqzkTFRbB00RcC`j
0.94 1.83 0.52 0.74 12.0 0.71 0.279
Table 3.27: The set of parameters used to calculate the band structure of Germanium, with the spin-orbit splitting are shown
Figure 3.27: The pseudopotentials Vlocal;screen (kx),Vlocal;screen (kx) +
Vnonlocal (kx)andVlocal;screen (kx) +
Vnonlocal (kx) +

Vspinorbit (kx)of Germanium are calculated
Figure 3.28: Calculated band structure of Germanium, a) using reference [5] and with the values of some important transitions
in eV, b) with the program, using the parameters of Table 3.27
T0T1T2T3T4T5T6T7T8T9T10 LB PRB32 PRB32
Experience 0.296 0.888 3.0 0.21 12.9 7.9 4.2 0.74 1.4 7.7 10.6 1.3 3.5 9.5
Calculated 0.296 0.882 2.6 0.23 11.7 8.1 3.8 0.76 1.3 6.9 9.8 1.0 2.9 8.1
Table 3.28: Important optical transitions of Germanium where calculated using the parameters of Table 3.27
We calculated a grid of 777points. The results are in Figure 3.29. We can compare this results with the
experiment from references [6] for D(E)and [7] for (!)andR(!), in the same Figure, and observe that in both
in the DOS, the peaks are in the same place, except in the calculated DOS in Figure 3.29 a) (and also the JDOS
in b)), we have the noise of the Gaussians, that as said before would be fixed if we did Ni!+1and!0. In
the real part of the dielectric function, the calculated peaks are slightly dislocated to higher energies of 0:5eV.
About the imaginary part, in Figure c), we indeed have one peak at 2:3eVand other at4:2eV, but they have
different intensities. The rest of the Figure c) hardly agrees with the experiment. About the reflectance, we can see
in Figure d), that the peak at 4:5eVis dislocated to a higher energy. The calculated Figures are very affected
by the small number of points in the grid. We could have a better quality of the figures if we used a grid with more
lxix

Figure 3.29: It is represented a) the calculated density if states (blue line),the photo emission spectroscopy and inverse photo
emission data obtained from reference [6] (yellow line) divided by a factor of 3, b) the calculated joint density of states, c)
calculated dielectric function, real part (blue) and imaginary (purple), d) calculated (blue line) and experimental (yellow line, [7])
1, e) calculated (blue line) and experimental (yellow line, [7]) 2f),g) calculated (blue line) and experimental (yellow line, [7])
reflectance for Germanium with the local pseudopotential of equation (3.7), the non-local projector for the pseudopotential of
equation (3.19) with the parameters written in Table 3.9 and the spin-orbit projectors using equations (2.74) and (3.25-3.27)
withlmax= 2
points, but in a software like M ATHEMATICA it increases significantly the computation time needed to do so.
lxx

4
Conclusions and Future Work
lxxi

In this research project we obtained a description of bulk Silicon, Carbon and Germanium, group IV elements,
using the pseudopotentials in the Empirical Pseudopotential Method. We obtained a successful description for
Silicon, since the band structure was pretty well adjusted, and we did that by using parameters in a pseudopotential
model that have physical meaning. The reflectance and the imaginary part of the dielectric function could be better,
but the calculation of this properties allowed us to diagnose a problem with the band structure that was drawn. As
a future work, the pseudopotential should also be adjusted to the experimental reflectance data, which seems to
promise a better success. Also the fact that we used a small number of points aggravated our predictions for the
optical properties of Silicon, because, since we used the M ATHEMATICA software, we could not use a big number
of points without costing a lot of time to do the calculations. For Germanium, the results are similar to Silicon. The
band structure and the density of states was pretty well adjusted to the experimental values, but in the future we
should also try to adjust to the optical properties such as reflectance, to obtain better results. To Carbon (Diamond),
we obtained the least successful description, since it was very difficult to adjust the band structure with just a small
number of parameters. This element requires some extra work and, in the future the research could go into finding
a pseudopotential model with more parameters to describe it.
All the tools where developed to continue this research and improve each of the descriptions. As a next step
we can introduce the pseudopotentials in a F ORTRAN program and increase the number of points in the cubic
grid, to obtain better results. Using the parameters found here, we already did that, as we can see in Figure 4.1
for Silicon, Figure 4.2 for Carbon and Figure 4.3 for Germanium. We compare the obtained band structure and
density of states with a DFT -MGGA (Density Functional Theory with a Meta Generalized Gradient Approximation).
DFT -MGGA is right now a state of the art ab initio theory but still doesn’t give the right gaps to semiconductors.
We can see in the figures that we have much better results only by adding a lot more points in the grid.
Figure 4.1: The figure shows the a) band structure and the b) density of states, calculated for Silicon using the pseudopotentials
obtained here and using a DFT -MGGA calculation
After having the best description possible of each of the elements, as a future research the pseudopotentials
will be used to be introduced in a supercell composed with these elements. This was already done for a Si29C
supercell, a supercell of Si with a C impurity. In Figure 4.4 is the calculated unfolded band structure.
The purpose is to find the best material possible with the best optoelectronic properties to be integrated in the
electronics and technology of today.
lxxii

Figure 4.2: The figure shows the a) band structure and the b) density of states, calculated for Carbon using the pseudopoten-
tials obtained here and using a DFT -MGGA calculation
Figure 4.3: The figure shows the a) band structure and the b) density of states, calculated for Germanium using the pseu-
dopotentials obtained here and using a DFT -MGGA calculation
Figure 4.4: The figure shows the a) band structure and the b) density of states, calculated for Germanium using the pseu-
dopotentials obtained here and using a DFT -MGGA calculation
lxxiii

lxxiv

Bibliography
[1] N. T. S Froyen, “A pseudopotential generation program.”
[2] H. M. P . Tserbak, C. and G. Theodorou, “Unified approach to the electronic structure of strained Si/Ge super-
lattices,” Physical Review B , 1993.
[3] H. J. Osten, “Band-gap changes and band offsets for ternary Si 1xyGexCyalloys on Si (001),” Journal of
applied physics , 1998.
[4] J. R. Chelikowsky and M. L. Cohen, “ Ab initio pseudopotentials and the structural properties of semiconduc-
tors,” Handbook on Semiconductors, ed. PT Landsberg, Elsevier, Amsterdam 1 , 1992.
[5] J. L. M. et all, “Program for calculation of the band structure using ab initio lda pseudopotentials.”
[6] W. J. Chelikowsky, Wagener, “Valence- and conduction-band densities of states for tetrahedral semiconduc-
tors: Theory and experiment,” Physical Review B , 1989.
[7] D. E. Aspnes and A. A. Studna, “Dielectric functions and optical parameters of si, ge, gap, gaas, gasb, inp,
inas, and insb from 1.5 to 6.0 ev,” Physical Review B , 1983.
[8] M. et al, “X-ray photoemission studies of diamond, graphite, and glassy carbon valence bands,” Physical
Review B , 1974.
[9] Roberts and Walker, “Optical Study of the Electronic Structure of Diamond,” Physical Review , 1967.
[10] W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Physical
Review , 1965.
[11] J. P . Perdew and Y . Wang, “Accurate and simple analytic representation of the electron-gas correlation energy,”
Physical Review B , 1992.
[12] F . Tran and P . Blaha, “Accurate band gaps of semiconductors and insulators with a semilocal exchange-
correlation potential,” Physical Review Letters , 2009.
[13] M. S. Aulbur, Wilfried G. and A. G ¨orling, “Exact-exchange-based quasiparticle calculations,” Physical Review
B, 2000.
[14] N. D. . S. D. J. Kouvetakis, J., “Synthesis and atomic and electronic structure of new Si-Ge-C alloys and
compounds,” Chemistry of materials , 1998.
[15] C. N. E. A. M. . C. M. Schmid, U., “Electronic and optical properties of strained Ge/Si superlattices,” Physical
Review B , 1991.
[16] . D. J. D. Newman, K. E., “Theory of deep impurities in silicon-germanium alloys,” Physical Review B , 1984.
lxxv

[17] L. J. W. C. T. . Z. A. d’Avezac, M., “Genetic-algorithm discovery of a direct-gap and optically allowed super-
structure from indirect-gap Si and Ge semiconductors,” Physical Review Letters , 2012.
[18] J. L. Martins, “Unpublished lecture notes in condensed matter physics.”
[19] Alloul, Henri, and Stephen Lyle, Introduction to the Physics of Electrons in Solids . Springer, 2010.
[20] M. Born and K. Huang, Dynamical Theory of Crystal Lattices . Clarendon Press, Oxford, 1988.
[21] N. Troullier and J. L. Martins, “Efficient pseudopotentials for plane-wave calculations,” Physical Review B ,
1991.
[22] J. R. Chelikowsky, “Empirical Pseudopotentials for Semiconductors,” KLUWER INTERNATIONAL SERIES IN
ENGINEERING AND COMPUTER SCIENCE , 1996.
[23] N. W. Ashcroft and N. D. Mermin., Solid State Physics . Saunders College, Philadelphia, 1976.
[24] L.-W. Wang and A. Zunger, “Local-density-derived semiempirical pseudopotentials,” Physical Review B , 1995.
[25] G. Theurich and N. A. Hill, “Self-consistent treatment of spin-orbit coupling in solids using relativistic fully
separable ab initio pseudopotentials,” Physical Review B , 2001.
[26] L. Kleinman and D. M. Bylander, “Efficacious form for model pseudopotentials,” Physical Review Letters , 1982.
[27] L. Kleinman, “Relativistic norm-conserving pseudopotential,” Physical Review B , 1980.
[28] Kittel, Charles, and Paul McEuen, Introduction to solid state physics . New Y ork: Wiley, 1976.
[29] Landau, L. D., and E. M. Lifshitz, Statistical physics, vol. 5 . Course of theoretical physics 30, 1980.
[30] Yu, Peter Y ., and Manuel Cardona, Fundamentals of semiconductors . Springer, 1996.
lxxvi

5
Appendix
.1 Spin-Orbit projectors for the pseudopotential
Here is explained the form of the “projectors” a`jmj. We start with the psedopotential on the Kleinman-Bylander
form considering the spin,
VKB=X
`;j;mj
VSO
`;jPP
`;j;mjED
PP
`;j;mj
VSO
`;j
D
PP
`;j;mj
VSO
`;jPP
`;j;mjE; (1)
withPP
`;j;mjE
=1p
2R`jY`mj1
2(;)
Y`mj+1
2(;)
and rewrite it in an other way:
VKB =X
`;j;mj
VSO
`;jPP
`;j;mjE
sgnD
PP
`;j;mj
VSO
`;jPP
`;j;mjED
PP
`;j;mj
VSO
`;j
D
PP
`;j;mj
VSO
`;jPP
`;j;mjE=
=X
`;j;mj
VSO
`;jPP
`;j;mjE
rD
PP
`;j;mj
VSO
`;jPP
`;j;mjEsgnD
PP
`;j;mj
VSO
`;jPP
`;j;mjED
PP
`;j;mj
VSO
`;j
rD
PP
`;j;mj
VSO
`;jPP
`;j;mjE:(2)
The value of b`jis
b`j=D
PP
`;j;mj
VSO
`;jPP
`;j;mjE
=Z2
0Z
0Z+1
0R
`j(r)1p
2
Y
`mj1
2(;)Y
`mj+1
2(;)

VSO
`j(r)R`j(r)1p
2Y`mj1
2(;)
Y`mj+1
2(;)
r2sindrdd =
=1
2Z2
0Z
0[Y
`mj1
2(;)Y`mj1
2(;) +Y
`mj+1
2(;)Y`mj+1
2(;)] sinddZ+1
0R
`j(r)
VSO
`j(r)R`j(r)r2dr=
=Z+1
0R
`j(r)
VSO
`j(r)R`j(r)r2dr; (3)
whereY`mare the spherical harmonic functions and we did m=mjms, wherems=1
2is the spin angular
momentum. It is independent of mj. The potential is, therefore,
1

VKB =X
`;j;mj
VSO
`;jPP
`;j;mjE
p
jb`jjsgn(b`j)D
PP
`;j;mj
VSO
`;j
p
jb`jj; (4)
and we define the “projector” spinor a`jmjas
a`;j;mj
=
VSO
`;jPP
`;j;mjE
p
jb`jj: (5)
and thus we have,
VKB =X
`;j;mja`;j;mj
sgn(b`j)
a`;j;mj: (6)
Therefore, the functions of the projectors are
a`jmj=1p
2jb`jj
VSO
`;j(r)R`j(r)Y`mj1
2(;)
Y`mj+1
2(;)
: (7)
.2 Clebsch-Gordon coefficients for mixing states lands=1
2
We want to calculate the Clebsh-Gordon coefficients for the total angular momentum jwhen we have a state
with orbital land spins=1
2angular momenta. To do so, we se the raising and lowering operators, J, thex;y;z
operatorsJx=1
2(J++J),Jx=1
2i(J+J)andJz, and also the J2operator.
When we apply the Joperators to a state jjmjithe results are:
J+jjmji=~
2q
j(j+ 1)mj(mj+ 1)jjmj+ 1i (8)
Jjjmji=~
2q
j(j+ 1)mj(mj1)jjmj1i: (9)
With theJx;yoperators:
Jxjjmji=1
2[J+jjmji+Jjjmji] =
=~
2q
j(j+ 1)mj(mj+ 1)jjmj+ 1i+q
j(j+ 1)mj(mj1)jjmj1i
(10)
Jyjjmji=1
2i[J+jjmjiJjjmji] =
=~
2iq
j(j+ 1)mj(mj+ 1)jjmj+ 1iq
j(j+ 1)mj(mj1)jjmj1i
: (11)
and theJzandJ2operators:
Jzjjmji=~mjjjmji (12)
J2jjmji=~2j(j+ 1)jjmji: (13)
Using the rules of the addition of angular momenta for ~J=~L+~S, the quantum numbers satisfy:
l1
2<j <l +1
2: (14)
2

Forl= 0, the only allowed value for jisj=1
2. For the other cases
l1
2<j <l +1
2; (15)
which means that the allowed values for jarej=l1
2andj=l+1
2. Thezprojections of the momenta always
satisfymj=ml+mswithms=1
2. So, the statejjmjican be described by a linear combination of the states
lmj1
21
21
2
andlmj+1
21
21
2
jjmji=Almj1
21
21
2
+Blmj+1
21
21
2
: (16)
First we apply the J2operator, for which:
J2=L2+S2+ 2~L
~S
~L=Lx~ ex+Ly~ ey+Lz~ ez
~S=Sx~ ex+Sy~ ey+Sz~ ez
~L
~S=LxSx+LySy+LzSz
J2=L2+S2+ 2(LxSx+LySy+LzSz): (17)
Applying eq. (17) to the jjmjistate, we have that:
J2jjmji=
L2+S2+ 2(LxSx+LySy+LzSz)
Almj1
21
21
2
+Blmj+1
21
21
2
=
=A
L2lmj1
21
21
2
+lmj1
2
S21
21
2
+
+2
Lxlmj1
2
Sx1
21
2
+
Lylmj1
2
Sy1
21
2
+
+
Lzlmj1
2
Sz1
21
2
+B
L2lmj+1
21
21
2
+lmj+1
2
S21
21
2
+
+2
Lxlmj+1
2
Sx1
21
2
+
Lylmj+1
2
Sy1
21
2
+
+
Lzlmj+1
2
Sz1
21
2
:
Using equations (10) to (13) and since the states1
2;3
2
do not exist (because msonly has values from sto
s) we have:
J2jjmji=A3
4~21
21
2lmj1
2
+~2l(l+ 1)1
21
2lmj1
2
+
+2"
~
21
21
2~
2 s
l(l+ 1)
mj1
2
mj+1
2lmj+1
2
+
+s
l(l+ 1)
mj1
2
mj3
2lmj3
2!
+
+i~
21
21
2~
2i s
l(l+ 1)
mj1
2
mj+1
2lmj+1
2
+
s
l(l+ 1)
mj1
2
mj3
2lmj3
2!
+
+~
21
21
2
~(mj1
2)lmj1
2
+
3

B3
4~21
21
2lmj+1
2
+~2l(l+ 1)1
21
2lmj+1
2
+
+2"
~
21
21
2~
2 s
l(l+ 1)
mj+1
2
mj+3
2lmj+3
2
+
+s
l(l+ 1)
mj+1
2
mj1
2lmj1
2!
+
+i~
21
21
2~
2i s
l(l+ 1)
mj+1
2
mj+3
2lmj+3
2
+
s
l(l+ 1)
mj+1
2
mj1
2lmj1
2!
+
~
21
21
2
~(mj+1
2)lmj+1
2
: (18)
As we said before, in eq. (16), the state jjmjiis only described bylmj1
21
21
2
andlmj+1
21
21
2
states, which means thelmj3
2
are not present:
J2jjmji=~2(
A3
4+l(l+ 1) +mj1
2
+Br
l(l+ 1)m2
j+1
4)1
21
2lmj1
2
+
~2(
B3
4+l(l+ 1)mj1
2
+Ar
l(l+ 1)m2
j+1
4)1
21
2lmj+1
2
; (19)
but also
J2jjmji=~2j(j+ 1)jjmji
=~2j(j+ 1)
Almj1
21
21
2
+Blmj+1
21
21
2
: (20)
So we have the system of equations:
8
<
:A
l(l+ 1) +1
4+mj
+Bq
l(l+ 1)m2
j+1
4=j(j+ 1)A
B
l(l+ 1) +1
4mj
+Aq
l(l+ 1)m2
j+1
4=j(j+ 1)B,
,8
<
:A
l(l+ 1)j(j+ 1) +1
4+mj
+Bq
l(l+ 1)m2
j+1
4= 0
B
l(l+ 1)j(j+ 1) +1
4mj
+Aq
l(l+ 1)m2
j+1
4= 0,
,(
A(a+mj) +Bb = 0
B(amj) +Ab = 0witha=l(l+ 1)j(j+ 1) +1
4
b=q
l(l+ 1)m2
j+1
4: (21)
Multiplying the first equation of the system by (amj)and the second by bwe have
A(a2m2
j) +Bb(amj) = 0
Bb(amj) +Ab2= 0: (22)
Subtracting the second equation from the first one:
A(a2m2
jb2) = 0,a2b2=m2
j,
,
l(l+ 1)j(j+ 1) +1
42
l(l+ 1) +m2
j1
4=m2
j,
,
l(l+ 1)j(j+ 1) +1
42
=l2+l+1
4=
l+1
22
,
,l(l+ 1)j(j+ 1) +1
4=
l+1
2
,
4

,j(j+ 1) =l(l+ 1)
l+1
2
+1
4,: (23)
Adding1
4to both sides of the equation:
j2+j+1
4=l(l+ 1)
l+1
2
+1
2,
j+1
22
=l2+ll1
2+1
2=l2
l2+l+l+1
2+1
2= (l+ 1)2,
,
j+1
2=l
j+1
2=(l+ 1),8
>><
>>:j =l1
2=l1
2
l1
2
j+1
2=(l+ 1)1
2=l+1
2
l3
2: (24)
But we know that j0so the only good values are j=1
2.
Having this result, we have that the values for aandbare:
a=l2+l
l1
2
l1
2+ 1
+1
4=
=l2+ll21
2ll1
2l1
41
2+1
4=
=l1
2=(l+1
2) (25)
b=s
l2+l+1
4
m2
j=s
l+1
22
m2
j=
=s
l+1
2+mj
l+1
2mj
: (26)
Putting this in eq. (21):
A

l+1
2
+mj
=A
l+1
2mj
=Bs
l+1
2+mj
l+1
2mj
,
,Ar
l+1
2mj=Br
l+1
2mj,jAj2
l+1
2mj
=jBj2
l+1
2mj
: (27)
ButjAj2+jBj2= 1
jAj2
l+1
2mj
= (1jAj2)
l+1
2mj
,jAj2l+1
2mj
l+1
2mj= 1jAj2,
,jAj2
1 +l+1
2mj
l+1
2mj
= 1,jAj2l+1
2mj+l+1
2mj
l+1
2mj=2l+ 1
l+1
2mj= 1: (28)
The Clebsch-Gordon coefficents are always real, and using the phase convention
l1
2;l1
2l+1
2l+1
2
0,A
has to be bigger than 0, so
A=s
l+1
2mj
2l+ 1=r
2l+ 12mj
4l+ 2: (29)
We can now calculate B:
Ar
l+1
2mj=Br
l+1
2mj,B=Aq
l+1
2mj
q
l+1
2mj=r
2l+ 1mj
2l+ 1; (30)
for which the upper and lower signs correspond to j=l1
2:
5

.3 Matrix elements of the momentum matrix operator
Considering a semiconductor crystal, a plane wave basis set is used for the wave functions inside the crystal
D
~ r~kn;nE
=1p
Vcrystalei~kn
~ r; (31)
where the position vector ~ rinside the crystal is defined by
~ r=y1N1~ a1+y2N2~ a2+y3N3~ a3; (32)
with0<yn<1and~ aiare the primitive lattice vectors, and in the reciprocal lattice,
~k=i1
N1~b1+i2
N2~b2+i3
N3~b3 (33)
withij= 0;:::;Nj1and~biare the primitive reciprocal lattice vectors. Vcell=~ a1
(~ a2~ a3)is the volume of the
primitive cell and the whole crystal has the size N1~ a1N2~ a2N3~ a3. The basis functions for the free electrons in
the conduction and valence bands, respectively are
D
~ r~kc;ncE
=1p
Vcrystalei~kc
~ randD
~ r~kv;nvE
=1p
Vcrystalei~kv
~ r: (34)
We now calculate the matrix element
D
~kc;nc~ e
~ p~kv;nvE
=iD
~kc;nc~ e
~r~kv;nvE
=i1
VcrystalZZZ
crystalei~kc
~ r(~ e
~r)ei~kv
~ rd3r
=i1
VcrystalZZZ
crystalei~kc
~ r(~ e
~kv)ei~kv
~ rd3r
=~ e
~kv
VcrystalZZZ
crystalei~kc
~ rei~kv
~ rd3r
=~ e
~kv
VcrystalZZZ
crystalei(~kv~kc)
~ rd3r: (35)
By making the change of variables ~ r=y1N1~ a1+y2N2~ a2+y3N3~ a3,
ZZZ
crystald3~ r=CZ1
0dy1Z1
0dy2Z1
0dy3,Vcrystal =C; (36)
and replacing ~kfor (33),
D
~kc;nc~ e
~ p~kv;nvE
=~ e
~kv
VcrystalVcrystalZ1
0Z1
0Z1
0ei(i1vi1c
N1~b1+i2vi2c
N2~b2+i3vi3c
N3~b3)
(y1N1~ a1+y2N2~ a2+y3N3~ a3)dy1dy2dy3:(37)
SinceD
~ ai~bjE
= 2ij,
D
~kc;nc~ e
~ p~kv;nvE
=~ e
~kv
VcrystalVcrystalZ1
0ei2(i1vi1c)y1dy1Z1
0ei2(i2vi2c)y2dy2Z1
0ei2(i3vi3c)y3dy3=
=~ e
~kvei2(i1vi1c)y1
i2(i1vi1c)1
0ei2(i2vi2c)y2
i2(i2vi2c)1
0ei2(i3vi3c)y3
i2(i3vi3c)1
0=
=~ e
~kvei2(i1vi1c)1
i2(i1vi1c)ei2(i2vi2c)1
i2(i2vi2c)ei2(i3vi3c)1
i2(i3vi3c)=
=~ e
~kvcos[2(i1vi1c)] +isin[2(i1vi1c)]1
i2(i1vi1c)
6

cos[2(i2vi2c)] +isin[2(i2vi2c)]1
i2(i2vi2c)
cos[2(i3vi3c)] +isin[2(i3vi3c)]1
i2(i3vi3c): (38)
Since 2(invinc)is an integer multiplied by 2, the value of the cosfunction with this argument is 1, while the
sinis always 0. Therefore, if ~kc6=~kv,D
~kc;nc~ e
~ p~kv;nvE
= 0. If~kc=~kv,
D
~kc;nc~ e
~ p~kv;nvE
=~ e
~kv
VcrystalZZZ
crystald3r=~ e
~kv
VcrystalVcrystal =~ e
~kv (39)
.3.1 Orthogonality of the basis functions
We have a 3-dimensional crystal with primitive lattice vectors ~ aj,Natom atoms in the primitive cell with atomic
positions~ i. The reciprocal lattice vectors are
~b1= 2~ a2~ a3
Vcell~b2= 2~ a3~ a1
Vcell~b3= 2~ a1~ a2
Vcell(40)
withVcell=~ a1
(~ a2~ a3)being the volume of the primitive cell. The whole crystal has the size N1~ a1N2~ a2N3~ a3.
The position vector ~ rinside the crystal is defined by
~ r=y1N1~ a1+y2N2~ a2+y3N3~ a3 (41)
with0<yn<1and in the reciprocal lattice:
~ki1i2i3=i1
N1~b1+i2
N2~b2+i3
N3~b3 (42)
withij= 0;:::;Nj1. The reciprocal lattice vectors ~G(the vectors that yield plane waves with the periodicity of
the Bravais lattice) are defined:
~G=I1~b1+I2~b2+I3~b3 (43)
A plane wave basis set is used for the wave functions inside the crystal:
D
~ r~ki1i2i3;~GI1I2I3E
=1pN1N2N3Vcellei(~ki1i2i3+~GI1I2I3)
~ r(44)
We can prove they are orthogonal for different vectors ~ki1i2i3and~GI1I2I3. Most derivations of this kind assume
normalizations. Here we show a lengthy but detailed derivation, using a normalization in the crystal.
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I0
3~ki1i2i3;~GI1I2I3E
=ZZZ
crystal1pN1N2N3Vcellei(~k0
i0
1i0
2i0
3+~G0
I0
1I0
2I0
3)
~ r 1pN1N2N3Vcellei(~ki1i2i3+~GI1I2I3)
~ rd3~ r=
=1
N1N2N3VcellZZZ
crystalei(~ki1i2i3~k0
i0
1i0
2i0
3)
~ rei(~GI1I2I3~G0
I0
1I0
2I0
3)d3~ r (45)
If~ki1i2i3=~k0
i0
1i0
2i0
3and~GI1I2I3=~G0
I0
1I0
2I0
3
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I0
3~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellZ
crystal1d3~ r
=1
N1N2N3VcellN1N2N3Vcell= 1 (46)
If not, we can make the change of variables ~ r=y1N1~ a1+y2N2~ a2+y3N3~ a3
ZZZ
crystald3~ r=CZ1
0dy1Z1
0dy2Z1
0dy3,N1N2N3Vcell=C (47)
7

D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I0
3~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN1N2N3Vcell
Z1
0Z1
0Z1
0ei(i1i0
1
N1~b1+i2i0
2
N2~b2+i3i0
3
N3~b3)
(y1N1~ a1+y2N2~ a2+y3N3~ a3)
ei((I1I0
1)~b1+(I2I0
2)~b2+(I3I0
3)~b3)
(y1N1~ a1+y2N2~ a2+y3N3~ a3)dy1dy2dy3=
=Z1
0ei2(i1i0
1+(I1I0
1)N1)y1dy1Z1
0ei2(i2i0
2+(I2I0
2)N2)y2dy2Z1
0ei2(i3i0
3+(I3I0
3)N3)y3dy3=
="
ei2(i1i0
1+(I1I0
1)N1)y1
i2(i1i0
1+ (I1I0
1)N1)#1
0"
ei2(i2i0
2+(I2I0
2)N2)y2
i2(i2i0
2+ (I2I0
2)N2)#1
0"
ei2(i3i0
3+(I3I0
3)N3)y3
i2(i3i0
3+ (I3I0
3)N3)#1
0=
=ei2(i1i0
1+(I1I0
1)N1)1
i2(i1i0
1+ (I1I0
1)N1)ei2(i2i0
2+(I2I0
2)N2)1
i2(i2i0
2+ (I2I0
2)N2)ei2(i3i0
3+(I3I0
3)N3)1
i2(i3i0
3+ (I3I0
3)N3)=
=cos [2(i1i0
1+ (I1I0
1)N1)] +isin [2(i1i0
1+ (I1I0
1)N1)]1
i2(i1i0
1+ (I1I0
1)N1)
cos [2(i2i0
2+ (I2I0
2)N2)] +isin [2(i2i0
2+ (I2I0
2)N2)]1
i2(i2i0
2+ (I2I0
2)N2)
cos [2(i3i0
3+ (I3I0
3)N3)] +isin [2(i2i0
2+ (I2I0
2)N2)]1
i2(i3i0
3+ (I3I0
3)N3)(48)
Since we have in the argument of the co-sinus and sinus functions, i2(ini0
n+InI0
n)Nn, which is an integer
number multiplied by 2, thecosfunction will be equal to 1, while the sinwill be zero. Therefore, if in6=i0
nand
In6=I0
nwe have thatD
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I0
3~ki1i2i3;~GI1I2I3E
= 0, so in conclusion:
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I0
3~ki1i2i3;~GI1I2I3E
=i1i0
1i2i0
2i3i0
3I1I0
1I2I0
2I3I0
3=~ki1i2i3~k0
i0
1i0
2i0
3~GI1I2I3~G0
I0
1I0
2I0
3(49)
.3.2 Local Pseudopotencial
Here we calculate the matrix elements of the local pseudopotential. As said before, most derivations use
no normalization. Here we show an accurate derivation, with normalization in the crystal. For a crystal with
primitive lattice vectors ~ aj,Natom atoms in the primitive cell with atomic positions i, and dimensions Nj~ aj, the
local pseudopotential is given by
V(L)(~ r) =X
J12ZX
J22ZX
J32ZN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1v(L)(~ r~tj1j2j3~TJ1J2J3~ i) (50)
where
~tj1j2j3=j1~ a1+j2~ a2+j3~ a3
~TJ1J2J3=J1N1~ a1+J2N2~ a2+J3N3~ a3 (51)
are the translation vectors for the finite crystal and for its periodic crystal images. The matrix elements of the
pseudopotential in a plane wave basis set
D
~ r~ki1i2i3;~GI1I2I3E
=1pN1N2N3Vcellei(~ki1i2i3+~GI1I2I3)
~ r(52)
with~ki1i2i3=i1
N1~b1+i2
N2~b2+i3
N3~b3, withij= 0;:::;Nj1, a vector in the primitive reciprocal unit cell and
~GI1I2I3=I1~b1+I2~b2+I3~b3withIj2Za reciprocal lattice vector, are
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=
=ZZZ
crystal1pN1N2N3Vcellei(~k0
i0
1i0
2i0
3+~G0
I0
1I0
2I3)
~ rV(L)(~ r)1pN1N2N3Vcellei(~ki1i2i3+~GI1I2I3)
~ rd3~ r=
=1
N1N1N3VcellZZZ
crystalei(~k0
i0
1i0
2i0
3+~G0
I0
1I0
2I3)
~ r
8

0
@X
J12ZX
J22ZX
J32ZN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1v(L)(~ r~tj1j2j3~TJ1J2J3~ i)1
Aei(~ki1i2i3+~GI1I2I3)
~ rd3~ r(53)
inverting the order the integral and the sums we have
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N1N3VcellX
J12ZX
J22ZX
J32ZN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ZZZ
crystalei(~k0
i0
1i0
2i0
3+~G0
I0
1I0
2I3)
~ rv(L)(~ r~tj1j2j3~TJ1J2J3~ i)ei(~ki1i2i3+~GI1I2I3)
~ rd3~ r=
=1
N1N1N3VcellX
J12ZX
J22ZX
J32ZN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ZZZ
crystalei(~ki1i2i3~k0
i0
1i0
2i0
3+~GI1I2I3~G0
I0
1I0
2I3)
~ rv(L)(~ r~tj1j2j3~TJ1J2J3~ i)d3~ r(54)
Making the change of variables ~ r0=~ r~TJ1J2J3:
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellX
J12ZX
J22ZX
J32ZN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ZZZ
crystal~Tei(~ki1i2i3~k0
i0
1i0
2i0
3+~GI1I2I3~G0
I0
1I0
2I3)
(~ r0+~TJ1J2J3)v(L)(~ r0~tj1j2j3~ i)d3~ r0(55)
But doing the sum in all the crystal translations Jn2Zand integrating in the crystal +~TJ1J2J3is the same as
integrating in all the space, so:
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ZZZ
ei(~ki1i2i3~k0
i0
1i0
2i0
3+~GI1I2I3~G0
I0
1I0
2I3)
(~ r0+~TJ1J2J3)v(L)(~ r0~tj1j2j3~ i)d3~ r0=
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ei(~ki1i2i3~k0
i0
1i0
2i0
3)
~TJ1J2J3eei(~GI1I2I3~G0
I0
1I0
2I0
3)
~TJ1J2J3
ZZZ
ei(~ki1i2i3~k0
i0
1i0
2i0
3+~GI1I2I3~G0
I0
1I0
2I3)
~ r0
v(L)(~ r0~tj1j2j3~ i)d3~ r0(56)
Since~TJ1J2J3is a translation of the crystal, we have that ei(~GI1I2I3~G0
I0
1I0
2I0
3)
~TJ1J2J3= 1, so
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ei(~ki1i2i3~k0
i0
1i0
2i0
3)
~TJ1J2J3
ZZZ
ei(~ki1i2i3~k0
i0
1i0
2i0
3+~GI1I2I3~G0
I0
1I0
2I3)
~ r0
v(L)(~ r0~tj1j2j3~ i)d3~ r0(57)
Making the change of variables ~ r00=~ r0~tj1j2j3~ i
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ei(~ki1i2i3~k0
i0
1i0
2i0
3)
~TJ1J2J3
ZZZ
ei(~ki1i2i3~k0
i0
1i0
2i0
3+~GI1I2I3~G0
I0
1I0
2I3)
(~ r00+~tj1j2j3+~ i)v(L)(~ r00)d3~ r00=
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ei(~ki1i2i3~k0
i0
1i0
2i0
3)
(~TJ1J2J3+~tj1j2j3+~ i)ei(~GI1I2I3~G0
I0
1I0
2I3)
(~tj1j2j3+~ i)
ZZZ
ei(~ki1i2i3~k0
i0
1i0
2i0
3+~GI1I2I3~G0
I0
1I0
2I3)
~ r00
v(L)(~ r00)d3~ r00(58)
But, since~tj1j2j3is a lattice translation in the crystal, we have that ei(~GI1I2I3~G0
I0
1I0
2I0
3)
~tj1j2j3= 1, so
9

D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ei(~ki1i2i3~k0
i0
1i0
2i0
3)
(~TJ1J2J3+~tj1j2j3+~ i)
ei(~GI1I2I3~G0
I0
1I0
2I3)
~ iZZZ
ei(~ki1i2i3~k0
i0
1i0
2i0
3+~GI1I2I3~G0
I0
1I0
2I3)
~ r00
v(L)(~ r00)d3~ r00(59)
Using the Fourier expansion of the potential v(L)
v(L)=X
~G00
I00
1I00
2I00
3a~G00
I00
1I00
2I00
3ei~G00
I00
1I00
2I00
3
~ r(60)
we have that
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ei(~ki1i2i3~k0
i0
1i0
2i0
3)
(~TJ1J2J3+~tj1j2j3+~ i)
ei(~GI1I2I3~G0
I0
1I0
2I3)
~ iZZZ
ei(~ki1i2i3~k0
i0
1i0
2i0
3+~GI1I2I3~G0
I0
1I0
2I3)
~ r000
BB@X
~G00
I00
1I00
2I00
3a~G00
I00
1I00
2I00
3ei~G00
I00
1I00
2I00
3
~ r001
CCAd3~ r00
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1X
~G00
I00
1I00
2I00
3a~G00
I00
1I00
2I00
3ei(~ki1i2i3~k0
i0
1i0
2i0
3)
(~TJ1J2J3+~tj1j2j3+~ i)
ei(~GI1I2I3~G0
I0
1I0
2I3)
~ iZZZ
ei(~ki1i2i3~k0
i0
1i0
2i0
3)ei(~GI1I2I3~G0
I0
1I0
2I3+~G00
I00
1I00
2I00
3)
~ r00
d3~ r00(61)
Using~ r00=y1N1~ a1+y2N2~ a2+y3N3~ a3, with 0<yn<1,~ki1i2i3=i1
N1~b1+i2
N2~b2+i3
N3~b3, withij= 0;:::;Nj1
and and~GI1I2I3=I1~b1+I2~b2+I3~b3, withIj2Z,
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1X
~G00
I00
1I00
2I00
3a~G00
I00
1I00
2I00
3ei(~ki1i2i3~k0
i0
1i0
2i0
3)
(~TJ1J2J3+~tj1j2j3+~ i)ei(~GI1I2I3~G0
I0
1I0
2I3)
~ i
ZZZ
ei(i1i0
1
N1~b1+i2i0
2
N2~b2+i3i0
3
N3~b3)
(y1N1~ a1+y2N2~ a2+y3N3~ a3)ei((I1I0
1+I00
1)~b1+ (I2I0
2+I00
2)~b2+ (I3I0
3+I00
3)~b3)
(y1N1~ a1+y2N2~ a2+y3N3~ a3)d3~ r00=
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1X
~G00
I00
1I00
2I00
3a~G00
I00
1I00
2I00
3ei(~ki1i2i3~k0
i0
1i0
2i0
3)
(~TJ1J2J3+~tj1j2j3+~ i)ei(~GI1I2I3~G0
I0
1I0
2I3)
~ i
ZZZ
ei2((i1i0
1)y1+(i2i0
2)y2+(i3i0
3)y3)ei2((I1I0
1+I00
1)N1y1+(I2I0
2+I00
2)N2y2+(I3I0
3+I00
3)N3y3)d3~ r00=(62)
As shown before, the integrals are zero unless i0
n=inand in this case also if I00
n=I0
nIn. This means that
the matrix elements of the local potential V(L)are different from zero only when ~k0
i0
1i0
2i0
3=~ki1i2i3and~G00
I00
1I00
2I00
3=
~G0
I0
1I0
2I0
3~GI1I2I3. With this we only not just proved that the matrix elements of the potential are zero for plane
waves with the same ~kvector but also the rules of the selection of the ~Gvectors.
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1X
~G00
I00
1I00
2I00
3a~G00
I00
1I00
2I00
3ei(~ki1i2i3~k0
i0
1i0
2i0
3)
(~TJ1J2J3+~tj1j2j3+~ i)ei(~GI1I2I3~G0
I0
1I0
2I3)
~ i~ki1i2i3~k0
i0
1i0
2i0
3~G0
I0
1I0
2I0
3~GI1I2I3;~G00
I00
1I00
2I00
3=
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1a~G0
I0
1I0
2I0
3~GI1I2I3ei(~G0
I0
1I0
2I0
3~GI1I2I3)
~ i(63)
10

Since the terms in the sum in jndo not depend on these indexes, it is equal to N1N2N3.
D
~k0
i0
1i0
2i0
3;~G0
I0
1I0
2I3V(L)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN1N2N3a~G0
I0
1I0
2I0
3~GI1I2I3NatomX
i=1ei(~G0
I0
1I0
2I0
3~GI1I2I3)
~ i
=1
Vcellv(L)(~G0
I0
1I0
2I0
3~GI1I2I3)S(~G0
I0
1I0
2I0
3~GI1I2I3) (64)
in whichv(L)(~G0
I0
1I0
2I0
3~GI1I2I3) =a~G0
I0
1I0
2I0
3~GI1I2I3in the Fourier term ~G0
I0
1I0
2I0
3~GI1I2I3ofv(L)
.3.3 Nonlocal Potencial
For a crystal with primitive lattice vectors ~ aj,Natom atoms in the primitive cell with atomic positions i, and
dimensions Nj~ aj, the non-local pseudopotential is given by
V(NL)(~ ra;~ rb) =X
J12ZX
J22ZX
J32ZN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1v(NL)
i(~ ra~tj1j2j3~TJ1J2J3~ i;~ rb~tj1j2j3~TJ1J2J3~ i)(65)
where
~tj1j2j3=j1~ a1+j2~ a2+j3~ a3
~TJ1J2J3=J1N1~ a1+J2N2~ a2+J3N3~ a3 (66)
The atomic non-local pseudopotential is
v(NL)
i(~ ra;~ rb) =lmaxX
l=0lX
m=ljailmisgn(bil)hailmj (67)
where
h~ rjailmi=fil(j~ rj)Ylm(br) (68)
and the empirical form of bilis
bil=Z1
0r2drfil(r) (69)
with the special property that fil(r)2Randfil(r) = 0 forr>rc, wherercis the core radius, typically smaller than
half the interatomic distance.
The matrix elements of the non-local pseudopotential in a plane wave basis set
D
~ r~ki1i2i3;~GI1I2I3E
=1pN1N2N3Vcellei(~ki1i2i3+~GI1I2I3)
~ r(70)
with~ki1i2i3=i1
N1~b1+i2
N2~b2+i3
N3~b3, withij= 0;:::;Nj1, a vector in the primitive reciprocal unit cell and
~GI1I2I3=I1~b1+I2~b2+I3~b3withIj2Za reciprocal lattice vector, are
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3V(NL)~ki1i2i3;~GI1I2I3E
=ZZZ
crystald3~ raZZZ
crystald3~ rb1pN1N2N3Vcellei(~k0
i0
1i0
2i3+~GI0
1I0
2I0
3)
~ ra
 X
J12ZX
J22ZX
J32ZN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1v(NL)
i(~ ra~tj1j2j3~TJ1J2J3~ i;~ rb~tj1j2j3~TJ1J2J3~ i)!
1pN1N2N3Vcellei(~ki1ei2i3+~GI1I2I3)
~ rb(71)
inverting the order of integrals and some sums we have
D
~ki0
1i0
2i0
3;~GI0
1I0
2I0
3V(NL)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ZZZ
crystald3~ raei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ raZZZ
crystald3~ rbei(~ki1i2i3+~GI1I2I3)
~ rb
11

X
J12ZX
J22ZX
J32Zv(NL)
i(~ ra~tj1j2j3~TJ1J2J3~ i;~ rb~tj1j2j3~TJ1J2J3~ i) (72)
making the change in variable ~ r0=~ r~tj1j2j3~ iwe have
D
~ki0
1i0
2i0
3;~GI0
1I0
2I0
3V(NL)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ZZZ
crystal+~tj1j2j3+~ id3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
(~r0a+~tj1j2j3+~ i)
ZZZ
crystal+~tj1j2j3+~ id3~ r0
bei(~ki1i2i3+~GI1I2I3)
(~ r0
b+~tj1j2j3+~ i)X
J12ZX
J22ZX
J32Zv(NL)
i(~ r0
a;~TJ1J2J3~ r0
b~TJ1J2J3)(73)
removing the exponentials terms that do not depend on the integrand from the integrals and reordering sums and
integrals we obtain
D
~ki0
1i0
2i0
3;~GI0
1I0
2I0
3V(NL)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1
ei(~ki1i2i3~ki0
1i0
2i0
3)
~tj1j2j3ei(~ki1i2i3~ki0
1i0
2i0
3)
~ iei(~GI1I2I3~GI0
1I0
2I0
3)
~tj1j2j3ei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
crystal+~tj1j2j3+~ id3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ r0
aZZZ
crystal+~tj1j2j3+~ id3~ r0
bei(~ki1i2i3+~GI1I2I3)
~ r0
b
X
J12ZX
J22ZX
J32Zv(NL)
i(~ r0
a~TJ1J2J3;~ r0
b~TJ1J2J3) (74)
The integrand is periodic in space with the periodicity of the crystal (special values of ~ki1i2i3) if we shift both
coordinates simultaneously, and ei~GI0
1I0
2I0
3
~tj1j2j3= 1, so we can simplify
D
~ki0
1i0
2i0
3;~GI0
1I0
2I0
3V(NL)~ki1i2i3;~GI1I2I3E
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1
ei(~ki1i2i3~ki0
1i0
2i0
3)
~tj1j2j3ei(~ki1i2i3~ki0
1i0
2i0
3)
~ iei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
crystald3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ r0
aZZZ
crystald3~ r0
bei(~ki1i2i3+~GI1I2I3)
~ r0
b
X
J12ZX
J22ZX
J32Zv(NL)
i(~ r0
a~TJ1J2J3;~ r0
b~TJ1J2J3) (75)
Grouping the terms that depend on ~tj1j2j3we have
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3V(NL)~ki1i2i3;~GI1I2I3E
=NatomX
i=11
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0ei(~ki1i2i3~ki0
1i0
2i0
3)
~tj1j2j3
ei(~ki1i2i3~ki0
1i0
2i0
3)
~ iei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
crystald3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ r0
aZZZ
crystald3~ r0
bei(~ki1i2i3+~GI1I2I3)
~ r0
b
X
J12ZX
J22ZX
J32Zv(NL)
i(~ r0
a~TJ1J2J3;~ r0
b~TJ1J2J3) (76)
The terms in the sum over the jnare zero unless ~ki0
1i0
2i0
3=~ki1i2i3. This results in the whole sum being equal to
N1N2N3i1i0
1i2i0
2i3i0
3.
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3V(NL)~ki1i2i3;~GI1I2I3E
=i1i0
1i2i0
2i3i0
31
VcellNatomX
i=1ei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
crystald3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ r0
aZZZ
crystald3~ r0
bei(~ki1i2i3+~GI1I2I3)
~ r0
b
X
J12ZX
J22ZX
J32Zv(NL)
i(~ r0
a~TJ1J2J3;~ r0
b~TJ1J2J3) (77)
Again the integrand has the periodicity of the crystal, so we can shift the domain so that the region ra< rc
andrb< rcare completely inside the crystal. In that case the integrals for all Jnare null except for the case
12

J1=J2=J3= 0, furthermore in that case integrating over the crystal or over all space gives the same result. Also
we have that ei~GI1I2I3
~TJ1J2J3= 1, so
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3V(NL)~ki1i2i3;~GI1I2I3E
=1
VcellNatomX
i=1ei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
d3~ raei(~ki1i2i3+~GI0
1I0
2I0
3)
~ raZZZ
d3~ rbei(~ki1i2i3+~GI1I2I3)
~ rbv(NL)
i(~ ra;~ rb) (78)
on which we did ~ r0!~ r. Replacing now the expression for v(NL)
iwe obtain, interchanging integrals and sums
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3V(NL)~ki1i2i3;~GI1I2I3E
=1
VcellNatomX
i=1lmaxX
l=0lX
m=lsgn(bil)ei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
d3~ raei(~ki1i2i3+~GI0
1I0
2I0
3)
~ raZZZ
d3~ rbei(~ki1i2i3+~GI1I2I3)
~ rbfil(ra)Ylm(bra)fil(rb)Ylm(brb)(79)
the integrals factorize and we have
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3V(NL)~ki1i2i3;~GI1I2I3E
=1
VcellNatomX
i=1lmaxX
l=0lX
m=lsgn(bil)

ei(~ki1i2i3+~GI0
1I0
2I0
3)
~ iZZZ
d3~ raei(~ki1i2i3+~GI0
1I0
2I0
3)
~ rafil(ra)Ylm(bra)

ei(~ki1i2i3+~GI1I2I3)
~ iZZZ
d3~ rbei(~ki1i2i3+~GI1I2I3)
~ rbfil(rb)Ylm(brb)
(80)
where the quantities in parenthesis are the complex conjugate of each other. Defining for each ~ki1i2i3the matrix
A(ilm) (I1I2I3)=ei(~ki1i2i3+~GI1I2I3)
~ iZZZ
d3~ rei(~ki1i2i3+~GI1I2I3)
~ rfil(r)Ylm(br) (81)
where the multi-indices are grouped in parenthesis, we have the structure
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3V(NL)~ki1i2i3;~GI1I2I3E
=1
VcellNatomX
i=1lmaxX
l=0lX
m=lsgn(bil)A(ilm) (I0
1I0
2I0
3)A(ilm) (I1I2I3) (82)
Using
ei~k
~ r= 41X
l=0lX
m=liljl(kr)Ylm(r;r)Ylm(k;k) (83)
we find that
ZZZ
d3rei(~k+~G)
~ rfil(r)Ylm(br) =Z1
0r2drZ
0sin(r)drZ2
0dr41X
l0=0l0X
m0=l0il0jl0(j~k+~Gjr)
Yl0m0(r;r)Yl0m0(k+G;k+G)fil(r)Ylm(r;r)
= 4ilYlm(k+G;k+G)Z1
0jl(j~k+~Gjr)fil(r)r2dr (84)
Where we used the fact that the Spherical Harmonic functions Ylmare orthonormal in a shpere of constant
radius, beingRR

YlmYl0m0d
=ll0mm0.
Therefore we have the final result
A(ilm) (I1I2I3)= 4ilei(~ki1i2i3+~GI1I2I3)
~ iYlm(k+G;k+G)Z1
0jl(j~k+~Gjr)fil(r)r2dr: (85)
Notice that we have two i, one is the sum index and the other one is the complex i. If you can’t tell which one
is which you shouldn’t be reading this thesis.
13

.3.4 Spin-Orbit contribution
For a crystal with primitive lattice vectors ~ aj,Natom atoms in the primitive cell with atomic positions i, and
dimensions Nj~ aj, the spin-orbit contribution to the pseudopotential is given by
V(SO)(~ ra;~ rb) =X
J12ZX
J22ZX
J32ZN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1v(SO)
i(~ ra~tj1j2j3~TJ1J2J3~ i;~ rb~tj1j2j3~TJ1J2J3~ i)(86)
where
~tj1j2j3=j1~ a1+j2~ a2+j3~ a3
~TJ1J2J3=J1N1~ a1+J2N2~ a2+J3N3~ a3 (87)
and the atomic spin-orbit pseudopotential is
v(SO)
i(~ ra;~ rb) =lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=jailjmj
sgn(bilj)
ailjmj (88)
where
ailjmj
=0
@ailjmj+1
2E
ailjmj1
2E1
A=ailjmj+1
2E
j"i+ailjmj1
2E
j#i (89)
,

~ railjmjms
=filj(j~ rj)Ylmjms(br) (90)
and
bilj=Z1
0r2drfil(r) (91)
with the special property that filj(r)2Randfilj(r) = 0 forr > rc, wherercis the core radius, typically smaller
than half the interatomic distance.
In this case we have to consider the spin in the plane wave basis set:
~ki1i2i3;~GI1I2I3E
=0
@~ki1i2i3;~GI1I2I3;+1
2E
~ki1i2i3;~GI1I2I3;1
2E1
A=~ki1i2i3;~GI1I2I3;+1
2
j"i+~ki1i2i3;~GI1I2I3;1
2
j#i (92)
where
~ r~ki1i2i3;~GI1I2I3;1
2
=1pN1N2N3Vcellei(~ki1i2i3+~GI1I2I3)
~ r(93)
with~ki1i2i3=i1
N1~b1+i2
N2~b2+i3
N3~b3, withij= 0;:::;Nj1, a vector in the primitive reciprocal unit cell and
~GI1I2I3=I1~b1+I2~b2+I3~b3withIj2Za reciprocal lattice vector, are, for two plane waves with spin up:

~k;~G;+1
2V(SO)~k0;~G0;+1
2
=X
~TX
~tNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=j
~k0;~G0;+1
2ailjmj
sgn(bilj)
ailjmj~k;~G;+1
2
=
=X
~TX
~tNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=j
~k0;~G0;+1
2ailjmj+1
2
h"j"i +
~k0;~G0;+1
2ailjmj1
2
h"j#i
sign(bilj)

ailjmj+1
2~k;~G;+1
2
h"j"i +
ailjmj1
2~k;~G;+1
2
h#j"i
=
=X
~TX
~tNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=j
~k0;~G0;+1
2ailjmj+1
2
sign(bilj)
ailjmj+1
2~k;~G;+1
2
(94)
The same way goes with two plane waves with spin down:
14

~k0;~G0;1
2V(SO)~k;~G;1
2
=X
~TX
~tNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=j
~k0;~G0;1
2ailjmj1
2
sign(bilj)
ailjmj1
2~k;~G;1
2
(95)
For the other matrix elements:

~k0;~G0;+1
2V(SO)~k;~G;1
2
=X
~TX
~tNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=j
~k0;~G0;+1
2ailjmj
sgn(bilj)
ailjmj~k;~G;1
2
=
=X
~TX
~tNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=j
~k0;~G0;+1
2ailjmj+1
2
h"j"i +
~k0;~G0;+1
2ailjmj1
2
h"j#i
sign(bilj)

ailjmj+1
2~k;~G;1
2
h"j#i +
ailjmj1
2~k;~G;1
2
h#j#i
=
=X
~TX
~tNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=j
~k0;~G0;+1
2ailjmj+1
2
sign(bilj)
ailjmj1
2~k;~G;1
2
(96)
and

~k0;~G0;1
2V(SO)~k;~G;+1
2
=X
~TX
~tNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=j
~k0;~G0;1
2ailjmj1
2
sign(bilj)
ailjmj+1
2~k;~G;+1
2
(97)
So, in conclusion:
D
~k0;~G0;m0
sV(SO)~k;~G;msE
=X
~TX
~tNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=jD
~k0;~G0;m0
sailjmjm0sE
sign(bilj)D
ailjmjms~k;~G;msE
(98)
We calculate the matrix elements like in the previous way, including the spin this time:
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=ZZZ
crystald3~ raZZZ
crystald3~ rb1pN1N2N3Vcellei(~k0
i0
1i0
2i3+~GI0
1I0
2I0
3)
~ ra
 X
J12ZX
J22ZX
J32ZN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1v(SO)
i(~ ra~tj1j2j3~TJ1J2J3~ i;~ rb~tj1j2j3~TJ1J2J3~ i)!
1pN1N2N3Vcellei(~ki1ei2i3+~GI1I2I3)
~ rb(99)
inverting the order of integrals and some sums we have
D
~ki0
1i0
2i0
3;~GI0
1I0
2I0
3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ZZZ
crystald3~ raei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ raZZZ
crystald3~ rbei(~ki1i2i3+~GI1I2I3)
~ rb
X
J12ZX
J22ZX
J32Zv(SO)
i(~ ra~tj1j2j3~TJ1J2J3~ i;~ rb~tj1j2j3~TJ1J2J3~ i)(100)
making the change in variable ~ r0=~ r~tj1j2j3~ iwe have
D
~ki0
1i0
2i0
3;~GI0
1I0
2I0
3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1ZZZ
crystal+~tj1j2j3+~ id3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
(~ r0
a+~tj1j2j3+~ i)
ZZZ
crystal+~tj1j2j3+~ id3~ r0
bei(~ki1i2i3+~GI1I2I3)
(~ r0
b+~tj1j2j3+~ i)X
J12ZX
J22ZX
J32Zv(SO)
i(~ r0
a;~TJ1J2J3~ r0
b~TJ1J2J3)(101)
15

removing the exponentials terms that do not depend on the integrand from the integrals and reordering sums and
integrals we obtain
D
~ki0
1i0
2i0
3;~GI0
1I0
2I0
3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1
ei(~ki1i2i3~ki0
1i0
2i0
3)
~tj1j2j3ei(~ki1i2i3~ki0
1i0
2i0
3)
~ iei(~GI1I2I3~GI0
1I0
2I0
3)
~tj1j2j3ei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
crystal+~tj1j2j3+~ id3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ r0
aZZZ
crystal+~tj1j2j3+~ id3~ r0
bei(~ki1i2i3+~GI1I2I3)
~ r0
b
X
J12ZX
J22ZX
J32Zv(SO)
i(~ r0
a~TJ1J2J3;~ r0
b~TJ1J2J3) (102)
The integrand is periodic in space with the periodicity of the crystal (special values of ~ki1i2i3) if we shift both
coordinates simultaneously, and ei~GI0
1I0
2I0
3
~tj1j2j3= 1, so we can simplify
D
~ki0
1i0
2i0
3;~GI0
1I0
2I0
3m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=1
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0NatomX
i=1
ei(~ki1i2i3~ki0
1i0
2i0
3)
~tj1j2j3ei(~ki1i2i3~ki0
1i0
2i0
3)
~ iei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
crystald3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ r0
aZZZ
crystald3~ r0
bei(~ki1i2i3+~GI1I2I3)
~ r0
b
X
J12ZX
J22ZX
J32Zv(SO)
i(~ r0
a~TJ1J2J3;~ r0
b~TJ1J2J3)(103)
Grouping the terms that depend on ~tj1j2j3we have
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=NatomX
i=11
N1N2N3VcellN11X
j1=0N21X
j2=0N31X
j3=0ei(~ki1i2i3~ki0
1i0
2i0
3)
~tj1j2j3
ei(~ki1i2i3~ki0
1i0
2i0
3)
~ iei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
crystald3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ r0
aZZZ
crystald3~ r0
bei(~ki1i2i3+~GI1I2I3)
~ r0
b
X
J12ZX
J22ZX
J32Zv(SO)
i(~ r0
a~TJ1J2J3;~ r0
b~TJ1J2J3) (104)
The terms in the sum over the jnare zero unless ~ki0
1i0
2i0
3=~ki1i2i3. This results in the whole sum being equal to
N1N2N3i1i0
1i2i0
2i3i0
3.
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=i1i0
1i2i0
2i3i0
31
VcellNatomX
i=1ei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
crystald3~ r0
aei(~ki0
1i0
2i0
3+~GI0
1I0
2I0
3)
~ r0
aZZZ
crystald3~ r0
bei(~ki1i2i3+~GI1I2I3)
~ r0
b
X
J12ZX
J22ZX
J32Zv(SO)
i(~ r0
a~TJ1J2J3;~ r0
b~TJ1J2J3)(105)
Again the integrand has the periodicity of the crystal, so we can shift the domain so that the region ra< rc
andrb< rcare completely inside the crystal. In that case the integrals for all Jnare null except for the case
J1=J2=J3= 0, furthermore in that case integrating over the crystal or over all space gives the same result. Also
we have that ei~GI1I2I3
~TJ1J2J3= 1, so
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=1
VcellNatomX
i=1ei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
ZZZ
d3~ raei(~ki1i2i3+~GI0
1I0
2I0
3)
~ raZZZ
d3~ rbei(~ki1i2i3+~GI1I2I3)
~ rbv(SO)
i(~ ra;~ rb) (106)
n which we did ~ r0>~ r. Replacing now the expression for v(SO)
iwe obtain, interchanging integrals and sums
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=1
VcellNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=jsgn(bilj)ei(~GI1I2I3~GI0
1I0
2I0
3)
~ i
16

ZZZ
d3~ raei(~ki1i2i3+~GI0
1I0
2I0
3)
~ raZZZ
d3~ rbei(~ki1i2i3+~GI1I2I3)
~ rbfilj(ra)Ylmjm0s(bra)filj(rb)Ylmjms(brb)(107)
the integrals factorize and we have
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=1
VcellNatomX
i=1lmaxX
l=0l+1
2X
j=jl1
2jjX
mj=jsgn(bilj)

ei(~ki1i2i3+~GI0
1I0
2I0
3)
~ iZZZ
d3~ raei(~ki1i2i3+~GI0
1I0
2I0
3)
~ rafilj(ra)Ylmjm0s(bra)

ei(~ki1i2i3+~GI1I2I3)
~ iZZZ
d3~ rbei(~ki1i2i3+~GI1I2I3)
~ rbfilj(rb)Ylmjms(brb)
(108)
where the quantities in parenthesis are the complex conjugate of each other. Defining for each ~ki1i2i3the matrix
A(iljmjms) (I1I2I3)=ei(~ki1i2i3+~GI1I2I3)
~ iZZZ
d3~ rei(~ki1i2i3+~GI1I2I3)
~ rfilj(r)Ylmjms(br) (109)
where the multi-indices are grouped in parenthesis, we have the structure
D
~ki0
1i0
2i0
3;~G0
I0
1I0
2I3;m0
sV(SO)~ki1i2i3;~GI1I2I3;msE
=1
VcellNatomX
i=1lmaxX
l=0lX
m=lsgn(bilj)A(iljmjms) (I0
1I0
2I0
3)A(iljmjms) (I1I2I3)
(110)
Using
ei~k
~ r= 41X
l=0lX
m=liljl(~k
r)Ylm(r;r)Ylm(k;k) (111)
we find that
ZZZ
d3~ rei(~k+~G)
~ rfilj(r)Ylm(br) =Z1
0r2drZ
0sin(r)drZ2
0dr41X
l0=0l0X
m0=l0il0jl0(j~k+~Gjr)
Yl0m0(r;r)Yl0m0(k;k)fijl(r)Ylm(r;r)
= 4ilYlm(k;k)Z1
0jl(j~k+~Gjr)filj(r)r2dr (112)
Where we used the fact that the Spherical Harmonic functions Ylmare orthonormal in a shpere of constant
radius, beingRR

YlmYl0m0d
=ll0mm0andm=ml=mjms.
Therefore we have that
A(iljmjms) (I1I2I3)= 4ilei(~ki1i2i3+~GI1I2I3)
~ iYlmjms(k;k)Z1
0jl(j~k+~Gjr)filj(r)r2dr (113)
17

18