Lattice dynamics : From fundamental research to practical applications

Lattice dynamics

Linköping Studies in Science and Technology Dissertation No. 2103

Johan Tidholm

Jo

han

T

idh

olm

La

ttic

e d

yn

am

ic

s: F

ro

m f

un

da

m

en

tal re

se

arc

h to p

ra

ctic

al a

pp

lic

atio

ns

20

FACULTY OF SCIENCE AND ENGINEERING

Linköping Studies in Science and Technology, Dissertation No. 2103, 2020 Department of Physics, Chemistry and Biology (IFM)

Linköping University SE-581 83 Linköping, Sweden

www.liu.se

From fundamental research

to practical applications

Linköping Studies in Science and Technology Disserta ons, No. 2103

La ce dynamics

From fundamental research to prac cal applica ons

Johan Tidholm

Linköping University

Department of Physics, Chemistry and Biology Division of theore cal physics

SE-581 83 Linköping, Sweden Linköping 2020

© Johan Tidholm, 2020 ISBN 978-91-7929-759-6 ISSN 0345-7524

URL http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-170800

Published articles have been reprinted with permission from the respective copyright holder.

Typeset using XƎTEX

POPULÄRVETENSKAPLIG SAMMANFATTNING

Den här avhandlingen handlar om beräkningar för material. När materialberäkningar utförs är det antingen för att förutsäga eller förklara egenskaper. De beräkningar som jag har gjort i denna avhandling är baserade på fundamentala fysiska lagar. Detta betyder att de är rent baserade på teori, och inte har anpassats efter resultat av experiment.

Jag har i mitt arbete använt mig mycket utav en teori som kallas gitter dynamik. Den är definierad för periodiska material, det vill säga att atomerna i dessa material upprepas i periodiska mönster. Vi kan då anta att det finns en jämviktspunkt för alla atomerna, som de vibrerar omkring. Dessa vibrationer kan beskrivas som om atomerna påverkar varandra med fiktiva fjädrar. Genom att beräkna styrkan för dessa fjädrar kan vi beskriva vibrationerna av atomerna. Dessa vibrationer i sin tur är avgörande för materialets egenskaper. För att beskriva ett material vid en specifik temperatur har jag använt mig utav olika metoder för att simulera det. En simulering kan ses som ett “dator experiment”. Problemet är dock hur vi ska mäta egenskaperna i simuleringen. Ju större och mera komplex en simulering är, desto svårare blir det att beräkna egenskaperna av det simulerade materialet. Vi hamnar i en situation likt den vi skulle befinna oss om vi hade gjort ett experiment i verkligheten, och tvingas använda förenklade modeler för att kunna tolka resultatet. Jag har därför använt mig utav metoder för att utvinna vibrationer av atomer, elektrontillstånd eller elastiska egenskaper, specifikt utvecklade för att användas på denna typ utav simuleringar. Mitt arbete har kretsat kring hur dessa egenskaper påverkas av extrema temperaturer och tryck.

De beräkningar jag har utfört vid höga tryck har varit för nyupptäckta nitrider och fa-ser av kiseldioxid. Nitriderna är porösa material som innehåller en stor mängd kväve. Det höga kväveinehållet gör så att det lagras en stor mängd kemisk energi i enkel- och dubbel-bindningar mellan kväveatomerna. De nya faserna av kiseldioxid har en betydelse för vår förståelse av jordens inre. Deras existens öppnar upp för att det kan finnas mera komplexa och ihoptryckbara flytande material, under jordens nedre mantel, än vad tidigare har varit antaget. Mina beräkningar har bekräftat strukturerna för dessa nyupptäckta material. Vid höga temperaturer har jag studerat för metallen niob hur vibrationerna av atomerna är relaterade till olika elektrontillstånd. För specifika vibrationer ökar frekvensen med ökad temperatur. Detta är något ovanligt eftersom vibrationernas frekvenser vanligtvis brukar minska med ökad temperatur. Mina simulering för denna metal överensstämmer med resul-tat från experiment. Orsaken till varför visa vibrationers frekvenser ökar kan jag förklara med att elektrontillståndens enskilda energier varierar över tid på grund av den ökade tem-peraturen. Jag har även använt mig av simuleringar för att beräkna elastiska egenskaper av legeringen Ti0.5Al0.5N. Ti1−xAlxN legeringar används som beläggningar på skärverktyg

som används för metall. För att öka effektiviteten av beläggningen, behövs det detaljerad kunskap av dess mekaniska egenskaper för den temperatur som de används vid. Jag be-räknade därför så noggrant som möjligt de elastiska egenskaperna för Ti0.5Al0.5N. Dessa beräkningar är avsedda för att användas som en referens för andra beräkningsmässigt bil-ligare metoder. Datan som genererades från mina simuleringar användes även för en sådan metod, baserad på maskininlärning.

The reason to perform calculations in material science usually falls into one of two categories: to predict or explain the origin of material properties. This thesis covers first-principle cal-culations for solids at extreme conditions, from both of the two mentioned categories. I primarily have studied the effects of high-pressure and high-temperature on lattice dynam-ics, mechanical and electronic properties. To treat the effects of temperature, ab initio molecular dynamics (AIMD) simulations and self-consistent phonon calculations, based on density functional theory, have been utilised. These approaches account for the temperature effects by considering thermally excited supercells as samples of a statistical ensemble. To extract properties from this representation, I have used methods which maps the supercell data to a unit cell representation or fits it to a simple model Hamiltonian.

The small displacement method was used to analyse the dynamical stability for nitrides and polymorphs of silica, synthesised at high-pressure in a diamond anvil cell. The nitride com-pounds consist of a high amount of nitrogen either as chains, forming a porous framework together with transition metal atoms or as dinitrogen molecules, occupying the channels of the framework. The nitrogen chains consist of single- or double-bonded nitrogen atoms, making these compounds highly energetic. Polymorphs of silica can be used to model deep Earth liquids. These new polymorphs, named coesite-IV and coesite-V, consist of four-, five-, and six-oriented silicon. Some of the octahedra of the six-oriented silicon atoms, of these new phases, are sharing faces, which according to Pauling’s third rule would make them highly unstable. My phonon calculations indicate these phases to be dynamically stable. Furthermore, my calculations predict higher compressibility for these new phases compared to the competing ones. By modelling silicate melts with IV and coesite-V, a more complex and compressible structure is expected, affecting the predicted seismic behaviour.

I studied Kohn anomalies for body-centered cubic niobium by simulating this material with self-consistent phonon calculations. The electronic structure was studied by using a band unfolding technique, for which I obtained an effective unit cell representation of the elec-tronic structure at elevated temperatures. Temperature primarily smeared the elecelec-tronic states but did not induce significant shifts of the bands. In parallel, the anharmonicity of this system was studied using the temperature dependent effective potential method. Even close to the melting temperature, this element is remarkably harmonic. The experimentally observed disappearance of the Kohn anomalies with increased temperature is predominantly dependent, according to my calculations, on the temperature-induced smearing of the elec-tronic states.

Using stress-strain relations, accurate high-temperature elastic properties were predicted for Ti0.5Al0.5N. The simulations were performed with AIMD. The stresses were fitted us-ing the least-squares method to a linear expression from which the elastic constants were derived. The results were compared with previously performed calculations that employed additional approximations. The results of the symmetry imposed force constant tempera-ture dependent effective potential (SIFC-TDEP) method agrees well with our results. I also compared my results with TiN calculations that employed a similar methodology. My and the SIFC-TDEP results are reporting lower values for the polycrystalline moduli than the calculations for TiN. The data I generated were also used for a machine learned interatomic potential method, where moment tensor potentials were trained and evaluated, using this data.

