First-Principles Study of Elastic Properties of Fe-Mg alloy at Earth’s core pressure

Department of Physics, Chemistry and Biology

Bachelor’s Thesis

First-Principles Study of Elastic Properties of

Fe-Mg alloy at Earth’s core pressure

Ulf Karg´

en

LiTH-IFM-G-EX-08/2011-SE

Department of Physics, Chemistry and Biology Link¨opings universitet, SE-581 83 Link¨oping, Sweden

Bachelor’s Thesis LiTH-IFM-G-EX-08/2011-SE

First-Principles Study of Elastic Properties of

Fe-Mg alloy at Earth’s core pressure

Ulf Karg´

en

Advisers: Igor Abrikosov IFM, Link¨oping Christian Asker

IFM, Link¨oping Examiner: Igor Abrikosov IFM, Link¨oping

Avdelning, Institution Division, Department

Theoretical physics

Department of Physics, Chemistry and Biology Link¨opings universitet, SE-581 83 Link¨oping, Sweden

Datum Date 2008-10-16 Spr˚ak Language  Svenska/Swedish  Engelska/English   Rapporttyp Report category  Licentiatavhandling  Examensarbete  C-uppsats  D-uppsats  ¨Ovrig rapport  ISBN ISRN

Serietitel och serienummer Title of series, numbering

ISSN

URL f¨or elektronisk version

Titel

Title First-Principles Study of Elastic Properties of Fe-Mg alloy at Earth’s core pressure

F¨orfattare Author

Ulf Karg´en

Sammanfattning Abstract

The purpose of this thesis has been to investigate the elastic properties of an fcc FeMg alloy with 10 at.% magnesium under high pressure. Recent research has shown that magnesium can be a possible candidate for light element impurities in the Earth’s inner core, something that was previously not considered possible because of the low miscibility of magnesium in iron at ambient pressure. Gaining knowledge about the composition of the Earth’s core can help us better understand such phenomena as seismic activity and the fluctuations of the Earth’s magnetic field.

The elastic constants of the FeMg alloy was calculated using ab-initio methods based on Density Functional Theory. The Exact Muffin-Tin Orbitals method was used in conjunction with the Coherent Potential Approximation.

The FeMg alloy was found to be overall considerably softer than pure iron, and the softening effect on the elastic constants was also found to increase with pressure. The results also showed that 10% Mg alloying increased the anisotropy with about 40% compared to pure iron.

Nyckelord Keywords

Density Functional Theory, EMTO, Earth’s core, High-Pressure Alloying, Fe-Mg alloy, Elastic Constants



http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-15742 —

LiTH-IFM-G-EX-08/2011-SE

Abstract

The purpose of this thesis has been to investigate the elastic properties of an fcc FeMg alloy with 10 at.% magnesium under high pressure. Recent research has shown that magnesium can be a possible candidate for light element impurities in the Earth’s inner core, something that was previously not considered possible because of the low miscibility of magnesium in iron at ambient pressure. Gaining knowledge about the composition of the Earth’s core can help us better understand such phenomena as seismic activity and the fluctuations of the Earth’s magnetic field.

The elastic constants of the FeMg alloy was calculated using ab-initio methods based on Density Functional Theory. The Exact Muffin-Tin Orbitals method was used in conjunction with the Coherent Potential Approximation.

The FeMg alloy was found to be overall considerably softer than pure iron, and the softening effect on the elastic constants was also found to increase with pressure. The results also showed that 10% Mg alloying increased the anisotropy with about 40% compared to pure iron.

Acknowledgements

First and foremost, I would like to thank my supervisor and examiner Igor Abrikosov for allowing me to do this diploma project, and for all his encouraging help. I would also like to thank my supervisor Christian Asker for all his help, especially with understanding the EMTO software. For providing useful comments and suggestions I thank my opponent Marcus Ekholm. I also thank the National Supercomputer Centre in Link¨oping for supplying the supercomputer resources for this work.

Contents

1 Introduction 1

1.1 Background . . . 1

1.2 Thesis outline . . . 2

2 Theory 5 2.1 Density functional theory (DFT) . . . 5

2.2 The Kohn-Sham equations . . . 6

2.3 The Exchange-Correlation energy functional . . . 7

2.3.1 The Local Density Approximation (LDA) . . . 7

2.3.2 The Generalized Gradient Approximation (GGA) . . . 8

3 Computational methods 9 3.1 Solving the Kohn-Sham equations numerically . . . 9

3.1.1 Finding solutions self-consistently . . . 10

3.2 The Exact Muffin-Tin Orbitals Method . . . 10

3.3 The Coherent Potential Approximation . . . 12

4 Elastic constants 13 4.1 Basic theory . . . 13

4.1.1 Elastic constants of cubic systems . . . 15

4.2 Calculating elastic constants . . . 16

4.3 The Face Centered Cubic lattice . . . 17

4.3.1 Polycrystalline elastic constants . . . 18

5 Computational details 19 5.1 The EMTO software . . . 19

5.2 Equation of State . . . 19

5.3 Elastic constants of the fcc lattice . . . 20

5.4 K-point convergence . . . 20

6 Results and discussion 23 6.1 Equation of State . . . 23

6.2 Elastic constants . . . 24

6.3 Polycrystalline shear moduli and anisotropy . . . 26 ix

7 Conclusions and outlook 29 7.1 Conclusions . . . 29 7.2 Outlook . . . 29

Chapter 1

Introduction

1.1

Background