Acknowledgments

To everyone, who has knowingly or unknowingly supported, inspired, or in any other way helped me, thank you. With that said, I feel like I need to especially mention some of you. To start with, I would like to thank my supervisor Igor Abrikosov, for his guidance, encouragement, and tireless feedback. To my first co-supervisor Olle Hellman and my favourite collaborator Nina Shulumba, thank you for making this time such an enjoyable adventure. I also have to thank Ferenc Tasnádi, my second co-supervisor, for his enthusiastic support and many valuable, informative, and amusing discussions and conversations. To the members of the Theory of Disordered Materials unit and its leader Björn Alling, thank you for your shared experience, knowledge, and feedback, especially during this last period of time. And to my remaining colleagues and friends from Linköping University, thank you for all the work related stuff and for things like after-works, travel adventures, culture exchange dinners, hitting the gym, and derailed gags.

Of all the comrades that e’er I had, I am also very thankful to. The ones from LRK, TLBC, and old times, thank you for all the time we spent together that have and have not become epic stories. To my family, I wish to express my sincere gratitude for everything they have done for me. And finally to my Marie, thank you, this would have been so much more difficult without you.

Contents

2 First principle calculations of electronic structure for solids 5

2.1 The many-body problem . . . 5

2.2 The Born-Oppenheimer approximation . . . 6

2.3 Density functional theory . . . 7

2.4 The Kohn-Sham auxiliary system . . . 8

2.5 Practicalities of performing Kohn-Sham DFT for solids . . . 10

3 Lattice dynamics 15 3.1 Force constants and the dynamical matrix . . . 15

3.2 Quantum mechanical description . . . 17

3.3 Vibrational free energy . . . 18

3.4 Quasi-harmonic approximation (QHA) . . . 19

3.5 The small displacement method . . . 19

3.6 Interacting phonons . . . 20

4 Simulating solids as thermally excited supercells 23 4.1 Ab initio molecular dynamics (AIMD) . . . 24

4.2 Self-consistent phonon (SCP) calculations . . . 25

4.3 Machine learned moment tensors potentials (ML-MTP) . . . 26

5 Calculating properties from thermally excited supercells 29 5.1 Electronic structure unfolding . . . 29

lations . . . 34

6 Results 37

6.1 Dynamical stability of nitrides synthesised at high-pressure . . 38

6.2 Theoretical investigation of high-pressure phases of silica . . . . 44

6.3 The temperature dependence of Kohn anomalies in bcc Nb . . 49

6.4 Temperature dependent elastic constants for cubic Ti0.5Al0.5N 54

7 Conclusion 61

Bibliography 63

Articles included in thesis 75

Paper I 79 Paper II 87 Paper III 97 Paper IV 105 Paper V 115 Paper VI 135

List of Figures

1.1 Phonon dispersion relation (a) and phonon spectral function (b) for

body-centered cubic cubic (bcc) niobium (Nb) at 2700 K. These figures represents the vibrations within a material. By explicitly treating the temperature, details are changed between these two representations. For example a “kink” between Γ and H on the

lower line is not as distinct in (b) as in (a). . . 2

1.2 Electron band structure (a) and Bloch spectral function (b) for

body-centered cubic (bcc) niobium (Nb) at 2700 K. These figures represents the electronic states of a material. By explicitly treating the temperature, featurs are changed between these two represen-tations. For example the three states at the Fermi energy (dashed

line) between Γ and N in (a) are indistinguishable in (b). . . . 3

2.1 Flowchart of the different steps in the iterative loop of Kohn-Sham

DFT. . . 11

5.1 Electron band structure (a) and Bloch spectral function (b) for

body-centered cubic (bcc) niobium (Nb). Note that the electron band structure in (a) is plotted in the first Brillouin zone (BZ) of the supercell. The smearing of the Bloch spectral function in (b)

corresponds to a temperature of 0 K. . . 30

5.2 Convergence test for the number of configurations needed in a set

to obtain consistent phonon dispersion relations for body-centered cubic (bcc) niobium (Nb). 50 calculations were performed for vari-ating sets with the same number of configurations. The sets were generated randomly from a pool of approximately 1000 configura-tions. The number of configurations in the sets were 1, 2, 4, 8, 16, 32, 64, and 128. The broadening for this figure is an estimate of the statistical error and should not be confused with the linewidths, in-verse phonon life time, of phonon spectral functions. Convergence is reached for this system when the sets consist of 8 configurations. 33

displacement method. In the crystal structure, the framework con-sists of rhenium atoms and nitrogen chains. Dinitrogen molecules are occupying the space inside the framework. No imaginary fre-quencies are observed in the dispersion relation, stating that this

crystal structure is dynamically stable. . . 38

6.2 Crystal structure and phonon dispersion relation for Os5N26⋅3N2

at 108 GPa. The phonon calculation was performed with the small displacement method. In the crystal structure, the framework con-sists of osmium atoms and nitrogen chains. Dinitrogen molecules are occupying the space inside the framework. No imaginary fre-quencies are observed in the dispersion relation, stating that this

crystal structure is dynamically stable. . . 39

6.3 Crystal structure and phonon dispersion relation for Hf4N20⋅N2at

92 GPa. The phonon calculation was performed with the small displacement method. In the crystal structure, the framework con-sists of hafnium atoms and nitrogen chains. Dinitrogen molecules are occupying the space inside the framework. Imaginary frequen-cies are observed in the dispersion relation, stating that this crystal

structure is dynamically unstable. . . 40

6.4 Phonon dispersion relations for WN8⋅N2 at the pressures 103, 44,

9, and 0 GPa. The crystal structure for this compound is the same

as for ReN8⋅N2, displayed in Fig. 6.1. These phonon calculations

were performed with the small displacement method. With de-creasing pressure a softening is observed, particularly noticeable for the acoustic branches. At 0 GPa small imaginary frequencies

are observed. . . 41

6.5 Electronic density of states (DOS) for WN8⋅N2 at the pressures

103, 44, 9, and 0 GPa. The electronic states were calculated with Kohn-Sham DFT. A growing peak at the Fermi energy (the dashed

vertical line at 0 eV) is observed at decreased pressure. . . 42

6.6 The ferromagnetic results of WN8⋅N2at ambient pressure for

elec-tronic density of states (DOS) (a) and phonon dispersion relation (b). The electronic structure and Hellman-Feynman forces were calculated with Kohn-Sham collinear spin DFT. The dispersion relation was calculated with small displacement method. In (a) a split of the peak originally observed in Fig. 6.5 for ambient pressure is seen. The corresponding acoustic phonon branch in (b), which had a small imaginary value in Fig. 6.4 is real. Stating dynamical

stability for this compound at ambient pressure. . . 43

6.7 Common fragment of the crystal structure for IV and

6.8 Crystal structure and phonon dispersion relation for coesite-IV at 40 GPA. The dispersion relation was calculated with the small dis-placement method. This polymorph of silica is dynamically stable

due to its phonon frequencies are all real. . . 45

6.9 Crystal structure and phonon dispersion relation for coesite-V at

57 GPA. The dispersion relation was calculated with the small dis-placement method. This polymorph of silica is dynamically stable

due to its phonon frequencies are all real. . . 47

6.10 The enthalpy differences (a) for seifertite, IV, and coesite-V in comparison to stishovite. A conspicuously large difference around 0.4 eV/ion in enthalpy is noticed for IV and coesite-V compared with stishovite. The enthalpy differences (b) only for coesite-IV and coesite-V compared with stishovite. Above

approx-imately 38 GPa, coesite-V has a lower enthalpy than coesite-IV. . 48

6.11 Temperature dependence of the lattice parameter, a, for bcc Nb. It is expressed as the percentage difference from a reference value,

(a − a0)/a0. The experimental results are from Miiller and

Cezair-liyan [83]. For both the temperature dependent effective potential

(TDEP) and the quasi-harmonic approximation (QHA) results, a0

was 3.33 Å at 300 K. The experimental a0 was taken at 293 K.

The second-order TDEP results agrees well with the experiments. In contrast to the QHA which underestimates the thermal

expan-sion and predicts a too linear dependence of the temperature. . . . 49

6.12 Phonon dispersion relations for bcc Nb at the temperatures 300 and 1030 K. The calculations were performed with the temper-ature dependent effective potential (TDEP) method, utilising a model truncated to the second-order, and the quasi-harmonic ap-proximation (QHA). The experimental values were obtained from inelastic neutron scattering by Powell, Martel, and Woods [79, 82].

The nesting vectors qa and qb are also drawn in the figure,

dis-playing an agreement between them and the “kinks”. At 300 K both of the computational approaches and the experiments agree.

At 1030 K, the QHA still predicts a kink at qa, which according

to experiments and the TDEP method should disappear. . . 50

6.13 Phonon spectral functions for bcc Nb at the temperatures 300 and 1200 K. With increasing temperature, the linewidth, inverse phonon lifetime, is not considerable affected. Nor are they notable shifted compared with the phonon dispersion relations in Fig. 6.12. This indicates a small anharmonic contribution, eliminating the phonon-phonon interactions as a significant contributor for the disappearance of the Kohn anomalies with increased temperatures. 51

perature. The linewidth of them though, is increased considerably with increased temperature. In Fig. 6.15, the states at the Fermi energy between the Γ and the N point are related to the nesting

vectors qa and qb. At 1200 K, these states are smeared beyond

recognition. Because they are smeared into one diffuse object, the effective Fermi surface is changed. Removing the parallel segments

witch was nested by qa and qb. Explaining why the Kohn

anoma-lies related to qa and qb disappears. . . 52

6.15 Phonon dispersion relations, electronic band structure, and(100)

cross section of the Fermi surface for bcc Nb at 0 K. In the figure

the nesting vectors, qa and qb are drawn both for the dispersion

relations and the Fermi surface cross section. The nesting vectors are related to segments of the electron band structure, labeled as “1”, “2”, and “3”. The “kinks” in the phonon dispersion relation coincide with the nesting vectors. The segments 1, 2, and 3 are de-noted both in the electronic band structure and the Fermi surface

cross section. . . 53

6.16 Temperature dependency of closet cubic symmetry projected

elas-tic constants, ¯C11, ¯C12, and ¯C44, for cubic Ti0.5Al0.5N. The solid

line is our results. The dashed line is the symmetry imposed force constant temperature dependent effective potential (SIFC-TDEP) results. The dash-dotted line is the result of the volume scaled (VS) calculations. Our mean values are denoted with an asterisk and the error bars displays a confidence interval of 95%. The SIFC-TDEP

results are shifted compared with ours for ¯C12 and ¯C44, but do

generally agree. . . 55

6.17 The temperature dependent trends of the closet cubic

symme-try projected elastic constants, ¯C11, ¯C12, and ¯C44, for cubic

Ti0.5Al0.5N. The solid line is our results. The dashed line is the

the symmetry imposed force constant temperature dependent ef-fective potential (SIFC-TDEP) results. The dash-dotted line is the result of the volume scaled (VS) calculations. Our mean values are denoted with an asterisk and the error bars displays a con-fidence interval of 95%. Again, the agreement between ours and the SIFC-TDEP results are shown. It also becomes more evident in this figure that the SIFC-TDEP method is underestimating the

temperature dependence of ¯C11 and overestimates it for ¯C12 and

¯

6.18 Temperature dependency of the Hill approximated bulk (BH),

shear (GH), and Young’s modulus (EH) for Ti0.5Al0.5N and TiN.

Our Ti0.5Al0.5N results are the solid blue lines with cubes. The the

symmetry imposed force constant temperature dependent

effec-tive potential (SIFC-TDEP) results for Ti0.5Al0.5N are the dashed

lines. The TiN results [87] are purple lines with cubes. For TiN, additional values at 1200 and 1800 K are marked out. A good agree-ment is shown between the SIFC-TDEP results and ours. Both of

the calculations for Ti0.5Al0.5N are displaying lower values for all

the polycrystalline moduli than TiN. . . 57

6.19 Temperature dependency of directional Young’s modulus for

Ti0.5Al0.5N. With increasing temperature, and increase in

anisotropy is predicted from our results the volume scaled (VS) calculations. The symmetry imposed force constant temperature dependent effective potential (SIFC-TDEP) results are display-ing a temperature independent anisotropy. From this, together

with the underestimation of the temperature dependence of ¯C11

and overestimation of ¯C12 and ¯C44, the conclusion that the

SIFC-TDEP method is averaging the temperature-induced effects on

List of Tables

6.1 Experimental reported (exp.) and DFT optimised (opt.) crystal

structures for coesite-IV and coesite-V. Both the experimental and optimised values agree except for coesite-IV at 49 GPa, which has during the optimisation become the crystal structure for coesite-V. 46

6.2 Bulk modulus and pressure derivate for stishovite, seifertite,

coesite-IV, and coesite-V. The bulk modulus for coesite-IV and

Chapter

1

Introduction

Theoretical calculations are invaluable resources for material science. They can be used either to predict or explain the behaviour of materials. This field is known as computational material science. The calculations I have per-formed for this field are on a microscopic scale, where quantum mechanical effects can not be ignored. The problem is that it is impossible to perform calculations where the entire system is fully treated quantum mechanically, except for the simplest of systems. To make these calculations possible, physi-cal motivated approximations are needed. In my work, related to this thesis, I have performed calculations both to predict and explain material behaviours at extreme conditions. This has either been for high temperatures or high pressures.

An approach to represent materials computationally at certain tempera-ture is to simulate them. A simulation can be seen as an in silico experiment, where, unlike a real experiment, all observables are accessible. But with grow-ing complexity, accessgrow-ing the observables becomes more challenggrow-ing. We find ourselves in a situation not too unfamiliar to an experimentalist, where we are forced to apply simplified models to interpret the results. The effect of explicitly treating the temperature in calculations on vibrational and elec-tronic properties is illustrated in Figs. 1.1 and 1.2. In (a) for both figures, the temperature is not explicitly treated, as it is in (b). These figures are for body-centered cubic (bcc) niobium (Nb) at 2700 K. They may seem to be similar, but important details differ between (a) and (b), changing the macro-scopic properties. This demonstrates the significance of how the temperature is treated in a calculation, affecting its result. Which is vital for when mate-rial properties are calculated for industmate-rial applications with high operating temperatures, like cutting tools.