Although we can observe stellar phenomena millions of light years away and send spacecraft to other planets we still know surprisingly little about the inner work-ings of our own planet. The main reason for this is the difficulty of doing direct observations of the interior of the Earth. No one has yet been able to drill beyond the thin layer of rock, known as the crust, on which we live. By analyzing the propagation of sound waves through the Earth, some insights in the structure of our planet have been gained. Below the crust is a solid region known as the man-tle extending approximately 2800 km down, which is about halfway to the center. Below this lies the core. The core is further divided into the liquid outer core and the solid inner core, which starts at a depth of about 5000 km. With state of the art experimental equipment, the conditions present about halfway through the outer core can be reproduced. The conditions in the inner core however, where the pressure is about 350 GPa (3.5 million atmospheres) and the temperature is over 6000 K, are still unattainable with experimental methods. Detailed knowl-edge of the chemical composition of the core, as well as crystal structures and other material properties are necessary in order to better understand and predict phenomena that could greatly affect us, such as Earth quakes and the fluctuation (and occasional flipping) of the Earth’s magnetic field.

With the rapid development of computer technology during the past few decades, advanced computer models have come to play an increasingly important role in understanding the interior of our planet. Early models used adjustable parameters that were fitted to experimental data. While these models provided many impor-tant insights, the entire approach has obvious shortcomings. Determining under which conditions an empirical model is valid can be very difficult. Experimentally obtained parameters may simply not be valid when the simulated conditions differ too much from the conditions at which the experimental data was obtained. An entirely different approach to the problem is so called ab-initio or first-principles methods. These methods are based on the fundamental and exact equations of quantum mechanics. In the ideal case, first-principles methods should be able

to exactly reproduce the properties of a wide range of materials under almost all kinds of conditions. In practice, of course, approximations must be made to reduce the computational effort, which in turn introduces constraints on the accuracy of the models.

From direct sampling of volcanic material and advanced experimental tech-niques combined with ab-initio calculations, the properties of the different regions of the mantle are relatively well known. The composition and material properties of the core on the other hand, is still rather poorly known. From observations of sound propagation through the Earth, cosmochemical data and studies of iron meteorites there is strong evidence that the core consist mainly of iron, with 5–15 % nickel content. The observed density of the core is however lower than it should be for the said iron-nickel alloy, which suggests the existence of small amounts of lighter elements in the core. The most commonly proposed candidates for lighter elements are sulfur, silicon and oxygen. Magnesium has previously not been con-sidered a candidate since it will not alloy with iron at ambient pressure. In 2005 however, Dubrovinskaia et al. [6] showed that under the pressure ranges present in the Earth’s core, iron and magnesium can form an alloy with about 10% mag-nesium.

Apart from the chemical composition, the crystal structure of the core is an-other very important property to investigate. It has long been assumed that iron should be in the hcp phase at the conditions present in the inner core. Recent research has however shown that both the bcc and fcc phases could be possible. Dubrovinsky et al. [7] has found that an iron-nickel alloy with 10% nickel under-goes a phase transformation to the bcc structure at about 225 GPa and 3400 K. These are about the same conditions as at the outer regions of the liquid core. Ab-initio studies in combination with experiment by Mikhaylushkin et al. [10] have also shown that at very high pressures and temperatures, similar to the ones in the inner core, the fcc phase of pure iron becomes as stable as the hcp phase. This indicates that the Earth’s inner core could have a complex “mixture” of crystal structures.

In light of these recent findings, there is clearly a need to further investigate the effects of small amounts of impurities on the material properties of iron. Elastic constants are some of the few properties of the Earth’s interior that can be directly measured, by analyzing sound waves traveling through the Earth. Since first-principles techniques have also been shown to work well for calculating elastic constants, the elastic properties are a good starting point for investigation.

The focus of this diploma project has been to use first-principles simulations to study the differences in elastic properties between a fcc Fe90Mg10 alloy and pure

iron, under extreme pressure.

1.2

Thesis outline

This diploma project has been carried out at the Theoretical Physics group within the Department of Physics and Measurement Technology (IFM), at the University of Link¨oping. The thesis consists of seven chapters. A brief description of all the

1.2 Thesis outline 3 chapters, except for this one, is given in the outline below.

Chapter 2: Theory

Here, the fundamental theory behind the computational methods used in this thesis is described. The Density Functional Theory and the Kohn-Sham equations are briefly described. The two most commonly used approxima-tions of the exchange-correlation energy are also treated.

Chapter 3: Computational methods

In this chapter, the basic concepts of solving the Kohn-Sham equations nu-merically are described. The basics of the Exact Muffin-Tin Orbitals Method and the Coherent Potential Approximation are explained.

Chapter 4: Elastic constants

This chapter treats the basic theory of elasticity, and the elastic constants for cubic systems are briefly described. The methods used for calculating the elastic constants of an fcc system using ab-inito software are also explained. Chapter 5: Computational details

This chapter explains some practical details on how the calculations were performed using the EMTO software.

Chapter 6: Results and discussion

Here, the results of the calculations are given and discussed. Chapter 7: Conclusions and outlook

The conclusions of the project are summarized and some suggestions of pos-sible continuations and extensions of the work are given.

Chapter 2

Theory

2.1

Density functional theory (DFT)

The material properties of condensed matter are governed by the interactions between electrons, and between electrons and nuclei. These interactions must be described quantum-mechanically, which in practice means that we need to solve the Schr¨odinger equation of the particle system. The general form of the (time-independent) Schr¨odinger equation is:

ˆ

HΨ = EΨ (2.1)

The Hamiltonian ( ˆH) will include the following terms: ˆ

H = ˆTnuc+ ˆT + ˆVnuc+ ˆVee+ ˆVext (2.2)

= −~ 2 2 X i ∇2 i Mi − ~ 2 2me X k ∇2 k+ 1 2 X i6=j ZiZj |Ri− Rj| +X k6=l e2 |rk− rl| −X i,k Zie |ri− Rk|

Where the first term ( ˆTnuc) is the kinetic energy of the nuclei, the second ( ˆT ) is the

kinetic energy of the electrons, the third ( ˆVnuc) is the Coulomb-interaction between

the nuclei, the fourth ( ˆVee) is the Coulomb-interaction between the electrons and

the last one ( ˆVext) is the external potential, the Coulomb-interaction between the

electrons and the nuclei. The wave function Ψ will have the general form1_:

Ψ = Ψ(r1, r2, r3, . . . , rN, R1, R2, R3, . . . , RM) (2.3)

That is, the wave function is a function of the position of all the electrons (r:s) and all the nuclei (R:s). Since the Schr¨odinger equation is only analytically solv-able for one-electron systems, we need to solve the equation numerically. Since a macroscopic sample of a material contains the order of 1023_{atoms, it is impossible}

to solve the system numerically without introducing some approximations or other descriptions of the system.

1_{Spin-coordinates have been left out.}

A first approximation is the so-called Born-Oppenheimer approximation, in which the nuclei are considered stationary. Since the mass of a single proton or neutron is approximately 2000 times larger than that of an electron, this approx-imation is usually valid. In our case this means that the first term in equation (2.2) disappears, and the second term becomes a constant.

While this approximation somewhat simplifies the problem, the number of variables in the equation is still far too large to allow a numerical solution. A clever reformulation of the problem is the Density Functional Theory (DFT). In DFT the three-dimensional electron (charge) density is considered a fundamental variable for describing many-electron systems. The degrees of freedom is thus reduced from 3N to 3 for a system of N electrons. The Density Functional Theory is based on two theorems known as the Hohnberg-Kohn theorems [3]. The two theorems state that:

1. The local external potential ˆVext(r) is within an additive constant

deter-mined by the ground-state electron density n(r).

2. The value of the charge density that minimizes the total energy functional2

is the ground state density.

2.2

The Kohn-Sham equations

While the theorems above are in themselves rather theoretical and abstract, the Kohn-Sham equations provide the means to use them to solve actual many-body problems. The main idea behind the Kohn-Sham equations is that the interacting many-body system could be replaced by an auxiliary non-interacting system with the same ground state density. The Hohnberg-Kohn theorems then ensures that we can also obtain the ground state energy.