1 2 3 4 5 6 Γ H N Γ P H 1 2 3 4 5 6 Fr eq ue nc y (T H Z) Fr eq ue nc y (T H Z) (a) (b) A “kink” No “kink”

Figure 1.1: Phonon dispersion relation (a) and phonon spectral function (b) for body-centered cubic cubic (bcc) niobium (Nb) at 2700 K. These figures represents the vibrations within a material. By explicitly treating the tem-perature, details are changed between these two representations. For example a “kink” between Γ and H on the lower line is not as distinct in (b) as in (a).

The calculations I have performed for materials at high pressures, have been for newly discovered nitrides and polymorphs of silica. For these mate-rials I have investigated the experimental reported crystal structure. Which essentially has been to evaluate if the reported structure is stable, or if it would fall apart. The implications of the existence of the nitrides would be the discovery of new nitrogen rich materials which could store a substantial amount of energy. For the new polymorphs of silica, the relevance of the study is that the silicate melts at the base of the Earth’s mantle, could have a structures and properties not previously considered.

For my work related to temperature, I have investigated a relationship between the vibrations in a solid and its electronic structure, namely Kohn anomalies [1]. This was done for bcc Nb. In (a) of Fig. 1.1, between Γ and H, a Kohn anomaly is seen as a “kink” on the lower line. Which is not present in (b) when the temperature is explicitly treated. The question to answer is why it disappeared.

To increase the hardness and wear resistance of cutting tools, they are

coated with Ti1−xAlxN. To optimise the performance for these tools,

accu-rate and detailed knowledge of their coatings mechanical properties at their working temperature is needed. I have therefore calculated temperature

previ--2 -1 0 1 2 Γ H N Γ P H -2 -1 0 1 2 En er gy (e V) En er gy (e V) (a) (b) Three clearly disctinct states Indistinguishable states

Figure 1.2: Electron band structure (a) and Bloch spectral function (b) for body-centered cubic (bcc) niobium (Nb) at 2700 K. These figures represents the electronic states of a material. By explicitly treating the temperature, fea-turs are changed between these two representations. For example the three states at the Fermi energy (dashed line) between Γ and N in (a) are indistin-guishable in (b).

ous calculations, which are computationally cheaper since they rely on more approximations. The data I generated was also used to machine learn in-teratomic potentials, which reduce the computational cost significant when they are used to simulate this particular system. If the cheaper approaches are accurate enough, the computational cost for calculating the temperature

dependent elastic constants for a wide span of compositions of Ti1−xAlxN,

would be drastically reduced.

This thesis is build on the foundation of first principles calculations. For which the theory is introduced in chapter 2. In chapter 3, I introduce the centrepiece of this thesis, lattice dynamics. Chapter 4 deals with the different approaches I have been in contact with to simulate materials. I continue in chapter 5 with the methods I have used to calculate material properties from these simulations. The results of my work, related to this thesis and the papers it includes, are presented in chapter 6. Chapter 7 is the conclusion of this work. In the following pages, I account for my contribution to the included papers of this thesis, followed by the papers themselves.

Chapter

2

First principle calculations

of electronic structure for

solids

It is my intention in this chapter to introduce the theory used throughout my work to calculate the electronic structure of solids. I will show in this chapter why we initially are describing the nuclei as static. The theory is said to be from first principles, equally referred to as ab initio, which means that it is based only upon basic physical laws. The inputs in these calculations are the atomic numbers and positions of the atoms. The focus of this chapter will be on how the electronic problem can be solved by describing the electrons by their density.

2.1 The many-body problem

To fully describe a system, we need to find a Hamiltonian for it. In the case

of solids, the system consist of Ne electrons and Nn nuclei that interacts. If

the system is considered to be non-relativistic, we can write the Hamiltonian operator as ˆ H = ˆTn+ ˆTe+ ˆVnn+ ˆVee+ ˆVne, (2.1) where ˆ Tn= − Nn ∑ i ̵h2 2Mi ∇2 i and (2.2) ˆ Te= − Ne ∑ k ̵h2 2me∇ 2 k (2.3)

are the kinetic energy of the nuclei and the electrons. Mi is the mass of the

nucleus indexed by i and me is the electron rest mass. The other terms in the

Hamiltonian are ˆ Vnn= 1 2 Nn ∑ i≠j ZiZje2 ∣Ri− Rj∣ , (2.4)

ˆ Vee= 1 2 Ne ∑ k≠l e2 ∣rk− rl∣ , and (2.5) ˆ Vne= − Nn ∑ i Ne ∑ k Zie2 ∣Ri− rk∣ (2.6) which corresponds to the potential energies due to nucleus-nucleus,

electron-electron, and nucleus-electron Coulomb interactions. Ziand Riare the atomic

number and the position, spatial coordinate, of nucleus i. rk is the position

of k-th electron and e is the elementary charge.

The time evolution of the system is given by the time-dependent Schrödinger equation [2],

ˆ

HΨ(R1, . . . , RNn; x1, . . . , xNe; t) = i̵h

∂Ψ(R1, . . . , RNn; x1, . . . , xNe; t)

∂t . (2.7)

Ψ is the many-body wave-function for Nnnuclei and Neelectrons. It depends

on the coordinates of all the particles. The coordinates for the nuclei are

spatial, Ri. For the electrons, both the spin, σk, and spatial coordinates, rk,

are assigned by xk. The problem encountered when trying to solve Eq. (2.7)

is the many degrees of freedom except for the simplest of systems [3]. The amount of data is too great to practically store. This is the reason why we cannot perform direct calculations where the solid is described by its many-body wave-function.

2.2 The Born-Oppenheimer approximation

The complexity of the many-body problem can be reduced by separating it into two parts, one for the nuclei and one for the electrons. This is done by using the Born-Oppenheimer approximation [4]. Its origin is the difference in masses between nuclei and electrons, from which the electrons are assumed to move much faster than the nuclei. Electrons will therefore reach equilibrium in a fraction of the time it would take for nuclei. From the perspective of the nuclei, the electrons will then always be at equilibrium and for the electrons the nuclei can be viewed as static. The motion of the nuclei will be dependent of the electronic system. By adopting to the perspective of the electrons and again using the fact that the masses of the nuclei are much larger, we can assume the kinetic energy of the nuclei to be small, to the point that we can ignore it [5]. The total Hamiltonian in Eq. 2.1 becomes for this perspective the following expression

ˆ

HBO= ˆTe+ ˆVee+ ˆVne+ Vnn, (2.8)

where Vnn is the classical potential interactions between nuclei. The positions

2.3. Density functional theory

means that nuclei-electron interactions, ˆVne, will be treated as an external

potential, specific for the static positions of the nuclei. The wave-function

for the system is reduced to be only for the electrons, Ψe, with a parametric

dependence of the the positions of the nuclei. This means if we could solve the time independent Schrödinger equation for the electronic structure,

ˆ

HeΨe(x1, . . . , xNe;{R1, . . . , RNn}) =

EeΨe(x1, . . . , xNe;{R1, . . . , RNn}) , (2.9)

with the electronic Hamiltonian ˆ

He= ˆTe+ ˆVee+ ˆVne, (2.10)

we would obtain the correct internal energy for the static system by only adding the classic Coulomb interactions of the nuclei to the electronic energy,

Ustat.= Ee+ Vnn. (2.11)

We would of course also obtain all the information of the electronic system if we could solve Eq. (2.9). The problem is that the degree of freedom still is to large for the electronic wave-function [3]. This is also true for the nuclei wave-function.

2.3 Density functional theory

We want to solve Eq. (2.9), which is the time independent Schrödinger equa-tion for the electronic structure. This is computaequa-tionally impossible except for systems consisting of a very small number of particles. However, reformu-lating it as an optimisation problem of the electronic density, this problem can be solved, in theory.

The electronic system is now solely considered. Therefore, I will adapt to a different notation for the terms of the Hamiltonian,

ˆ

He= ˆTe+ ˆVee+ ˆVne= ˆT+ ˆVint.+ ˆVext.. (2.12)

The change in notation reflects that we are using the “electronic perspective” of the Born-Oppenheimer approximation, presented in the previous section. The interactions between the electrons and the nuclei are viewed as the elec-trons are interacting with an external potential. The interactions among the electrons are considered to be internal for the system. The subscript for the electronic kinetic energy operator can be removed because the kinetic energy of the nuclei is discarded.

Hohenberg and Kohn introduced density functional theory (DFT) in 1964 by suggesting that the electronic system could be described by its density, replacing the electronic wave-function [6], significantly reducing the number of degrees of freedom. DFT can be summed up in two theorems:

1. The external potential, Vext.(r), of a system consisting of charged inter-acting particles is uniquely defined, up to an additive constant, by the

ground state density of the charged particles, n0(r).

2. For the density of the charged particles, n(r), a universal energy

func-tional exist for any external potential, F[n(r)]. With it, together with

the energy contribution of the external potential the following functional can be defined:

EVext.[n(r)] = F [n(r)] + ∫ Vext.(r)n(r)dr (2.13)

The ground state energy and density is found by minimising this func-tional,

E0= min

˜

n EVext.[˜n(r)] = EVext.[n0(r)]. (2.14)

The problem now is to first find the functional F[n(r)] and then the ground

state density n0(r). The issue is that the theory of Hohenberg and Kohn does

not provide any tools to accomplish this.

2.4 The Kohn-Sham auxiliary system

Almost exactly a year after Hohenberg and Kohn published their new formu-lation of the electron structure problem, Kohn and Sham published a practical methodology to find the ground state density [7]. The core of their approach is to replace the electrons with non-interacting particles. An effective poten-tial which will yield the same ground state density as the original problem is replacing the external potential. We can therefore describe the system with a

set of single-particle orbitals, ψk(r) fulfilling the Schrödinger-like equation:

ˆ

HKSψk(r) = ϵkψk(r). (2.15)

In (2.15) the Hamiltonian is defined as ˆ HKS= − ̵h2 2me ∇2_{+ V} eff.(r) (2.16) where

Veff.(r) = Vext.(r) + VHar.(r) + Vxc(r). (2.17)

Here, the Hartree potential, VHar.(r), is the classical potential energy per

electron for the electron-electron interactions and it is the classical counterpart

to ˆVint. in Eq. (2.12). Vxc is the modification of the potential to compensate

for changing from interacting particles to non-interacting. The density is obtained from the single-particle orbitals as

n(r) = Ne ∑ k ∣ψk(r)∣ 2 . (2.18)

2.4. The Kohn-Sham auxiliary system We can rewrite the ground state energy functional, Eq. (2.13), as

EKS= Ts[n(r)] + ∫ n(r)Vext.(r)dr + EHar.[n(r)] + Exc[n(r)], (2.19)

were Ts[n(r)] is the non-interacting kinetic energy of the particles. The

Hartree energy, which is the classical repulsion energy of a electron gas with

the density n(r) is written as

EHar.[n(r)] = ∫

e2n(r)n(r′)

∣r − r′_∣ drdr′. (2.20)

From this expression, we can define the Hartree potential by taking the func-tional derivative of it,

VHar.(r) =

δEHar.[n(r)]

δn(r) = ∫

e2n(r′)

∣r − r′_∣dr′. (2.21)

Eq. (2.19) is exact and we can define the exchange-correlation energy by setting it equal to Eq. (2.13):

Exc[n(r)] = F [n(r)] − (Ts[n(r)] + EHar.[n(r)]) . (2.22)

Rewriting F[n(r)] as the sum of the expectation values of the electronic

kinetic energy, ˆT , and the internal potential, ˆVint., one obtains a more

com-prehensible form:

Exc[n(r)] = ⟨ ˆT⟩ − Ts[n(r)] + ⟨ ˆVint.⟩ − EHar.[n(r)]. (2.23)

Vxc(r) is obtained from the functional derivative of Exc[n(r)],

Vxc(r) =

δExc[n(r)]

δn(r) . (2.24)

Eq. (2.23) is good for telling us more clearly what exchange-correlation is: the difference in kinetic energy between the electrons and the non-interacting particles and the difference in the repulsion energy between the system of electrons treated quantum mechanically and classically. But Eq. (2.23) is of no use to actually calculate the exchange-correlation term. An exact numerical expression for it does not exist. Therefore this term is approximated. A large number of approximations are available to us, but the simplest and still useful is the local density approximation (LDA), suggested by Kohn and Sham [7]. Within this approximation the exchange-correlation energy per particle at

point r, ϵhom.

xc (n(r)), is assumed to be as for an homogenous electron gas with

the density n(r). The exchange-correlation energy for the entire system is

obtained by integrating ϵhom.

xc (n(r)) multiplied by n(r) at r over the entire

spatial space,

To refine the approximation for the exchange-correlation term, a dependency

of both the density n(r) and its gradient ∇n(r) is assumed in the generalised

gradient approximation (GGA),

EGGAxc [n(r)] = ∫ n(r)ϵxc(n(r), ∇n(r))dr. (2.26)

It is possible to formulate a GGA functional in numerous ways, which has lead to many versions of this approximation. But the most commonly used for solid state calculations is the version suggested by Perdew, Burke, and Ernzerhof (PBE) [8].

The ground state density is calculated in Kohn-Sham DFT by an iterative

process. From an initial guess of the particle density, n0(r), the effective

potential Veff.(r) is constructed. For this potential the Kohn-Sham equation

can be solved, Eq. (2.16), and a new density is generated from the single-particle orbitals, Eq. (2.18). If the difference in densities between the latest and the second latest iteration are smaller than a pre-decided convergence value, self-consistency is reached. The iteration process is stopped and the ground state density is found. A graphical representation of this process is presented in Fig. 2.1.

2.5 Practicalities of performing Kohn-Sham DFT for

solids

For all the solids treated in this thesis, there is a periodic crystal structure. Since the nuclei are positioned in a periodic pattern, the effective potential

Veff.(r) will have the same periodicity as the crystal structure. Therefore,

ac-cording to Bloch’s theorem, can the Kohn-Sham orbitals, ψk(r), be expressed

as a plane-wave multiplied by a periodic function with the same periodicity as the crystal structure,

ψkn(r) = ukn(r)eik⋅r. (2.27)

The orbitals are now indexed by their wave-vector k and a band index n.

ukn(r) is the periodic function fulfilling the expression

ukn(r) = ukn(r + R), (2.28)