The total ground state energy functional could be written as [3]: E[n(r)] = Ts[n] + Z Vext(r)n(r)dr + Exc[n] + 1 2 Z Z _n(r)n(r0₎ |r − r0_| drdr 0 _(2.4)

Where Ts[n] is the kinetic energy of non-interacting particles, and Exc[n] is the

correlation energy for a system of interacting electrons. The exchange-correlation term contains all the manybody-effects of the system. This term can not be determined exactly, but it can often be approximated with sufficient accu-racy (see section 2.3)

To transform the interacting system in equation (2.4) to a non-interacting equivalent system, a new effective potential, Vef f is introduced. It can be shown

that this potential must be:

Vef f(r) = Vext(r) + Vxc(r) +

Z _n(r0₎

|r − r0_|dr

0 _(2.5)

2_{A functional is a mapping from a vector space, usually of functions, to real space. That is,} a functional takes functions as its argument, and outputs a scalar value.

2.3 The Exchange-Correlation energy functional 7

The Schr¨odinger equation for the system of non-interacting particles is: n − ~ 2 2m∇ 2_{+ V} ef f(r) o ψi(r) = εiψi(r) (2.6)

When we have the solutions to equation (2.6) we can obtain the electron charge density by the equation

n(r) =

N

X

i=1

|ψi(r)|2 (2.7)

The three equations (2.5), (2.6) and (2.7) are known as the Kohn-Sham equations and can be solved iteratively until convergence is reached. (See chapter 3.) It should be noted that the solutions ψi to equation (2.6) do not represent real

elec-trons, but non-interacting quasi-particles used to obtain the total charge density. The one-electron eigenvalues εi has thus no direct physical meaning.

One can show that the total energy functional of equation (2.4) can be rewritten as [2]: E[n(r)] = N X i εi− 1 2 Z Z _n(r)n(r0₎ |r − r0_| drdr 0₋Z _V xc(r)n(r0)dr + Exc[n(r)] (2.8)

This equation can be used to obtain the total energy once the Kohn-Sham equa-tions have been solved.

2.3

The Exchange-Correlation energy functional

The exchange3 _{and correlation}4_{effects of the many-body system are contained in}

the exchange-correlation term Exc in equation (2.4). As mentioned earlier, this

term can not be determined exactly, but there exists several approximations of it.

2.3.1

The Local Density Approximation (LDA)

The LDA is the simplest approximation to Exc but has proven to be sufficiently

accurate for a wide range of problems. In the LDA the exchange-correlation energy of the system is calculated as:

E_xcLDA[n] = Z

εxc[n(r)]n(r)dr (2.9)

where εxc[n(r)] is the exchange-correlation energy per particle in a uniform

elec-tron gas. This means that in every point being summed over in the integral, the exchange-correlation energy is the same as in a uniform electron gas with the same density as the density at that specific point in the real system. This approximation works well in many cases, but the LDA is known to produce wrong results for some systems. One of the most common improvements of the LDA is the generalized gradient approximation.

3_{Due to the Pauli repulsion.}

2.3.2

The Generalized Gradient Approximation (GGA)

The GGA is a modification of the LDA, which also takes into account the rate at which the electron density changes in every point. This is done by including the gradient of the electron density in the equation:

E_xcGGA[n] = Z

f [n(r), ∇n(r)]n(r)dr (2.10)

where the functional f [n(r), ∇n(r)] is determined by parametrization and fitting to calculations. Several different methods for doing this exist.

Chapter 3

Computational methods

3.1

Solving the Kohn-Sham equations numerically

In order to solve the Kohn-Sham equations numerically, equation (2.6) must be transformed from a differential equation into a linear-algebraic equation, which is more suitable for solving on a computer. This is done by first expanding the one-electron wave functions in some basis set:

|ψii =

X

j

cj|ϕji (3.1)

Choosing a good basis set for the expansion is of great importance when doing the actual calculations. A simple basis set, such as ordinary plane-waves, will require more terms to accurately describe the wave function than a more complex one. Which way to go depends on the problem at hand. The basis set used for the calculations in this project is described in section 3.2.

By inserting the expansion (3.1) into the one-particle Schr¨odinger equation [eq. (2.6)] and multiplying from the left with the bra-vector hϕj| we get the

fol-lowing equation:

X

j

cj{Hkj− εiOkj} = 0 (3.2)

where Hkj are the matrix elements of the Hamiltonian, Okj is the overlap matrix

and εi is the one-particle eigenvalue for the i :th orbital. The equation will only

have (non-trivial) solutions if

detnHkj− εiOkj

o

= 0 (3.3)

which is the linear-algebraic variant of the differential Schr¨odinger equation.

3.1.1

Finding solutions self-consistently

The Kohn-Sham equations must be solved self-consistently. This process could be briefly described as:

1. Make a qualified guess for the electron density.

2. Feed that guess into the equation for the effective potential. [eq. (2.5)]. 3. Use the obtained potential for solving the Schr¨odinger equations for the

non-interacting particles [eq. (2.6)].

4. Calculate the total charge density with equation (2.7).

5. Repeat the process from step 2, using the recently calculated charge density as a new (better) guess.

At the end of each iteration the newly obtained charge density is compared to the one from the previous iteration. If the difference is sufficiently small, the calculations are said to be converged and the solution has been found. The required precision in the calculations of course determines what is “sufficiently small”. When feeding the output density directly back into the loop there can sometimes be problems with oscillations near the minimum, which will make the calculations take longer or perhaps not finish at all. One way of avoiding this is to use so called mixing, which means that the newly obtained density is mixed to some extent with the previous densities. There are several different schemes for doing this.

3.2

The Exact Muffin-Tin Orbitals Method

The choice of basis set will usually depend heavily on the potential of the specific system being treated. A basis set that more correctly describes the potential will require fewer terms in the expansion. In the computer codes used for this diploma project the so called Exact Muffin Tin Orbitals Method (EMTO) is used. The EMTO theory was developed by O. K. Andersen [1] and is built around the so called muffin-tin approximation. In this approximation the actual Kohn-Sham potential is approximated with a spherically symmetric one, centered around the atom core. (These are the “muffin-tins”.) Beyond the muffin-tin radius, in the interstitial region, the potential is approximated as constant. The EMTO method improves upon this simple model by allowing the muffin-tins to overlap in a manner that best represents the real Kohn-Sham potentials [12]. When using the EMTO method, the Kohn-Sham wave function ψi(r) is expanded in terms of the exact

muffin-tin orbitals: ψi(r) = X RL ¯ ΨaRL(εi, rR)uaRL,i (3.4)

3.2 The Exact Muffin-Tin Orbitals Method 11 Different basis sets are used in the different regions of the muffin-tin potentials [3]:

ψi(r) = X RL φ(εi, rR)Θ(rR)uaRL,i (3.5) +X RL Ψa_RL,i(κi, rR) h 1 − Θ(rR) i v_RL,ia where Θ(rR) = 

1, if rR inside Muffin-Tin sphere

0, if rR outside Muffin-Tin sphere

(3.6) Here εiis the one-particle energies, rRis the distance from the center of the nucleus

R, L denotes the l and m quantum numbers and κ is defined as κ2= ε − VM T Z,

where VM T Z is the constant potential between the muffin-tins. The coefficients

ua_RL,iand va_RL,iare determined from the condition that the wavefunction ψ(r) and its first derivatives must be continuous at the potential (muffin-tin) sphere.

The basis functions φ(ε, rR) have the form:

φ(ε, rR) = φRL(ε, rR) · YL(ˆr) (3.7)

where the φRL(ε, rR) are solutions to the radial Schr¨odinger equation and YL(ˆr)

are spherical harmonics.

The basis set used outside the muffin-tin spheres consists of so-called screened spherical waves, that are solutions to the Helmholtz equation:

 1 2∇

2_{+ κ}2 _Ψa

RL,i(κi, rR) = 0 (3.8)

The boundary conditions for this equation are defined together with non-overlapping hard-spheres with radii aRl. The screened spherical waves behave

like pure real spherical harmonics on their own hard-spheres, while the projections of the spherical harmonics on all other hard-spheres vanish. The screened spherical waves form a complete basis set in the region between the hard spheres.

The connection between the hard-spheres and the muffin-tin spheres is done with the aid of another wavefunction, ϕRl(ε, rR)YL(ˆrR). At the muffin-tin sphere,

this additional wavefunction joins continuously and differentiably to the partial wave, while at the hard-sphere it joins continuously to the screened spherical wave [3].

Finally, the exact muffin-tin orbitals can now be written: ¯ Ψa_RL(εi, rR) = h φRL(ε, rR) − ϕRl(ε, rR) i YL( ˆrR) + ΨaRL,i(κi, rR) (3.9)

where the radial part of φRL is truncated outside the muffin-tin sphere, and the

radial part of ϕRlis truncated outside the hard-spheres. The ¯ΨaRL:s are continuous

and differentiable in whole space, except at the hard-sphere boundary, where the wave function will have kinks. The requirement that these kinks must vanish leads to the so-called screened Korringa-Kohn-Rostocker equations, whose solutions will yield the one-electron eigenvalues and eigenfunctions.

3.3

The Coherent Potential Approximation

One of the major obstacles to overcome when doing ab-initio calculations on alloys is how to represent random alloys. In a perfectly random alloy of elements A and B, where the concentration of A-atoms is c, the probability to find an A-atom at a site is c and the probability to find a B-atom is (1 − c). In order to completely accurately describe this alloy one needs to do calculations over an, in the ideal case, infinite lattice with randomly placed atoms. Since this is clearly not feasible it’s necessary to use some kind of approximation.

One common way to approximate a random alloy is the Coherent Potential Approximation (CPA), originally proposed in 1966 by Paul Soven [11]. The CPA is a so-called effective medium approximation. This means that the real medium is replaced with an effective one, in which each of the atoms in the alloy is inserted one at a time. The effective medium should represent the real environment as accurately as possible. This is done by replacing the actual randomly distributed potentials at each site with an “average” one. A schematic figure of this is shown in Fig. 3.1. Creating the effective medium is not a trivial task. In the CPA the effective medium is constructed from the requirement that an electron traveling through the effective medium should, on the average, not be scattered further if one of the effective potentials was replaced by a real one [2].

Figure 3.1. Schematic drawing of an effective medium. Black and white circles represent the real atoms, while the gray circles represent ”average” effective-medium atoms. The variable c is the concentration of ”white” atoms. The picture is based on figure 3.4 in Ref. [3].

Since each of the real potentials are inserted into the effective medium one at a time, the effective medium approaches are single site approximations. The main drawback with single site approximations is that local effects can not be treated, since there is no information about the local environment.

Chapter 4

Elastic constants

This chapter will give a brief overview of the theory of elastic constants and how they can be calculated using first principles techniques. The first section in this chapter follows to the most part the treatment of elastic properties in Introduction to Solid State Physics by C. Kittel [9].

4.1

Basic theory

The two most fundamental quantities when discussing elastic properties of solids are stress and strain. The stress is defined as the force acting on a unit area in the solid [9]. The deformation that is the result thereof is the strain. The basic theory of elasticity in solid materials is based upon Hooke’s law, which states that in an elastic solid the strain is directly proportional to the stress, provided that the strains are small.

We start by considering an elastic solid that is subject to small strain so that the three orthogonal vectors ˆx, ˆy, ˆz have been slightly distorted. The new axes can then be written as

x0 = (1 + xx)ˆx + xyy + ˆ xzzˆ

y0 = yxx + (1 + ˆ yy)ˆy + yzˆz (4.1)

z0 = zxx + ˆ zyy + (1 + ˆ zz)ˆz

where the coefficients αβ define the deformation. The displacement R of the

deformation can be defined as:

R ≡ x(x0− ˆx) + y(y0− ˆy) + z(z0− ˆz) (4.2)

= (xxx+ yyx+ zzx)ˆx + (xxy+ yyy+ zzy)ˆy + (xxz+ yyz+ zzz)ˆz

We introduce u, v, w such that the displacement can be written:

R = u(r)ˆx + v(r)ˆy + w(r)ˆz (4.3)

By Taylor expansion of R around R(0) = 0 we get xxx≈ x ∂u ∂x ; yyx≈ x ∂u ∂y ; · · · (4.4)

The strain components exx, eyy, ezz are defined by the relations:

exx≡ xx= ∂u ∂x ; eyy ≡ yy= ∂v ∂y ; ezz≡ zz= ∂w ∂z (4.5)

The other strain components eαβ are symmetrical (eαβ= eβα) and are defined

as: exy≡ x0· y0∼= yx+ xy= ∂u ∂y + ∂v ∂x eyz ≡ y0· z0∼= zy+ yz= ∂v ∂z+ ∂w ∂y (4.6) ezx≡ z0· x0 ∼= zx+ xz= ∂u ∂x+ ∂w ∂x

Since the coefficients  are small, terms of order 2_{and higher have been neglected.}

It’s clear that the strain components in equation (4.5) are simply defined as the change of length per unit length of the three spatial axes. The other strain com-ponents eαβ are defined as the change in angle between the axes α and β [8].

We now turn to the stress components, who are here denoted as Xx, Xyetc. In

this notation Xxis a force acting in the x direction on a unit plane perpendicular

to the x direction while Xy is a shearing force acting in the x direction on a unit

plane perpendicular to the y direction, and so on. The nine stress components are reduced to six by the requirement that the net torque of a system at rest must be zero, which means that some of the stress components must be equal. The independent stress components could then be chosen as Xx, Yy, Zz, Yz, Zx, Xy.

According to Hooke’s law the stress components are linear functions of the strain components: Xx= C11exx+ C12eyy+ C13ezz+ C14eyz+ C15ezx+ C16exy Yy= C21exx+ C22eyy+ C23ezz+ C24eyz+ C25ezx+ C26exy Zz= C31exx+ C32eyy+ C33ezz+ C34eyz+ C35ezx+ C36exy (4.7) Yz= C41exx+ C42eyy+ C43ezz+ C44eyz+ C45ezx+ C46exy Zx= C51exx+ C52eyy+ C53ezz+ C54eyz+ C55ezx+ C56exy Xy= C61exx+ C62eyy+ C63ezz+ C64eyz+ C65ezx+ C66exy

The quantities C11, C12, . . . are called elastic stiffness constants or just elastic

constants and have the dimensions force per unit area1_.

4.1 Basic theory 15

In the approximation of Hooke’s law, the elastic energy density is a quadratic function of the strains:

U = 1 2 6 X λ=1 6 X µ=1 ˜ Cλµeλeµ (4.8)

The indices 1–6 are defined as:

1 ≡ xx 2 ≡ yy 3 ≡ zz 4 ≡ yz 5 ≡ zx 6 ≡ xy Note that ˜Cλµ is not the same as the elastic constants Cλµ.

By the definition of potential energy we have that the stress components are found as the derivative of U with respect to the corresponding strain component. Ap-plying this for the stress component Xxyields:

Xx= ∂U ∂exx ≡ ∂U ∂e1 = ˜C11e1+ 1 2 6 X β=2 ( ˜C1β+ ˜Cβ1)eβ (4.9)

Only combinations of the form1₂( ˜Cαβ+ ˜Cβα) appear in the expression for the stress

component. It is easily verified that the same holds for all the stress components. This shows that the elastic constants are symmetrical:

Cαβ=

1

2( ˜Cαβ+ ˜Cβα) = Cβα (4.10) The thirty-six elastic constants are thus reduced to twenty-one.

Since we have ( ˜Cαβ+ ˜Cβα) = 2Cαβwe can get rid of the ˜C:s in equation (4.8) and

replace them with the ordinary elastic constants C by introducing a new notation for the strain components2:

  u1 u6 u5 u6 u2 u4 u5 u4 u3  =   e1 2e6 2e5 2e6 e2 2e4 2e5 2e4 e3  =   exx 2exy 2ezx 2exy eyy 2eyz 2ezx 2eyz ezz   (4.11)

Now the energy density U can be expressed as:

U =1 2 6 X α=1 6 X β=1 Cαβuαuβ (4.12)

4.1.1

Elastic constants of cubic systems

Cubic systems have four symmetry axes across the four diagonals through the cube. Rotating a cubic system about these axes is equivalent to interchanging the x, y, z axes according to four different schemes, one for each symmetry axis. As some of the interchanged axes are negative, some of the quadratic terms of strain 2_{Note that the u:s here have nothing to do with the u in equation (4.3). The choice of symbols} have been made to keep notations consistent with the sources used.

components will change sign under rotation. The requirement that the elastic energy density should be invariant under rotation implies that these terms must disappear. For this to be satisfied, it can be shown3 _{that there can only be three}

independent elastic constants in a cubic system, namely C11, C12and C44. All the

elastic constants in equation (4.7) for a cubic system are presented below.

exx eyy ezz eyz ezx exy Xx C11 C12 C12 0 0 0 Yy C12 C11 C12 0 0 0 Zz C12 C12 C11 0 0 0 Yz 0 0 0 C44 0 0 Zx 0 0 0 0 C44 0 Xy 0 0 0 0 0 C44

As can be seen from the table, the elastic constants C11and C12determines how

the material responds to stress perpendicular to one face, while C44 determines

how much shearing that will result from stress parallel to one face. C11determines

the strain in the parallel direction of the stress, while C12 determines the strain

in the perpendicular direction of the stress. The geometrical interpretation of the elastic constants is depicted in figure 4.1.

4.2

Calculating elastic constants

When using ab-initio software to calculate elastic constants, the crystal is distorted in such a way that the strain is governed completely by the specific elastic constant under consideration. The program output is usually the total energy per atom in the system. Equation (4.12) is therefore multiplied by the volume of a unit cell to get the total energy instead of the energy density. Also, the energy difference between the distorted and undistorted crystal is used in the final calculations, rather than the total energy. For many deformations, the relation between the energy difference and the elastic constant can be written [3]

∆E ≈ A · V · C · δ2 (4.13)

where A is a constant that depends on the specific crystal type and deformation, V is the volume of the system (unit cell), C is the elastic constant under consideration and δ represents the corresponding strain component(s).

The actual value of the elastic constant is obtained by calculating the energy for several small deformations and plotting ∆E as a function of δ2. By performing linear fitting one can then obtain the slope of the curve, which can be divided by A and V to obtain C.

4.3 The Face Centered Cubic lattice 17 Stress Strain in parallel direction, determined by C11 Strain in perpendicular directions, determined by C12 Stress Shear strain, determined by C₄₄ Undistorted Crystal

Figure 4.1. Geometrical interpretation of the elastic constants C11, C12and C44 for a cubic system. The figure shows the cases of no stress (left), stress perpendicular to one face (middle) and stress parallel to one face (right).

4.3

The Face Centered Cubic lattice

When doing the actual calculations there are often many different deformations that can be used to obtain a specific elastic constant. A good choice of deformation can often reduce calculation time and work effort. In order to be able to calculate an elastic constant as a function of volume (or pressure), one must use volume conserving deformations. (See chapter 5.) By using the elastic constant C0, defined as:

C0= 1

2(C11− C12) (4.14)

one can calculate both C11 and C12 at once by applying the following volume

conserving deformation to the fcc lattice:

u1= δ , u2= −δ , u3=

δ2

1 − δ2 (4.15)

C11and C12 can then be obtained by the following relations [3]:

C11= B + 4 3C 0 _(4.16) C12= B − 2 3C 0 _(4.17)

where B is the bulk modulus, which can be obtained directly from the total energy as a function of volume. See chapter 5.

To calculate C44 one can use a body centered orthorhombic (bco) structure

with c/a =√2, which is equivalent to an fcc crystal. The deformation used is:

u3=

δ2

1 − δ2 , u6= 2δ (4.18)

which is also volume conserving.

4.3.1

Polycrystalline elastic constants

Real macroscopic samples of materials are often polycrystalline, meaning that they are not composed of one single, continuous lattice, but of several smaller crystal grains, oriented in different directions. This in turn means that the elastic properties of a macroscopic sample of a material is usually direction-dependent or anisotropic. In this work, the Voigt-Reuss-Hill definition of anisotropy is used:

A = GV − GR GV + GR

(4.19) where GV and GR are the Voigt and Reuss shear moduli respectively. The Voigt

and Reuss shear moduli are defined as follows [3]: GV = C11− C12+ 3C44 5 = 2C0+ 3C44 5 (4.20) GR= 5(C11− C12) · C44 4C44+ 3(C11− C12) = 5C 0_C 44 2C44+ 3C0 (4.21) When C0 and C44 have been calculated using the methods described earlier,

equation (4.19) can be used to calculate a measure of the anisotropy of the material. The anisotropy is an important property, since it’s one of the few properties of the Earth interior that can be directly measured by analyzing the propagation of sound waves traveling in different directions through the Earth.

Chapter 5

Computational details

5.1

The EMTO software

The ab-initio calculations in this work were carried out with the EMTO 5.6 soft-ware package written by L. Vitos et al. As the EMTO softsoft-ware is very large and complex, the purpose of this chapter is not to go into detail in describing how to set up the calculations, but rather to give an overview of how the software was used. A thorough description of the theories behind the EMTO software is given in Ref. [12], written by the author of the software.

The core purpose of the EMTO computer code is to solve the Kohn-Sham equations [equations (2.5)–(2.7)] iteratively, as outlined in section 3.1.1. This is done by expansion into the EMTO basis set described in section 3.2. As the GGA sometimes introduces problems with convergence, the EMTO software uses the LDA to approximate the exchange-correlation term while solving the Kohn-Sham equations, and the GGA when calculating the total energy [equation (2.8)].

All calculations were done without magnetism taken into account, since iron is nonmagnetic at high pressures. The nuclei were approximated as stationary (see section 2.1), which effectively means that the calculations were carried out for a temperature of 0 K.

5.2

Equation of State

An equation of state (EOS) is a thermodynamic equation, describing the relation between two state variables, such as temperature, pressure, volume or internal en-ergy. Many (approximate) equations of state exist, describing the thermodynamic properties of different classes of materials at different conditions.

When calculating the EOS using the EMTO software, the (average) total ground state energy per atom is calculated for several volumes. This yields a ”plot” of the relation between volume and energy for a few selected volumes. The parameter used for adjusting the volume is the so called Wigner-Seitz radius, RW S, which is

defined as the radius of a sphere with the same volume as the Wigner-Seitz cell of 19

the lattice1 _{The calculations were performed for Wigner-Seitz radii of 2.10 – 2.65}

au2_{, in steps of 0.05 au.}

With the data from the ab-inito calculations, a more precise EOS can be inter-polated. The program used for this was written by I.A. Abrikosov. The program allows fitting of the plot data in a number of ways, such as polynomial fit, cubic splines or fitting against different idealized equations of state. In this work, the data was fitted against the Birch-Murnaghan EOS [5], which is an EOS frequently used to model condensed matter under high pressure.

With knowledge of the energy and volume, the pressure and bulk modulus can easily be calculated. This allows for some basic elastic properties to be obtained directly from the EOS, such as the bulk modulus as a function of pressure.

5.3

Elastic constants of the fcc lattice

The elastic constants were calculated with the method described in section 4.2, for the same set of volumes (Wigner-Seitz radii) as the EOS.

For each volume the distortion (δ) was varied between 0 and 0.05 in steps of 0.01. A simple shell script was used to extract the energy difference as a function of δ2 and the data was fitted against a linear function using Octave. By dividing the slope of the curve with the constants A and V in equation (4.13) the elastic constant could be obtained. The volume V is simply the volume of the sphere corresponding to the Wigner-Seitz radius, and the constant A was equal to 2 (two) for both deformations (C0 and C44).

From the EOS [B = B(V )] and the plot of C0 as a function of volume, C11

and C12 as functions of volume could then be obtained by applying equations

(4.16) and (4.17) ”point-by-point” for each volume. From this, plots of the elastic constants as functions of pressure could also be constructed.

5.4

K-point convergence

To ensure accurate results when using the EMTO software, one must make sure that the resolution in k-space is sufficient. This is done by performing the entire set of calculations at least twice when calculating some property, and letting one of the sets use a higher number of ”k-points” to build the k-space mesh. If the end results from the two sets differ only negligible, then the number of k-points is sufficient. Otherwise the process needs to be repeated with a higher number of k-points, until convergence is reached.

Figure 5.1 shows an example of k-convergence testing for the FeMg alloy. The results for 25x25x25 and 21x21x21 k-points are indistinguishable, showing that the calculations are properly converged with respect to the number of points used to build the reciprocal space.

1_{In other words, the Wigner-Seitz radius is the radius of a sphere with a volume equal to the} volume per atom in the lattice.

5.4 K-point convergence 21 2.1 2.2 2.3 2.4 2.5 2.6 2.7 R ws (au) 0 100 200 300 400 500 Pressure (GPa) 25x25x25 21x21x21

Figure 5.1. An example of k-convergence for the FeMg alloy, showing P as a function of RW Sfor both 25x25x25 and 21x21x21 k-points

Chapter 6

Results and discussion

6.1

Equation of State

Figure 6.1. Volume as a function of pressure for fcc Fe90Mg10and pure Fe.

The equation of state [V = V (P )] for Fe90Mg10is shown in figure 6.1. The EOS

for Pure Fe has also been included for reference. The results for pure iron have been calculated the same way as with FeMg. All the plots have been interpolated using cubic splines in Grace. The point at RW S = 2.10 au has not been included,

since the EOS-program had problems fitting the data to the Birch-Murnhagan EOS when that point was included.

From the plot in Fig. 6.1, it can be seen that at the same pressure, the volume of Fe90Mg10 is larger than that of Fe, and that the volume converges towards that

of iron with increasing pressure. This seems reasonable, since it is the large size mismatch between the Fe and Mg atoms that prevents the FeMg-alloy to form at ambient pressure [6]. With increasing pressure, the Mg atoms are ”compressed”, leading to smaller volume of the alloy and increased miscibility.

6.2

Elastic constants

Bulk modulus (GPa)

Pure Fe Fe₉₀Mg₁₀

Figure 6.2. Bulk modulus as a function of pressure for fcc Fe90Mg10and pure Fe.

Figure 6.2 shows the bulk modulus as a function of pressure. The bulk modulus of Fe90Mg10is about 55 GPa lower than that of pure Fe in the whole pressure range,

showing that the FeMg-alloy is overall softer than pure iron.

Figure 6.3 shows the results for the elastic constants C11, C12, C44and C0. The

two elastic constants calculated directly, C44 and C0, are shown in the two upper

graphs of Fig. 6.3. As can be seen from the figure, both are considerably lower than for pure Fe, and the difference increases with pressure. At a pressure of 300 GPa, corresponding to Earth core conditions, C44is 147 GPa or 17% smaller than

for pure Fe. The relative change of C0 is even bigger, with a decrease of 87 GPa or 37%.

To verify that no mistake had been made in the setup of the calculations, an entire set of identical calculations for C44 and C0 were made, with the

concen-trations changed to 100% Fe, 0%Mg. Each value from the test run are displayed with an ”×” in Fig. 6.3. (The point at RW S = 2.50 had problems converging and

has thus been left out.) As can be seen from the figure, the test reproduces the results for pure Fe from Ref. [3] to a good degree. For C0 the results are almost identical, while for C44there is a small difference. The reason for this is unknown,

6.2 Elastic constants 25 C₄₄ 0 250 500 750 1000

Elastic constant (GPa)

Fe Fe 90Mg10 Fe - test C’ 0 100 200 300 0 100 200 300 400 C 11 0 500 1000 1500 2000 Pressure (GPa) 0 100 200 300 400 C 12 0 500 1000 1500

Figure 6.3. C11, C12, C44and C0 as a function of pressure. The results for pure Fe are taken from Ref. [3]. The values obtained in the test run with pure Fe are marked with an ”×”.

one possible explanation is that the reference results for Fe were calculated with an older version of the EMTO software.

The lower two graphs shows the two elastic constants C11 and C12. An

inter-esting property of the FeMg alloy is that C11 is significantly lower than for pure

Fe, also showing an increasing difference with pressure as with C44and C0, while

C12 remains almost unchanged and is even a little higher than for pure Fe at high

pressures. Since C11 and C12 can not be independently measured by experiment,

it is difficult to say how this change in the ratio between C11 and C12 affects the

directly observable elastic properties.

6.3

Polycrystalline shear moduli and anisotropy

The averaged shear moduli GV and GRare shown in figure 6.4. The results follows

the same trend as C0 and C44with substantially lower values for FeMg compared

to pure iron, and an increasing difference with pressure. At a pressure of 300 GPa the difference is 123 GPa or 20% and 134 GPa or 31% for GV and GRrespectively.

Finally, figure 6.5 shows the anisotropy of the FeMg alloy. As can be seen from the graph, the addition of 10% magnesium gives a huge increase in anisotropy. At 300 GPa the anisotropy is 0.073 higher than for pure iron, corresponding to an increase of 40%. The increase of anisotropy with pressure is also considerably higher for the FeMg alloy compared to pure Fe.

The anisotropy-increasing effect of magnesium on iron is interesting, since the Earth is known to be elastically anisotropic; sound waves travels 3–4% faster along the spin axis of the Earth than in the direction along the equatorial plane [4]. Of course, these calculations were performed for a temperature of 0 K, which is very far from the temperature of the Earth core (about 6000–7000 K). One can thus not draw any direct conclusions on possible effects of magnesium on the anisotropy of the Earth core from these results, but they can perhaps serve as a starting point for further research.

6.3 Polycrystalline shear moduli and anisotropy 27 0 100 200 300 400 Pressure (GPa) 0 100 200 300 400 500 600 700 800 Shear moduli Pure Fe: GV Pure Fe: GR Fe90Mg10: GV Fe90Mg10: GR

Figure 6.4. Voigt and Reuss shear moduli as a function of pressure for fcc Fe90Mg10 and pure Fe.

0 100 200 300 400 Pressure (GPa) 0 0,05 0,1 0,15 0,2 0,25 0,3 Anisotropy Fe₉₀Mg₁₀ Pure Fe

Chapter 7

Conclusions and outlook

7.1

Conclusions

From the results in chapter 6 one can draw the conclusion that magnesium has a softening effect on fcc iron. For all elastic constants except the bulk modulus, the softening effect also increases with pressure.

The overall conclusion of this work is that even small amounts of magnesium can have a large effect on the elastic properties of iron, and that these effects also appear to increase with pressure.

7.2

Outlook

One obvious shortcoming of these results is that the method used requires the cal-culations to be performed for a temperature of zero Kelvin. To be able to further investigate the possibility of magnesium in the Earth core, one has to perform calculations with realistic Earth-core temperatures.

The reasons for the increasing softening effect of magnesium with pressure and the large jump in anisotropy, as well as the reason why the C12 elastic constant

remains almost unchanged, are not known and may be investigated further. Finally, it is of course necessary to investigate the properties of FeMg alloys in other crystal structures than fcc, and with different concentrations of magnesium.

Bibliography

[1] O.K. Andersen, O. Jepsen, and G. Krier. Exact muffin-tin orbital theory. World Science Publishing Co., 1994.

[2] Christian Asker. First principles calculation of core-level shifts in random metallic alloys: The transition state approach. Master’s thesis, Link¨oping University, 2004.

[3] Christian Asker. Spectroscopic and Elastic Properties in metallic systems From First-Principles methods. Licentiate thesis, Link¨oping University, 2007. [4] Anatoly B. Belonoshko, Natalia V. Skorodumova, Anders Rosengren, and B¨orje Johansson. Elastic anisotropy of earth’s inner core. Science, 319:797– 800, 2008.

[5] Francis Birch. Finite elastic strain of cubic crystals. Physical Review, 71:809– 824, 1947.

[6] N. Dubrovinskaia, L. Dubrovinsky, I. Kantor, W. A. Crichton, V. Dmitriev, V. Prakapenka, G. Shen, L. Vitos, R. Ahuja, B. Johansson, and I. A. Abrikosov. Beating the miscibility barrier between iron group elements and magnesium by high-pressure alloying. Physical Review Letters, 95(245502), 2005.

[7] L. Dubrovinsky, N. Dubrovinskaia, O. Narygina, I. Kantor, A. Kuznetzov, V. B. Prakapenka, L. Vitos, B. Johansson, A. S. Mikhaylushkin, S. I. Simak, and I. A. Abrikosov. Body-centered cubic iron-nickel alloy in earth’s core. Science, 316:1880–1883, 2007.

[8] Y. C. Fung and P. Tong. Classical and computational solid mechanics. World Scientific, 2001.

[9] C. Kittel. Introduction to solid state physics. John Wiley & Sons, 8th edition, 2005.

[10] A. S. Mikhaylushkin, S. I. Simak L. Dubrovinsky, N. Dubrovinskaia, B. Jo-hansson, and I. A. Abrikosov. Pure iron compressed and heated to extreme conditions. Physical Review Letters, 99(165505), 2007.

[11] Paul Soven. Coherent-potential model of substitutional disordered alloys. Physical Review, 156:809–813, 1967.

[12] L. Vitos. Computational Quantum Mechanics for Materials Engineers. Springer, 2007.