where R is a lattice vector in any direction. All the Bloch states will be inside the primitive unit cell of the reciprocal lattice, which is defined as the first Brillouin zone (BZ). The number of allowed k-states is equal to the number of primitive cells in the crystal. The allowed states are decided by periodic boundary conditions. For a piece of a solid which we are used to handle in our macroscopic world, there will be a continuum of k-states, which would demand infinite computer resources to be calculated with this theory. Although there is a continuum in k-states, a finite mesh can be used, and beyond a certain

2.5. Practicalities of performing Kohn-Sham DFT for solids

Generate effective potential

Initial guess of density

Solve Kohn-Sham equations

Calculate new density

Reached self-consistency?

Yes

Done!

No

Figure 2.1: Flowchart of the different steps in the iterative loop of Kohn-Sham DFT.

density in k-points, the results should converge. A parameter the user must handle therefore is the k-point density. This is the number of wave-vectors included in the calculation, representing the periodicity of the simulation cell. The user should ensure that the calculations are converged with respect to the number of k-points used to probe the first BZ.

Numerous methods have been developed to solve the Kohn-Sham equa-tions. The approach I have used throughout my work is to expand the

Kohn-Sham orbitals in a basis set of plane-waves. The periodic function ukn(r) of

the Kohn-Sham orbitals is then expanded as

ukn(r) = ∑

G

Ckn(G)eiG⋅r. (2.29)

The orbitals can therefore be expressed as

ψkn(r) = ∑

G

Ckn(G)e(k+G)⋅r. (2.30)

In Eqs. (2.29) and (2.30) Ckn(G) is the plane-wave expansion coefficient and

G is a reciprocal lattice vector. The number of plane-waves used to expand the periodic function for the Bloch functions is also a parameter that the user need to define. It is usually set as a limit of the plane-wave kinetic energy

̵h2

2me∣k + G∣

2

< Ecut, (2.31)

where the only plane-waves included in the calculation are those whose wave-vector corresponds with the G-wave-vectors fulfilling Eq. (2.31). Again, the user

needs to ensure that the calculations are converged with respect to Ecut.

In my work, I have used exclusively the projector augmented-wave (PAW) method [9] as implemented in the Vienna Ab initio Simulation Package (VASP) [10–12]. This method works under the premise that a linear

op-erator can be constructed that relates the pseudo wave-function,∣ ˜ψ⟩, with the

all-electron wave-function, ∣ψ⟩,

∣ψ⟩ = ˆT ∣ ˜ψ⟩ . (2.32)

The pseudo wave-functions used in VASP consist of plane-waves, but with a modified potential. The modified potential is designed to produce smooth wave-functions, even close to the nuclei where the all-electron wave-functions for the valence electrons are rapidly oscillating. The plane-wave expansion for the pseudo wave-functions will therefore need a smaller number of plane-waves compared with the all-electron wave-function.

The operator in (2.32) can be written as ˆ

T = 1 + ∑

i

2.5. Practicalities of performing Kohn-Sham DFT for solids

where ∣ϕi⟩ − ∣ ˜ϕi⟩ is the difference between partial all-electron wave-functions

and pseudo wave-functions and∣pi⟩ are projection functions. These functions

are precomputed. The difference between the all-electron wave-functions and the pseudo wave-functions are by design zero outside an augmented sphere centred around the the nuclei.

Chapter

3

Lattice dynamics

The atoms in a solid are under constant motion. Even at absolute zero, they are still moving. Lattice dynamics is a theory that allows us to study the vibrational part of these motions. In my work, lattice dynamics has been my weapon of choice. It has been used to study mechanical and thermal properties, dynamical stability of crystal structures and how the lattice vi-brations can be related to the electronic structure. I will initially introduce it as a classical theory and thereafter a quantum mechanical, where the lattice vibrations are described as quasi-particles we know as phonons. I will also in this chapter include how the vibrational Helmholtz energy is calculated, the quasi-harmonic approximation, and how the inclusion of three-phonon interactions affect the phonon calculations.

3.1 Force constants and the dynamical matrix

Let me start with assuming that we have an infinite periodic structure. This structure consist of an infinite repetition of a unit cell. The position of each unit cell is described by

rl= l1a1+ l2a2+ l3a3. (3.1)

In (3.1) l1, l2, and l3are integers whose combination is specific and individual

to the index l and a1, a2, and a3are primitive translation vectors. Within this

unit cell different atoms will be positioned in a basis described by the vectors

rµ, we assume this basis to be the equilibrium positions. The positions for

the atoms when their displacement, ulµ, are included are then expressed by

A natural approach is to assume that the atoms do not venture too far from their equilibrium positions and that the energy of the system can be described by a Taylor expansion in the displacements [13]:

U = U0+ ∑ l Na ∑ µ ∑ α Φα_lµuα_lµ +1 2∑_lm Na ∑ µν ∑ αβ Φαβ_lmµνuα_lµuβ_mν+ 1 3!_lmn∑ Na ∑ µνξ ∑ αβγ Φαβγ_lmnµνξuα_lµuβ_mνuγ_nξ+ ⋯. (3.3)

In Eq. (3.3) lmn, µνξ, and αβγ are indices to unit cells, atoms within the

basis, and cartesian components. Nais the number of atoms in the basis. The

Taylor expansion coefficients

Φαlµ= ∂U ∂uα lµ RRRRR RRRRR R{u}=0 , (3.4a) Φαβ_lmµν= ∂ 2_U ∂uα lµ∂u β mν RRRRR RRRRR RR{u}=0 , and (3.4b) Φαβγ_lmnµνξ= ∂ 3_U ∂uα lµ∂u β mν∂uγ_nξ RRRRR RRRRR RR{u}=0 (3.4c)

are the atomic force constants. The subscript {u} = 0 emphasises that the

evaluations of the partial derivatives are executed at the atoms’ equilibrium positions.

The first order force constants, Eq. (3.4a), give the α-component of the negative of the force on the atom µ in the unit cell l when it and all the other atoms are in their equilibrium positions. By definition, these should be zero due to, as recently stated, the evaluation of the partial derivatives are assumed to be at the atoms equilibrium positions. The second order force constants, expressed in Eq. (3.4b), give the force acting on atom µ in the unit cell l when atom ν in unit cell m is moved in the β direction. All the remaining higher order terms can, in a similar way, be physically interpreted. If we truncate the Taylor expansion of the potential energy, Eq. (3.3), to the second order and are taking the assumed equilibrium positions to be true, the expression for the potential energy then becomes

U = U0+ 1 2∑_lm Na ∑ µν ∑ αβ Φαβ_lmµνuαlµu β mν. (3.5)

This is the harmonic approximation. Its Hamiltonian for the vibrational sys-tem is Hvib.class.= U0+ ∑ l Na ∑ µ ∑ α Pα lµ 2Mµ +1 2∑_lm Na ∑ µν ∑ αβ Φαβ_lmµνuα_lµuβ_mν (3.6)

3.2. Quantum mechanical description for which we can write the equation of motion (EOM) as

Mµu¨αlµ= − ∑ m Na ∑ ν ∑ β Φαβ_lmµνuβ_jν. (3.7)

The solution for this EOM can be written as a set of plane-waves [13, 14],

uα_lµ∣qs= √1 Mµ

ϵα_µ∣qsei(q⋅rlµ−ωqst)_{. (3.8)}

In Eq. (3.8) ϵα

µ∣qs is the cartesian α-component for the polarisation vector for

all the µ atom in the unit cell related to the frequency ωqs. The polarisation

vector describes the direction of the displacements. q is a wave-vector in the first Brillouin zone (BZ) like k in Eq. (2.27). s is the polarisation, which

there will be a total of three times the number of atoms in the unit cell, 3Na.

A plane wave described by a q-vector and the index s is called a normal mode of the crystal structure. All the normal modes are independent of each other, and for each atom µ in the unit cell there will be three individual modes, two transverse and one longitudinal.

ωqs and ϵqs are obtained as the eigenvalues and eigenvectors of the

dy-namical matrix

ω_qs2 ϵqs= D(q)ϵqs. (3.9)

The dynamical matrix is of the size 3Na× 3Na and is constructed as

D(q) =⎛⎜ ⎝ D11(q) ⋯ D1Na(q) ⋮ ⋱ ⋮ DNa1(q) ⋯ DNaNa(q) ⎞ ⎟ ⎠, (3.10) where Dµν(q) = ∑ l Φ0lµν √ MµMν eiq⋅(rlν−r0µ)_{. (3.11)}

Note that in Eq. (3.11) Φ0lµν is a 3× 3 matrix. The index 0 denotes an

arbitrary reference unit cell in the crystal. The displacements from the equi-librium positions should result in restoring forces, acting to counteract the displacements. This corresponds to real positive eigenvalues of the dynamical matrix. The crystal structure is dynamically stable. Dynamical instability is when the forces are not restoring: the forces acting on the atoms are moving them away from the assumed equilibrium positions. The real matrices for the second order force constants will not be positive definite, resulting in nega-tive eigenvalues of the dynamical matrix. The frequencies from a dynamically unstable solid will therefore be imaginary.

3.2 Quantum mechanical description

We have so far treated the vibrations of atoms in solids classically, by solving Newton’s equations for the atoms, which we have viewed as classical particles.

We noted that the solutions are sets of uncoupled normal modes, Eq. (3.8). A quantum mechanical representation of the normal modes is possible by

using creation, ˆaqs, and annihilation, ˆa

†

qs, operators of quantum harmonic

oscillators, ˆ uαlµ= ¿ Á Á À ̵h 2NucMµ∑q ∑ s ϵα_µ∣qs √_ω qs eiq⋅rlµ(ˆa qs+ ˆa†qs) . (3.12)

Nuc is the number of unit cells in the crystal. The vibrational Hamiltonian

can be expressed by ˆaqsand ˆa

† qsas ˆ H_vib.qm = ∑ q 3Na ∑ s ̵hωqs(ˆa†qsˆaqs+ 1 2) . (3.13)

The vibrations are now treated to consist of non-interacting quasiparticles

known as phonons. ˆa†qsˆaqs is the number operator giving the number of

phonons occupying the state qs. Phonons are bosons and will therefore occu-pying states according to Bose-Einstein statistics,

nqs=

1

eβT̵hωqs− 1. (3.14)

Eq. (3.14) is the Bose-Einstein occupation number and βT is the inverse of

the temperature, βT = 1/kBT .

3.3 Vibrational free energy

Within the harmonic approximation, the system consists of uncoupled har-monic oscillators. The partition function of such a system is expressed by

Z= ∏

qs

Zqs=

e−βT ̵hωqs2

1− e−βT̵hωqs. (3.15)

From which Helmholtz free energy can be expressed as

Fvib.= −kBT ln Z= ∑ q 3Na ∑ s (̵hωqs 2 + kBT ln(1 − e −β̵hωqs)) = ∫₀∞g(ω) (̵hω 2 + kBT ln(1 − e −β̵hω_{)) dω,} (3.16)

where the second equality is the expression rewritten as dependent of the

phonon density of states, g(ω). Which can be expressed as

g(ω) = ∑

s

Vuc

(2π)3∫_BZδ(ω − ωqs)dq, (3.17)

3.4. Quasi-harmonic approximation (QHA)

3.4 Quasi-harmonic approximation (QHA)

A first step to include the effects of temperature in our simulations is to include it indirectly, as its effect on a parameter. This parameter is usually taken to be the volume of the simulation cell. Phonon frequencies are dependent on the temperature, but if we stay in the harmonic approximation we cannot explain why they should change nor why the simulated material should expand with increased temperature. But by assuming that the harmonic approximation is valid for each value of the volume [15], we could trace the change in the harmonic potential. This approximation is described by the expersion

ω(V, T ) ≈ ω(V (T )). (3.18)

Within this approximation, we can define the Helmholtz free energy as

FQHA= Ustat.(V ) + ∑

qs

(̵hωqs(V )

2 + kBT ln(1 − e

−β̵hωqs(V ))) , _(3.19)

where Ustat.(V ) is the internal energy of the static lattice, as expressed in

Eq. (2.11), at the volume V . By finding the volume minimising FQHA for a

specific temperature, the thermal expansion can be calculated, and changes in the corresponding phonon frequencies can be traced.

3.5 The small displacement method

Methods for calculating the dynamical matrix can be divided into two types. One is the linear response method [16–18] where the dynamical matrix is constructed from the inverse dielectric matrix. The other one is the direct

method, which is divided into two subtypes. These are the frozen-phonon method [19–21] and the small displacement method [21, 22]. In the

frozen-phonon metod the frozen-phonon energy is calculated as the difference in energy between the undistorted and a distorted supercell. The phonon energy is calculated as a function of the amplitude of the displacements, which the force constants are approximated from. The small displacement method, which I have used in my work, uses the Hellmann-Feynman forces [23],

Flµ= − ⟨Ψe∣ ∇Rlµ( ˆHe+ Vnn) ∣Ψe⟩ = −∇RlµUstat({R}), (3.20)

to approximate the force constants [24]:

Φαβ_lmµν= ∂F α lµ ∂uβmν ∣ {u}=0 ≈ F α lµ ∆uβmν ∣ {∣u∣}=∆u . (3.21)

In Eq. (3.21),{∣u∣} = ∆u emphasises that the derivatives are evaluated at the

displacements. The core of this approximation is written in its name, that the displacement should be very small, they are usually around 0.01 Å.

3.6 Interacting phonons

Hitherto, we have considered phonons to be non-interacting quasiparticles with infinite lifetime. In reality, they are interacting and have a finite lifetime. From the third order force constants, the effect of including three-phonon interactions can be evaluated. This is done by calculating the self-energy from the third order force constants. The self-energy is frequency dependent and consists of a real and imaginary part,

Σqs(Ω) = ∆qs(Ω) + iΓqs(Ω). (3.22)

The imaginary part of the self-energy for three-phonon interactions can be expressed as [25] Γqs(Ω) = ̵hπ 16 _q∑′q′′ 3Na ∑ s′s′′ ∣Φqq′q′′ss′s′′∣2 ×{(nq′s′+ nq′′s′′+ 1)δ(Ω − ωq′s′− ωq′′s′′) +(nq′s′− nq′′s′′)[δ(Ω − ωq′s′+ ωq′′s′′) −δ(Ω + ωq′s′− ωq′′s′′)]}, (3.23)

where Φqq′q′′ss′s′′is the three-phonon matrix element, which can be calculated

as [14] Φqq′q′′ss′s′′= ∑ lm Na ∑ µνξ ∑ αβγ ϵα µ_∣qsϵ β ν∣q′s′ϵ γ ξ∣q′′s′′ √_ω qsωq′s′ωq′′s′′ Φαβγ_0lm √ MµMνMξ ei(q⋅r0µ+q′⋅rlν+q′′⋅rmξ)_. (3.24)

In Eq. (3.23) nqs is the Bose-Einstein occupation number, Eq. (3.14). The

real part of the self-energy can be obtained from a Kramers-Kronig transfor-mation of the imaginary part,

∆qs(Ω) = 1 π ∫ ∞ −∞ Γqs(ω) ω− Ω dω. (3.25)

To compare our calculated results with the ones obtained from inelastic neutron scattering, we can calculate the one-phonon scattering cross section [25], σ(q, Ω) ∝ ∑ s 2ωqsΓqs(Ω) [Ω − ωqs2− 2ωqs∆qs(Ω)] 2 + 4ω2 qsΓ2qs(Ω) . (3.26)

An example of a one-phonon scattering cross section, in this work usually referred to as a phonon spectral function but not uncommonly also called

3.6. Interacting phonons lineshape, was presented in the introduction, see Fig. 1.1 (b). We can inter-pret Eq. (3.26) and Fig. 1.1 (b) as the probability to excite a phonon with the wave-vector q and the energy ̵hΩ. For a specific q-vector there will be a probability distribution for the excitation energy. If we assume the real and imaginary part to be small compared with the square root of the eigenvalues

of the dynamical matrix, ωqs ≫ ∆qs, Γqs, the spectral function will be

re-duced to a sum of Lorentzian functions. The real part of the self-energy shifts

the peak of the Lorentzian functions from the frequencies ωqs, centring them

at ωqs+ ∆qs. The imaginary part, Γqs, is the half width at half maximum

of the probability distribution. Giving the inverse of the phonon lifetime as

Chapter

4

Simulating solids as

thermally excited supercells

With DFT, we can calculate the internal energy of the static system, Ustat.,

for a specific spatial configuration of the nuclei, which I hereon will refer to as a “configuration”. However, to simulate a solid at a certain temperature, there is not only one correct configuration, there are many of them. A simulation of a solid consists of a set of configurations. These configurations are thermally excited supercells. A supercell is generated by repeating a unit cell along its primitive translation vectors. We use supercells to include the interactions among the nuclei, which are temperature dependent. The size of the supercell also dictates from the boundary conditions which normal modes propagates throughout it.

If the simulation is at a specific volume and temperature, these configu-rations would be representatives of a canonical ensemble. The task is now, to make a set of configurations that mimics the distribution of representa-tives for the canonical ensemble. I will introduce two classical approaches, which can solve this problem. In these two different methods, the nuclei are treated as classical particles. By solving Newton’s equation of motion for the nuclei, their trajectories can be obtained and are the configurations added to the set. I will also introduce a method in this chapter that will conserve the Bose-Einstein statistics of the phonons for the configurations.

4.1 Ab initio molecular dynamics (AIMD)

If we view the nuclei to be classical particles, their motion can be described by Newtonian mechanics, ˙ Ri= Pi Mi = Vi (4.1a) ¨ Ri= ˙ Pi Mi = Fi Mi . (4.1b)

In Eq. (4.1), ˙R and Vi are the velocity, Pi is the momentum, and Mi is the

mass for the nucleus positioned at Ri. ˙P and Fi are the force acting on this

nucleus, obtained by using the Hellmann–Feynman theorem [23],

Fi= − ⟨Ψe∣ ∇Ri( ˆHe+ Vnn) ∣Ψe⟩ = −∇RiUstat({R}). (4.2)

The forces obtained account for the quantum mechanical behaviour of the electrons, but the nuclei are treated as classical particles and will therefore still move like classical particles. The trajectories are obtained by a numerical integration of the Newtonian equation of motion. The integration itself can be performed with different algorithms [26]. One of the simplest but also usually best working is the velocity Verlet algorithm [27],

Ri(t + δt) = Ri(t) + ˙Riδt+ 1 2 ¨ Riδt2+ O(t4) (4.3a) ˙ Ri(t + δt) = ˙Ri(t) + ¨ Ri(t) + ¨Ri(t + δt) 2 δt+ O(t 4_{). (4.3b)}

A physical observable is the ensemble average. By letting the system evolve and tracing the trajectories, we are sampling representatives of the ensemble. But collecting them all for the entire ensemble, to be used as a set of configu-rations, is not possible. We do not have the computational resources to do so, and never will have. But we can replace the ensemble average by the average of the set of configurations, if the system has no memory of its initial state, according to the ergodic principle:

⟨A⟩ = lim Nc→∞ 1 Nc Nc ∑ c Ac. (4.4)

The type of ensemble considered to this point, have a constant number of particles, volume, and energy. It is a microcanonical ensemble (NVE). In my work, I have studied the effects of temperature. Therefore, the representation of the system needs to be for a specific temperature. An ensemble with con-stant number of particles, volumes, and temperature is a canonical ensemble (NVT). We can change to this type of ensemble by introducing a thermostat

4.2. Self-consistent phonon (SCP) calculations in our simulation. If the Nosé-Hoover thermostat is used [28, 29], Eq. (4.1b) is rewritten as ¨ Ri= Fi Mi − γPi Mi = Fi Mi − γ ˙Ri= Fi Mi − γVi, (4.5)

were γ is a friction term that evolves in time as

˙γ= 1

Q(∑_i Mi

˙

R2i − χkBTset) . (4.6)

In Eq. (4.6) Q is the Nosé-mass, which couples the system with an external

heat-bath and χ is the number of degrees of freedom. Tsetis the temperature

of the external heat-bath, which the system is converging towards.

4.2 Self-consistent phonon (SCP) calculations

A major drawback of using AIMD to generate a set of configurations is the computational cost. To make matter even worse, tracing the trajectories is an ineffective way to sample configurations because only a small volume of the phase space is probed for a large number of sequential configurations. An additional disadvantage is that the nuclei in the generated configurations will be distributed according to Maxwell–Boltzmann statistics, due to the classi-cal treatment of them, excluding quantum mechaniclassi-cal effects like zero-point motion. A method allowing a more effective sampling of the configurations is self-consistent phonon (SCP) calculations [30, 31]. An additional advan-tage for this method is that the nuclei in the generated configurations could be distributed according to Bose–Einstein statistics, which includes zero-point motion. From an initial guess of the phonons, configurations can be generated by sampling the phase space,

uα_i = ∑ q ∑s ⟨Ai_∣qs⟩ ϵαi∣qs √ −2 ln ξ1sin 2πξ2 (4.7a) ˙ uα_i = ∑ q ∑ s ⟨Ai∣qs⟩ ωi∣qsϵαi_∣qs √ −2 ln ξ1cos 2πξ2. (4.7b)

In Eqs. (4.7a) and (4.7b) uα

i and ˙u

α

i are the displacement from the

equilib-rium position and velocity of the i-th nucleus for the cartesian component

α. The wave-vectors, q, included in the sums are the ones allowed by the

periodic boundary conditions for the simulation cell. ωi∣qs and ϵi∣qs are the

phonon frequencies and polarisation vectors related to nucleus i. ξ1 and ξ2

are two-dimensional normal distributed scalars between 0 and 1. The normal distribution in two dimensions is ensured by a Box-Muller transformation [32]. The thermal average of the amplitude for the normal mode qs is given by [33]

⟨Ai∣qs⟩ = ¿ Á Á À ̵h Miωqs (nqs+ 1 2) ≈ 1 ωqs √ kBT Mi . (4.8)