Document [original]

Magnetic Systems studied by First-Principles Thermodynamics

Fritz Körmann

Paderborn 2011

Magnetic Systems studied by First-Principles

Thermodynamics

Dissertation

zur Erlangung des akademischen Grades

Doktor der Naturwissenschaften (Dr. rer. nat.)

vorgelegt dem

Department Physik der Fakultät für Naturwissenschaften

an der Universität Paderborn

Fritz Körmann

Promotionskommission

Vorsitzender Prof. Dr. Klaus Lischka

Gutachter Prof. Dr. Jörg Neugebauer

Gutachter Prof. Dr. Arno Schindlmayr

Gutachter Prof. Dr. Igor A. Abrikosov

Tag der Einreichung: 14. Februar, 2011

Tag der mündlichen Prüfung: 03. Mai, 2011

TEXT IN WHITE in order to keep this page empty

To my family...

Danksagung

Zuallererst danke ich Herrn Prof. Dr. Jörg Neugebauer nicht nur für das spannende Thema

dieser Dissertation, sondern auch dafür, dass ich viel von ihm lernen durfte. Er hat die Arbeit durch

gezielte Fragestellungen in den entscheidenden Phasen auf Kurs gehalten. Mein besonderer Dank

gilt auch Tilmann Hickel und Alexey Dick. Durch viele Gespräche haben sie die Arbeit entscheidend

mitgeprägt. Ich möchte Beiden für die unbeschreibliche Unterstützung während der gesamten Zeit

(nicht nur bei physikalischen Fragestellungen) als auch für detaillierte Verbesserungsvorschläge zu

dieser Arbeit danken.

Blazej Grabowski danke ich für unzählige Diskussionen rund um die Physik, insbesondere für

die Einführung in die Berechnung vibronischer und elektronischer Entropiebeiträge. Nicht uner-

wähnt möchte ich Ugur Aydin, Sixten Boeck, Björn Lange, Gernot Pfanner, Niko Sandschneider

und Alexander Udyansky lassen (in alphabetischer Reihenfolge), denen ich viele, viele interessante

Gespräche zu verdanken habe. Ein lieber Dank gilt Ugur, dass er mir den größeren Teil der Kaﬀee-

Verwaltung unserer Abteilung abgenommen hat. Frau Helemann danke ich für die Bewältigung

vieler organisatorischer Aufgaben sowie aller Fragen rund um das Thema Fussball. Der gesamten

Abteilung verdanke ich ein ausgezeichnetes Arbeitsklima, welches den entscheidenden Grundstein

für die Arbeit gelegt hat. Albert, Alexander, Björn, Johann, Niko und Ugur danke ich für das Kor-

rekturlesen. Der Max-Planck-Gesellschaft sowie ICAMS bin ich für die ﬁnanzielle Unterstützung

während der Promotion zu Dank verpﬂichtet.

Ausserhalb des Instituts gilt mein Dank vor allem meinen Eltern und Geschwistern, für ihre

stete Unterstützung in allen Zeiten. Mein letzter Dank gilt meiner kleinen Familie. Dank ihrer

kann ich jeden Tag neue Kraft tanken und mich jeden Morgen mit einem Lächeln auf den Lippen

erwischen.

"Zu einer friedlichen Familie kommt das Glück von selbst."

Sprichwort aus China

Abstract

A major challenge in applying state-of-the-art density functional theory (DFT) simulation tools

on real materials is that they have been originally designed to predict ground-state properties.

Since ﬁnite temperature eﬀects are crucial for practically all real-world applications, extending DFT

capabilities towards ﬁnite temperature predictions is indispensable. Within the last few years it has

been shown that a careful consideration of diﬀerent temperature induced excitation processes such as

lattice and electronic contributions can provide extremely accurate ﬁnite-temperature predictions

for nonmagnetic materials. However, for the vast variety of materials, such as metallic alloys,

magnetic shape memory alloys, or diluted magnetic semiconductors, where magnetic interactions

and excitation processes are critical and cannot be neglected, the nonmagnetic approaches are not

suﬃcient. So far practically no theoretical concepts to include magnetic excitations and to bridge

between (i) the complexity of real structural materials, (ii) magnetic theories which are designed to

describe very particular model systems, and (iii) DFT calculations exist.

In this work we have closed the gap between these three ﬁelds. A key challenge is that highly

accurate techniques such as Spin Quantum Monte Carlo work best for frustration-free spin systems,

whereas practically all real world materials show an oscillatory and long-range magnetic interac-

tion resulting into spin frustration when being mapped onto model Hamiltonians. Therefore, ﬁnite

temperature magnetism of real systems is commonly described using classical approaches such as

classical Spin Monte Carlo. These methods work well for temperatures well above the critical (Curie,

Néel) temperature, but fail for low temperatures, where spin-quantization becomes crucial. A key

concern in this work is the correct incorporation of quantum eﬀects into our models which turned

out to be mandatory even up to the critical temperature. We, therefore, developed a hierarchy

of numerically exact quantum Monte Carlo based methods and analytical (Green’s functions) ap-

proaches to treat the magnetic free energy of complex magnetic materials. The proposed approach

allowed us to describe the free energies with a hitherto not achievable accuracy for a wide range

of materials. Successful applications include key materials in steel manufacturing like ferrite and

cementite as well as various magnetic metals (e.g., nickel and cobalt). A careful analysis of these

results allowed us for the ﬁrst time to evaluate the relevance of spin quantum eﬀects on thermody-

namic properties. In contrast to common belief, we could show that spin quantum eﬀects have a

dramatic impact on various quantities such as free energies, heat capacities, and magnetizations, all

the way up to the critical (Curie, Néel) temperature. The presented methodologies also provide a

clear and systematic separation of diﬀerent physically relevant contributions and are even applicable

to phases which are not or only hardly accessible by experiment (such as, e.g., metastable phases),

or systems where experimental input simply does not exist.

Contents

1 Introduction 1

2 Theoretical Background 5

2.1 Many-Body Problem and Density Functional Theory . . . . . . . . . . . . . . . . . . 8

2.1.1 The Many-Body Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.1.2 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2 Extension to Finite Temperatures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.2.1 Magnetic Excitations from DFT . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.2.2 The Heisenberg Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.2.3 Thermodynamic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.2.4 Vibronic and Electronic Free Energy . . . . . . . . . . . . . . . . . . . . . . . 33

3 Results 37

3.1 Determination of TC(p)in bcc iron . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.2 Thermodynamics of Magnetic Metals: Free Energy and Speciﬁc Heat Capacity . . . 46

3.2.1 Analytic Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.2.2 Numerical Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3.3 Structural Phase Transitions in Iron . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

3.3.1 Thermodynamics of γ-Iron . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

3.3.2 Impact of Vibronic, Electronic, and Magnetic Contributions . . . . . . . . . . 77

3.4 Nearest-Neighbor Approach: Application to Fe, Co, and Ni . . . . . . . . . . . . . . 83

3.5 Impact of Magnetism on the Phase Stability of Cementite . . . . . . . . . . . . . . . 89

3.5.1 Thermodynamic Properties of Fe3C....................... 89

3.5.2 Phase stability of Fe3C............................... 93

4 Summary and Outlook 95

A Appendix 101

A.1 Mean Field Approximation for Ferromagnets . . . . . . . . . . . . . . . . . . . . . . 101

A.2 Random Phase Approximation for Ferromagnets . . . . . . . . . . . . . . . . . . . . 104

A.3 Random Phase Approximation for Noncollinear Magnets . . . . . . . . . . . . . . . . 107

Bibliography 113

TEXT IN WHITE in order to keep this page empty

Chapter 1

Introduction

The development and improvement of functional materials, such as, e.g., steels, is a key concern since

thousands of years. To give an example, the earliest known production of steel (i.e., iron with a small

amount of carbon) is about 4000 years old [1]. Today, steel is one of the most important materials in

our daily life (buildings, ships, tools, cars, etc.) with more than 1300 million tons produced annually

worldwide, see Fig. 1.1. Along with the production, the demand on new structural materials with

1900 1925 1950 1975 2000

Year

500

1000

1500

World Steel Production [million tpa]

2008: 1329 million tons

Annual growth rate

of Steel production in [%]

1950 - 60: 6,2

1960 - 70: 5,5

1970 - 75: 1,6

1975 - 80: 2,2

1980 - 85: 0,1

1985 - 90: 1,4

1995 - 00: 2,4

2000 - 05: 5,9

2005 - 10: 3,8

Figure 1.1: World Steel production for 1900-2008. Data taken from Refs. [2] and [3].

well deﬁned properties (e.g., ultra-light weight steels) has also signiﬁcantly increased. Every new

steel generation goes along with new and more complex physical mechanisms often increasing the

involved number of parameters (constituents) as well. A typical example of a modern manganese

steel is given by [4]

Fe77.98Mn20.98C0.71Si0.17Ni0.04N0.02Cr0.02Nb0.02Mo0.01Al0.01Cu0.01Co0.01... (in wt.%)

with 16 ingredients in total. Reﬁnement of such compositions by pure experimental approaches

becomes not only time-consuming, but more and more economically expensive. To cope with this

Figure 1.2: The temperature driven structural phase transitions of pure iron. Details are given in

the text.

large conﬁgurational space, approaches such as calphad (calculation of phase diagrams) have

been developed, which are based on experimental input and sophisticated interpolation schemes.

However, in absence of experimental input, such approaches can fail dramatically.

A very recent approach has attracted a lot of attention starting from the most elementary

microscopic level - quantum mechanics. This approach is only based on the fundamental physical

constants such as nuclei and electron mass and charge, and allows completely parameter-free and

unbiased materials properties prediction without any experimental input. This method is, therefore,

called ﬁrst principles or, using the Latin term, ab initio technique.

A major challenge in applying state-of-the-art density functional theory (DFT) simulation tools

to real materials is that they have been originally designed to predict ground-state properties (T=

0K). Since ﬁnite-temperature eﬀects are crucial for practically all real-world applications, extending

DFT capabilities towards ﬁnite-temperature predictions is indispensable. Fortunately, nowadays

sophisticated combined concepts exist (e.g., the so-called quasiharmonic approximation and ﬁnite

temperature DFT), which allow a highly accurate prediction of properties of nonmagnetic materials

even up to the melting temperature [5]. As shown recently [6], such an approach is, however, not

suﬃcient for the main ingredient of steels and the most famous magnetic material, namely iron.

Iron shows a complex and untypical order of structural phase transitions, illustrated in Fig. 1.2.

At low temperatures iron crystalizes in the bcc structure called α-iron or ferrite. Above the Curie

temperature TC= 1043 K, the paramagnetic bcc Fe (β-iron) transforms at 1184 K into the more

close-packed fcc lattice, also called γ-iron or austenite. At 1673 K Fe undergoes a second structural

transition and becomes bcc again (δ-iron) until it melts at about 1800 K. It is known that in

particular the magnetic excitations are of fundamental importance for the structural transitions

and phase stability [7, 8]. However, the theoretical concepts to predict magnetic contributions to

thermodynamic properties based on DFT input is still in its infancy.

The reason for this is that so far practically no theoretical concepts to include magnetic exci-

tations and to bridge between (i) the complexity of real structural materials, (ii) magnetic theories

which are designed to describe very particular model systems, and (iii) DFT calculations exist. Start-

ing from the microscopic theory, the extension of DFT towards the spin-polarized DFT (SDFT),

ﬁrst introduced by von Barth and Hedin [9], is nowadays well established for the computation of

ground-state properties of magnetic systems [10]. To investigate ﬁnite-temperature properties it is,

however, indispensable to recourse to model considerations. Besides the Ising model and eﬀective

theories such as the Weiss mean ﬁeld theory for ferromagnets, more realistic theoretical descriptions

have become possible by the quantum-mechanical Heisenberg model [11]. In the last two decades a

number of techniques (such as the frozen-magnon approach [12]) emerged, which allow to extract

realistic Heisenberg model parameters from SDFT methods. A key challenge is that highly accurate

techniques to solve the model such as Spin Quantum Monte Carlo work best for frustration-free spin

systems, whereas practically all real world materials show an oscillatory and long-range magnetic

interaction resulting into spin frustration when being mapped onto model Hamiltonians. Therefore,

ﬁnite-temperature magnetism of real systems is commonly described using classical approaches such

as classical Spin Monte Carlo [13–17]. These methods work well for temperatures well above the

critical (Curie, Néel) temperature, but fail for low temperatures, where spin-quantization becomes

crucial.

A key concern in this work is the correct incorporation of quantum eﬀects into our models. We,

therefore, developed a hierarchy of numerically exact quantum Monte Carlo based methods and

analytical (Green’s functions) approaches to treat the magnetic free energy of complex magnetic

materials. The proposed approach allowed us to describe the free energies with a hitherto not

achievable accuracy for a wide range of materials. Successful applications include key materials in

steel manufacturing like ferrite and cementite as well as various magnetic metals (e.g., nickel and

cobalt). A careful analysis of these results allowed us for the ﬁrst time to evaluate the relevance

of spin quantum eﬀects on thermodynamic properties. The presented methodologies also provide a

clear and systematic separation of diﬀerent physically relevant contributions and are even applicable

to phases which are not or only hardly accessible by experiment (such as, e.g., metastable phases),

or systems where experimental input simply does not exist.

This work is organized as follows: After introducing the theoretical fundamentals (Chap. 2), we

focus in Sec. 3.1 on the pressure dependence of the critical temperature TC(p)of pure iron. We

show that the experimentally observed weak dependence on pressure for iron (Texp

C(p)≈0) can

be understood by a competition between two mechanisms. These investigations revealed that in

particular the RPA provides a promising approach to study magnetic properties. The theory, which

was originally designed for predicting the magnetization as well as TC[18], is therefore extended

to compute also quantities such as free energies and heat capacity capacities. The results of the

derived theory are ﬁnally combined with vibronic and electronic contributions (Sec. 3.2.1) in order

to achieve a fully ab initio framework for the thermodynamic properties of bcc iron.

A shortcoming of the analytic RPA treatment is, by construction, the neglect of magnetic short-

range order. A numerically exact Quantum Monte Carlo (QMC) simulation appears as the ideal

candidate for this purpose. We show, however, that this technique is not applicable for a realistic

system such as bcc iron due the negative sign problem. In an extensive comparison between the

numerically exact quantum and classical results of the Heisenberg model for various thermodynamic

properties, the advantages and disadvantages of a classical approach to the quantum Heisenberg

model become apparent (Sec. 3.2.2). Based on the results of the model calculations, a semi-empirical

rescaling function is derived. The rescaling method allows an ad hoc inclusion of quantum eﬀects in

the classical description and its application to bcc iron is discussed in comparison with experimental

data.

The fcc phase of iron is in the focus of Sec. 3.3. Fcc iron constitutes a challenging task for the

theoretical modeling: Recent theoretical calculations reveal a non-collinear spin-spiral ground state

[19], which is a particular challenge for analytic treatments such as the RPA, typically requiring

a global quantization axis (as given for ferro- or antiferromagnets). Working in a locally rotated

coordinate system allows us to generalize the RPA theory, which is developed for bcc iron, to non-

collinear magnetic spin system. We apply both, the generalized RPA as well as the classical Monte

Carlo method to fcc iron, and combine the results with vibronic and electronic contributions.

Particularly, a model study (Sec. 3.2.2) shows that the thermodynamic properties are robust

with respect to the speciﬁc magnetic conﬁguration or the range of magnetic exchange interaction,

if considered on a reduced t=T/TCtemperature scale. This encourages a decoupling of TCfrom

other thermodynamic properties such as the speciﬁc heat capacity. Based on this insight we propose

a nearest-neighbor model, which can be solved quantum-mechanically and numerically exact using

QMC (Sec. 3.4). The shortcomings and advantages of this ansatz are critically evaluated on the

three prototype ferromagnets Fe, Co, and Ni.

Finally, the thermodynamic properties of cementite, the most common precipitate phase in Fe-C

steels (see Fig. 1.2), are investigated in Sec. 3.5. The impact of magnetic contributions for selected

properties such as the speciﬁc heat capacity or formation enthalpy are discussed. This example

highlights that the developed methods in this work not only allow insights into fundamental aspects

of pure elementary magnets such as iron, cobalt, or nickel, but could be employed to address even

properties of complex magnetic materials.

Chapter 2

Theoretical Background

The ﬁrst-principles thermodynamic description of magnetic materials is the central aim of this work.

A key quantity in the thermodynamic description is the thermodynamic potential. In principle its

knowledge provides access to all thermodynamic equilibrium properties. For the methodologies in

this work it is convenient to start with the Helmholtz free energy (Helmholtz potential) F(T, V ),

with temperature Tand volume Vas natural variables. It is completely determined by the partition

sum Z(T, V ) = Tr e−βHas

F(T, V ) = −kBTlnZ(T, V ).(2.1)

The key challenge in order to access F(T, V )is the computation of the trace in Z. This requires

the solution of the Schrödinger equation

HΨν=EνΨν,(2.2)

where Ψν≡Ψν({riσi},{RI})represents the many-particle wave function, {riσi}and {RI}denote

the complete set of coordinates for the electrons with spin σi= (↑,↓)and nuclei respectively, and

Eνrepresents the ν’s energy levels of the system. However, although the microscopic eigenvalue

problem (2.2) can be readily formulated, a direct solution is in general impractical. Describing a

solid involves a many-body problem of an immense number of particles (the number of nuclei Nn

and electrons Neare in the order of ∼1023 particles per cm3). A direct treatment of Eq. (2.2) is

therefore not feasible and simpliﬁcations are mandatory.

One of the most important approximations is the adiabatic decoupling of the ionic, electronic,

and magnetic degree of freedom. This is justiﬁed by considering the diﬀerent time scales of the

underlying physical excitations. For the d-metals considered in the present study, the electronic

excitations are inherently connected to the d-band width Wdrevealing the underlying excitations

in the order of ~/Wd∼10−15 s. The characteristic time of spin-waves are inverse proportional to

their excitation energy, 1/ωSW ∼10−13 s, and vibronic excitations are characterized by their inverse

Debye temperature ∼10−12 s.1As sketched in Fig. 2.1 the adiabatic decoupling for the Helmholtz

free energy results in

F(T, V ) = Fel(T, V ) + Fvib(T, V ) + Fmag(T, V ) + . . . , (2.3)

where the “...” symbolize higher order terms (such as, e.g., coupling terms or vacancy contributions)

and Fel,Fvib, and Fmag denote the electronic, vibronic, and magnetic contribution, respectively.

In a similar way the adiabatic decoupling can be employed on the microscopic level for the

Schrödinger equation (2.2). Regarding the energy scale, the fastest degree of freedom is realized

for the electronic subsystem. Its solution involves still a huge number of particles Ne. A very

elegant and powerful method to solve the electronic subsystem in its ground state was introduced in

1963 by W. Kohn, the reformulation of the many-body problem Eq. (2.2) in terms of the electronic

density, namely the development of the density functional theory (DFT). Although the DFT is

in principle an exact theory, for a practical implementation an approximation for the so called

exchange-correlation (xc) functional is required. For the xc-functional well established and carefully

tested approximations exist nowadays so that DFT has become the most common applied ﬁrst-

principles approach in theoretical solid state physics for the prediction of properties at T= 0 K.

The extension towards ﬁnite temperatures is a challenging and nontrivial task. In order to

address T > 0K properties, and in particular to determine the free energy contributions, the

knowledge of ground-state properties alone is not suﬃcient. At T > 0K the number of possible

microstates increases dramatically making the usage of statistical concepts in combination with DFT

mandatory. In the last decades, the quasiharmonic approximation and the ﬁnite temperature DFT

have been proven to provide a very accurate and realistic approach for Fel and Fvib respectively,

see Fig. 2.1. Although the principle concepts behind were known since long time, the practical

implementation in combination with DFT turns out to be a highly nontrivial and challenging task

itself [5]. The inclusion of magnetic excitations into a full thermodynamic description is, however,

still in its infancy. This is inter alia caused by the fact that even model Hamiltonians capturing the

main underlying excitations are in general not exactly solvable.

This chapter is structured as follows: We start with the formulation of the microscopic many-

body problem Eq. (2.2) and the introduction of the spin-polarized DFT. This provides the funda-

mentals for the further analysis at ﬁnite temperatures. The concepts to determine the electronic

and vibronic contributions based on the DFT input are shortly discussed in Sec. 2.2.4.

In Sec. 2.2 an introduction into the diﬀerent magnetic treatments at ﬁnite temperatures is

given. In particular the Heisenberg model, chosen for the further analysis, is discussed in detail.

Diﬀerent analytic and numerical solution techniques are introduced. The section is concluded with

the adiabatic spin dynamics. It is shown how spin-wave energies and eﬀective exchange coeﬃcients

for the Heisenberg model can be determined by DFT.

1The decoupling of the electronic motion and that of the nuclei (Born-Oppenheimer approximation) is a priori

backed up by the different masses of both particles. Note that a similar mass argument for the decoupling of the

magnetic degree of freedom does not exist.

Thermodynamics

F(T, V ) = −1/β ln Tre−βH

Fvib Fel Fmag

Quasiharmonic

approximation

Finite-temperature

DFT

Quantum

Heisenberg model

DFT for different

ionic configurations

DFT for different

electronic temperatures

DFT for different

magnetic structures

microscopic Schrödinger equation

HΨ=EΨ

Adiabatic decoupling

Dynamical

matrix

Single-particle

energies

Effective

exchange coefficients

Adiabatic decoupling

Figure 2.1: Sketch of the integrated ab initio approach to access thermodynamic properties of mag-

netic materials. The adiabatically decoupled free energy contributions as well as the adiabatically

decoupled treatment via spin density functional theory are shaded in light green. The statistical

concepts to capture the vibronic and electronic contributions are shaded in yellow, whereas the

Heisenberg model, accounting for the magnetic degree of freedom, is highlighted in orange.

82.1 Many-Body Problem and Density Functional Theory

2.1 Many-Body Problem and Density Functional Theory

2.1.1 The Many-Body Problem

In principle, the many-body problem of interacting electrons and nuclei is completely determined

by Eq. (2.2) as

HΨν({riσi},{RI}) = EνΨν({riσi},{RI}).(2.4)

Deﬁning the solid as a system of localized ions (nuclei) and a speciﬁc number of electrons, the

full Hamiltonian2

H=He+HI+HeI +Hz(2.5)

is given by the following contributions:

The Hamiltonian of the electronic subsystem Heis split into the kinetic energy Teand the

electron-electron repulsion term Vee:

He=Te+Vee =−

2me∇2

i+1

j6=i

4πǫ0|riσi−rjσj|.(2.6)

Here ~denotes the reduced Planck constant, methe electron mass, ethe elementary charge, ǫ0the

electric constant (vacuum permittivity) and Nethe total number of electrons.

Analogously the operator of the ionic subsystem HIcan be decomposed into a kinetic Tnand a

repulsion term Vnn as

HI=Tn+Vnn =−

2MI∇2

I+1

J6=I

ZIZJe2

4πǫ0|RI−RJ|,(2.7)

where MIdenotes the mass and ZIthe number of protons of the I-th of the Nnnuclei.

The interaction of the electrons with the nuclei HeI(repulsion term) and the (possible) interaction

with an (external) magnetic ﬁeld Hz(Zeeman term) are deﬁned as

HeI +Hz=−

ZIe2

4πǫ0|riσi−RI|+ 2µB

B(r)σi≡υext({riσi}),(2.8)

where µBdenotes the Bohr magneton, Ban external ﬁeld, σi= (σix, σiy, σiz)the spin operator,

and σi,(xyz)being the Pauli matrices. As already mentioned, the huge number of involved particles

rules out a direct solution of the eigenvalue problem and several simpliﬁcations have to be accepted,

which will be discussed in the following.

The Born-Oppenheimer Approximation

In 1927 Born and Oppenheimer suggested to decouple the electronic and ionic motion. Such an

approximation is justiﬁed due to the huge unbalance between the mass of the electron and of the

nucleus, me≪mI, yield diﬀerent time scales of electronic and nuclei motion. Consequently, the

2This is the non-relativistic formulation of the Hamiltonian. It does not account for relativistic effects such as

spin-orbit coupling. These contributions are, however, not relevant for this work and are neglected in the following.

2.1 Many-Body Problem and Density Functional Theory 9

electronic systems could be solved for a fixed (given) set of nuclei coordinates {˜

RI}. Physically, this

corresponds to the situation where the electrons follow the nuclei motion almost instantaneously

or saying diﬀerently: the nuclei move on the energy surface EBOS

νdeﬁned by the energy of the

electronic subsystem Ee

ν({RI})and the nuclei-nuclei-repulsion Vnn({RI}). The energy surface on

which the nuclei move,

EBOS

ν=Ee

ν({RI}) + Vnn({RI}),(2.9)

is called Born-Oppenheimer surface. It determines together with the kinetic energy the ionic con-

tribution of Eq. (2.4),

HI({RI})Θµν =Tn({RI}) + EBOS

ν({RI})Θµν =EI

µν({RI})Θµν ,(2.10)

providing the EI

µν, where µis running over all eigenenergies for a given potential surface ν. In prac-

tice, one usually considers out of the various potential energy surfaces EBOS

νonly the energetically

lowest one corresponding to the electronic ground state (T= 0 K), i.e., EBOS

ν≡EBOS

In principle, however, an analogous adiabatic decoupling could also be derived for ﬁnite tempera-

tures using the same mass argument as for the derivation of the equation above. It is straightforward

to show [20] that for T > 0K the Born-Oppenheimer approximation leads to

Tn({RI}) + FBOS({RI}, T)˜

Θµ=˜

µ({RI}, T )˜

Θµ,(2.11)

where the temperature dependent Born-Oppenheimer surface is deﬁned as

FBOS =−kBTln X

e−βEBOS

ν.(2.12)

It turns out, however, that ˜

µ({RI})≈EI

µ0({RI}), resulting in a negligible impact of applying

Eq. (2.11) compared to Eq. (2.10) for the vibronic free energy.3In this work Eq. (2.10) is used in

combination with the T= 0 K energy surface EBOS

ν0.

Employing the Born-Oppenheimer approximation, the electronic eigenvalue problem for a given

set of nuclei coordinates ({RI}),

HeΦ{RI}

ν({riσi}) = He+υextΦ{RI}

ν({riσi}) = Ee

ν({RI})Φ{RI}

ν({riσi}),(2.13)

remains to be solved, involving a gigantic number of particles though.

The density functional theory, introduced next, allows an enormous reduction of the involved

degrees of freedom.

2.1.2 Density Functional Theory

In this section we brieﬂy discuss the fundamental aspects of the density functional theory and the

most important involved simpliﬁcations and techniques.

The introduction of the density functional theory (DFT) in 1965 by W. Kohn, is certainly one

3In Ref. [5] it is shown on the specific example aluminum that incorporation of the electronic temperature for the

derivation of Fvib results in energy corrections which are in magnitude smaller compared to the contribution of the

quasiharmonic approximation or its higher order corrections (explicit anharmonic contributions).

10 2.1 Many-Body Problem and Density Functional Theory

of the fundamental breakthroughs in condensed matter physics in the last century. Since magnetic

systems are in the focus of the present study we directly discuss the extension of DFT to spin-

polarized systems, the so called spin-polarized DFT (SDFT), introduced by von Barth and Hedin

[9] in 1972. The main concepts of both versions are, however, very similar.

The basic idea of the SDFT is to express the eigenvalue problem Eq. (2.13) not in terms of the

many-body wave function Φ=Φ({riσi}), but instead in terms of the electronic density

n(r) = hΦ|

δ(r−ri)|Φi,(2.14)

and the magnetization density

m(r) = −2µBhΦ|

σiδ(r−ri)|Φi.(2.15)

Both are for convenience often combined in the compact 2×2spin-density matrix as

n(r) = 1

2(n(r)+σm(r)) .(2.16)

The original version of the First Hohenberg Kohn theorem states that for a nonmagnetic, nonde-

generate electronic system in a given external potential vext (in absence of external ﬁelds), this

potential is uniquely deﬁned (within a trivial constant) by the ground-state electron density n0(r).

Since by deﬁnition the kinetic and electron-electron interaction operators in He(see Eq. (2.6)) are

fully determined by the total number of electrons Ne(which can be derived from the electronic den-

sity nas well), the complete eigenvalue problem Eq. (2.13) can be expressed as a functional of n. A

similar unique correspondence in the SDFT between the spin density matrix nand a given external

potential does, however, not exist. This was already recognized a long time ago by von Barth and

Hedin [9], and formally shown independently by Capelle and Vignale [21] and Eschrig and Pickett

[22]. The nonuniqeness is a generic feature of multidensity DFTs (as the spin-polarized DFT or the

current-DFT) and the consequences are still under investigation [23–27]. There are indications that

the non-invertibility is restricted to exceptional and rather artiﬁcial cases [24]. This assumption

is empirical substantiated by the fact that the SDFT is the most widely used formulation of DFT

in practical applications. If provided that the ﬁrst theorem holds for SDFT as well, the external

potential vext, and hence the total ground-state energy and all ground-state properties could be

expressed as a functional of the ground-state spin density n0. Evidently, this yields a substantial

reduction of the degrees of freedom (from the many-body wave function Φtowards the 2×2spin

density matrix nhaving 12 degrees of freedom).

The Second Hohenberg-Kohn theorem states that the total energy for a given density nis always

greater or equal to the true ground-state energy deﬁned by the corresponding ground-state density

n0. It is based on the Rayleigh-Ritz variational principle [28], according to which for any given

2.1 Many-Body Problem and Density Functional Theory 11

Hamiltonian H(and in particular ˜

Hein Eq. (2.13)) and any trial wave function Φtrial, the expectation

value of hHiΦtrial is always greater than the ground-state energy E0:

Etot[Φtrial]≡hΦtrial|He|Φtriali

hΦtrial|Φtriali≥Etot

0.(2.17)

The theorem allows by virtue of the minimization above a straightforward way to ﬁnd the ground-

state energy by minimizing the functional Etot[Φtrial]with respect to Φtrial. Since the ground-state

density n0uniquely determines the ground-state wave function n0→Φ0,Etot

0is completely deﬁned

by n0as well,

Etot[Φ0]↔Etot[n0]≡Etot

0,(2.18)

and can be determined by minimizing (2.17).

The energy functional

Etot[n(r)] = F[n(r)] + Zn(r)υext(r)dr(2.19)

contains, however, an unknown (yet universal, i.e., not depending on υext though) functional

F[n(r)] = min

Φtrial→n(r)hΦtrial|(Te+Vee)|Φtriali,(2.20)

including Teand Vee given by the electronic subsystem Eq. (2.6).

The Hohenberg-Kohn theorems provide a major simpliﬁcation by transferring the problem of

dealing with the many-body wave function to the much simpler ground-state density. A (practical)

route to treat the functional F[n]was proposed by Kohn and Sham in 1965.

Kohn-Sham Equations

In order to determine the functional F[n], Kohn and Sham introduced a ﬁctitious auxiliary system

of non-interacting electrons. The “trick” is to map the complex problem of interacting electrons in

a given external potential υext onto a system of non-interacting electrons in an eﬀective, so called

Kohn-Sham potential VKS [29]. The potential VKS of the auxiliary system is constructed in such

a way that the electronic density of the non-interacting electron system equals the density of the

interacting electronic system, i.e., n≡nKS.4This allows to reformulate the functional F[n]as

F[n] = TKS[n] + J[n] + Exc[n],(2.21)

where the kinetic energy functional TKS[n]of the non-interacting electrons deﬁned by the single-

particle wave functions ϕKS

iα could be expressed in a comparable simple form5

TKS[n]↔TKS[ϕKS

iα ] = −X

α=↑,↓X

2hϕKS

iα |∇2

i|ϕKS

iα i,(2.22)

4There is no exact proof that this mapping is always possible, see e.g. [10].

5Provided that the ˘ϕKS

iα ¯are orthonormalized.

12 2.1 Many-Body Problem and Density Functional Theory

where the sum runs over all Neoccupied eigenstates. The Hartree term J[n](classical electrostatic

electron-electron interaction) in Eq. (2.21) is deﬁned as

J[n(r)] = J[n(r)] = e2

4πǫ0

2Z Z n(r)n(r′)

|r−r′|drdr′.(2.23)

Apart from J[n]and TKS[n], which are explicitly known [see Eqs. (2.22) and (2.23)], the (unknown)

correction including the kinetic energy diﬀerence of the non-interacting and interacting electronic

system in Eq. (2.21) and the nonclassical contribution of the electron-electron repulsion is deﬁned

as the exchange-correlation energy

Exc[n]≡Te[n]−TKS[n] + Vee[n]−J[n].(2.24)

Although the explicit formulation for the functional dependence of Exc[n]is unknown, successful

approximations exist (some are discussed in the next section). On the other hand, we gain (due to

the mapping) a much more convenient situation: The many-body problem in Eq. (2.13) is reduced

to the problem of solving the single-particle Schrödinger equation6

−1

2∇2+Zn(r′)

|r−r′|dr′+VKS(r) ϕKS

i↑(r)

ϕKS

i↓(r)!=ǫi ϕKS

i↑(r)

ϕKS

i↓(r)!.(2.25)

The eﬀective Kohn-Sham potential VKS of the one-particle electron system is given by the 2×2

matrix

VKS

αβ (r) = υext

αβ (r) + Vxc

αβ(r),(2.26)

including the known contributions υext and α, β representing the spins ↑,↓. The unknown exchange-

correlation potential is derived from the corresponding energy, Eq. (2.24), as

Vxc

αβ =δExc[n(r)]

δnαβ(r).(2.27)

Due to the construction of the auxiliary system, the electronic density matrix Eq. (2.16) can be

expressed in terms of the introduced Kohn-Sham single-particle wave functions as

nαβ(r) =

ϕKS

iα (r)ϕ∗KS

iβ (r).(2.28)

Based on the variational principle Eq. (2.17), the Kohn-Sham formalism [Eqs. (2.21)-(2.28)] pro-

vides a practical recipe to compute the ground-state density and the total energy of the electronic

subsystem for a given set of nuclei positions {˜

RI}.

In principle the whole formalism is very general and applicable even to non-collinear spin sys-

tems, i.e., systems where the spin density matrix n(r)in Eq. (2.28) is not diagonal for all rin space.

For these situations, the evaluation of Eq. (2.25) including Eq. (2.27) is not convenient for practical

implementation due to the non-diagonal elements in n. More important, the approximations for

6For the sake of clarity we use atomic units here.

2.1 Many-Body Problem and Density Functional Theory 13

Exc, which enter (2.27), are commonly based on diagonal n. Following Sandratskii et al. [30], the

spin matrix can be locally diagonalized using the unitary transformation

αβ

UiαnαβU+

βj =δijni,(2.29)

where Uis given by the S= 1/2-rotation matrix

U= exp 1

2iθσyexp 1

2iφσz,(2.30)

and θ, φ represent the spherical angles. To express Eq. (2.27) in terms of the (locally) diagonal ni

one uses [31]

∂ni

∂nαβ

=UiαU+

βi.(2.31)

Employing the above relation, Vxc

αβ in Eq. (2.27) can be expressed in terms of the locally diagonal

nias

Vxc

αβ =δExc

δnαβ

i=↑,↓

δExc

δni

δnαβ

i=↑,↓

δExc

δni

UiαU+

βi.(2.32)

The single-particle Schrödinger equation, Eq. (2.25), can be rewritten in terms of the diagonal and

non-diagonal part of the spin density. Therefore it is convenient to express the eﬀective KS potential

VKS in terms of the rotated Pauli-matrix,

σz=UσzU+= cos θ e−iφ sin θ

eiφ sin θ−cos θ!,(2.33)

VKS =V0+∆V ˜

σz(2.34)

with diagonal terms

V0=υext +Zn(r′)

|r−r′|dr′+X

i=↑,↓

δExc

δni

,(2.35)

and mixing terms

∆V =δExc

δn↓−δExc

δn↑.(2.36)

This results in the general expression for the single-particle Schrödinger equation for non-uniform

magnetic systems7

(−1

2∇2+V0)+∆V ˜

σz ϕKS

i↑(r)

ϕKS

i↓(r)!=ǫi ϕKS

i↑(r)

ϕKS

i↓(r)!.(2.37)

7Practically, one usually coarse grain the space rin a discrete mesh containing atomic spheres ν, i.e., θ(r), φ(r)→

θν, φν. This is justified in most cases, where the spin density is almost collinear within the atomic sphere νand is

vanishing outside.

14 2.1 Many-Body Problem and Density Functional Theory

Since V0as well as ∆V0depend only on the (locally) diagonal spin density matrices Eq. (2.29), the

solution is much more straightforward [30] than solving Eqs. (2.25)-(2.27) directly.

Exchange-correlation Functionals

As pointed out above, the only remaining unknown part is the exchange-correlation (xc) functional

Exc, Eq. (2.24), for which several approximations have been proposed.

A well established approximation is the so called local-spin-density approximation (LSDA) which

is based on the homogeneous electron gas model. It is assumed that the exchange-correlation energy

per electron ǫxc at every coordinate ris a function of n↑(r)and n↓(r)at that point only, i.e.,

ǫxc

LSDA ≡ǫxc

LSDA(n↑(r), n↓(r)) which results for the Exc in the LSDA in8

Exc

LSDA[n↑(r), n↓(r)] = Zn(r)ǫxc

LSDA(n↑(r), n↓(r))dr,(2.38)

where n↑,↓represent the spin density of the spin-up/spin-down electrons and n=n↑+n↓.

The xc-energy ǫxc

LSDA can be further split into the exchange and correlation energy:

ǫxc

LSDA =ǫx

LSDA +ǫc

LSDA.(2.39)

For the homogeneous electron gas model the ǫx

LSDA can be derived analytically via Hartree-Fock

theory [20] whereas for the remaining correlation part in ǫc

LSDA eﬃcient parameterizations [32] based

on highly accurate Quantum Monte Carlo calculations of the homogeneous electron gas exist [33].9

Despite the simple character of Eq. (2.38), the LSDA provides (surprisingly) extremely accurate

results.

However, as it turned out, one of the most famous deﬁciencies of the LSDA is that it predicts

a nonmagnetic fcc ground state for iron (which is one of the key materials in the present study)

instead of the experimentally observed ferromagnetic bcc state [34]. This issue was resolved by the

introduction of the generalized gradient approximation (GGA) which accounts in addition for an

explicit dependence of ǫxc on ∇n(r):

Exc

GGA[n(r)] = Zn(r)ǫxc

GGA n↑(r), n↓(r),∇n↑(r),∇n↓(r)dr.(2.40)

Among the diﬀerent proposed GGA parameterizations (PBE, PBE0, rPBE, PW91, etc.), the PBE-

GGA type is one of the most established ones and is employed in the present work.

Although the GGA turned out to improve various properties with respect to the LSDA, Eq. (2.40)

does not allow to quantify the introduced error as well it is not a priori clear if it improves LSDA

in principle. For selected applications in the present study, both functionals were used to estimate

the intrinsic error within the approximations of the xc-functional.

8For the sake of clearness we present here the representation of the LSDA for collinear spin polarization (e.g.

ferro-/anti-ferromagnets). The generalization to non-collinear spin polarization can be achieved by a local rotation

of the coordinate system similar to Eq. (2.29) [see e.g. Ref. [31]].

9The parameterization performed in [32] are performed for the polarized system ζ= (n↑−n↓)/n = 1, the

unpolarized system (ζ= 0) and an interpolation in between as proposed by Barth and Hedin [9].

2.1 Many-Body Problem and Density Functional Theory 15

Bloch Theorem, Periodic Boundary Conditions, and Plane Wave Approach

The DFT formalism discussed above maps the problem of solving a system of interacting electrons

onto an equivalent system of non-interacting electrons in an eﬀective Kohn-Sham potential. Nev-

ertheless, the amount of the eﬀective KS wave functions ϕKS for the description of solids (>1023

per cubic cm) is much too high to be computationally handled. A drastic reduction of the required

wave functions is obtained by taking the crystal symmetries into account. According to Bloch’s

theorem, the periodicity of a (external) potential Vext(r+R) = Vext(r)with respect to a lattice

translation Rn=Piniaiin terms of multiple of the lattice unit vectors aiis reﬂected in the wave

functions of the electrons (also called Bloch functions):

ϕKS

ikα(r) = eik·ruikα(r).(2.41)

Here idenotes the band index, ka vector in the ﬁrst Brillouin zone ΩBZ and uikα(r) = uikα(r+R)

the lattice periodic part of the wave function. Since the single-particle Hamiltonian, the wave

functions and hence all observables inhibit the translational invariance, solving the Schrödinger

equation in the unit cell spanned by the aiis suﬃcient to describe the whole translational periodic

system.

Formulated diﬀerently, the problem of solving an inﬁnite number of electrons occupying a ﬁnite

number of states in an inﬁnite crystal is transformed into the problem of solving a ﬁnite number of

electrons in a ﬁnite cell, but occupying an inﬁnite number of states (k-points). This is reﬂected in

the computation of the density

n(r) = X

α=↑,↓

occ

i|ϕKS

iα (r)|2=1

ΩBZ X

α=↑,↓

occ

iZBZ |ϕKS

ikα(r)|2dk.(2.42)

In practice, the integration over the inﬁnite number of k-points in ΩBZ is replaced by a discrete

grid, and typically crystal symmetries are employed to reduce the computational eﬀort. The explicit

chosen number of k-points for the numerical integration in Eq. (2.42) has to ensure the desired

convergence criteria and depends on the system under investigation.10

Although it is in principle possible to solve the unknown part of the wave function uikα(r)

directly numerically, it is more convenient to expand the latter ﬁrst in an appropriate basis set.

One of the most common approaches in solid state physics is the so called plane-wave approach. It

is based on a plane-wave expansion of uikα(r)as

uikα(r) = X

ci,k+G,α exp[iG·r],(2.43)

where Gdenotes the reciprocal lattice vectors. The expansion coeﬃcients ci,k+G,α can be calculated

by combining Eqs. (2.37), (2.41), and (2.43) as

G′−1

2|k+G|2+V0(G,G′)δGG′+∆V (G,G′)˜

σz ci,k+G,↑

ci,k+G,↓!=ǫik ci,k+G,↑

ci,k+G,↓!.(2.44)

10For each material system under investigation, the explicit chosen parameters are given in the corresponding

section.

16 2.1 Many-Body Problem and Density Functional Theory

The coeﬃcients with small kinetic energy ≤1

2|k+G|2are typically more important than those with

large kinetic energy [35]. In practice, the basis set Eq. (2.43) is therefore truncated at a certain

cutoﬀ value 1

2|k+G|2≤Ecut. The increase of Ecut yields to a consistently better description of the

ground-state density. Therewith Ecut represents a very convenient convergence parameter to control

the basis set and the error introduced due to the truncation within the plane wave approach. The

controllability of the convergence parameters (i.e., in particular the cutoﬀ energy) is one of the big

advantages of the plane-wave basis set.

Summary

Let us brieﬂy summarize this section. For a feasible solution of the many-body problem for the bulk

system, Eq. (2.4), two major steps are required. First, the decoupling of the nuclei and magnetic

degrees of freedom from the electronic motion. In order to compute the ground-state energy for a

given (ﬁxed) set of nuclei {˜

RI}and/or a given magnetic conﬁguration (direction of spin density ˜

m),

the electronic Hamiltonian ˜

Hein an eﬀective potential, Eq. (2.13), has to be solved. This is possible

by the introduction of the DFT reducing the problem to a minimization procedure for the ground-

state spin density n0. The mapping onto an eﬀective, auxiliary system of non-interacting electrons,

the Kohn-Sham approach, reveals a practical recipe to perform this energy minimization, yielding

n0and the ground-state energy. The only remaining approximation is given by the treatment of

the (unknown) exchange-correlation functional Exc.

Ultimately, we are interested in a complete description of the Helmholtz free energy (2.3). In

the next sections we discuss how the DFT computed properties for given constraints can be used

as input data for ﬁnite temperature concepts.

2.2 Extension to Finite Temperatures 17

2.2 Extension to Finite Temperatures

The methodology discussed so far provides an eﬃcient ﬁrst-principles scheme to compute ground-

state properties at zero temperature. For the description of finite temperature properties, two chal-

lenges have to be accomplished:

(a) The description of ﬁnite temperature excitations within the Kohn-Sham DFT formalism.

(b) The contribution of these excitations to the thermodynamic properties.

Both tasks (a) and (b) represent a wide ﬁeld of research on their own. We focus in the following

in particular on the magnetic excitations. In the beginning, a brief overview of the diﬀerent types

of magnetic excitations (a) is given and several techniques to capture these excitations within the

framework of DFT are introduced. For this purpose we discuss in particular the so called adiabatic

spin dynamics.

The treatment of thermodynamic properties (b) is discussed afterwards. For this purpose the

Heisenberg model is introduced and discussed in detail followed by diﬀerent solution strategies.

Finally, the concepts for the determination of vibronic and electronic excitations are shortly sum-

marized.

2.2.1 Magnetic Excitations from DFT

The dominant magnetic excitation in ferromagnets are the so called (single-particle) Stoner ex-

citations and the (collective) spin-wave excitations (see also Fig. 2.2). At low temperatures the

Figure 2.2: Sketch of Stoner continuum (shaded region) and spin-wave excitation spectrum (red

line) for a strong ferromagnet.

spin-waves dominate the magnetic contribution to the thermodynamic properties. They are for

instance responsible for the famous T3/2Bloch law for the magnetization in ferromagnets as well

as for their T3/2contribution to the speciﬁc heat capacity at low T.

Going to higher temperatures, in addition single-particle spin-ﬂip processes take place which also

contribute to the thermodynamics of the system. These Stoner excitations are electron transitions

between bands of opposite spin. For q→0the energy of such excitation corresponds to the

exchange splitting energy (e.g., ∆ex = 1.5,1.1, and 0.3eV for Fe, Co, and Ni, respectively [36]). At

18 2.2 Extension to Finite Temperatures

high temperatures, the crossover of the spin-wave spectrum and Stoner continuum results in the

damping of spin-waves.

A rigorous way to describe both types of excitations starting from the microscopic many-body

problem is the determination of the dynamical spin susceptibility,

χαβ

ij (ri, ti;rj, tj) = gµB

∂mα(ri, ti)

∂Bβ(rj, tj),(2.45)

which describes the (time-dependent) response of the magnetization

mα(r, t) = gµBhσαi(r, t)(2.46)

under variation of a magnetic ﬁeld Br(t). In the above equation, α, β = (x, y, z, +,−)denote

the index of the Pauli matrices including spin annihilation “−“and creation “ + “, and ri,rj

represent the real space coordinates at times ti, tj, respectively. The collective excitations (including

the interference with the Stoner continuum) can be readily obtained by the poles of the Fourier

transformed transversal spin susceptibility χ+−(q, ω).

A direct determination of the time-dependent response function Eq. (2.45) is, however, not

possible within the framework of the time-independent DFT formalism introduced above. Two

extensions to DFT turned out to allow a very accurate description of the dynamical susceptibility,

namely the linear-response DFT [37] based on the time-dependent DFT (TDDFT) formalism [38]

as well as the application of many-body perturbation theory (MBPT) (see e.g. [39]). Both methods

are, compared to the pure DFT formalism, conceptually and computationally much more complex

and expensive, respectively. We brieﬂy expose in the following one of these concepts, namely the

MBPT, to illustrate how χαβ

ij can be in principle addressed within the DFT framework. For a

detailed review we refer to [40]. Starting point is the local spin density in Eq. (2.46),

hσα

rm(tm)i=−iX

σα

ijGji(rm, tm;rm, t+

m),(2.47)

where t+

mindicates an inﬁnitesimal increase of the time variable (to ensure proper time ordering).

The single-particle Green’s function is given by

Gij(r1, t1;r2, t2) = hhai(r1, t1); a+

j(r2, t2)ii ≡ Gij(1,2),(2.48)

where a±

i(rm, tm)denotes the creation (“+”) or annihilation (“−”) of an electron with spin iat time

tmand space rm. Using the short-hand notation introduced above, i.e., 1 = (r1, t1), the single-

particle Green’s function of the full system Gij(1,2) is determined by the Dyson equation (see e.g.

[41]) as

Gij(1,2) = G0

i(1,2)δij +X

kZ Z G0

i(1,3)Σik(3,4)Gkj(4,2)d3d4,(2.49)

where G0

idenotes the single-particle Green’s function of the free (non-interacting) system and Σ

the non-local and energy-dependent self-energy. The Kohn-Sham formalism directly provides via

construction as a byproduct the eigenenergies ǫiand eigenfunctions φKS

iof the auxiliary (non-

2.2 Extension to Finite Temperatures 19

0 0.1 0.2 0.3

ω (eV)

Im χ+−(q,ω) (arbitrary units)

x=0.05

x=0.1

q=(x,x,x)

x=0.2

collective spin wave

excitation

Figure 2.3: Transversal spin susceptibility χ+−(q, ω)of fcc Ni for diﬀerent q-vectors along the [1 1 1]-

direction. For small q-vectors the collective excitations are clearly dominating (sharp peaks). Due

to the interference with the Stoner continuum (see also sketch in Fig. 2.2) the spin-wave excitation

energies are smeared out at higher q. The data are taken from [39].

interacting) Kohn-Sham system. Although they are strictly speaking no physical observables, as a

ﬁrst approximation a non-interacting DFT-KS Green’s function GKS can be deﬁned. Therewith the

single-particle excitations can directly be studied. To determine the collective excitations, access

to the Green’s function Gij of the full, interacting system is, however, required. One of the key

ingredients within a many-body-perturbation approach is therefore the self-energy Σin Eq.(2.49).

An established approach within the framework of DFT is the so called GW approximation [42]

Σij(1,2) ≈iGij(1,2)W(1,2),(2.50)

where Wdenotes the screened interaction. First-principles GW approaches have been successfully

carried out since thirty years [40]. The second ingredient within a MBPT is the treatment of the

electron-hole interaction which is well described by the Bethe-Salpeter equation via a functional

derivative of Σ[40]. As an example we show in Fig. 2.3 the computed χ+−(q, ω)for fcc Ni and

diﬀerent qvectors. For small q, the collective excitations are clearly pronounced by the sharp

peak in χ+−(q, ω). For higher qand higher energies, the interference with the Stoner continuum

(see Fig. 2.2) results in a broadening of the excitation spectrum emphasizing the importance of

single-spin ﬂip excitations in fcc Ni.

Considering the elementary ferromagnets Fe, Co, and Ni, very recent MBPT results by Şaşıoğlu

et al. [39] reveal that the single-particle excitations in Fe and Co are in a large energy window

negligible compared to the collective spin-wave excitations. These ﬁndings are in agreement with

employment of the TDDFT by Savrasov in Ref. [37].

On the other hand, Stoner excitations are found to be much more signiﬁcant in case of fcc Ni

due to the much lower band-splitting [37, 39]. Focusing on the collective spin-wave excitations in

Fe and Co, it is found that both approaches (TDDFT and MBPT) [37, 39] nicely agree with the so

called adiabatic approximation [12, 43] which we discuss next.

20 2.2 Extension to Finite Temperatures

Adiabatic Spin Dynamics from DFT

In the last decade most of the theoretical studies of magnetic excitations were based on an adiabatic

treatment of the spin degree of freedom. The much slower motion of magnetic moments is separated

from the faster electronic motion by mapping the many-body problem onto an eﬀective Heisenberg

model. In such an ansatz, the single-particle spin-ﬂip excitations are neglected.

The derivation provided in the following is based on the adiabatic decoupling as given in [12,

43, 44]. We start by substituting ∆V ˜

σzwith σBeff in Eq. (2.37),

(−1

2∇2+V0)+σBeff  ϕKS

i↑(r)

ϕKS

i↓(r)!=ǫi ϕKS

i↑(r)

ϕKS

i↓(r)!,(2.51)

which emphasizes an important point: The term σBeff indicates that the spin-dependent part of

the exchange-potential11 (i.e. ∆V ) acts like an effective magnetic ﬁeld which is given by

Beff

r=Beff

r





cos φrsin θr

sin φrsin θr

cos θr





,(2.52)

with Beff taking the role of ∆V i.e.

Beff

r:= ∆Vr.(2.53)

The spin-waves are caused by the time-dependent motion of the magnetization mr(t)at each

point in space ras function of time t. To solve the equation of motion for the time-dependent

magnetization we have in principle to solve the full time-dependent Schrödinger equation to obtain

nr(t)and hence mr(t) = Tr (σnr(t)). A full time-dependent DFT would render the current study

to be infeasible.

Since the motion of mr(t)(spin-waves ∼meV) is much slower than the electronic motion (∼eV)

it is justiﬁed to use an adiabatic approach similar to the decoupling of the nuclei and electronic

motion (Born-Oppenheimer approximation in Sec. 2.1.1). This implies that the electrons follow the

time-dependent change of mr(t)instantaneously.

The main assumption is that the time-dependent Schrödinger equation holds also in the frame-

work of DFT in the adiabatic approximation [38], i.e.,

i∂

∂t ϕKS

i↑

ϕKS

i↓!=(−1

2∇2+V0)+σBeff ϕKS

i↑

ϕKS

i↓!.(2.54)

The equation above describes the time evolution of the spinor components ϕKS

i↑and ϕKS

i↓. Solving

this equation yields the spin density nr(t)and hence the time dependence of the magnetization

density via ∂

∂tmr= Trσ∂

∂tnr.(2.55)

11We remind that V0includes the spin-independent part of the exchange-correlation potential, see Eq. (2.35).

2.2 Extension to Finite Temperatures 21

Using Eq. (2.54) and some algebra (see e.g. [31]) one gets a rather simple form12,

∂

∂tmr=−2mr×Beff

r,(2.56)

for the motion of the magnetization density.

Above equation could be identiﬁed with the famous (classical) Landau-Lifshitz Equation which

was introduced by Landau and Lifshitz already in 1935 to model the precessional motion of the

magnetization under an (eﬀective) magnetic ﬁeld Beff [45].13 In our case, the eﬀective ﬁeld Beff is

caused by the spin-dependent part of the exchange-correlation potential, Eqs. (2.52) and (2.53).

If we assume the z-axis as a global quantization axis θr≡θ(which is in absence of magnetic

anisotropy possible), an ansatz fulﬁlling the diﬀerential equation (2.56) is given by

mr(t) = mr





cos φr(t) sin θ

sin φr(t) sin θ

cos θ



,(2.57)

where θdenotes the time- and location-independent polar (cone) angle and φr(t)the time-dependent

azimuthal angle at every point in the crystall r.

In principle, the magnetization direction mr(t)/|mr(t)|as well as the eﬀective ﬁeld Beff

rmay

vary in every point in space r. However, for most practical calculations mr(t)is in a good ap-

proximation collinear within the atomic spheres Ωνof each atom ν[31]. One can therefore proceed

by coarse graining mr(t)and introducing an effective point magnetic moment Mνlocated at the

position Rνof every atom ν,

Mν=ZΩν

mrdr.(2.58)

The above mapping of the magnetization onto eﬀective local (point) moments allows one to deﬁne

the spin quantum number via

|Mν|=Mν=gµBSν(2.59)

with the Landé factor g≈2and the Bohr magneton µB.

By taking the lattice symmetry into account we can deﬁne φr(t)in Eq. (2.57) as

φr(t)→φν(t) = ωt +qRν.(2.60)

The so called spin spiral characterized by Eqs. (2.57)-(2.60) is illustrated in Fig. 2.4. All moments

Mνrotate with a ﬁxed magnitude (projection) around the z-axis (determined by a ﬁxed cone angle

θ). At each time t, the azimuthal angle φν(t)at atom νis determined by a wave vector q, Eq. (2.60).

Since the situation sketched in Fig. 2.4 appears like a snapshot of a frozen spin-wave, this approach

is sometimes also referred to as frozen magnon approach.

The total energies of the such constructed spin spirals can be associated with the spin-wave

energy. To identify the total energy calculations with the spin-wave energy we assume a ferromag-

12Terms which can be related to the change of the orbital contribution of the magnetization are neglected since we

are interested in the spin part only [31].

13Actually, Landau and Lifshitz treated a more general case of Eq. (2.56) including a damping factor.

22 2.2 Extension to Finite Temperatures

qRj

RkRl

Figure 2.4: Sketch of a spin spiral Eq. (2.57).

netic system, which is perturbed by a change of the direction of the magnetic moments Mν. The

resulting energy change ∆E can be expressed as

∆E =X

αβ

JCL

αβ MαMβ(2.61)

which is exactly the classical solution of the Heisenberg model (the model will be discussed in detail

afterwards). The perturbation ∆E gives rise to an eﬀective ﬁeld Bνin every atomic sphere ν,

Bν=−∆E

∂Mν

=−2X

JCL

αβ Mβ.(2.62)

Assuming a small perturbation (θ→0) and combining Eqs. (2.56) and (2.57)-(2.62), one ends up

with

ωDFT

q= 4 X

JCL

αβ Mβeiq(Rα−Rβ)−1=4

M(˜

Jq−˜

J0),(2.63)

where we assumed in the last step a homogeneous magnetization M≡Mβand absorbed Min the

eﬀective exchange coeﬃcients JCLM2=˜

J.14

For ﬁnite cone angles θ,∆E is given as

∆E(q, θ) = sin2θ(˜

Jq−˜

J0).(2.64)

Thus the dispersion ωDFT

qEq. (2.63) reads as [12, 31]

ωDFT

q=4

Mlim

θ→0

∆E(q, θ)

sin2θ.(2.65)

This is one of the key expressions which links the total ground state energies of diﬀerent spin spirals

(obtained employing the DFT formalism) with the magnetic spin-wave excitations. Furthermore,

Eqs. (2.64)-(2.65) allow one to extract directly eﬀective exchange coeﬃcients ˜

J, which could be

further used to derive other (thermodynamic) magnetic properties.

Before concluding this section we note one important point for the practical evaluation of

Eq. (2.65). This concerns the choice of the cone angle θ. From a physical point of view, it should

14For the sake of simplicity the derivation is shown for single-species systems. The extension to multi-species

systems is, however, straightforward, see e.g. [12, 31].

2.2 Extension to Finite Temperatures 23

be as small as possible since spin-waves are small perturbations from the ground state. Because for

small cone angles the total energy ∆E goes to zero as sin2θ[Eq. (2.64)], calculations are carried

out at a ﬁnite θ, which provides suﬃcient large energy diﬀerences ∆E within a reasonable choice of

technical convergence parameters (such as cutoﬀ energy, k-point grid, etc). Apart from this rather

technical aspect, the “best” choice of the cone angle for some applications is still debated [46, 47].

For the determination of high temperature properties such as TC, some works suggest [47] to employ

larger cone angles since they may better mimic the situation at high T(paramagnetic regime). We

will discuss this point in more detail for a realistic material system (namely bcc iron) and also

the impact of the choice of θon TC. However, we would like to note already here that for strong

ferromagnets (large exchange splitting) such as bcc iron, the actual choice of θis rather uncritical.

This has been shown in Refs. [12, 48] and very recently by Kübler [44]. In fact, already the use of

planar spin spirals (the limit case θ=π/2), i.e.

ωDFT

q=4

M∆E(q, π),(2.66)

allows reliable predictions [43, 44].

Another point concerns the computation of incommensurable (wrt. the unit cell) spin spiral

structures. In this case, the “conventional” translational invariance (Bloch theorem) is broken for

translations non-orthogonal to q. However, the atoms separated by a lattice translation are still

equivalent apart from the rotated local moments (if we neglect the small lattice anisotropy due to

spin orbit coupling). Transformations combining a translation (by Rn)and a rotation (by q·Rn)

leave the spiral structure invariant and provide a very elegant way to treat spin spirals. Following

Sandratskii [49], this can formally be expressed by a generalized translation Tnwhich includes the

translation Rnand the spin rotation U[Eq. (2.30)], and commutes with the Hamiltonian:

TnH(r)ϕ(r) = H(r)U(−qRn)ϕ(r+Rn).(2.67)

In analogy to Bloch’s Theorem the eigenstates can be chosen as

Tnϕik(r) = U(−qRn)ϕik(r+Rn) = eik·Rnϕik(r).(2.68)

The lattice periodic part of the eigenstates above follows the chemical lattice (via Rn). This

implies that spin spirals can be computed in the chemical unit cell, which drastically reduces the

computational eﬀort, or even makes the computation of incommensurable spiral structures possible

in ﬁrst place. Further reduction of the computational eﬀort is possible if the perturbation is small

(small qvectors and cone angles θ). In this limit, the magnetic force theorem [15] can be applied

allowing to reduce the computational costs as well.

The eﬀective spin Sν(2.59) and the eﬀective exchange coeﬃcients ˜

J(2.63)-(2.65) completely

determine the Heisenberg model which we discuss next. Therefore, we have all ingredients in hands

to study the thermodynamic properties of magnetic materials based on ﬁrst-principles results.

24 2.2 Extension to Finite Temperatures

Figure 2.5: Sketch of the Heisenberg model Eq. (2.70). Exchange interaction JQM

ij between spins

Siand Sjlocalized at lattice sites Riand Rj.

2.2.2 The Heisenberg Model

The frozen magnon approach based on the adiabatic spin dynamics discussed in the previous section

allows us to map the complex many-body problem onto a more simple, effective magnetic Hamil-

tonian. Explicitly, we have introduced in Eq. (2.61) the classical formulation of the Heisenberg

model.15

Originally, the model was proposed to study the magnetism caused by localized moments (e.g.

due to partially ﬁlled electron shells, such as e.g. EuO or EuS). The quantum-mechanical Hamilto-

nian is given by

H=−X

JQM

ij SiSj,(2.70)

with exchange coeﬃcients JQM

ij , which mediate the magnetic exchange between spins Silocalized

at lattice sites iand j. The Sidenote the spin operators obeying the usual commutator rules for

orbital moments and having discretized eigenvalues. The Heisenberg model should be interpreted

as an effective operator simulating the Coulomb interaction via the spin-spin interaction Si·Sj.

Although originally designed to describe the situation of ideal localized moments, the model has

been successfully16 applied to magnetic metals showing strongly (but not ideally) localized magnetic

moments (such as Fe and Co) see e.g. [15, 16, 43, 48, 54].

15For the sake of completeness we mention here also two other famous magnetic model Hamiltonians. The first one

is the so called Hubbard model which is supposed to describe in particular itinerant magnetic materials. These are

systems where both phenomena, conductivity and magnetism, are caused by the same group of (itinerant) electrons.

It is therefore in particular well suited to account for the single-particle Stoner excitations. The Hubbard model was

proposed in 1963 by Hubbard [50], Gutzwiller [51] and Kanamori [52]:

H=X

ijσ

Tij a+

iσajσ +UHub

iσ

niσni−σ.(2.69)

Here, a+

iσ (aiσ) denotes the creation-(annihilation-) operator of an electron having spin σ=↑,↓.Tij is the hopping

integral between lattice sites Riand Rj. The operator niσ =a+

iσaiσ denotes the occupation number operator and

UHub (also called Hubbard-U) represents a simplified version of the Coulomb interaction energy between the electrons,

which is reduced to a local, intra-atomic interaction. The lattice information is effectively included in the Tij . The

Hubbard model allows hence to study the interplay between kinetic energy, Coulomb interaction and lattice structure

under the Pauli principle.

Another famous magnetic model was proposed for systems where two different groups of electrons are responsible

for the magnetism and the conductivity. A prototype is the 4f-system Gd. Here, the itinerant electrons can interact

with the localized spin system. The appropriate model is called Kondo-lattice model [53].

16We note that this holds only true if one is interested in the magnetic properties of these materials and not in the

interplay of, e.g., conductivity and magnetism.

2.2 Extension to Finite Temperatures 25

The subscript “QM” underlies that the above representation is a quantum model. In fact, we

will show in Sec. 3.2.2 explicitly that the quantum eﬀects, neglected by the classical formulation

(2.61), are mandatory for an accurate description of thermodynamic properties. The solution of the

quantum model is, however, not that straightforward.

Before discussing diﬀerent solution techniques we use the opportunity to deepen the concept

of spin-waves. Up to now they were phenomenologically introduced as small perturbations (“spin

spirals”) from the ferromagnetic ground state. This results ﬁnally in the expression (2.65), which

combines the total energy calculations ∆E(q, θ)with the excitation energy ωDFT

q. We provide in

the following a more formal approach to the spin-wave concept.

Spin-waves

To show that spin-waves are the natural eigensolutions of the Heisenberg model it is convenient to

express Eq. (2.70) in terms of raising (lowering) operators S±

i=Sx

i±iSy

ias

Si·Sj=Sx

iSx

j+Sy

iSy

j+Sz

iSz

j=1

2(S+

iS−

j+S−

iS+

j) + Sz

iSz

j,(2.71)

where the z-axis is chosen as the quantization axis.17

After Fourier transformation,

Sσ

q=X

e−iqRiSσ

i,(σ= +,−, z)(2.72)

JQM

q=1

JQM

ij eiq(Ri−Rj),(2.73)

the Hamiltonian can be expressed in q-space as

H=−1

JQM

qS+

qS−

−q+Sz

qSz

−q−1

~gµBBz

0Sz

0.(2.74)

Denoting the ground state (ferromagnetic solution) as |0ione can show (see e.g. [56]) that

|qi=1

~√2SN S−

q|0i(2.75)

is a (normalized) eigenstate of (2.74) with an eigenenergy given by

Eq=E0+Eq,(2.76)

where

E0=−NJ0~2S2−NgµBBz

0(2.77)

17In high symmetric 3D systems the magnetic lattice anisotropy due to spin-orbit coupling as well as the dipolar

anisotropy is at least one order of magnitude smaller than the exchange interaction Jand is therefore neglected in

the present study. We note, however, that these anisotropies are crucial in low dimension, in particular to realize

finite magnetization in 2D films (Mermin Wagner theorem [55]).

26 2.2 Extension to Finite Temperatures

denotes the ground-state energy.

The energy

Eq=gµBBz

0+ 2S~2JQM

0−JQM

q(2.78)

is interpreted as the excitation energy of one magnon. From the Zeeman term one directly recognizes

that the total moment of the system is reduced by gµB. The quasiparticle has therefore the spin 1,

i.e. magnons are bosons.

The expectation value of the lattice site dependent magnetization,

hq|Sz

i|qi=hSz

ii=~S−1

N,(2.79)

reveals an astonishing outcome: the spin deviation is lattice site independent! This directly implies

the concept of a spin-wave: The spin deviation Eq. (2.79) is uniformly distributed over the whole

lattice and characterized by a wave vector q. Expressed diﬀerently, the local spins at each lattice

site Siprocess around the z-axis with a fixed projection onto the z-axis. In the thermodynamic

limit N→ ∞, the cone angle θ, determining the deviation from the z-axis in Eq. (2.65), goes to

zero.

Since at T= 0 K the classical and quantum solution of the ferromagnetic Heisenberg model are

identical, one can formally derive Eq. (2.65) by postulating JQM S2≡˜

Jand S≡2Mi.e. g≈2in

Eq. (2.59).

2.2.3 Thermodynamic Properties

Having introduced the magnetic model Hamiltonian (2.70) and underlying DFT methodology to

extract the required model parameters, we can proceed in deriving ﬁnite temperature properties.

For this purpose an important quantity is the internal energy U(T, V ) = hHi(T, V )of the system.

It is directly connected to the speciﬁc heat capacity via

CV(T) = dU

dT V

,(2.80)

which is in turn related to the entropy via the integration

SV(T) = ZT

CV(T′)/T′dT′.(2.81)

Combining (2.80) and (2.81), one can derive an expression for the Helmholtz free energy solely18

depending on U(T, V )[56]

Fmag(T, V ) = U(0, V )−TZT′

dT′U(T′, V )−U(0, V )

T′2.(2.82)

However, a numerical or analytical exact solution for hHi(T, V )is in general not feasible. Pro-

viding eﬃcient approaches for hHi(T, V )is therefore one of the key concerns in the present study.

18Eq. (2.82) is, strictly speaking, only valid for systems providing a unique ground state, e.g., a perfect crystal

lattice.

2.2 Extension to Finite Temperatures 27

To prepare this we discuss next diﬀerent solution techniques which are commonly used to determine

thermodynamic properties.

Analytic Approaches

We start by introducing two analytical approaches which are commonly used to determine the Curie

temperature TC.

For an analytical treatment it is convenient to start from the Heisenberg model in the formulation

of raising/lowering operators as in Eq. (2.74), i.e.

H=−X

JQM

ij (S+

iS−

j+Sz

iSz

j).(2.83)

To obtain a solution for TC, we investigate the average magnetization, i.e. hSzi(T), which

vanishes as the temperature approaches the Curie temperature, hSzi(TC)→0.

A rather simple but often successful approach is the mean ﬁeld (MF) or often called molecular

ﬁeld approximation. Essentially, ﬂuctuations of the observables A, B around their mean values

hAi,hBiare neglected in MF, i.e. explicitly

(A−hAi)(B−hBi)→0.(2.84)

Therewith Eq. (2.83) becomes:

HMF =JQM

0hSzi2

MF(T)−BMF(T)X

i.(2.85)

Physically, spin operator products are now decoupled reducing the Hamiltonian to the one of an

ideal paramagnet [56] under an “eﬀective” (mean) ﬁeld

BMF(T) := 2JQM

0hSziMF(T),(2.86)

caused by the exchange interaction JQM

ij between all spins, JQM

0= 1/N Pij JQM

ij according to

(2.73).

Based on Eq. (2.85) it is straightforward to derive the reduced magnetization (see A.1 for a

detailed derivation):

mMF(T) = hSziMF/S =BS 2S2JQM

kBThSziMF

S!,(2.87)

where the Brillouin function BS(x)is deﬁned as

BS(x) := 2S+ 1

2Scoth 2S+ 1

2Sx−1

2Scoth x

2S.(2.88)

The limit T→TCimplies x∼ hSzi → 0in Eq. (2.88) and we can expand the Brillouin function:

BS(x)x→0

−→ S+ 1

3Sx+O(x3).(2.89)

28 2.2 Extension to Finite Temperatures

Combining Eqs. (2.87) and (2.89) we obtain the expression for the mean ﬁeld critical temperature:

kBTMF

C=2

3S(S+ 1)JQM

0=2

3˜

J0.(2.90)

In Eq. (2.90) we have introduced the eﬀective exchange parameter for the quantum model, i.e.

S(S+ 1)JQM =˜

J. (2.91)

One has to note that there is no “strict” mapping between the classical and quantum model [57].

Another choice would be S2JQM =˜

Jsimilar as done for the classical formulation which does,

however, not account for the quantum nature of the spins. In contrast to this, the deﬁnition above

is rooted in the quantum nature of the spin operators ˆ

S, i.e. due to the eigenvalue relation of angular

momentum operators ˆ

S2|Si=S(S+ 1)|Si. Such a deﬁnition of eﬀective exchange coeﬃcients is

useful in particular when results of the classical and quantum solution are compared [57, 58], since

thermodynamic quantities approximately scale in the quantum case as S(S+ 1).19

Obviously, the simple character of the MF approximation reveals some shortcomings. Having

spin conservation at each lattice site in mind, i.e. hSi+i=hSi−i ≡ 0, the spin ﬂip terms in the

Hamiltonian (2.83) were neglected by the mean ﬁeld ansatz

hS+

iS−

jiMF

−→ 0.(2.92)

This means that also the fundamental spin-wave excitations are suppressed in the MF ansatz. In

addition short-range order is fully neglected since in the paramagnetic regime T≥TC(where the

total magnetization goes to zero hSzi → 0), the MF Hamiltonian Eq. (2.85) vanishes, HMF →0.

Although it is well known that the mean ﬁeld solution Eq. (2.90) overestimates TC, it usually

provides (at least qualitatively) a good prediction of some physical trends [17, 60].

A more sophisticated approach, which includes the collective (spin-wave) excitations, is given

by the following Green’s function (GF) based formalism. The main idea is to use appropriate GFs

to access correlation functions which determine the magnetization. We now discuss a ﬁrst order GF

decoupling procedure which is known to provide a reliable prediction of TC, even on a quantitative

level. We start with the following operator identity

hS−

iS+

ii=~2S(S+ 1) + ~hSz

ii−h(Sz

i)2i.(2.93)

In order to compute the magnetization hSz

iion the right hand side of this equation, we need to

determine the correlation function hS−

iS+

iias well as the higher order moment h(Sz

i)2i. Following

the formalism proposed by Callen in 1963 [61], both can be obtained by the following GF

G(a)

ij =hhS+

i;B−(a)

jii := hhS+

i;eaSjz S−

jii.(2.94)

19The critical temperature within the quantum case does not strictly follow TC(S)∼JQMS(S+ 1) as was, e.g.,

recently shown using high temperature series expansion [59]. The complex dependence of TC(S)is, however, beyond

the scope of this work and we use (2.91) instead.

2.2 Extension to Finite Temperatures 29

The above introduced parameter aallows in a very elegant way to determine h(Sz

i)nifor various n.

The equation of motion for G(a)

ij is given by

EG(a)

ij (E) = ~h[S−

i;B+(a)

j]−i+hh[S−

i,H]−;B+(a)

jiiE.(2.95)

A direct solution of G(a)

ij is impractical due to the higher order GFs on the right hand side, which

appear due to the commutator [S−

i,H]−. To overcome this, the higher order GFs are typically

decoupled on a certain level. The Tyablikov decoupling, as proposed by Tyablikov in 1959 [18], or

also called random phase approximation (RPA), is given by the ﬁrst order decoupling of the higher

order Green’s function on the right hand side. Explicitly, the GFs are decoupled such as

hhS+

mSz

i;S−

jiiERP A

−→ hSz

iihhS+

m;S−

jiiE(2.96)

hhS+

iSz

m;S−

jiiERP A

−→ hSz

mihhS+

i;S−

jiiE.(2.97)

Physically, this corresponds to the assumption that the z-component of the magnetization hSz

iiis

not strongly correlated to the surrounding and hence its ﬂuctuations are small. This holds true

at very low temperatures (close to the ground state) as well as at very high temperatures where

the thermal induced randomization dominates over the correlation due to the exchange interaction.

An important outcome is that the results of the RPA agree with non-interacting spin-wave theory

at low temperatures as well as statistical theory at high temperatures (series expansion technique)

[18, 62]. For this reason, the RPA is assumed to provide a reasonable description also between these

limits. Indeed it is found that despite the simple character of the approximation, the RPA procedure

describes the transversal correlation function hS+

iS−

jisatisfactory over the entire temperature range

[63].

Using above decouplings [Eqs. (2.96) and (2.97)] it is possible to solve Eq. (2.95) and to derive

the magnetization (a full derivation is given in A.2) as

hSziRPA =(1 + ϕ)2S+1(S−ϕ) + ϕ2S+1(S+ 1 + ϕ)

(1 + ϕ)2S+1 −ϕ2S+1 ,(2.98)

with the eﬀective magnon occupation number

ϕ=1

eβE(q)−1(2.99)

and eﬀective magnon energies

Eq= 2~hSziRPA(JQM

0−JQM

q).(2.100)

The critical temperature is obtained by the limit of vanishing magnetization, i.e. for hSziRPA →0

in Eq. (2.98) resulting in [18]

kBTRPA

C=2

3S(S+ 1) X

J0−Jq!−1

3 X

J0−˜

Jq!−1

.(2.101)

30 2.2 Extension to Finite Temperatures

Together with the MF result (3.2), the above expression for TCis standardly used in combination

with DFT derived input parameters [13–17, 58].

Numerical Approaches

In principle, there are diﬀerent methodologies to solve the Heisenberg model exactly by using

numerical approaches. Generally, the most straightforward way would be an exact diagonalization

of the Hamiltonian (2.83) in its matrix representation. However, due to the fast increase of the

Hilbert space, which scales as ∼(2S+1)Nwith the number of spins N, these methods are restricted

to small systems (N∼10 −100) and are usually applied to study small magnetic clusters [64]. The

magnetic interaction in bulk systems such as in the 3d ferromagnets Fe, Co, or Ni are typically very

long-ranged (usually not negligible below at least twenty neighboring shells) and show RKKY-type20

oscillating (ferromagnetic/antiferromagnetic) interactions.

The Quantum Monte Carlo (QMC) technique allows one to study much larger system sizes

(N∼104) [65]. For spin systems as the Heisenberg model,21 the standard methods are employing

so called loop algorithms (loop [67] or direct loop [68]) based on diﬀerent representations of the

partition function (e.g. stochastic series expansion (SSE) [69] or path integral representation).22 A

detailed discussion of these methods goes beyond the scope of this work. We brieﬂy discuss the main

idea behind the QMC method and the main inherent limitations and refer for a detailed overview

to Ref. [67].

For calculating the phase space average of a quantity Oin a classical Monte Carlo (cMC)

simulation, the direct evaluation of the sum

hOi=1

c∈Ω

O(c)p(c),(2.102)

where Z=Pc∈Ωp(c)denotes the partition sum, over the high-dimensional phase space Ωis circum-

vented by choosing a subset of conﬁgurations {ci}from Ωaccording to the probability distribution

p(ci). The average is afterwards approximated as the mean value

hOi ≈ O=1

Nsteps

O(ci).(2.103)

Typically in statistical physics, p(c) = e−βE(c)is the Boltzmann weight and E(c)is the energy of

the conﬁguration c.

20RKKY stands for Ruderman-Kittel-Kasuya-Yosida and was originally proposed to describe the indirect exchange

coupling between nuclear spins via the conducting electrons. The ferromagnetic/antiferromagnetic oscillation of the

exchange coupling depending on distance is one prediction of the RKKY theory.

21In presence of very strong magnetic fields, the so called Worm algorithm provides an alternative treatment [66].

22A very recent methodology is the so called Quantum Wang-Landau method (QWL) [70]. It is the analogon of the

classical Wang-Landau technique [71] and was very recently implemented in the ALPS package [72] for the special

case of a nearest-neighbor S= 1/2model. Since our attention is to study also the spin-dependence S≥1/2and the

impact of long-range interactions we do not discuss the method in the present work.

2.2 Extension to Finite Temperatures 31

For a quantum system described by a quantum-mechanical Hamiltonian H, a thermal average

requires the solution of an operator expression

hOi=Tr Oe−βH

Tre−βH.(2.104)

The main problem for solving expressions as the one above is the exponential e−βHin the trace.

A direct calculation requires the complete diagonalization of H. As pointed out above, this is

computationally not feasible for larger systems due to the huge conﬁguration space. Monte Carlo

techniques can again be applied, but require the mapping of the quantum problem onto a classical

one. A common approach is to rewrite the exponent in a Taylor expansion23

ZQM = Tre−βH=

∞

n=0

(−β)n

n!TrHn=

∞

n=0 X

{α}

(−β)n

n!hα|Hn|αi(2.105)

in the basis {α}={|Sz

1, Sz

2,...,Sz

Ni}.

This procedure is called stochastic series expansion (SSE) technique [69]. The effective weights

p(c)are deﬁned by the product of matrix elements Hnand the factor (−β)n/(n!). The reformulation

of Eq. (2.104) in terms of p(c)allows one to employ (classical) Monte Carlo schemes. The main

limitation is, however, that by construction the eﬀective weights p(c)for some conﬁgurations ccan

be in principle negative, i.e. p(c)<0. This is the famous “negative sign problem” in QMC methods

[73, 74]. These situations occur in the Heisenberg model when magnetic frustration is present (e.g.

in the antiferromagnetic fcc lattice [75]). A prototype example to illustrate magnetic frustration is

the antiferromagnetic triangular Ising model

HIsing

Triangular =J(σ1σ2+σ2σ3+σ3σ1)(2.106)

with σ=↑,↓as sketched in Fig. 2.6. Explicitly, negative weights arise due to matrix products such

p(c)∼(−1)3h↑↓↓ |H| ↓↓↑ih↓↓↑ |H| ↓↑↓ih↓↑↓ |H| ↑↓↓i <0(2.107)

where three spin ﬂips occur, yielding a factor (−J)3.

There is up to now no general solution for the sign problem.24 In principle, simulations can still

be performed by taking the absolute value of p(c)

hOi=PcO(c)p(c)

Pcp(c)=PcO(c)sgn(p(c))|p(c)|/Pc|p(c)|

Pcsgn(p(c))|p(c)|/Pc|p(c)|=hOsgni|p|

hsgni|p|

.(2.108)

Unfortunately, the statistical errors increase exponentially with increasing particle number Nand

inverse temperature β. To show this we consider the mean value of the sign hsgni|p|=ZQM/ZQM

|p|in

above equation being the ratio of the partition function ZQM =Pcp(c)and the auxiliary function

23Alternatively the exponent can also be discretized in imaginary time τ(inverse temperature) β=M∆τ and

e−βH= (1 −∆τH)M+O(∆τ). This representation is usually referred to as world line representation.

24In fact one can show that a generic solution of the sign problem would solve all problems in the complexity class

NP (nondeterministic polynomial) in polynomial time [76].

32 2.2 Extension to Finite Temperatures

J<0

Figure 2.6: Sketch for the antiferromagnetic (J < 0) triangular. There is no conﬁguration which

minimizes the energy of all three bonds, causing magnetic frustration.

ZQM

|p|=Pc|p(c)|used for the sampling. Since the partition functions are exponentials of the free

energies, this ratio is exponential proportional to the free energy diﬀerence ∆F

hsgni|p|=ZQM/ZQM

|p|= exp(−βN∆F).(2.109)

Consequently, the relative error ∆sgn/hsgniincreases exponentially with increasing particle number

Nand inverse temperature β:

∆sgn

hsgni=p(hsgn2i−hsgni2)/Nsteps

hsgni=p(1 −hsgni2)

pNstepshsgni∼exp(βN∆F )

pNsteps

.(2.110)

Similar the error in the numerator in (2.108) increases exponentially [76] and hence the computa-

tional time needed to achieve a certain relative error scales exponentially with inverse temperature

βand spins N.

As mentioned in the beginning of this section, in real materials such as Fe, Co or Ni, the interac-

tions are typically long-ranged and show RKKY-type oscillating (ferromagnetic/antiferromagnetic)

interactions. As we will show explicitly in Sec. 3.2.2, these type of interactions cause a negative sign

problem and limit therewith the application of QMC to these materials. Considering the shortcom-

ings of the quantum methods, it is not surprising that for properties of magnetic systems addressed

in the present work (e.g. ab initio derived TC), apart from the analytic methods [13–17], only

classical cMC calculations are commonly employed [13, 43, 48, 54].

2.2 Extension to Finite Temperatures 33

2.2.4 Vibronic and Electronic Free Energy

So far, the magnetic contribution to the free energy is considered. In order to accurately describe the

transversal spin degree of freedom we recourse to the Heisenberg model, where the single-particle

excitations (spin-ﬂips) were neglected. Increasing the temperature does not only induce single-

particle excitations between bands of opposite spins. It goes without saying that due to the Fermi

distribution the electrons also occupy states above the Fermi energy ǫFin the same bands (without

spin-ﬂips). In the framework of ﬁnite temperature DFT this contribution is usually referred to as

electronic contribution. The other main mechanism contribution to the free energy is driven by the

motion of atoms, the so called vibronic contributions. In the following we brieﬂy outline how both

contributions can be eﬀectively captured within the DFT framework.

Electronic Free Energy

For convenience we separate the electronic free energy Fel into the T= 0 K ground-state energy

Etot and the ﬁnite temperature entropy term Sel,

Fel(T, V ) = Etot(T= 0 K, V )−T Sel(T, V ).(2.111)

A very common and accurate parameterization of Etot(T= 0 K, V )is given by the Murnaghan

equation of state which will be used for the further evaluation [77].

To derive Sel(T, V )we employ the ﬁnite temperature formulation of DFT given by Mermin [78].

The key idea to simulate the electronic excitations is to combine the Fermi function25 fFermi(ǫi) =

[1 + exp (βǫi)]−1with the single-particle electron energies ǫi[Eq. (2.25)]. The approach employs the

analytically known (see e.g. [79]) expression for the electronic entropy of an ideal, non-interacting

Fermion gas (e.g., the auxiliary Kohn-Sham particles)

Sel(T, V ) = −kBX

i{[1 −fFermi(T, ǫi(V))] ln[1 −fFermi(T, ǫi(V))]

+fFermi(T, ǫi(V)) ln fFermi(T, ǫi(V))}.(2.112)

Physically, the ﬁrst expression on the right-hand side of above equation originates from the con-

tribution of unoccupied states (holes) whereas the second term originates from the occupied states

(particles). In practice, both contributions were computed on a ﬁnite mesh26 of temperatures T

and volumes V. Eq. (2.112) is widely used to compute DFT based electronic free energies and has

been successfully applied to various materials (see e.g. [80, 81]).

Vibronic Free Energy

A major contribution to the free energy in solids is given by the thermally driven nuclei motion, the

vibronic excitations. Thermal energy induces collective lattice vibrations which can be described by

Bose quasiparticles, called phonons. It turns out that the so called quasiharmonic approximation

25For the ideal Fermi gas (of the electrons) at T= 0 K, the chemical potential µis equal to the Fermi-energy EF.

For convenience the single-particle energies are shifted such that µ≡0.

26The explicit technical parameters ensuring a convergence of less than 1meV per atom is system dependent. The

chosen parameters are therefore explicitly given in the text at the corresponding place.

34 2.2 Extension to Finite Temperatures

(qh) is able to predict the vibronic contribution to a large extent whereas explicit27 anharmonic

contributions (which go beyond the qh) are in magnitude smaller [80, 81].

To derive the qh contribution to Fvib one can employ constrained DFT calculations for a ﬁxed

(given) set of nuclei positions {RI}. Within the adiabatic approximation the nuclei move in an

eﬀective potential given by the Born Oppenheimer surface EBOS, Eq. (2.9). The harmonic approx-

imation is derived by a Taylor expansion of EBOS. Deﬁning the ionic equilibrium positions {R0

and ionic displacements {∆I}we can expand EBOS({R0

I+∆I})for small {∆I}:

EBOS({R0

I+∆I}) = EBOS({R0

I}) + X

Iα

ΦIα∆Iα +1

Iα,Jβ

ΦIα,Jβ∆Iα∆Jβ +O(∆3).(2.113)

Here,

ΦIα =∂EBOS

∂RI,α R0

I≡0,(2.114)

denotes the derivative of EBOS with respect to the displacement of nuclei Iin direction αat R0

which vanishes by deﬁnition (equilibrium position). Analogously, one deﬁnes the so called force

constant matrix as

ΦIα,Jβ =∂2EBOS

∂RIα∂RJβ

=−∂F HF

Iα

∂RJβ

.(2.115)

FHF

Iα represent the Hellmann-Feynman forces which can be very conveniently computed within a

plane-wave approach using the Hellmann-Feynman theorem [82].

Neglecting the higher order terms O(∆3)in Eq. (2.113) is called harmonic approximation.

Eq. (2.113) is therewith reduced to the problem of a harmonic oscillator with energies ωifor which

the known analytic expression for the free energy is given by

Fvib(T, V ) = 1

3Nn

i1

2~ωi+kBTln 1−exp −~ωi

kBT,(2.116)

where the sum runs over the 3Nnphonon states of the Nnatoms in the unit cell.28 The ﬁrst

summand corresponds to the zero-point vibrations and the second represents the entropy of the ideal,

non-interacting Bose-gas (here phonon gas) [79]. There are in principle two, nowadays standard

approaches to obtain the phonon frequencies ωi. Starting point is the construction of the dynamical

matrix DIα,Jβ =1

MΦIα,Jβ which can be done in real space (direct-force method) or in reciprocal

space employing typically perturbation theory (linear-response methods).29 For a detailed overview

we refer to [83]. The eigenvalues of Dprovide eventually the squared frequencies ω2

i. In practice,

crystal symmetries are usually employed to reduce the number of required atom displacements {∆I}

to a minimum set.30

27Anharmonic contributions are (partially) implicitly included in the quasiharmonic approximation due to the

volume dependence of the frequencies. The remaining part, which is not included in the qh approximation, is usually

called explicit anharmonic contribution.

28For convenience we outline here the approach for one-species unit cell. The extension to multi-species systems is

straightforward.

29Both are complementary and should yield in principle to the same results. We restrict ourselves in the present

study, therefore, to the direct-force method.

30In the present work, the SxPhonons add-on of s/phi/nx[84] has been employed to fully automatically determine

2.2 Extension to Finite Temperatures 35

In principle the phonon frequencies can depend on temperature and volume respectively, i.e.

ωi=ωi(T, V ). However, as it turned out, the explicit dependence on temperature is rather small

[81] whereas the major contribution is included in the volume dependence. Therefore, one conven-

tionally assumes ωi≡ωi(V), which is known as the quasiharmonic (qh) approximation. The qh

approximation has been proven to provide reliable results for a large range of metals [80] and will

be also used in this work whenever the vibronic contribution is computed.

Summary

In this section it was shown, how ﬁnite temperature properties can be obtained by ﬁrst principles

by combining statistical concepts with DFT input data. It was discussed how electronic and vi-

bronic contributions to the Helmholtz free energy are commonly calculated based on ﬁrst-principles

methods. To account for magnetic excitations, the concept of adiabatic spin dynamics and the

Heisenberg model were introduced. Theoretical concepts to derive the input data from ﬁrst princi-

ples and diﬀerent ways to solve the magnetic model were presented.

the set of minimum displacements and to compute the frequencies based on the ab initio calculated Hellmann-Feynman

forces.

36 2.2 Extension to Finite Temperatures

Chapter 3

Results

In the previous chapter the magnetic Heisenberg model was derived based on the adiabatic spin

dynamics. We remind that strictly speaking this magnetic model is only well deﬁned for a localized

spin system while single spin-ﬂip transitions are neglected. Despite many successful applications

(see e.g. [15, 16, 48]), employing the Heisenberg model for itinerant magnetic systems is, therefore,

in particular for iron, still critically discussed in literature [48, 85].

In addition, as already indicated in the previous section, there is no numerically exact solution

to the quantum Heisenberg Hamiltonian, Eq. (2.70), with realistic exchange parameters available. A

more detailed discussion in Sec. 3.2.2 shows that the most promising approach for an exact solution,

the QMC, suﬀers from the negative sign problem caused by the long-ranged magnetic interactions.

Therefore, one has to choose between

(a) an approximate solution to the full Hamiltonian, or

(b) an exact solution to a simplified Hamiltonian.

Based on the general concepts explained in the previous chapter, we have developed for both

choices promising and reliable methods. The ultimate aim of these eﬀorts is to describe and predict

thermodynamic properties such as, e.g., the speciﬁc heat capacity or the free energy for magnetic

materials. In order to achieve a full description, the magnetic contribution will be combined with

other excitation processes such as vibronic contributions. To large extend the magnetic free energy

of a metal depends on its magnetic state. Therefore, we ﬁrst focus on the determination of TC.

This also gives us the chance to prove the applicability of our chosen magnetic model Hamiltonian.

3.1 Determination of TC(p)for Bcc Iron:

Performance and Robustness of Different Approaches

We introduce the approximate solution (a) on a question of high relevance for applications: The

pressure dependence of the Curie temperature TCof iron. Apart from the impact for technological

applications, we will later identify TCas the crucial quantity for building up a simpliﬁed model

Hamiltonian ((b) above). This will be discussed in more detail in Sec. 3.4.

In Chapter 1 we pointed out the need of parameter free approaches in particular for supplement-

ing/complementary data within the well known and heavily used calphad (calculation of phase

38 3.1 Determination of TC(p)in bcc iron

diagram) approach. We remind that the calphad methodology allows to compute phase diagrams

of multi-component alloys on the basis of ﬁtting thermodynamic data. Indeed, TCis a crucial quan-

tity within the calphad ansatz since it is one of the few ﬁtting parameters entering the empirical

free energy expression [86, 87]. Since in particular the description of pressure dependent quantities

is still limited in calphad (most of experimental input is given for ambient pressure), predicting

pressure dependent TC(p)is a very important task for the proposed ab initio ansatz.

For one of the key materials in the present study, bcc iron, the experimentally known weak de-

pendence of TC(p)has by now not been understood from ﬁrst principles and a key question remains

still unanswered: What is the physical origin of the experimentally observed negligible change of

TCwith applied pressure (∆TC(p)≈ ±2Kfor pup to 20 kbar [88])?

In a recent work [89] it is claimed that the Heisenberg model itself (which we have chosen for

our analysis) is inapplicable for a meaningful description of TC(p)in iron. This is despite the fact

that the model performs well even for more complex magnetic systems such as, e.g., Fe-Co-alloys

[54], Heusler alloys [90], or diluted magnetic semiconductors [91]. This raises the question: what is

the reason behind this selective failure?

The following section is devoted to these questions. At the same time the analysis oﬀers an

ideal opportunity to test the introduced methods. We would like to address explicitly the following

points:

(1) How accurate is the DFT method including the necessary approximations such as for the xc-

functional (LDA/GGA) to describe the spin-wave spectra (and hence the exchange coeﬃcients)

of magnetic systems?

(2) Is the combined method (DFT based Heisenberg model+analytic approximation) able to

reproduce the experimental ﬁndings for bcc iron, in particular the pressure independence

of the Curie temperature TC(p)≈0?

(3) What is the physical mechanism behind this pressure independence?

In order to make sure that the performance of our ab initio approach is not limited by the

numerical convergence, the total energy diﬀerences are converged within less than 1meV per atom

which corresponds to uncertainties for the calculated TCof less than 10 K. This accuracy limit

has been achieved by using the following parameters: a density of the k-points generated with the

Monkhorst-Pack scheme of ≈5000 per atom, a plane-wave energy cutoﬀ of Ecut = 450 eV, and

≈3000 q-points for the magnetic integration in Eqs. (3.2) and (3.3). All calculations are performed

employing the VASP [92, 93] package and using the projector augmented wave method (PAW) [94]

within LDA and GGA (Perdew-Burke-Ernzerhof parameterization [95]).

We start the analysis with total energy calculations of spin spirals (2.57) for diﬀerent volumes

V, where the energy diﬀerence between the unperturbed (collinear) ground state (q≡0) and the

spin spiral with wave vector qis denoted as ∆EV(q, θ), according to Eq. (2.64). We remind that

the ∆EV(q, θ)are directly connected to the spin-wave energy ωV(q)via Eq. 2.65:

ωq(V) = 4

Mlim

θ→0

∆E(q, θ)

sin2θ.(3.1)

3.1 Determination of TC(p)in bcc iron 39

Since all calculations were performed using fully relaxed non-collinear calculations1, the choice of

the θangle is restricted to the two local minima at θ=π/2and θ= 0. We, therefore, determine

the magnon energies at θ=π/2but we will discuss the impact of the cone angle by comparing

the magnon spectra with previous calculations which employ a perturbation theoretical ansatz for

θ→0as well as by employing a “cutoﬀ-analysis” for ωq.

In order to check the stability of our approach with respect to the inherent approximation of the

DFT formalism (i.e. the treatment of the xc-functional), the magnon spectrum ωDFT

q(V)of bcc iron

is computed for both, the LDA and the GGA-PBE functional. All calculations are carried out for

diﬀerent volumes V∈(11.39Å3,12.45Å3)to obtain eventually the volume dependence TC(V). The

experimental volume at the Curie temperature under ambient pressure is Vexp(TC)≈12.1Å3. The

results for ωDFT

q(V)at two diﬀerent volumes, V= 11.21 Å3and V= 12.45 Å3, are shown in Fig. 3.1

in comparison to experimental and previously published theoretical data from Refs. [15] and [89].

The equilibrium lattice constants at T= 0 K are given by VGGA

0= 11.39 Å3,VLDA

0= 10.31 Å3,

and is extrapolated to Vexp

0= 11.69 Å3from experiment, respectively.

The agreement between the experiment and the calculated frequencies is remarkable for both

LDA and GGA functionals, indicating that the choice of the xc-functional is probably not the

limiting factor. In fact both functionals provide a similar accurate description of the magnetic

excitations in the system. The agreement between the previously and presently calculated magnon

spectra for experimentally accessible magnon energies is also good. For magnon energies higher than

≈0.25 eV there exist qualitative diﬀerences between these results, particular around the H-point.

As was shown by Shallcross et al. in Ref. [46], these diﬀerences can be traced back to diﬀerent spin-

wave cone angles θ(see Eqs. (2.64) and (2.65)) used for the calculation of the magnon spectrum.

Since the diﬀerences for bcc iron are only at relative high magnon energies (0.25 eV corresponds to

≈3000 K) they can be safely neglected. We will show this explicitly by introducing a cutoﬀ-analysis

for the magnon energy.

It should be mentioned that, as discussed in Sec. 2.2.2, spin-waves are only small perturbations

from the ground state. From a principle point of view, a cone angle θ→0seems to provide the most

reasonable description [12]. For the extracted exchange coeﬃcients (which are directly connected to

ω(q)via Eq. (2.64)), such a principle physical argument does not hold and the appropriate choice

of θfor other systems (such as Nickel) at high temperatures is still debated in literature [46, 47]. In

Ref. [47] it is argued that a larger cone angle may better mimic the situation at high temperatures

(paramagnetic regime). In particular in case of bcc iron the choice of θis, however, not crucial.

This is consistent with previous calculations of planar spirals θ=π/2(as used in this work) which

yield satisfactory results for the magnon spectra and the critical temperature [43, 44] of bcc Fe.

Considering the ﬁrst point (1) we raised in the beginning, we can conclude that the DFT

formalism in combination with the adiabatic spin dynamics, the so called frozen magnon approach,

indeed allows a reliable and accurate description of the collective magnetic excitations (spin-waves).

1In principle vasp offers the possibility to add a penalty contribution to the total energy expression in order

to force the direction of the local magnetic moments into a given direction via the flag i_constrained_m. This

contribution is controlled via a parameter lambda.

40 3.1 Determination of TC(p)in bcc iron

We now closer investigate the critical temperature. Combining Eqs. (2.63)-(2.65) with Eqs. (2.90)

and (2.101), the analytic solutions for mean ﬁeld are given as

kBTMF

C(V) = M(V)

ωq(V),(3.2)

and for RPA as

kBTRPA

C(V) = M(V)

6 X

ωq(V)!−1

.(3.3)

To compare the theoretical values for TC(V)with the experimental Texp

C(p)we transformed the

experimentally measured Texp

C(p)to Texp

C(V)employing experimental data for Vexp(p, T =TC).2

In order to elucidate the sensitivity of the Heisenberg model with respect to the quality of the

model input data, we investigate the volume dependence of the Curie temperature employing three

diﬀerent DFT data sets. Speciﬁcally, we compare the qualitative behavior of TC(V)obtained on the

basis of the DFT calculations (i) with the PAW basis set and LDA functional, (ii) with the PAW

basis and the GGA-PBE xc-functional, and (iii) with the tight-binding linear-muﬃn-tin-orbital

(TB-LMTO) basis set within the atomic sphere approximation (ASA) and LDA xc-functional with

the experimental results (see Fig. 3.2). Comparison of the solutions (i) and (ii) allows to ﬁgure

out the eﬀect of the xc-functional, while comparison of solutions (i) and (iii) allows us to identify

possible uncertainties related to the use of diﬀerent DFT basis sets and methodologies for extracting

the magnetic exchange constants (frozen-magnon vs frozen-magnon-torque method [97]).

Figure 3.2 highlights that TC(V)of the Heisenberg model is, apart from a constant shift, very

insensitive with respect to the choice of the xc-functionals for both, the RPA and the MF solution,

when the same basis set is used. Even more noticeable, the use of the same exchange-correlation

functional but diﬀerent basis sets in the DFT calculations might lead to dramatic changes in the

theoretical predictions for TC(V), as can be seen, e.g., on Fig. 3.2(b). The comparison with experi-

ment reveals that all theoretical approaches evaluated in this work, i.e., GGA-PBE/LDA input with

RPA/MF, are in good qualitative agreement with experiment showing a very weak dependence of

the Curie temperature on the volume.

In order to ﬁgure out whether the use of diﬀerent spin-wave cone angles θ(discussed above) might

be responsible for the observed changes in the calculated volume-dependent Curie temperature, we

have calculated TCaccording to Eqs. (3.2)-(3.3) by introducing a cutoﬀ magnon energy ωcut =

0.25 eV 3for the magnon frequencies, i.e. we set ωq≡ωcut for ωq> ωcut. This treatment is

motivated by the pronounced changes of the magnon spectrum at high energies (not accessible by

2In order to derive Vexp(p, T =TC)we have parameterized the experimental dependence of the total energy

E(V, T =TC)on the crystal volume Vemploying the Murnaghan equation of state. [77] The parameters of the

Murnaghan equation at TC(the bulk modulus Bexp

0= 119 GPa, its derivative B′exp

0= 5.3, and the equilibrium

lattice constant Vexp

0= 12.1Å3) are taken from experiment performed by Zhang and Guyot. [96] The minimum

of the enthalpy H(p, V ) = E(V) + pV at the Curie temperature and a given pressure pdetermines the equilibrium

volume Vexp(p, T =TC)and the corresponding lattice constant Vexp(p, T =TC).

3The threshold value ωcut is an energy above which the magnon spectra calculated employing θ=π/2and θ→0

start to significantly (by more than 45 meV) deviate. The value ωcut = 0.25 eV is deduced from the comparison of

our results with the corresponding data of Ref. [15] (see Fig. 3.1). A similar threshold value applies for the different

spin spirals calculated in Ref. [47].

3.1 Determination of TC(p)in bcc iron 41

ΓN P ΓH N

200

400

600

800

1000

ωq (meV)

GGA 11.21 Å³

GGA 12.45 Å³

LDA 11.21 Å³

LDA 12.45 Å³

LDA 11.21 Å³ Ref. [89]

LDA 12.45 Å³ Ref. [89]

GGA 11.82 Å³ Ref. [15]

11,2 11,6 12 12,4

V(Å3)

200

240

280

D (meV A2)

Figure 3.1: Calculated magnon spectra using planar spin spirals for two lattice constants, a= 2.82 Å

(dashed lines) and a= 2.92 Å (solid lines), within GGA-PBE (black lines) and LDA (blue lines).

For comparison the spectra obtained using a small cone angle taken from Ref. [89] (orange lines)

and Ref. [15] (red line) are shown. The spectra taken from Refs. [15] and [89] were calculated using

TB-LMTO employing GGA-PBE and TB-LMTO+ASA within LDA respectively. Experiments:

Ref. [98], measured at 10 K (open squares), and Ref. [99], Fe(12-at.%Si), measured at 300 K (open

circles).

Inset: spin-wave stiﬀness Dobtained by using a least square ﬁt of the low energy branch up to

|q|= 0.4(2π/a0).

experiment) due to the diﬀerent cone angles and allows to estimate the resulting uncertainty of the

calculated Curie temperatures.

The comparison of the Curie temperatures evaluated from the magnon spectrum with and with-

out the cutoﬀ energy are shown in Fig. 3.2. The introduction of the cutoﬀ ωcut results in an

essentially volume-independent shift of the Curie temperature to smaller values over the complete

range of the studied volumes. This already justiﬁes that the choice of the spin-wave cone angle θ

is not essential for a qualitative analysis of the volume dependence of the Curie temperature in bcc

Fe. We notice, however, that quantitatively such shifts diﬀer for the MF and RPA solutions of the

Heisenberg Hamiltonian. The most pronounced change of the calculated Curie temperatures occurs

for the MF solution (Eq. (3.2)), independent on whether LDA or GGA exchange-correlation func-

tionals are used. Within the RPA theory, however, an accurate description of the high-temperature

magnons is not decisive for the evaluation of the Curie temperature since those energies entering the

approximation for the Curie temperature have a lower weight [Eq. (3.3)]. Indeed, the maximum er-

ror introduced by the approximative treatment of the spin spirals for the RPA solution in this work

42 3.1 Determination of TC(p)in bcc iron

11,1 11,4 11,7 12 12,3 400

600

800

1000

1200

1400

1600

11,1 11,4 11,7 12 12,3

400

600

800

1000

1200

1400

1600

GGA

GGA+ωcut

LDA

LDA+ωcut

LDA Ref. [89]

Experiment

V (Å³)

TC (K)

RPA

Figure 3.2: The Curie temperatures within the RPA (left panel) and the MF (right panel) approx-

imation, for both LDA (blue lines) and GGA (black lines) exchange-correlation functionals. For

comparison the results of Ref. [89] (LDA, orange curves), and the experimental Curie tempera-

tures Texp

C(a) for diﬀerent lattice constants a(open green circles) are shown. The transformation

of Texp

C(p)[88] to Texp

C(a)is described in the text.

is less than 10 K for LDA, and 20 K for GGA. The assumption of the planar spin spirals within

RPA only marginally aﬀects the overall accuracy of our calculations. Since the frozen-magnon-

torque method [97] should be comparable with the frozen-magnon-approach used in this work, the

remaining diﬀerences between our results and Ref. [89] are likely related to the diﬀerent basis sets

(plane-wave vs. LMTO basis) used within the DFT calculations.

Since the PAW basis set used in this work allows for consistent convergence with respect to the

quality of the basis set, in contrast to the LMTO basis employed in Ref. [89], we focus on the DFT

calculations obtained with the PAW basis set in the following. As the scatter of the theoretical

predictions for TC(V)is still signiﬁcant, we identify the best theoretical model for the prediction of

the Curie temperature in bcc iron by comparing the theoretical data with the experiment. Regarding

the qualitative behavior we ﬁnd that all combinations (GGA-PBE/LDA input with RPA/MF)

provide a reasonable agreement with experiment by showing a very weak dependence of the Curie

temperature on the volume. The combination GGA-PBE and RPA provides moreover a quantitative

agreement with experiment. The open debate in Ref. [89] if the model is applicable for bcc Fe is

therewith resolved. The diﬀerent results obtained in [89] can be traced back either to the employed

basis set (tight-binding LMTO within ASA) or the frozen-magnon-torque method [97].4

Our analytic solutions Eqs. (3.2) and (3.3) allow us to obtain a deeper physical insight into

the dependence of the critical temperature on pressure. The two physical quantities determining

4In fact, in the description of the method [97] the authors clearly state that their approximation is “reasonable

fast and semiquantitative” and “certainly superior to model descriptions like the Debye model”. Our PAW approach

together with a full self-consistent treatment of the spin-waves goes one step further and allows to systematically

overcome some of the approximations necessary in Ref. [89].

3.1 Determination of TC(p)in bcc iron 43

Table 3.1: Derivatives of the Curie temperature TCwith respect to lattice constant aand pressure

p.∂TC/∂a is obtained by using a linear ﬁt over the entire volume range from a= 2.82 Å up to

a= 2.92 Å. The values in brackets are obtained by the linear ﬁt in a volume range from a= 2.88 Å

to a= 2.9Å, which corresponds to the pressure range applied in experiment (≈20 kbar). The

procedure of transforming TC(V)to TC(p)is described in the text.

GGA LDA

∂TC/∂V ∂TC/∂p ∂TC/∂V ∂TC/∂p

K/10−2Å3(K/kbar) K/10−2Å3(K/kbar) Ref.

MF 1.63(1.25) -0.14(-0.11) 5.35(3.04) -0.60(-0.34) this work

MF+ωcut 1.79(1.58) -0.13(-0.11) 5.26(2.93) -0.48(-0.27) this work

MF+LMTO - - -24.1 1.6 Ref. [89]

RPA -3.81(-3.30) 0.32(0.28) -5.60(-11.51) 0.61(1.26) this work

RPA+ωcut -4.57(-4.79) 0.31(0.33) -6.96(-14.31) 0.62(1.27) this work

RPA+LMTO - - -27.4 1.8 Ref. [89]

Experiment 0.0±0.04 0.0±0.03 0.0±0.04 0.0±0.03 Refs. [88], [96]

its behavior are the spin-wave energies ωqas well as the local magnetic moments M. To analyze

in the following the reason for the weak volume (pressure) dependence of the calculated Curie

temperature we consider ﬁrst the corresponding changes of the magnon spectra. The results are

shown in Fig. 3.1 and we observe that the magnon branches become softer with increasing volume

of the atomic cell, i.e., the energies of the corresponding magnons become smaller. This can be

interpreted as a decrease of the magnetic coupling between the localized magnetic moments. To

quantify this eﬀect we show in the inset of Fig. 3.1 (upper panel) the volume-dependent spin-wave

stiﬀness D, derived in the low wavelength limit ωDFT

q→0≈Dq2. According to Eqs. (3.2) and (3.3) the

decrease of the magnon energies leads to a decrease of the Curie temperature. For both functionals

GGA as well as for LDA the decrease is apparently compensated by another mechanism, since the

Curie temperature remains essentially constant with pressure (see Tab. 3.1).

The reason for the stability of TCversus the applied pressure can be understood recalling that

the increase of the atomic volume does not only weaken the magnon excitation energies, but also

changes the magnitude of the local moments. We show the local magnetic spin moment M(V)

at diﬀerent atomic volumes in Fig. 3.3 for both DFT functionals. At the theoretical ground-state

volumes, i.e. VLDA

0and VGGA

0, LDA underestimates the local spin moment while GGA slightly

overestimates the experimental ﬁndings Mexp(T= 0 K) = 2.14µB[100].5(The small contribution

to the local magnetic moment due to spin orbit coupling (<0.1µB) is neglected.) Both function-

als show a smoothly increasing local magnetic spin moment with increase of the atomic volume,

which is due to the increasing distance between the atoms and a consequent enhancement of the

localization of the 3d-electrons that form the local magnetic moments. According to Eqs. (3.2) and

(3.3) the Curie temperature is directly proportional to the magnitude of the local magnetic moment

M(V), i.e., the increase of the local moments with increase of the volume should lead to a corre-

5At the experimental volume the situation is reversed and MLDA(Vexp)comes closer to the corresponding exper-

imental value.

44 3.1 Determination of TC(p)in bcc iron

11.4 11.7 12 12.3

V (Å³)

2.1

2.2

2.3

2.4

m (µB)

GGA

LDA

Figure 3.3: The local magnetic spin moment of ferromagnetic bcc iron as function of the volume

calculated within GGA (solid orange line) and LDA (dashed blue line) xc-functionals.

sponding increase of the Curie temperature. Therefore, two coupled volume-dependent magnetic

eﬀects determine the behavior of the Curie temperature: The increase of the volume (decrease of

the pressure) leads to a weakening of the spin-wave energy (reﬂecting the strength of the magnetic

coupling) and simultaneously leads to an increase of the local magnetic moments. The sensitive

interplay between these two eﬀects determines the actual dependence of the Curie temperature on

the volume (pressure) for a particular material. All calculated volume and pressure dependencies6

are summarized in Tab. 3.1. The volume dependence is obtained by using a linear ﬁt over the

entire atomic volume range from V= 11.21 Å3up to V= 12.45 Å3. Since the pressure applied in

experiment (≈20 kbar) corresponds to a change of the volume of only ≈0.25 Å3(see Refs. [96]

and [101]), we also ﬁt the volume dependence of TCin a volume range starting from V= 11.94 Å3

to V= 12.19 Å3. The corresponding results are also shown in Tab. 3.1 (shown in brackets). Us-

ing the smaller ﬁt-range the volume and pressure dependencies are only slightly changed, except

for the combination of LDA-RPA where the volume dependence is more pronounced. Regarding

the quantitative agreement between theory and experiment, GGA appears to be closer to exper-

iment with respect to the Curie temperature. The constant shift between LDA and GGA is a

direct consequence of the higher calculated magnetic moment since the calculated magnon spec-

tra within both GGA and LDA are almost similar. In overall the Curie temperatures calculated

with both functionals do not signiﬁcantly depend on the volume (see Tab. 3.1), and we get within

GGA (∂T RPA

C/∂V )≈4.7K/10−2Å3. This corresponds to a change of the Curie temperature of

∆TC≈5K for the pressure range investigated in the experiment [88] (up to 17.5 kbar).

6The pressure dependence of the Curie temperature TC(p)can be obtained from the corresponding volume de-

pendence TC(V)via ∂TC/∂p = (∂TC/∂V )(∂V/∂p)|V=V(TC), where (∂V/∂p)|V=V(TC)=−[V/B0(TC)] is obtained by

calculating7the temperature dependent bulk modulus B0, as described in Refs. [16, 80] in detail.

3.1 Determination of TC(p)in bcc iron 45

Conclusion

Based on DFT input, the critical temperature can be accurately determined by approximative

solutions of the Heisenberg model. This has important consequences: As will be discussed in more

detail in Sec. 3.4, TCcan be identiﬁed as a key quantity for constructing an alternative, simpliﬁed

model. A main advantage of the employed analytic access is that the diﬀerent physical mechanisms

which determine TCcan be discriminated. Both analytic approaches provide a promising way for

the next step which is consequently the extension of the approach for thermodynamic quantities

such as the free energy or the speciﬁc heat capacity.

46 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific

Heat Capacity

A central concern of this work is the study of temperature dependent quantities such as the free

energy or the speciﬁc heat capacity. These quantities are not only determined by magnetic con-

tributions but require also the incorporation of other excitation mechanisms such as vibronic and

electronic excitations.

For this purpose, the Helmholtz free energy surface F(T, V )as a function of the temperature T

and the crystal volume Vis considered in the adiabatic approximation, see Eq. (2.3):

F(T, V ) = Fvib(T, V ) + Fel(T, V ) + Fmag(T, V ),(3.4)

that treats the vibrational Fvib, the electronic Fel (including the energy of the static lattice) and the

magnetic Fmag contribution separately (see also Fig. 2.1). For bcc iron, the adiabatic approximation

is well justiﬁed since the underlying excitation mechanisms, i.e. phonon, electron and magnon exci-

tations, reside on diﬀerent time scales [102]. The validity of the adiabatic approximation is further

supported by recent explicit ab initio calculations of the magnon-phonon coupling in iron [103].

The contributions Fvib and Fel are calculated within the quasiharmonic approximation and ﬁnite

temperature DFT [78], respectively, as discussed in Sec. 2.2.4. As the ﬁrst step, we have therefore

calculated the phonon dispersion of ferromagnetic bcc iron which serves as an input to the partition

function needed for computing Fvib, see Eq. (2.116). The comparison with the experimental phonon

spectrum shows a good agreement (Fig. 3.4). The small diﬀerence between experimental and our

theoretical phonon energies of 5 meV is due to the xc-functional. Using a Born-von Kármán ﬁt to

the experimental as well as to the theoretical phonon dispersions we estimate the diﬀerence in the

free energy to be around 5 meV/atom at 1300 K. Such diﬀerences are consistent with the ﬁndings

in [80] where the uncertainty due to the treatment of the xc-functional on the vibronic free energy

was empirically found for a large range of materials to be in the order of 10 meV up to 1000 K.

The electronic contribution Fel, Eq. (2.111), is obtained by computing the total energies Etot

and the band occupations determined by the Fermi distribution function at a given temperature T,

see Eq. (2.112). The resulting combined vibronic and electronic free energy (Fvib +Fel) is shown

in Fig. 3.5 (orange line) and compared to calphad data. Clearly, the deviation between ab initio

(Fvib +Fel) and calphad data increases rapidly with temperature. For instance, at 1200 K, which

corresponds to the experimental bcc to fcc transition temperature, the diﬀerence is ≈45 meV.

An even more sensitive physical quantity that allows to check the predictive accuracy of the T

dependent ab initio calculations is the heat capacity, CP, which is the second derivative of the

free energy surface with respect to temperature. Figure 3.6 shows that the combined vibronic and

electronic heat capacity agrees well with experimental data only up to ≈300 K, in agreement with

previous results [6]. At higher temperatures, however, Fig. 3.6 clearly reveals that it cannot account

for the rapid increase of the experimental heat capacity. It is therefore obvious that an accurate

determination of the free energy contribution from magnetic excitations is crucial.

For this purpose we ﬁrst employ the analytic approximations, which already showed accurate

results for the critical temperature in the previous section. We discuss the known MF solution for

the magnetic free energy Fmag which we employ in an integrated approach (including Fvib and Fel).

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 47

Phonon energy (meV)

200

400

600

800

1000

1200

Magnon energy (meV)

Figure 3.4: Ab initio phonon (red dashed line, left axis) and magnon (black line, right axis) spectrum

of ferromagnetic bcc iron. The results are respectively compared with available neutron scattering

data [98, 99, 104].

In contrast to MF, no consistent RPA solution for Fmag existed previous to our work for the case

of a general spin system. Consequently, apart from TCpredictions, the RPA has not been applied

up to now for Fmag predictions of realistic materials systems. The extension and application of the

RPA theory is therefore a major subject of the next section.

3.2.1 Analytic Approaches

The key quantity for deriving the magnetic free energy

Fmag(T) = U(0) −TZT′

dT′U(T′)−U(0)

T′2(3.5)

is the internal energy Uof the system. Therefore, the expectation value of the magnetic Hamiltonian

Eq. (2.83), i.e. hHi ≡ U, has to be calculated.

We start with the MF ansatz: HMF is given by Eq. (2.85) from which we can directly derive the

internal energy UMF(T). It is proportional to the square of the magnetization,

UMF(T) = JQM

0hSzi2

MF =S

S+ 1 ˜

J0m2

MF(S, T ),(3.6)

where the reduced MF magnetization mMF is given by Eq. (2.87). Inserting Eq. (3.6) into (2.82)

yields straightforwardly the MF magnetic free energy Fmag

MF .

The exchange coeﬃcients (discussed in the previous section) and the spin quantum number S

which enter (3.6) can be derived employing SDFT. The spin number Sis connected to the local

magnetic moment Mby (2.59). For iron, Mis almost entirely determined by the strongly localized

48 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

200 400 600 800 1000 1200

T(K)

-700

-600

-500

-400

-300

-200

-100

F (meV/atom)

vib+el

vib+el+MF

vib+el+iMF

vib+el+RPA

vib+el+iRPA

CALPHAD

200 400 600 800 1000 1200 1400 1600

-90

-60

-30

120

Figure 3.5: Ab initio free energy Fat zero pressure as a function of temperature Tincluding various

combinations of the vibrational (vib), electronic (el), and magnetic contribution. The latter has been

calculated within the mean ﬁeld (MF), improved mean ﬁeld (iMF), random phase approximation

(RPA), and improved RPA (iRPA). The calphad data has been obtained with the thermocalc

program and the SGTE unary database [105]. Inset: The calphad data is taken as the reference

at each temperature.

3d-electrons and is therefore given by the total magnetization scaled by the number of atoms and a

Landé factor of g≈2. Our T= 0 K PBE result of the total magnetic moment gives8

M≈2.2µB→S≈1.1.(3.7)

Using the computed spin quantum number and the corresponding exchange coeﬃcients, we calculate

the internal energy UMF(T)and obtain Fmag by integration, Eq. (2.82). The MF free energy and

speciﬁc heat capacity are shown in Figs. 3.5 and 3.6 (violett dashed lines), respectively. Up to

≈700 K both quantities are reasonably improved with respect to the pure vibrational and electronic

contributions, but MF fails to reproduce the qualitative behavior of the heat capacity at the Curie

temperature. In addition, as also shown in Fig. 3.2, the calculated TMF

Cis ∼1.5times higher than

the experimental Texp

C= 1044 K.

As discussed in Sec. 3.1, the RPA is known to give more accurate Curie temperatures than

MF. We remind that within RPA we obtained a value much closer to experiment: TRPA

2hPq1/(˜

J0−˜

Jq)i−1/(3kB) = 1082 K (see also Fig. 3.2). However, in contrast to MF there is

no RPA Hamiltonian which can directly be used. Instead the expectation value of the full Heisen-

berg model Eq. (2.83) needs to be calculated in order to derive the expression for the internal

energy URPA =hHi. We need two correlation functions, namely the transversal hSi+Sj−iand

the longitudinal hSizSjzione. Via the corresponding Green’s function Gij, an expression for the

8We remind that at the theoretical ground-state volumes LDA underestimates the local spin moment while GGA

slightly overestimates the experimental findings Mexp(T= 0 K) = 2.14µB

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 49

0 200 400 600 800 1000 1200

T(K)

Heat capacity CP(kB)

vib+el

vib+el+MF

vib+el+iMF

vib+el+RPA

vib+el+iRPA

Experiment

Figure 3.6: Heat capacity CPin units of the Boltzmann constant kBat zero pressure pas a function

of temperature T. The experimental data is from [106, 107].

transversal correlation function hSi+Sj−ihas already been derived (Sec. 2.2.3) and has been proven

[18, 62, 63, 108] to yield a reasonable description for low (limit of spin-wave theory) as well as high

temperatures. However, an expression for hSizSjziis not directly given within the original theory

proposed by Tyablikov [18]. In fact, the few years later derived expressions for hSizSjzibased on

similar decoupling techniques [108] reveals inconsistencies with respect to rotational invariance.9

The underlying reason is that established Green’s function truncation techniques (as e.g. the RPA

decoupling) treat the transversal and longitudinal correlation function on a diﬀerent level of approx-

imation [108]. Therefore, the accuracy of any given approximation can be in general diﬀerent for

the spin-correlation function of the transverse and longitudinal components [108]. Many diﬀerent

schemes were proposed to obtain an improved description for hSizSjzi(e.g. [62, 108–110]) yielding,

however, a huge increase of complexity in the formalism.

However, as we will show in the following, it is still possible to obtain an accurate estimation of

the longitudinal correlation with reasonable eﬀort! In the special case of spin S= 1/2, exact rela-

tions to express the challenging correlation function hSizSjziin terms of hSi+Sj−ican be employed

and a S= 1/2solution for the internal energy can be derived as was, e.g., shown by Tyablikov

[63].10 Therefore the RPA solution for TCis available for S≥1/2, whereas the internal energy

(mandatory for deriving the free energy) is only derived for S= 1/2(see also Fig. 3.7). However,

since we are interested in systems such as bcc iron, a pure S= 1/2theory is of course insuﬃcient

(see Eq. (3.7)). To overcome this issue a rescaling factor in the S= 1/2theory can be introduced

aﬀecting the exchange coeﬃcients JQM [16].

To obtain the rescaling factor we compare the known expressions for TC. Regarding the spin

dependence within the RPA, the spin dependence within TC(S)can be transferred to the exchange

9For example, above TCthe xx,yy and zz correlation functions are not equivalent.

10We present a more general derivation in the Appendix A.3.2.

50 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

RPA

General Spin

(S≥1/2)

Special Spin

S= 1/2

m(T)TC(S)TC(S)m(T)

U(T, S)U(T, S)

Figure 3.7: Sketch of the RPA solution for S= 1/2and general spin S≥1/2. The unknown

S≥1/2solution for the inner energy (grey) is obtained by rescaling the S= 1/2solution for U(T).

The rescaling factor is determined by the known scaling of TCwith S.

parameters as ˜

J=S(S+ 1)JQM, see Eq. (2.101). Employing the TCformula, Eq. (2.101), for

S= 1/2yields:

kBTS=1/2

C=S X

JS=1/2

0−JS=1/2

q!−1

.(3.8)

Comparing Eq. (3.8) with the S≥1/2expression

kBTRPA

C=2

3S(S+ 1) X

JS≥1/2

0−JS≥1/2

q!−1

(3.9)

we ﬁnd the following expression for the eﬀective exchange coeﬃcients

J=S(S+ 1)JQM

S≥1/2=SJQM

S=1/2/γ (3.10)

with the scaling factor γ= 2/3. Using Eq. (3.10) within the S= 1/2solution for the internal

energy [63] we ﬁnally end up with

URPA

mag (T) = Umag

T=0 +M0

2m(T)[m(T) + 1]γh˜

J−˜

Jqi

exp h2βm(T)γ˜

J−˜

Jqi−1

(3.11)

where the magnetization is given by Eq. (2.98) for S= 1/2:

m(T) = (1 + 2

exp[2βm(T)γ(˜

J−˜

Jq)])−1

(3.12)

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 51

0 0,2 0,4 0,6 0,8 1

Reduced temperature T/TC

0,5

1,5

Magnetization M (µB)

improved

RPA

Langevin (classical)

Monte Carlo (classical)

Experiment

Figure 3.8: Comparison of several theoretical approximations for the magnetization M(T) =

M0m(T)in comparison with experimental data [100] and Monte Carlo calculations [43]. The results

are plotted as a function of the reduced temperature T/TC.

and the T= 0 contribution is given by

Umag

T=0 =M0/2γJ0.(3.13)

Employing (3.11), the description of the heat capacity is quantitatively and qualitatively im-

proved as compared to MF indicating that RPA allows a more accurate prediction of the disconti-

nuity of the heat capacity (second order phase transition) at TC(see Fig. 3.6). Being proportional

to the magnetization dUmag/dT vanishes within both theories, MF and RPA, above TC, which, in

contrast to experiment, yields a vanishing magnetic contribution to the heat capacity CP. Both

approximations, MF and RPA, are not constructed for describing the region above TCbecause these

approximations do not capture the local magnetic order which is still present in this temperature

region. This restriction does not, however, aﬀect the absolute free energy as strongly as it aﬀects its

second derivative, the heat capacity. Figure 3.5 indeed reveals that the integrated RPA free energy

together with Fvib and Fel deviates from the calphad data by less than 15 meV at 1200 K.

Since the magnetization shape m(T)is the central quantity entering the internal energies of

both theories, the description of the speciﬁc heat is inherently connected with the description

of m(T). Therefore, we investigate next whether the remaining diﬀerences between ab initio and

experimental data can be decreased by using an improved description of the temperature-dependent

magnetization m(T). The RPA relation between the magnetic internal energy URPA

mag (T)and the

magnetization m(T)follows from Eq. (3.11). The fact that m(T)is the central quantitiy appearing

in Eq. (3.11) is related to the important role of hSzi ∼ m(T)within the RPA decoupling, see also

Eqs. (2.96) and (2.97). In Fig. 3.8, we show the magnetization curves obtained with MF and RPA

and compare them with experimental results. The classical solution found by numerically exact

Monte-Carlo simulations [43] and a Langevin function are also shown for comparison. Clearly, both

52 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

approximations, MF and RPA, are not able to accurately reproduce the experimental data. Aiming

at a better description of the free energy and speciﬁc heat capacity for fm bcc iron, we, therefore,

use an improved magnetization curve

mimp(T) = 1−s(T/TC)3

2−(1 −s)(T/TC)41

3,(3.14)

that was empirically derived from a detailed analysis of the experimental m(T)data for several

ferromagnetic systems [111], and fulﬁlls Bloch’s famous T3/2-law for low temperatures. The only

free parameter sin Eq. (3.14) can be calculated using ab initio input [112].11 The obtained shape

of the magnetization describes perfectly the experimental data, as shown in Fig. 3.8 (black line).

The improved magnetization given by Eq. (3.14) is now used as input to Umag in MF and RPA

(labeled iMF/iRPA).12 For the improved MF, the calculated speciﬁc heat capacity (violet solid line

in Fig. 3.6) shows the correct qualitative behavior and comes signiﬁcantly closer to experiment.

Please note that for the improved MF we not only correct for m(T)but also for the critical temper-

ature TCbeing signiﬁcantly overestimated in MF. For the free energy, however, the deviation from

experiment becomes more pronounced and even undergoes a sign reversal (Fig. 3.5). Therefore,

we conclude that even an application of the improved temperature-dependent magnetization is not

suﬃcient to account for the magnetic contribution to the free energy within MF approximation.

In contrast, a noticeable improvement of both, CPand F, is obtained for the improved RPA. The

experimental heat capacity is now nearly perfectly predicted up to TC. The diﬀerence between

the theoretical integrated free energy with a magnetic contribution calculated with iRPA and the

corresponding calphad data decreases to less than 10 meV at 1200 K.

Conclusion

The big advantage of analytic approaches such as MF or RPA is that they allow physical insights into

the diﬀerent mechanism contributing to the magnetic contribution to the thermodynamic properties.

We were, e.g., able to quantify the relation between the magnetization shape and the free energy

and speciﬁc heat capacity.

However, it is also found that RPA as well as MF do not capture the short-range order above

TCproperly. Reminding that the structural transition of pure iron at about 1200 K takes place just

above the critical temperature TC= 1043 K of bcc Fe, a careful treatment of the short-range order

in iron is of particular importance. The inclusion of the short-range order in our concepts will be

the content of the next sections.

11The calculated s=0.41 is close to the empirical value s=0.35 found by fitting Eq. (3.14) to experimental data [111].

The ab initio values used to calculate the theoretical sare: The ground state magnetization M0= 2.23 µB; the spin-

wave stiffness D= 249 meVÅ2obtained from the low excitation spectrum ω(q→0) ≈Dq2(which agrees with other

theoretical works [15]); TRPA

C=1082 K.

12We note that replacing m(T)within MF and RPA yields to internal inconsistency of the theories. However, such

an ansatz can be useful to figure out routes for improved theoretical descriptions.

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 53

3.2.2 Numerical Approaches

In the previous section two diﬀerent analytic approximate solutions for the Heisenberg model were

discussed. It was thus possible to get a combined ab initio concept, which allows a very accurate

prediction of thermodynamic properties of magnetic systems up to the critical temperature. In this

section we extend the methods and focus on the thermodynamic description in particular above

the critical temperature. In the paramagnetic regime, i.e. T > TC, the analytic approximations

discussed so far (MF/RPA) are limited due to the neglect of short-range order. Consequently,

one important point we would like to address in this section is to obtain a reliable prediction in

particular at high temperatures. A second issue is fundamentally rooted in the approximate nature

of the analytic approaches: The error due to the approximations are not known a priori. Clearly, a

numerically exact solution to overcome these issues would be highly desirable.

The QMC introduced in Sec. 2.2.3 is a promising candidate for this purpose. However, as already

discussed, its main limitation is the negative sign problem accompanied by an increase of statistical

error occurring if magnetic frustration is present. As shown in Sec. 2.2.3, the eﬀective weights p(c)

of these conﬁgurations cbecome negative. QMC simulations can still be performed by employing

the absolute value |p(c)|to obtain the expectation value, e.g. for the heat capacity (2.108):

CV≡ hCVi=PcCV(c)p(c)

Pcp(c)=PcCV(c)sgn(p(c))|p(c)|/Pc|p(c)|

Pcsgn(p(c))|p(c)|/Pc|p(c)|=hCVsgni|p|

hsgni|p|

.(3.15)

The statistical error ∆sgn|p|on hsgni|p|increases, however, exponentially with decreasing temper-

ature T, Eq. (2.110). To illustrate this, we employ QMC calculations using realistic exchange

parameters taken from Ref. [15]. Indeed, as shown in Fig. 3.9 (a), we observe a negative sign

problem reﬂected in hsgni|p|approaching zero for decreasing temperature, if we include more and

more interaction shells. This goes along with a dramatic increase of the statistical error on hsgni|p|

resulting consequently in a large statistical error on CV, Eq. (3.15), shown in Fig. 3.9 (b).

Employing (2.110)

∆sgn|p|

hsgni|p|

=qhsgn2i|p|−hsgni2

|p|

pNsteps ≈1

pNstepshsgni|p|

,(3.16)

we can estimate how much computational time a “brute force” method would take (i.e. by simply

running the simulation longer and increasing the number of iteration steps Nsteps) to achieve conver-

gence within a reasonable statistical error. Estimating the computational time needed to reduce the

statistical error given in Eq. (3.16) to less than 1% (typically needed for an accuracy of <1meV of

free energy calculations), we end up with infeasible long calculation times in the order of thousands

of years on a conventional PC (see Fig. 3.9 (c)).13 We therefore conclude: A numerically exact

quantum solution of the Heisenberg Hamiltonian using realistic exchange parameters (e.g. for bcc

iron) is not feasible with the present techniques!

The question arising from the fact is, of course: If we do not have access to the numerically

13Conventional PC means here explicitly a computer having 2×AMD Opteron DP 250 2.4 GHz and 4 GB RAM.

54 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

0 1 2 3

T/TC

10-4

10-3

10-2

10-1

100

<sgn>|p|

0 1 2 3

T/TC

0.0

0.5

1.0

1.5

2.0

CV(kB)

2 shells (FM)

3 shells

5 shells

10 shells

0 1 2 3

T/TC

100

103

106

109

cpu time (hours)

1 week

1 year

1000 years

(a)

(c)

(b)

Figure 3.9: (a) Expectation value of the sign hsgni|p|for increasing number of interaction shells. (b)

Corresponding heat capacity CV. (c) Estimated CPU time on a conventional PC for 1% statistical

error on hsgni|p|. The exchange parameters are taken from [15].

exact quantum solution, how much do we gain from deriving a numerically exact classical solution

instead? In the classical description, the total energy is directly given by

ECL =−X

Jijeiej(3.17)

where ˜

Jij =JCL

ij S2are the eﬀective exchange parameters including the length S=|S|of the

(continues) spin vectors Si=Seion lattice site iand eidenoting the unit vectors. Hence, for

each conﬁguration of spins the energy (3.17) can be straightforwardly computed. The classical

expectation values can be readily obtained via

hOiCL =PiOie−βECL

Pie−βECL

,(3.18)

where the sum runs over all ipossible conﬁgurations. To tackle the above equation, cMC techniques

are nowadays well established. Typically, the conﬁguration space is sampled using Metropolis

importance sampling i.e. a series of conﬁgurations is generated following usually the Boltzmann

statistics.

The introduced error for a classical treatment is, however, not known in advance. The critical

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 55

discussion of the limits of a classical treatment as well as the impact of quantum eﬀects on various

thermodynamic properties are in the focus of the next section.

Rescaled Classical Monte Carlo Approach

In the following we propose an alternative scheme to calculate the thermodynamic properties of

magnetic Heisenberg systems where magnetic frustrations are weak or absent. The essence of the

approach is a numerically exact classical Monte Carlo (cMC) solution of the (classical) Heisenberg

model [Eq. (3.17)]. The quantum-mechanical eﬀects which are neglected in the classical treatment of

the originally quantum-mechanical Heisenberg model are afterwards approximately included. This

inclusion is divided into two parts: On the one hand, we ensure the quantum-mechanically correct

high temperature limit of the free energy by aligning the calculated classical entropy to the analyt-

ically known quantum solution. On the other hand, the impact of the neglected low temperature

quantum eﬀects is estimated by a careful comparison of model systems for which the numerically

exact quantum solution via QMC calculations is also accessible (no negative sign problem). The

comparison allows to formulate a semi-empirical rescaling approach (in the following denoted as

rMC), which modiﬁes the classical solution towards the quantum-mechanical counterpart. To show

the performance of the method we ﬁnally combine the rMC approach with the ab initio computed

exchange coeﬃcients of Sec. 3.1 as well as ab initio vibronic and electronic free energy contributions

and compare with the results obtained in the previous section.

The use of DFT based input parameters in Eq. (2.70) and the subsequent solution of the model by

means of classical cMC is nowadays a standard and successful approach [13, 43, 48, 54]. Nevertheless,

having in mind that quantum eﬀects are not accounted for within the cMC approach, a number

of crucial questions arises: (i) How signiﬁcant are quantum eﬀects for thermodynamic properties

like, e.g. the heat capacity CV(T), the internal energy U(T)or the magnetic free energy Fmag(T)?

(ii) Above which temperatures Tand for which spin numbers Sare classical calculations suﬃcient?

(iii) Is there a practical scheme to estimate the magnitude and inﬂuence of quantum eﬀects without

having the quantum solution available, i.e. based solely on the classical results?

In order to tackle these questions we have evaluated the thermodynamic properties for a set

of carefully selected magnetic systems using the exact QMC data as a reference and comparing

them with the corresponding data obtained employing the cMC approach. This comparison pro-

vides direct insight into the importance of quantum eﬀects as function of spin quantum number and

temperature. The choice was limited to magnetically unfrustrated systems to avoid the negative

sign problem in the QMC approach. We included both kinds of magnetic interactions (ferromag-

netic/antiferromagnetic) for various 3D bulk lattice types (sc, bcc, fcc) and used diﬀerent numbers

of interaction shells (one and two shells).

To allow for a convenient comparison of the results for various spins S, and in particular for the

classical limit S→ ∞, we use common eﬀective exchange parameters ˜

J=S(S+ 1)JQM =S2JCL,

where JQM and JCL are the nearest-neighbor exchange integrals respectively in the quantum and

classical case. In the following we express all energies in units of ˜

Jand temperatures in units of

J/kB.

All cMC and QMC calculations were performed using the ALPS code [72]. For both calculations,

cMC and QMC (the Loop algorithm [67] in the stochastic series expansion method is used), cluster

56 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

12 3 4

0.5

1.5

rMC

QMC

cMC

12 3 4

0.5

1.5

12 3 4

0.5

1.5

12 3 4

0.5

1.5

Temperature (J/kB)

CV (kB)

S=1/2

S=7/2

S=3/2

S=1

Figure 3.10: Heat capacity vs. temperature for a 8x8x8 bcc lattice with ferromagnetic nearest-

neighbor coupling solved using cMC (spin-independent) and QMC for various spin quantum numbers

S. For comparison the rMC data (explained in the text) is shown.

updates were employed to calculate at least 2·105steps for equilibrating the system and 2·106

steps for thermal averaging. All supercells contained 512 atoms. For each conﬁguration ≈400

temperatures were evaluated. These settings ensure statistical errors of less than 0.01kBon CVin

the considered temperature range except for the region close to the critical temperature. Increasing

the supercell to 1728 atoms results in a more pronounced peak in CVat the critical temperature,

but has only a marginal eﬀect on the formalism discussed below. In the following we illustrate

our ﬁndings on the example of a bcc lattice with ferromagnetic nearest-neighbor coupling. The

other investigated systems yield very similar dependencies and conclusions and will be therefore not

discussed in detail here.

Within the ALPS code, the speciﬁc heat capacity CVfor cMC and QMC is calculated by virtue

of the ﬂuctuation-dissipation theorem (see, e.g. Ref. [113]) and shown in Fig. 3.10. The critical

temperature TCwas determined by the position of the peak in CV. For all considered spin quantum

numbers (S= 1/2,1,3/2,7/2), the classical and quantum solutions for CVagree well for suﬃciently

high temperatures (T&3˜

J/kB). Apart from the lower spin quantum numbers, i.e. for S > 1, the

classical approach provides even a good description of the quantum solution for temperatures above

TcMC

C≈2˜

J/kBimplying that the short-range order (dominating in the paramagnetic regime) is

well captured by the cMC approach and quantum eﬀects have a minor impact in this temperature

regime. At the same time the quantum solution of TC(S)is rapidly approaching the classical one

(indicated by a dashed line in Fig. 3.10) with increasing spin: TC(S→ ∞) = TcMC

C≈2˜

J/kB. For

a detailed discussion of the dependence of TC(S)we refer to Ref. [114]. The fact that the classical

solution for the heat capacity does not vanish for T→0(giving rise to a diverging magnetic

entropy), is a direct consequence of the classically non-discretized spectrum of magnetic excitations

in the Heisenberg model.

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 57

12 3 4

-4

-3

-2

-1

QMC

rMC

cMC

12 3 4

-4

-3

-2

-1

12 3 4

-4

-3

-2

-1

12 3 4

-4

-3

-2

-1

Temperature (J/kB)

Internal Energy (J)

S=1/2

S=7/2

S=3/2

S=1

Figure 3.11: Internal energy vs. temperature for a 8x8x8 bcc lattice with ferromagnetic nearest-

neighbor coupling solved using cMC and QMC for various spin quantum numbers S. For comparison

the rMC data is shown.

From the numerical calculations we further obtain the internal energy shown in Fig. 3.11. The

cMC and QMC results agree again for temperatures larger than T&3˜

J/kB. For lower temperatures

(T.2˜

J/kB) there are dramatic diﬀerences between the classical and quantum solutions for the

internal energy, which increase with decreasing spin quantum numbers. The largest diﬀerence occurs

at T= 0 and is equal to

UQM

T=0(S)−UCL

T=0 =−(˜

J0S/(S+ 1) −˜

J0) = ˜

J0/(S+ 1) (3.19)

with ˜

J0=Pi,j ˜

Jij.

The magnetic entropy Smcan be directly calculated from the heat capacity via integration either

starting from low temperatures (T= 0), or, if the maximum entropy Smax

m=Sm(T→ ∞)is known,

from the high temperature limit:

Sm(T) = ZT

CV(T′)/T′dT′(3.20)

=Smax

m−Z∞

CV(T′)/T′dT′.(3.21)

For the quantum spin system represented by Eq. (2.70), the maximum entropy is spin dependent,

Smax

m(S) = kBln(2S+ 1), reﬂecting that in the high temperature limit all 2S+ 1 spin states are

equally occupied.

Due to the above-mentioned non-vanishing classical heat capacity at T= 0, the corresponding

magnetic entropy is not well deﬁned for T→0. To compare the classical results with the quantum

results, we therefore align the classical maximum entropy to the quantum one using (3.21). The

corresponding calculated entropies are shown in Fig. 3.12. With decreasing temperature the entropy

58 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

12 3 4

0.25

0.5

0.75

12 3 4

0.25

0.5

0.75

12 3 4

0.5

1.5

12 3 4

0.5

1.5

QMC

rMC

cMC

Temperature (J/kB)

Entropy (kB)

S=1/2

S=7/2

S=3/2

S=1

Figure 3.12: Magnetic entropy vs. temperature for a 8x8x8 bcc lattice with ferromagnetic nearest-

neighbor coupling solved using cMC and QMC for various spin quantum numbers S. For comparison

the rMC data is shown.

decreases from the paramagnetic limit given by Smax

m(S)for fully disordered spins. Due to the treat-

ment of spins as continues vectors, the classical solution overestimates the change in the magnetic

entropy with temperature and yields a negative entropy at temperatures below a spin-dependent

oﬀset value Toff(S). The range of temperatures T < Toff(S)where cMC predicts a qualitatively

wrong physical behavior of the system is shaded in Fig. 3.13 with light-grey color. Nevertheless,

the alignment of the classical entropy to the quantum limit ensures that the free energy obtained

with cMC for high temperatures (T > 2˜

J/kB) accurately reproduces the quantum solution.

Despite the aligning scheme, the cMC approach still suﬀers from the neglect of the quantum

eﬀects, which turn out to signiﬁcantly inﬂuence the free energy data even at moderate temperatures

(up to TC). Having the importance of the classical approach for realistic materials science problems

in mind, one of the central points of the present study is, therefore, the improvement of the cMC

scheme in this temperature regime.

To achieve this goal we develop a mean ﬁeld-based semi-empirical approach, which provides on

the one hand a practical and eﬃcient way to account for quantum eﬀects and on the other hand

keeps the full ﬂexibility of the cMC scheme. We concentrate on CVsince all other thermodynamic

quantities can be easily calculated once the heat capacity is known [115]. The two main deviations

between the classical and quantum solution are a) the shape of CVand b) the spin-dependence

of the critical temperature. Since the latter would increase the complexity of our approach, we

introduce a normalized temperature t=T/TCand focus on the deviations regarding the shape as

a function of temperature and spin quantum number. These can be expressed in a compact form

by introducing

f(t, S, σ) := CQMC

V(t, S, σ)/CcMC

V(t, σ),(3.22)

where σincludes other dependencies apart from Sand t, such as, e.g. the lattice and/or magnetic

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 59

12 3 4

-3

-2

-1

12 3 4

-5

-4

-3

-2

12 3 4

-6

-5

-4

-3

-2

12 3 4

-8

-6

-4

QMC

rMC

cMC

Temperature (J/kB)

Free Energy (J)

S=1/2

S=7/2

S=3/2

S=1

Figure 3.13: Free energy vs. temperature for a 8x8x8 bcc lattice with ferromagnetic nearest-neighbor

coupling solved using cMC and QMC for various spin quantum numbers S. For comparison the

rMC data is shown. The range where the magnetic entropy of the aligned cMC and the rMC data

becomes negative is highlighted by light and dark grey shading, respectively.

conﬁguration (e.g. number of interaction shells). In principle, if fwere known, CQMC

Vcould be

computed solely based on the classical CcMC

Vsolution. A numerical estimate of fcan be obtained

for the model systems such as those considered in Fig. 3.10, but is not available for conﬁgurations

present in realistic material systems such as bcc iron.

To analyze the shape and functional dependence of the scaling function fwe compute this

ratio between the quantum-mechanical and the classical heat capacity for an extensive set of model

systems for which QMC results are accessible. The set of chosen conﬁgurations σincludes diﬀerent

lattice structures (sc, bcc, fcc) with nearest-neighbor ferromagnetic conﬁgurations (J1>0) as well

as antiferromagnetic conﬁgurations (J1<0) for bcc and sc [116]. To check the robustness of

the approach, we additionally included ferromagnetic conﬁgurations with second-nearest-neighbor

interactions taking J2= 2.5J1>0and J2= 5J1>0. The numerical results of Eq. (3.22) for two

diﬀerent spins, S= 1 and S= 7/2, are shown in Fig. 3.14. They show a signiﬁcant scatter for high

temperatures, where the quantum and classical solution should coincide, f(t≫1) →1. This is due

to CV(t≫1) →0[see Eq. (3.22)] and the connected increase of the statistical errors. For the more

important regime below t= 1, however, the inﬂuence of σon the function fis surprisingly weak

and signiﬁcantly smaller as compared to that of Sand t.

The observed weak dependence of the function fon the speciﬁc conﬁguration σthat includes

besides the crystal structure and magnetic state, also the range of interaction, is a central insight

of this study and allows to neglect the σdependence in ﬁrst order. For a practical realization

of this scheme we determine the averaged scaling function f(t, S)that averages over all considered

conﬁgurations. This function is expected to provide a good approximation also for realistic material

systems, which have not been included in the current set of structures (e.g. due to the not accessible

QM solution). A practical example will be given in Sec. 3.2.2 for the material system bcc iron. For

60 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

0.5 11.5 2

0.2

0.4

0.6

0.8

f(t,S,σ)

bcc fm

bcc afm

fcc fm

sc fm

sc afm

bcc J2=2.5J1 (fm)

fcc J2=2.5J1 (fm)

sc J2=2.5J1 (fm)

bcc J2=5J1 (fm)

fcc J2=5J1 (fm)

sc J2=5J1 (fm)

S=7/2

S=1

Figure 3.14: Numerically evaluated f(t, S, σ)ratio (symbols) [see Eq. (3.22)] vs. temperature for

two spins S= 1 and S= 7/2and various σ(lattice/magnetic conﬁgurations). Symbols: (circles,

squares, diamonds)=(sc, bcc, fcc); colorcode: (black, blue, orange, green)=(fm, afm, J2= 2.5J1

(fm), J2= 5J1(fm)).

an easy implementation of this approach an analytical expression for f(t, S)is desirable. In the

following we discuss an empirical ansatz for this function, and test its accuracy/limits by comparing

the rescaled cMC data based on the empirical ansatz with numerically exact results.

We start this discussion by considering the mean-ﬁeld approximation of the Heisenberg model,

which is known to describe thermodynamic properties at low temperatures reasonably well (see,

e.g. Figs. 3.5- 3.8). An advantage of this approximation is the availability of analytical expressions

for the magnetization in the quantum (Brillouin function) and in the classical (Langevin function)

case. The resulting ratio of heat capacities, fMF(t, S), can therefore straightforwardly be computed.

The shape of fMF(t)is in good agreement with that of fin the low temperature regime except of a

scaling factor due to the well known overestimation of TCin the MF approximation. To elaborate

this in more detail we analytically expand the low temperature MF ratio (see Appendix A.1.3) and

obtain

fMF(t≪1, S)≈2tS/t

exp(tS/t)2

(3.23)

with the spin-dependent temperature

1/tS=α(S+ 1).(3.24)

αis a proportionality factor. Its MF value is αMF ≡2/3.

For temperatures above the critical temperature (t > 1) MF predicts CV≡0due to the neglect of

short-range order. Therefore, the analytic function Eq. (3.23) needs to be extended by an additional

term, such that the resulting empirical function fapprox fulﬁlls the following boundary conditions:

First, for high temperatures both solutions (quantum and classical) should coincide, i.e. fapprox(t→

∞, S)→1. Second, the same should hold for an inﬁnite spin quantum number, i.e. fapprox(t, S →

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 61

0.5 11.5

0.5

f(t,S)

S=7/2

S=3/2

S=1

S=1/2

(a)

0.5 11.5 22.5 33.5

0.5

1.5

2.5

3.5

1/tS

(b)

fit

Figure 3.15: (a) Averaged data from Eq. (3.22) for diﬀerent spin quantum numbers Sof a bcc

nearest-neighbor ferromagnet calculated using ≈400 temperature steps (solid lines) and corre-

sponding ﬁt fapprox(t, S)obtained using Eq. (3.25) (dashed lines). (b) Linear relation of 1/tS

versus spin S. Symbols: (circles, squares, diamonds)=(sc, bcc, fcc); colorcode: (black, blue, orange,

green)=(fm, afm, J2= 2.5J1(fm), J2= 5J1(fm)).

∞)→1, reﬂecting that the classical model converges to the QM model in the limit of an inﬁnite

spin S. The choice

fapprox(t, S) := 2tS/t

exp(tS/t)−exp(−tS/t)2

(3.25)

obeys both conditions and provides a good ﬁt of the data as will be discussed next.

In a second step we ﬁt the numerically exact scaling data Eq. (3.22) to the functional dependence

of Eq. (3.25), using tSas ﬁt parameter. This procedure has ﬁrst been performed for each of the

investigated spin quantum numbers S, each structure (sc, bcc, fcc) and magnetic conﬁguration

(J > 0,J < 0,J2= 2.5J1>0,J2= 5J1>0) individually. As can be seen from Fig. 3.15 (b), the

used ﬁtting function again appears to be remarkably independent of the conﬁguration σ. Regarding

the dependence on S, the parameter 1/tSturns out to be nearly linearly proportional to the spin

quantum number, similarly to MF in Eq. (3.24). This is consistent with the observation that the

characteristic temperature above which the cMC and the QMC solution agree is inverse proportional

to S. Performing a linear regression over all calculated datasets for diﬀerent systems yields

1/tS≈0.54(1)S+ 0.54(2),(3.26)

close to the mean ﬁeld value as given by Eq. (3.24). The mean deviation on the last digit is shown

in brackets.

The quality of the analytical expressions for the rescaling functions Eqs. (3.25) and (3.26), is

demonstrated in Fig. 3.15 (a). Here, we observe for all considered parameters a fair agreement

between f(t, S)and fapprox(t, S), with the smallest deviations for low temperatures tand high

spin values S. In order to test the accuracy of this approach, we apply the rescaling scheme to the

previously discussed cMC data of model systems. This allows a direct comparison of the temperature

62 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

dependences obtained by the two-step procedure (ﬁtting+linearizing 1/tS) resulting in fapprox(t, S)

with the numerically exact QMC results. In contrast to Fig. 3.15 (a), it has been performed for all σ

individually and would only be trivial if the rescaling function f(t, S, σ)were used. Since in practice

the quantum solution is usually not known a priori, we also do not align the classical and quantum

solution with respect to the temperature, i.e. we do not enforce TrMC

C=TQMC

C. As an example we

show the heat capacity for the bcc, ferromagnetic nearest-neighbor case (black line in Fig. 3.10).

For all spin quantum numbers the classical results are found to be signiﬁcantly improved. A similar

accuracy is found for all other considered systems.

Based on the rescaled heat capacity we compute, similar to Eq. (3.21), the internal energy

UrMC(T) = UcMC(TH)−ZTH

CrMC

V(T)dT. (3.27)

Here we have used that UcMC(TH)≈UQMC(TH)and CV(TH)≈0for suﬃciently high cutoﬀ

temperatures TH> TC. The rescaled internal energy Eq. (3.27) and the rescaled entropy Eq. (3.21)

are compared in Figs. 3.11 and 3.12 with the cMC and QMC results. Altogether, the classical de-

scription is again improved for all considered spin quantum numbers and thermodynamic properties.

As expected, deviations are larger for small spin quantum numbers (e.g. S= 1/2) (see Fig. 3.10).

A major source for the remaining discrepancies are the deviations between the Curie temperatures

for small spins and in the classical limit S→ ∞ [114]. In general, with increasing spin number the

predictive accuracy of rMC improves and already for S≥3/2the results are very close to the full

QMC calculations. A similar behavior is found for the internal energy and magnetic entropy.

Generally we ﬁnd that the rescaling approach consistently improves the problematic low tem-

perature regime. For integrated quantities such as the entropy, internal, and the free energy, the

remaining errors accumulate in the low temperature region due to the backward integration from

high to low temperatures by virtue of Eqs. (3.21) and (3.27). A consequence is the presence of

unphysical negative entropies at low temperatures (dark shaded temperature range in Fig. 3.13)

similarly to the purely classical approach (light shaded temperature range in Fig. 3.13), despite the

good description of the rescaled heat capacity both at low and high temperatures with the largest

deviations close to the critical temperature.

To quantitatively evaluate the accuracy achievable by the rescaling approach we summarize

in Tab. 3.2 the mean deviation in percent (standard deviation of the various structures shown

in brackets) between the rMC (ﬁrst rows) and the cMC (second rows) with the QMC data for

diﬀerent temperature ranges and averaged over all model structures σconsidered in this study.

For low temperatures (below the critical temperature) the aligned cMC approach is not suﬃcient

to describe the quantum system whereas the rMC approach provides already for S≥1a good

approximate prediction of the quantum-mechanical free energy with an error of less than 4%. For

higher temperatures (above the critical temperature) the rMC reduces the deviations between the

cMC and the QMC data for all thermodynamic quantities by almost a factor of two. The fact

that the scatter between the diﬀerent structures, as reﬂected in the standard deviation shown in

Tab. 3.2, is almost always smaller than the mean deviation and that it consistently decreases with S,

indicates that the rescaling approach works equally well for all investigated lattice types, magnetic

conﬁgurations, and number of included interaction shells. We, therefore, expect a similar accuracy

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 63

Table 3.2: Mean value of the mean deviation between the rMC/cMC and the QMC calculations for

the heat capacity, internal energy, entropy, and free energy in percent using ∆A = 1/NσPNσ∆Aσ

with ∆Aσ= 1/N Pti(ArMC/cMC

σ(ti)−AQMC

σ(ti))/AQMC

σ(ti)with N= 100 in two temperature

regimes for all structures σ(sc, bcc, fcc; fm, afm, J2= 2.5J1(fm), J2= 5J1(fm)). The standard

deviation of the diﬀerent structures is shown in brackets.

0.1< T/TcMC

C<1

∆A S = 1/2S= 1 S= 3/2S= 7/2

A=CVrMC 41.2(7.1) 21.9(4.3) 15.8(2.9) 5.5(2.0)

cMC 1897.3(951.9) 885.5(413.8) 533.3(171.4) 112.3(12.1)

A=UrMC 32.2(11.9) 9.8(5.4) 4.3(3.4) 1.0(1.2)

cMC 125.5(17.9) 53.9(6.2) 31.5(3.4) 9.(1.0)

A=SmrMC 226.2(81.0) 48.9(42.8) 28.2(32.7) 9.3(10.3)

cMC 3630.(1614.4) 1810.3(646.4) 1087.4(343.0) 223.1(57.0)

A=FrMC 12.8(5.0) 3.9(2.2) 1.7(1.4) 0.7(0.7)

cMC 49.3(6.2) 20.8(2.3) 11.9(1.3) 3.4(0.6)

1< T/TcMC

C<2

∆A S = 1/2S= 1 S= 3/2S= 7/2

A=CVrMC 49.8(25.9) 15.2(9.7) 7.1(5.2) 2.2(3.4)

cMC 99.7(34.8) 33.4(11.3) 16.9(5.6) 3.7(1.4)

A=UrMC 14.1(9.5) 5.2(4.2) 2.7(2.6) 1.1(1.8)

cMC 25.7(10.4) 10.7(4.2) 5.9(2.3) 1.5(0.6)

A=SmrMC 3.7(3.1) 1.3(1.3) 0.8(0.8) 0.7(0.5)

cMC 6.5(3.3) 2.2(1.4) 1.1(0.8) 0.6(0.5)

A=FrMC 1.1(0.8) 0.5(0.3) 0.5(0.3) 0.6(0.4)

cMC 1.5(0.9) 0.7(0.4) 0.6(0.4) 0.7(0.4)

of the proposed rescaling approach for other lattice geometries and if interaction shells beyond

second-nearest-neighbors are included.

Application of the rMC Approach to FM Bcc Iron

In the following we employ the rMC approach to compute the magnetic contribution of bcc iron as

an example of a realistic magnetic material system. Despite its predominant ferromagnetic coupling,

positive (i.e. antiferromagnetic) exchange integrals occur in higher-order shells. These exchange

integrals cannot be neglected due to the very long-ranged character of the interactions [15] and lead

to weak magnetic frustrations which rule out the use of the QMC approach (see also Fig. 3.9) [117].

To compute the thermodynamic properties we use the adiabatic approximation (3.4). In Fig. 3.16

we show once more the sum of Fvib and Fel and compare it with calphad [86, 87] data. The strong

deviations between the calculated and the experimental free energy and heat capacity (Fig. 3.17)

above room temperature highlight the signiﬁcance of magnetic contributions, as we already discussed

in Sec. 3.2.1 [16].

In order to evaluate the magnetic contribution, we start with the comparison of the classical

64 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

200 400 600 800 1000 1200

T (K)

-700

-600

-500

-400

-300

-200

-100

F (meV/atom)

vib+el

vib+el+MF

vib+el+RPA

vib+el+cMC

vib+el+rMC

CALPHAD

300 600 900 1200

-50

Figure 3.16: The free energy of bcc iron as a function of temperature. Inset: Diﬀerence between

the diﬀerent theoretical and the calphad data is shown. For comparison the theoretical MF and

RPA results of Sec. 3.2.1 are also shown.

Monte Carlo calculations and the rMC data for the speciﬁc heat capacity (Fig. 3.17). Classical Monte

Carlo calculations were performed using a 12x12x12 supercell (containing 1728 atoms) including

magnetic interactions of up to 25 neighbor shells. The eﬀective real-space exchange coeﬃcients ˜

Jij

entering the Heisenberg model have been obtained by performing a Fourier transformation of the

reciprocal-space exchange integrals, which we have already extensively discussed in Sec. 3.1. The

magnetic heat capacity is calculated assuming that the magnetic moments of Fe do not signiﬁcantly

change due to the thermal volume expansion and temperature, resulting in CV(T)≈CP(T). We

now add the magnetic contribution to CP(T)(Fig. 3.17). The peak in the heat capacity (Fig. 3.17)

reﬂects the magnetic transition from the ferromagnetic to the paramagnetic state. The resulting

critical temperature is TcMC

C≈1060 K which is close to the experimental value Texp

C= 1044 K

and similar to the analytic result TRPA

C= 1083 K (Sec. 3.1). In the paramagnetic regime the cMC

approach reproduces the experimental behavior very well, conﬁrming our ﬁnding in the previous

section that magnetic short-range order is well captured within the cMC calculations (for S > 1)

and that quantum eﬀects are not critical in this temperature regime (see also Fig. 3.10).

To obtain FcMC

mag we use the calculated cMC internal energy UcMC and the classical entropy ScMC

aligned to the maximum entropy Smag

m(S)for the quantum spin system at high temperatures [see

Eq. (3.21)]. The eﬀective spin quantum number S= 1.1, which corresponds to the ab initio calcu-

lated ground-state magnetic moment 2.2 µB, is used. In order to obtain the magnetic excitation

energies (the T= 0 K-energy is already included in Fel), the T= 0 K contribution UQM

T=0(S)is sub-

tracted from the calculated FcMC

mag after the alignment of Fand the calphad data.14 The resulting

free energy FcMC is shown in Fig. 3.16. The obtained agreement at high temperatures is excellent.

14As discussed above, the classical treatment of the quantum Heisenberg Hamiltonian fails in particular at low

temperatures where, however, the computed free energies are typically aligned to calphad data (i.e. 200 K) for

comparison. After the alignment we, therefore, subtract the deviation caused by the neglected quantum effects

employing ∆FDFT−calphad (200 K) = FDFT

cMC (200 K)−Fcalphad(200 K) := UQMC(200 K)−UcMC(200 K)≈˜

J0/(S+1),

3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity 65

0 200 400 600 800 1000 1200

T(K)

CP(kB)

vib+el

vib+el+MF

vib+el+RPA

vib+el+cMC

vib+el+rMC

Experiment

bcc fcc

Figure 3.17: Heat capacity vs. temperature for the diﬀerent theoretical data vs. experiment. The

experimental transition temperature (bcc →fcc)Tαγ ≈1184 K is indicated by a dashed line. For

comparison the theoretical MF and RPA results of Sec. 3.2.1 are also shown.

Above the critical temperature, the aligned cMC free energy provides a very good description of the

magnetic contribution to the free energy of bcc iron. At lower temperatures, however, the deviations

are clearly visible and increase in magnitude with decreasing the temperature. Below ≈700 K the

aligned cMC entropy becomes negative (highlighted by light grey color), which leads to a wrong

asymptotic behavior of the calculated free energy. This region is, therefore, not properly described

within the cMC ansatz, although quantitatively the deviations of FcMC from the experimental data

remain moderate even at temperatures below 700 K.

We use now the rescaling framework and multiply the cMC data for CP(T)with the rescaling

function given by Eqs. (3.25)-(3.26). The results are shown in Fig. 3.17 (black diamonds) and

reveal a signiﬁcant improvement over the cMC results below the room temperature, yielding an

overall excellent agreement of the calculated heat capacity with the available experimental data.

The largest deviations occur close to the critical temperature where the ﬁnite size eﬀects in the

MC calculations hinder an accurate description of the heat capacity in the thermodynamic limit

of inﬁnite system size L→ ∞ [118]. Applying our rescaling ansatz to the cMC entropy and

internal energy we derive the rMC free energy (black straight line in Fig. 3.16). As can be seen,

our rescaling approach provides an accurate agreement with the experiment with an error bar of

less than 10 meV/atom over the whole temperature range of thermodynamically stable bcc iron.

We note that the temperature below which the rMC entropy becomes negative (indicating the

breakdown of this approach) is signiﬁcantly shifted by ≈300 K to lower temperatures (marked by

the dark shaded region) when compared to the aligned cMC. The results show that the rescaling

see (3.19). In the last step UQM

T=0 K≈UQM

T=200 Kis assumed due to the negligible contribution of CP(T < 200K)≈

0 = dU/dT .

66 3.2 Thermodynamics of Magnetic Metals: Free Energy and Specific Heat Capacity

procedure is indeed able to correct for the cMC deﬁciencies in the low temperature (T < TC) regime

and lead to a signiﬁcant decrease of the errors in the calculated thermodynamic properties.

Conclusion

Classical MC calculations employing realistic magnetic exchange coeﬃcients can be used to accu-

rately determine the Curie temperature TCas well as to account for magnetic short-range order

above TC. The numerical treatment provides therewith a complementary approach in particular to

the analytic RPA method (accurate description in the ferromagnetic regime below TC) enabling us

to describe bcc iron over the whole temperature range. Below TC, the cMC method is not suﬃcient

to describe the quantum Heisenberg model. This implies that quantum eﬀects are important even

up to TC, which is in case of bcc Fe already above 1000 K!

The almost universal behavior of the derived rescaling function for the diﬀerent tested model

systems indicates that the overall shape of thermodynamic quantities, such as the speciﬁc heat

capacity, are largely independent on the speciﬁc magnetic conﬁguration if considered on a reduced

temperature scale. This insight provides the basis for an eﬀective model discussed in Sec. 3.4.

Having tested and veriﬁed the developed methods for bcc iron we proceed in the next section

with the investigation of fcc iron. A major aim is to resolve the fundamental question of the

underlying driving forces of the bcc - fcc phase transitions.

3.3 Structural Phase Transitions in Iron 67

3.3 Structural Phase Transitions in Iron

So far we have concentrated our discussion in particular on the magnetic transition in bcc iron,

i.e., the transition from the ferromagnetic to the paramagnetic regime at the Curie temperature TC

(lower part in Fig. 3.18). We have developed highly eﬃcient methods to incorporate the magnetic

entropy contribution into a combined ab initio framework for a parameter-free prediction of the free

energy F(T). Since this procedure can be applied to several phases, it now provides the possibility

to investigate even structural phase transitions.

400 600 800 1000 1200 1400 1600 1800 2000

ferromagnetic paramagnetic

-bcc

-fcc liquid

-bcc

1184 K

1673 K

T (K)

1044 K

< 1 meV

∆

Ftot (meV)

-5 0 5 10

Figure 3.18: The calphad [119] free energy diﬀerence ∆Ftot between the bcc and fcc phase is

shown (red line). For ∆Ftot <0the bcc phase is more stable than the fcc phase and for ∆Ftot >0

vice versa.

The structural transitions in pure iron have been intensively investigated from experimental

as well as theoretical sides for more than half a century [7, 8, 120–123]. However, due to the

complexity involved in the magnetic degrees of freedom, a full theoretical understanding is still

lacking. One of the key issues is the detailed balance between the various entropy contributions

and their impact on the structural stability of the bcc and fcc phase, resulting in questions such

as: are phonon excitations the driving force of the bcc to fcc transition [120] or are magnetic

contributions playing the more signiﬁcant role [8]? Based on pure geometry one could argue that

phonon excitations in a more open structure are softer than in a more close packed lattice structure

and that consequently the vibronic entropy should stabilize the bcc structure at high T. If we

accordingly assume that vibronic contributions are the driving force for the second transition (γ→

δ), the magnetic contributions are consequently responsible for the ﬁrst one (α→γ). However,

a solely geometric argumentation (“more open structure”) does not necessarily apply to magnetic

systems, where for instance the amplitude of phonons sensitively depends on the magnitude of the

local magnetic moments [124].

As mentioned in the beginning, the bcc-fcc phase transitions in iron were investigated by many

authors over the last ﬁfty years: In 1956, within a calphad-type assessment of experimental data,

Weiss and Tauer [120] demonstrated the importance of vibronic contributions for the stabilization

of the fcc phase. Later, in 1963, Kaufman et al. [7] suggested the existence of two spin states in

fcc Fe (called ferromagnetic high-spin and low-spin state), which yield an additional entropy con-

tribution which is responsible for the stabilization of the fcc phase (sometimes also denoted as the

68 3.3 Structural Phase Transitions in Iron

phenomenological two-state-hypothesis).15 Twenty years later, Hasegawa and Pettifor [8] proposed

a purely theoretical ansatz based on a mean ﬁeld single-band Hubbard model. In this way they

were able to show that both phase transitions could be explained qualitatively solely by magnetic

contributions. However, the proposed theory included a number of approximations (no spin-wave

contributions, ignoring short-range order, etc.) and a ﬁnal answer regarding the impact of magnetic

contributions compared to the vibronic ones was not possible. Instead, a few years later, Zarestky

and Stassis reanimated the discussion of the impact of vibronic contributions by performing neutron

scattering experiments for annealed fcc iron [122]. Based on these experiments, the Debye temper-

ature of fcc Fe was estimated to be about 25%below the bcc Fe Debye temperature strengthening

the initial calphad evaluation of Weiss and Tauer. In 2010, Lavrentiev et al. [125] employed a

newly developed magnetic cluster expansion [126] to investigate the phase transitions in iron. The

authors combined the computed magnetic entropy with empirically derived vibronic contributions16

and were able to reproduce the experimentally observed phase transitions. Their approach relies,

however, on a number of assumptions as well as experimental input.17 Furthermore, important

contributions such as the electronic entropy were neglected. A fully theoretical understanding of

the impact of diﬀerent contributions (electronic, vibronic, and magnetic) on the phase transitions

is obviously not yet established.

The following section is divided into two parts: The ﬁrst part is devoted to the treatment of

the fcc phase. Previous DFT calculations reveal a non-collinear spin spiral ground state for fcc iron

(see e.g. Ref. [19]). In this case, the developed RPA theory in Sec. 3.2.1 needs to be extended

and generalized to non-collinear magnetic conﬁgurations. We apply the more general analytic

expression as well as the rescaled Monte Carlo approach from Sec. 3.2.2 to fcc Fe and compare

them with each other. The results are further combined with vibronic and electronic entropy

contributions and compared with the results obtained by the calphad approach. Having obtained

a reliable prediction of the fcc phase, we discuss in the second part in detail the temperature

dependent free energy diﬀerence of both phases, ∆Fbcc−fcc, as well as the impact of the various

entropy contributions.

3.3.1 Thermodynamics of γ-Iron

The experimental investigation of the fcc phase of iron is complicated due the fact that it is ther-

modynamically stable only in a rather high and narrow temperature window between 1184 K and

1673 K, see Fig. 3.18. Experimentally fcc iron is stabilized mostly in combination with fcc Cu,

e.g., by annealing Fe alloys in Cu (forming Fe precipitates of the parent (fcc Cu) lattice structure)

[128–130] or by epitaxial growth of thin Fe layers on Cu(001) (up to a critical layer thickness)

[131–138].

15We note that although the existence of the two spin states is nowadays well established by electronic structure

calculations, see e.g. Ref. [121], it is usually not considered as an additional source of entropy.

16The contribution of bcc Fe is derived from experimental elastic constant measurements at 1173 K [127]. For fcc

Fe, a force-constants model is employed to fit experimental inelastic neutron scattering data taken at 1428 K [122].

17Concerning the magnetic model, the long-ranged exchange interactions are ignored and only the first two (three)

nearest-neighbor exchange coefficients for bcc (fcc) Fe are included. The error introduced due to the classical treatment

of the magnetic Hamiltonian is not known and also not discussed by the authors. The vibronic contribution is

determined by fitting experimental data at fixed temperature which is even not the same for the bcc and fcc phase.

The dependence of the vibronic contribution on volume and temperature is not included. Finally, the electronic

entropy is completely neglected.

3.3 Structural Phase Transitions in Iron 69

For precipitates sizes up to 150 ˚

A, a simple type-I antiferromagnetic (afm) structure is found

[139]. In the ultrathin ﬁlms, a nonmonotonous variation of the magnetization was found, which

was assigned to the formation of a helical spin-wave (characterized by Qexp = 2π/a(0.1,0,1)) in

the deeper layers being close to an afm conﬁguration Qafm = 2π/a(0,0,1) in (001)-direction. The

experimental results (for the precipitates as well as the layered structures) show a strong tendency

of antiferromagnetic order accompanied with potential distortions (modulations) depending in a

crucial way on the external constraints (e.g. the experimental ﬁlm growing technique or precipitate

size).

To investigate the magnetic ground state of fcc iron, several DFT studies have been performed,

see e.g. Ref. [19] and references therein. Regarding collinear magnetic conﬁgurations [121], ferro-

magnetic high-spin and low-spin phases are found to be energetically more stable than the non-

magnetic (nm) phase. The lowest energies are found for antiferromagnetic phases, in particular for

the antiferromagnetic single-layer (afm) and antiferromagnetic double-layer (afmd) conﬁgurations

in (001)-direction. In a previous work [19], which includes diﬀerent spin spiral conﬁgurations in the

total energy calculations, two local energy minima are found for spin spirals Q1= 2π/a(0.2,0,1)

(close to the experimental one found in ultrathin ﬁlms, see above) and Q2= 2π/a(0,0,0.6). At

rather low volumes (about 7%lower than the experimental volume) the spin spiral state Q1turns

out to be more stable than Q2, whereas for larger volumes the latter is found to be energetically

more favorable.

To compute the magnetic entropy we employ in the following a theory similar to bcc Fe, which

covers the full Q-space, i.e., the energies corresponding to all possible spin spiral conﬁgurations.

Therefore, similar calculations as in [19] are performed. Regarding the vibronic entropy, phonon

calculations for diﬀerent spin spirals may become computational very expensive18. We, therefore,

perform phonon calculations in two representative collinear phases, i.e., the afm and afmd phase

discussed above.

Ground-state Properties

We ﬁrst consider the T= 0 K total energy calculations for various magnetic phases (including the

above introduced afm and afmd phase) of fcc Fe in Fig. 3.19. For the considered collinear magnetic

phases, the afmd phase represents the energetically lowest magnetic structure. We performed in

addition total energy calculations for diﬀerent non-collinear spin spiral structures Qand ﬁnd that the

Q2structure is lower in energy than the collinear phases in the considered volume range. We note,

however, that all magnetic phases are very close in energy indicating that the speciﬁc magnetic

arrangements of the magnetic moments in fcc Fe have only small inﬂuence on the total energy

compared to bcc Fe, see Fig. 3.19. Since the energy diﬀerence ∆E between the fm and afm phase

is in ﬁrst order approximation19 proportional to the nearest-neighbor exchange coeﬃcient ˜

J12, the

magnetic exchange in fcc Fe (∆E ≈50 meV≈580 K·kB) is about one order of magnitude smaller

than in bcc Fe (∆E ≈450 meV≈5200 K·kB).

18For the spiral state Q2= 2π/a(0,0,0.6) at least 10 atomic layers in the (0,0,1)-direction are required to construct

the unit cell. In addition, the symmetry breaking in non-collinear calculations renders such calculations to be at least

one order of magnitude more expensive than common collinear phonon calculations.

19For a classical Heisenberg system defined by the energy E= 1/N Pij Jij~ei~ej, the energy difference at T= 0 K

between the fm and afm configuration is ∆E =PjJ1,2j≈J1,2.

70 3.3 Structural Phase Transitions in Iron

65 70 75 80 85 90

Volume / atom (Bohr3)

100

200

300

400

500

600

Energy / atom (meV)

bcc afm

bcc fm

bcc nm

fcc afm

fcc afmd

fcc fm

fcc nm

fct afmd

fcc Q=0.6XΓ

Figure 3.19: The total energy ground-state phase diagram of Fe in the fcc as well as the bcc structure

for various magnetic phases. The energy diﬀerence between the diﬀerent magnetic conﬁgurations is

much larger in bcc than in fcc Fe indicating a smaller magnetic exchange in the latter, see text for

details.

We included in the ground-state search the tetragonal distorted fct structure. Allowing in

addition relaxations (internal coordinates as well as cell shape) in the total energy calculations, the

afmd fct structure (breaking the axial symmetry) shown in Fig. 3.19 is pronounced as the lowest

magnetic phase resulting in c/a ≈1.09.

To illustrate that in the ideal fcc structure the Q2state provides the lowest spin spiral structure

of all possible q-spirals, we show in Fig. 3.20 (b) the computed spin spiral spectra at two volumes

V= 10.61 Å3and V= 11.19 Å3. As mentioned above, a second local minimum occurs at Q1=

2π/a(0.2,0,1) close to the experimental as well as previous theoretical ﬁndings [19]. Comparing

the spin spiral spectra for fcc Fe (Fig. 3.20 (b)) with the previously obtained spectra for bcc Fe

(Fig. 3.20 (a)) we directly observe that the eﬀective magnetic exchange coeﬃcients ˜

Jqare indeed

about one order of magnitude smaller in the fcc than in the bcc phase. This is consistent with the

various magnetic phases in Fig. 3.19 being close in energy. Also the Γ-point (ferromagnetic phase)

represents no longer the energetic minimum, which reﬂects the fact that the system is magnetically

unstable in the ferromagnetic phase. At both considered volumes the spin spiral with the wave

vector Q2represents the magnetically most stable conﬁguration.

3.3 Structural Phase Transitions in Iron 71

ΓN P ΓH

100

200

300

400

500

∆Eq (meV)

ΓX W L Γ

-100

-75

-50

-25

∆Eq (meV)

for V=V0

afm

for V=V0

afmd

Figure 3.20: Total energy of diﬀerent spin spirals qas given by Eq. (2.61). Above: bcc Fe for

V= 11.21Å3. Below: fcc Fe for two diﬀerent volumes Vfcc,afm

0= 10.61Å3and Vfcc,afmd

0= 11.19Å3,

see text for details.

Analytic Approximations

We start the investigation of the magnetic entropy and the magnetic critical temperature of fcc Fe

by considering analytic approaches, i.e., the RPA and the MF solution. In the context of the RPA,

the main diﬃculty for a system showing a non-collinear spin spiral ground state arises from the fact

that there is no global quantization axis available anymore (see Fig. 3.21). A common approach

for analytic approaches in such a case is to express the Hamiltonian in a locally rotated coordinate

system which allows eventually again to work with one global quantization axis. As proposed in

Ref. [140], we rotate the local coordinate system Σiof each spin Siinto a new one Σ′

iparallel to

the z-axis. The spins in the unrotated system Σican be expressed in the rotated coordinates via







Six

Siy

Siz





=





cos(QRi)S′

ix + sin(QRi)S′

S′

−sin(QRi)S′

ix + cos(QRi)S′





,(3.28)

where S′

ik is the k-component of the spin S′

iat lattice site iin the rotated system and Riis the

lattice vector. We remark that such a general solution naturally includes the ferromagnetic solution

Q= (0,0,0) which we already applied in Sec. 3.2.1. The Heisenberg Hamiltonian (2.83) in this

72 3.3 Structural Phase Transitions in Iron

Figure 3.21: Sketch of non-collinear spin spiral state.

(locally rotated) coordinates reads as [140]

H′=−X

ij KijS′

ixS′

jx +LijS′

iyS′

jy +MijS′

izS′

jz+ 2 X

NijS′

ixS′

jz,(3.29)

with anisotropic exchange constants

Kij =Mij =Jij cos [Q(Ri−Rj)] ,(3.30)

Nij =Jij sin [Q(Ri−Rj)] ,(3.31)

Lij =Jij.(3.32)

Analogously to the case of the ferromagnet, the free energy Fmag we are interested in can be obtained

by deriving the internal energy Uand integrating it afterwards, see Eq. (2.82).

The key idea is exactly the same as for the ferromagnet, i.e., to rescale a S= 1/2RPA solution.

For the sake of clarity we will discuss the main results for the non-collinear RPA solution. The

intermediate steps are discussed in the appendix in detail. We proceed in a similar way as in

Sec. 3.2.1:

In the ﬁrst step we derive a rescaling factor γwhich allows to rescale the RPA internal energy

solution for S= 1/2into a general spin solution. For this purpose we consider ﬁrst the expression

for the critical temperature TNof a non-collinear magnet for S= 1/2as well as S≥1/2and obtain

(see Sec. A.3.1 for a detailed derivation):

S≥1/2kBTN=4

3S(S+ 1) X

q1

K−Kq

K−Jq!−1

,(3.33)

S= 1/2kBTN= 2S X

q1

K−Kq

K−Jq!−1

.(3.34)

Comparing the solutions of TN, the rescaling factor γfor the normalized exchange coeﬃcients ˜

Jis

deduced as

S(S+ 1)J≡˜

J⇒SJ =: γ˜

J, (3.35)

3.3 Structural Phase Transitions in Iron 73

with γ= 1/(S+ 1) = 2/3for S= 1/2. We want to point out that the known expression for the

critical temperature above in terms of commonly employed eﬀective exchange coeﬃcients ˜

kBTRPA

N´eel =4

3 X

q"1

K−˜

Jq#!−1

,(3.36)

has already been derived and successfully applied to europium [140, 141] and is consistent with

the conventional solution Eq. (2.101) of the ferromagnet (Q= (0,0,0), i.e., ˜

Kq≡˜

Jq), which we

employed to investigate the pressure dependence of TC(p)in the bcc phase in Sec. 3.1. The major

task for the next step is to derive a S= 1/2expression for the internal energy and replace all SJ

by γ˜

To compute the magnetic entropy, we need in the second step an expression for the RPA internal

energy for S= 1/2within the rotated coordinate system (3.29), which could be afterwards combined

with Eq. (3.35). We want to stress that such a theory is not available in the literature. In principle

our derivation is similar to the one for the ferromagnet, but involves a more lengthy derivation due

to the rotated coordinate system. For the sake of readability we present here the ﬁnal results only

and provide the full derivation in the appendix A.3.2. For the internal energy Uwe obtain

URPA(T) = U0+M

2m(T)X

qA+

q+A−

q+B+

q+B−

q,(3.37)

where A±

qand B±

qare given by

A±

q=χq±±ωq+ 2γ(˜

K−˜

Jq)

exp(±βωq−1) ,(3.38)

B±

q=±2γm(T)(±ωq−2γ˜

Jq)( ˜

K−˜

Jq)

4Eqexp(±βωq−1) ,(3.39)

and with anisotropic exchange coeﬃcients

Kq= 1/2( ˜

JQ+q+˜

JQ−q).(3.40)

Here, the reduced magnetization is given by

mS=1/2

RPA (T) = 1

1 + 2ϕ,(3.41)

with the eﬀective magnon occupation

ϕ(T) = 1

qχq+

eβωq−1+χq−

e−βωq−1.(3.42)

The eﬀective magnon excitation energies are given by

ωq= 2γmS=1/2

RPA r˜

K−˜

Kq˜

K−˜

Jq(3.43)

74 3.3 Structural Phase Transitions in Iron

and

χq±=1

21±aq

ωq(3.44)

represent the weights for magnon creation (“+”) and annihilation (“-”) respectively with

aq= 2γm(T)( ˜

K−1/2( ˜

Kq+˜

Jq)).(3.45)

Eqs. (3.37)-(3.45) form a closed system of equations which can be solved self-consistently. It rep-

resents the generalization of the theory proposed in Sec. 3.2.1 and is equivalent to the latter for

Q= (0,0,0) (ferromagnetic ground state).

For comparison we also employ a known MF theory for non-collinear spin systems, where the

internal energy is given as [31]

UMF(T) = M0hS′

zi2=S

S+ 1 ˜

M0m2

MF(T) = S

S+ 1 ˜

JQm2

MF(T).(3.46)

The reduced magnetization mMF(T)can be determined in full analogy to (2.87), resulting in a

Brillouin function which is characterized by the MF critical temperature

kBT′MF

N´eel =2

3˜

JQ.(3.47)

The Néel Temperature of Fcc Iron

Before discussing the diﬀerent entropy contributions (vibronic, electronic, and magnetic) we consider

ﬁrst the Néel temperature by employing the derived expressions within MF and RPA, Eq. (3.47) and

(3.36). Similar to fm bcc Fe, we employ the eﬀective exchange coeﬃcients ˜

Jobtained at the ground-

state volume. As discussed above, there are many, nearly degenerate, energetically metastable

magnetic phases. We consider, therefore, in the following two representative structures, i.e., the

afm and the afmd case and perform the spin spiral calculations at the corresponding ground-state

volumes, which are Vafm

0= 10.61 Å3and Vafmd

0= 11.19 Å3. Within MF, TMF

N´eel is only determined by

the lowest spin spiral energy, Eq. (3.47), and we obtain TMF,afm

N´eel = 183 K and TMF,afmd

N´eel = 235 K. We

recall that an experimental measurement of TN´eel for pure bulk fcc iron is not possible because it is

not stable at low temperatures. Fcc iron is often stabilized by alloying it with manganese (Mn). For

this reason, we compare in Fig. 3.22 the MF results with the phenomenological calphad prediction

as well as experimental data for fcc Fe-Mn alloys and diﬀerent Mn-concentrations. Consistent with

our ﬁndings for ferromagnetic (fm) bcc iron, the MF approximation provides an upper limit for

TN´eel being slightly above the calphad estimation.

The RPA approximation provided a very reasonable estimation of TCfor the ferromagnetic

systems. Employing (3.36), we obtain TRPA

N´eel = 141 K/199 K for the exchange coeﬃcients extracted

at the ground-state volumes of the afm and afmd phase, respectively. Due to the weaker ˜

Jij in fcc

iron compared to the bcc phase, both approximations (MF/RPA) provide very similar predictions.

The much lower critical temperature in the austenite compared to the ferrite is a direct consequence

of the much ﬂatter spin spiral dispersion in Fig. 3.20. The results for TRPA

N´eel are smaller than

the corresponding MF values and also closer to the calphad prediction, see Fig. 3.22. Apart

from the analytic methods, the classical Monte Carlo approach provided for bcc iron also a very

3.3 Structural Phase Transitions in Iron 75

0 0.1 0.2 0.3 0.4 0.5

Mole-fraction Mn

100

200

300

400

500

600

TNeel (K)

Experiment

calphad

MF for afm

MF for afmd

RPA for afm

RPA for afmd

cMC for afm

cMC for afmd

Figure 3.22: Computed Néel temperature in fcc iron in comparison with experimental data [142]

for diﬀerent Mn-concentrations. The empirical calphad estimation [142] is shown for comparison.

reasonable estimation of the critical temperature. As a crosscheck, cMC calculations were also

carried out for fcc20 based on the computed exchange coeﬃcients in real space (obtained by a

Fourier transformation). The critical temperatures we obtain are, as in the bcc case, similar to the

ones obtained using the analytic RPA approach. Explicitly, we obtain TcMC

N´eel = 122 K/159 K for the

afm/afmd phase. To estimate the theoretical uncertainty on the TN´eel prediction we average over

the three approximations and the two magnetic structures and obtain

Tγ,theo

N´eel ≈173 ±41 K,(3.48)

where 41 K denotes the standard deviation.

The Helmholtz Free Energy of Fcc Iron

After achieving a reliable description of the magnetic phase transition temperature we pass over

to the Helmholtz free energy, Eq. (3.4). For this purpose we identify ﬁrst the afmd tetragonal

distorted fct structure as the energetically lowest collinear magnetic phase from the the total energy

calculations in Fig. 3.19. To estimate the uncertainty due to the chosen magnetic ground-state

conﬁguration, we performed (similarly to the determination of the critical temperature) each of the

steps for the afm phase as a crosscheck as well.

To derive the vibronic contribution to the total free energy within the quasiharmonic approx-

imation (Sec. 2.2.4), phonon calculations were carried out at diﬀerent volumes for both magnetic

phases (afm/afmd). Since there are diﬀerent atom types in the unit cell (spin-up and spin-down),

optical branches appear in addition to the acoustic ones.21 We show in Fig. 3.23 the phonon density

20The cMC calculations involve 4×106steps including thermalization and thermal averaging. The calculations

are performed for a system containing N= 1728 spins and taking into account up to 25 interaction shells.

21For the afm structure, a small part of the phonons were imaginary (less than 3%), ruling out the direct use of the

76 3.3 Structural Phase Transitions in Iron

10 20 30 40

ωPh(meV)

500

1000

1500

phonon DOS (arb. units)

bcc, fm

fcc, afmd

fcc, afm

0 10 20 30 40

ωPh(meV)

integrated phonon DOS (arb. units)

300 K

Figure 3.23: Phonon DOS of fcc iron deduced from the antiferromagnetic single- (red shaded areas)

as well as double-layer phase (blue shaded area). The DOS of the fm bcc phase is shown for

comparison. The vertical dotted line corresponds to the thermal energy at room temperature

(300 KkB≈25 meV).

of states (DOS) deduced from the fcc afm and afmd phases (inset) as well as the integrated phonon

DOS. We directly observe that the DOS of the afm phase are shifted to lower energies compared

to the afmd phase indicating a higher vibronic contribution. In order to quantify the impact of

the chosen magnetic phase (afm or afmd) on the vibronic free energy contribution, we compare in

Fig. 3.24 the results of the qh approximation with calphad data. Within the qh approximation

the impact of the chosen magnetic phase (solid and dashed green line) is in the order of 10 meV.

This indicates the robustness of the approach with respect to the chosen magnetic ground-state

conﬁguration. Regarding the absolute agreement strong deviations (∼150 meV at 1800 K) between

the vibrational analysis and the empirical calphad data are found. Based on our previous inves-

tigations for bcc Fe, a considerable amount of entropy stems from the electronic degree of freedom

which has to be taken into account for free energy considerations.

After calculating the electronic contribution, we also compare the free energy of both (phonon

+ electronic) contributions in Fig. 3.24 with the calphad data. It is worth mentioning that with

respect to the starting conﬁguration (afm or afmd) the diﬀerence of the total free energy between

the two considered phases is less than 30 meV at 1800 K which corresponds to a relative deviation

of less than 3%. However, in contrast to the bcc phase, the missing magnetic contribution to the

total free energy compared to the calphad data is already visible at relatively low temperatures.

This is consistent with the computed Néel temperature, Eq. (3.48), being in the order of the room

temperature.

To compute the magnetic contribution to the free energy we apply the derived expressions for the

qh approximation (negative argument of logarithm in Eq. (2.116)), and indicating that the afm structure is unstable

for certain distortions. Since only a small part of the computed phonon frequencies are affected, we apply a common

first order approximation and ignore such phonon modes for the evaluation of the vibronic free energy.

3.3 Structural Phase Transitions in Iron 77

500 1000 1500 2000

200 400 600 800 1000 1200 1400 1600 1800 2000

T (K)

-1200

-1000

-800

-600

-400

-200

F (meV)

Vib (afm)

Vib (afmd)

Vib+El (afm)

Vib+El (afmd)

RPA (afm)

RPA (afmd)

rMC (afmd) +Vib+El

rMC (afm)

MF (afm)

MF (afmd)

calphad

Figure 3.24: Free energy of fcc iron vs. calphad data for diﬀerent treatments of the magnetic

contribution. Inset: The calphad data is taken as the reference at each temperature.

MF and RPA, Eqs. (3.47) and (3.36), and employ in addition the developed rMC ansatz (Sec. 3.2.2)

based on the cMC data (which has been employed to derive TcMC

N´eel ). For the latter, the spin quantum

numbers derived from the ground-state moments Mfcc,afm

0= 1.39µBand Mfcc,afmd

0= 1.67µBare

used for the rescaling function f(S, t)provided in Eq. (3.25). The computed total free energies

including the magnetic contributions are also shown in Fig. 3.24. For both considered phases the free

energy signiﬁcantly improves with respect to the experimental data. It is interesting to note that,

in contrast to bcc iron, even the MF approximation allows an excellent estimation of the magnetic

entropy contribution. This is related to the much lower critical temperature of fcc compared to bcc

iron and the fact, that MF reproduces the high temperature limit of fully disordered spins.

We can conclude that the developed methodologies allow a promising (even quantitative) pre-

diction of the full Helmholtz free energy of γ-iron. Starting the analysis from two representatives

(afm and afmd phase), the total computed free energies provide a conﬁdence interval for the true

experimental energy with an energy window of about ∼20 meV per atom. Starting from other

(magnetic) phases (e.g. paramagnetic phase) should therefore not change the obtained results cru-

cially. We proceed in the next section with the discussion of the relative diﬀerence between the bcc

and the fcc phase and the discussion about the diﬀerent impact of vibronic, magnetic, and electronic

contributions on the phase stability of these phases, respectively.

3.3.2 Impact of Vibronic, Electronic, and Magnetic Contributions on the Struc-

tural Phase Transition in Iron

In thermodynamic equilibrium the phase transitions between bcc and fcc correspond to a double

crossing of the underlying Helmholtz free energies of Fe, i.e., of the bcc phase, Fbcc(T), and the fcc

phase, Ffcc(T). Deﬁning the free energy diﬀerence ∆Ftot =Fbcc −Ffcc, a positive value ∆Ftot >0

corresponds to a stabilization of the fcc phase whereas a negative value ∆Ftot <0indicates a stable

78 3.3 Structural Phase Transitions in Iron

0500 1000 1500 2000

T (K)

∆Fvib (meV)

vib: bcc - fcc (afm)

vib: bcc - fcc (afmd)

vib+el: bcc - fcc (afm)

vib+el: bcc - fcc (afmd)

Figure 3.25: Impact of electronic and vibronic contributions on the stabilization of the bcc and fcc

phase of iron. The dotted lines denote the experimental observed transition temperatures between

the bcc and fcc structure.

bcc structure. We show in Figs. 3.18 the free energy diﬀerence ∆Ftot obtained by the calphad

method [119]. It provides very important information for our further analysis: the observed diﬀer-

ence clearly shows that the theoretical accuracy needed to reproduce the experimental ﬁndings is

in the order of ∼1meV for temperatures up to 2000 K.

Regarding the thermodynamic methods to capture the magnetic entropy contributions we have

developed so far, a numerical accuracy of ∼1meV is indeed feasible (accompanied with extremely

accurate and expensive calculations though). With the term "numerical accuracy" we speciﬁcally

denote all errors which can be controlled ensuring suﬃciently high convergence parameters, for

instance for the DFT calculations the chosen k-point grid, cutoﬀ energy or supercell, or for the

Monte Carlo simulations the total number of measurement steps, number of included interaction

shells and chosen supercell. All these parameters can be chosen such that the required convergence

criterion (∼1meV) is fulﬁlled.

We denote the remaining source of uncertainty as methodological errors which can not be de-

creased by choosing "better" convergence parameters. Regarding the treatment of the magnetic

entropy, this includes the error due to the simpliﬁcation of the magnetic system (mapping onto an

eﬀective Heisenberg Hamiltonian) as well as the necessary approximations to solve the model (rMC,

MF or RPA). In addition, on the DFT level, the potential error due to the approximative treatment

of the xc-functional. For bcc and fcc iron we found empirically that the methodological errors are

at the considered high temperatures in the order of ∼10 −20 meV. Regarding the DFT input,

similar methodological limits are found for the quasiharmonic and electronic entropy considerations

(see e.g. [5, 80]). We, therefore, do not aim at a quantitative agreement with the calphad data in

the following but rather concentrate on balance and impact of the diﬀerent entropy contributions

which is still a open question in the literature.

Diﬀerent authors [120, 122, 143] suggested that the vibronic contribution of the fcc phase has

an important impact on the stabilization of the fcc phase and hence on the bcc-fcc transition (at

3.3 Structural Phase Transitions in Iron 79

0500 1000 1500 2000

T (K)

∆F mag (meV)

MF (afm)

MF (afmd)

rMC (afm)

rMC (afmd)

RPA (afm)

RPA (afmd)

0500 1000 1500

T(K)

0.2

0.4

0.6

0.8

1.2

Smag(kB)

critical temperature

fcc bcc

Figure 3.26: Left side: Impact of magnetism on the stabilization of the bcc and fcc phase of iron for

diﬀerent treatments of the magnetic contribution. Right side: Magnetic entropy obtained within

MF for the fcc afm and the bcc phase.

1200 K). On the other hand, as mentioned above, from purely geometric arguments, the atoms in

the more open bcc structure should have more space to vibrate and hence it should be easier to

excite phonons. However, as shown in Fig. 3.23, the bcc phonon DOS is in fact altogether higher

in energy than the DOS of the fcc phase. This ﬁnding is consistent with recent ﬁxed spin moment

phonon calculations for bcc iron [124]. It is shown that the magnetic moment of bcc iron has an

considerable impact on the phonon spectra. As a consequence, the computed phonon modes for bcc

iron decrease in energy with decreasing moment [124]. Since the bcc phase has a higher magnetic

moment than the fcc phase, one can conclude that the smaller vibronic contribution to the total

free energy of the bcc phase can be attributed to its larger magnetic moments.

In order to elucidate its quantitative impact on the phase transition, we show in Fig. 3.25 the

diﬀerence of the vibronic contribution ∆Fvib =Fbcc

vib −Ffcc

vib. For both considered magnetic phases of

the fcc structure, the vibronic contribution is indeed larger than for the bcc phase, i.e., it stabilize the

fcc phase. This is consistent with the obtained phonon densities shown in Fig. 3.23 and strengthens

again the fact that a geometric argumentation is not suﬃcient in case of magnetic systems. The

impact of the electronic contribution on the phase transition is similar as found for the vibronic one,

stabilizing the fcc phase with increasing temperature (also shown in Fig. 3.25). Considering the

impact of both contributions on the phase stability leads us to a very important conclusion: Vibronic

and electronic contributions deﬁnitely drive the transition bcc→fcc, but neither the vibronic nor

the electronic contribution can explain the observed second transition (fcc →bcc). This clearly

shows the importance of magnetic contributions which have to be taken into account.

To clarify their role, we show in Fig. 3.26 the diﬀerence of the magnetic contribution of the bcc

and fcc phase, ∆Fmag =Fbcc

mag −Ffcc

mag. All theoretical treatments show qualitatively the same be-

havior, i.e., both transitions could be explained on qualitative level solely by taking into account the

magnetic contributions, which is in agreement with the results of Ref. [8]. The qualitative behavior

80 3.3 Structural Phase Transitions in Iron

200 400 600 800 1000 1200 1400 1600 1800

T (K)

-50

∆Ftot (meV)

rMC (afm)

rMC (afmd)

RPA (afm)

RPA (afmd)

MF (afm)

MF (afmd)

800 1200 1600

-10

Figure 3.27: Diﬀerence of the Helmholtz free energies of the fcc and the bcc phase of pure iron

in comparison with calphad data (shown as black diamonds) [119]. The theoretical results are

aligned to the calphad data at 200 K.

can be understood as follows: In general, the magnetic entropy of both phases increases rapidly

when approaching the critical temperature (maximum entropy for full disorder). However, since

the critical temperature of the fcc phase is much lower than that of the bcc phase (TN´eel ≪TCurie),

the contribution of the fcc phase dominates the diﬀerence up to ≈TC. At these high temperatures

the fcc phase already reached its maximum entropy while in the bcc phase a substantial amount of

magnetic entropy is still conserved in the magnetic order (shown in Fig. 3.26 for the MF result).

At very high temperatures, however, the bcc phase ﬁnally overcomes the magnetic contribution

of the fcc phase due to the higher magnetic moment, i.e., higher maximum entropy. Therefore,

the magnetic contributions tend to stabilize the fcc phase at low temperatures (due to the lower

critical temperature) and hence support (in combination with the vibronic contribution) the ﬁrst

phase transition (bcc→fcc). For higher temperatures the stronger local magnetic moments of the bcc

phase are responsible for the, in the end, higher magnetic contribution as well as the weaker vibronic

contribution of the bcc phase, both being responsible for the second phase transition (fcc→bcc).

We note that the quantitative impact of these eﬀects diﬀers for the diﬀerent theoretical approaches

(MF, RPA and rMC), since the latter provide slightly diﬀerent magnetic entropy contributions.

We ﬁnally compare in Fig. 3.27 the absolute free energy diﬀerences with the calphad data.

Since errors in the theoretical prediction accumulate with increasing temperature, the spread be-

tween diﬀerent methods increases at high temperatures. Overall, the theoretical calculations are

in fair agreement with experiment and correctly describe the ﬁrst transition from the bcc to the

fcc phase (apart from the fcc RPA afm result). However, regarding the second transition, a ﬁnal

conclusion whether the transition back to the bcc phase takes place or not, is - solely based on the

theoretical predictions - not possible. In fact, the error bar provided by the diﬀerent methods allows

either a transition back to the bcc phase as well as a stabilization of the fcc phase within an energy

window consistent with the methodological accuracy limits of the employed schemes.

3.3 Structural Phase Transitions in Iron 81

Conclusion

The scatter between the results of the diﬀerent theoretical treatments for the magnetic entropy

contribution of fcc Fe is much smaller compared to bcc Fe. The reason is that the critical temperature

and eﬀective spin moment of iron in the fcc lattice is much lower than in the bcc lattice. As a

consequence, a mean ﬁeld treatment already provided a reasonable description of the magnetic free

energy of the austenite above room temperature.

Regarding the impact of magnetism on the structural transitions, the two most important factors

are (again) the eﬀective spin quantum number (determining the high temperature entropy limit)

as well as the critical temperature. Both quantities were suﬃcient to understand qualitatively the

phase transitions in iron.

82 3.3 Structural Phase Transitions in Iron

3.4 Nearest-Neighbor Approach: Application to Fe, Co, and Ni 83

3.4 Nearest-Neighbor Approach: Application to Fe, Co, and Ni

In the beginning of Chapter 3 we discussed two possibilities for solving the Heisenberg model for

realistic systems namely (a) approximative solutions for a full (long-ranged) Hamiltonian and (b)

an exact solution of an eﬀective (simpliﬁed) one. In the previous sections we developed and applied

diﬀerent approaches for (a). In doing so we learned that, apart from the crystal structure and the

spin quantum number, also the critical temperature critically aﬀects many aspects of the underlying

magnetic model Hamiltonian. In addition, we observed a universal behavior of the thermodynamic

properties of the Heisenberg model in Sec. 3.2.2, if the critical temperature is kept ﬁxed. In the

following we employ this insight to propose an appropriate ansatz (b), which is more convenient

than solving the full Hamiltonian and provides a similar accuracy. For this purpose we introduce an

auxiliary, eﬀective nearest-neighbor Hamiltonian. The model is characterized by a nearest-neighbor

exchange parameter J, as well as the crystal structure and spin Sof the system, and can be

solved numerically exact by QMC. The previously discussed negative sign problem does not apply

here, since magnetic frustration does not occur for a nearest-neighbor ferromagnet. The eﬀective

Jparameter is constructed such that it contains the long-ranged interactions22 Jij in a mean ﬁeld

manner. Our previous investigations in Sec. 3.2.2 have revealed that the critical temperature TCis

the decisive quantity in this context and it needs to be ensured that it does not alter during the

mapping. Such an approach has the big advantage of decoupling the evaluation of TC(for which

the approximative methods such as RPA and the cMC work very well) from the computation of the

other thermodynamic properties, see Fig. 3.28.

If the critical temperatures are provided, such an approach is computational much less expensive

(due to the reduction to just one J-parameter) compared to a full cMC simulation (including

long-ranged exchange parameters). To test the transferability of the approach we include in the

following all 3d ferromagnets, namely Fe, Co, and Ni. Finally, the computed magnetic contribution is

combined with ab initio derived vibronic, and electronic contributions to access the full free energy

and heat capacities. To account for the full free energy we employ the adiabatic approximation

and separate the electronic, vibronic and magnetic contributions as in the previous sections. As

for bcc iron, electronic and vibronic contributions have been determined using the quasiharmonic

approximation and ﬁnite temperature DFT respectively, as in [6, 16, 58]. All QMC calculations

are done employing the direct-loop algorithm in the stochastic series expansion as implemented in

the ALPS code [72]. Monte Carlo calculations involve 2.5×106steps including thermalization and

statistical averaging.

The theoretical ground-state magnetic moments23 M0are used to deﬁne the spin quantum

number ˜

S=M0/(gµB)where g≈2denotes the Landé factor, see Eq. (3.7). We note that in

systems containing itinerant electrons the value of M0and hence of ˜

Sis not restricted to those of

the (ideally localized) quantum spin model (S= 1/2,1,3/2,···), i.e., ˜

Scould be any real number

S∈ ℜ. In practice we compute the exact solutions for both neighboring spin numbers of a given ˜

22Since by construction the exchange interactions are reduced using a single parameter (i.e. TC), the extension

to multi-magnetic alloys (systems containing different magnetic species in the unit cell and hence several effective

parameters) is not straightforward.

23The computed ground-state moments are: 2.21 µB(Fe), 1.64 µB(Co), and 0.63 µB(Ni).

84 3.4 Nearest-Neighbor Approach: Application to Fe, Co, and Ni

Ab initio:

Exchange parameters Jij ,

Spin S, crystal structure

Experiment:

Critical temperature TC

Full Hamiltonian:

H=−Pij Jij ~

Si~

Critical temperature TC

Effective Hamiltonian:

H=−JPij ~

Si~

Critical temperature of

effective model

Thermodynamic properties:

Free energy, Heat Capacity, ...

Input Input

Mapping Keeping spin and

crystal structure

Approx.

solution

Adjustment of

effective exchange

parameter J

Exact

solution

Exact solution

Figure 3.28: Sketch of the eﬀective nearest-neighbor approach. The DFT data is used as input

for the full Hamiltonian which is solved approximately to obtain TC. Alternatively, the critical

temperature can be taken from experiment. The eﬀective nearest-neighbor exchange parameter J

is chosen as such that it reproduces the provided TC. Based on this, the eﬀective Hamiltonian is

solved exactly providing the thermodynamic properties.

S1≤˜

S≤S2(S2−S1= 1/2). The physical quantities ˜

f(T, ˜

S)(e.g. heat capacity or magnetization)

are computed employing a linear interpolation afterwards, i.e.,

f(T, ˜

S) = αf(T, S2) + (1 −α)f(T, S1),(3.49)

with α= 2( ˜

S−S1). In practice we performed calculations for S= 1 and S= 3/2for bcc Fe, and

for S= 1/2and S= 1 for fcc Co and fcc Ni. For the latter a linear extrapolation instead of an

interpolation is used in order to estimate the solution for ˜

S= 0.3.

Apart from ˜

Sthe remaining input parameters for Eq. (2.83) are the crystal structure and the

critical temperature TC. As explained above, the latter is used to determine the effective nearest-

neighbor exchange coupling J.TCcould be taken either from a theoretical estimation [15, 48] or

from experiment. Here, the experimental24 TCis chosen (Texp

C→J). In practice, we ﬁrst perform

the model calculations for an arbitrary Jstart parameter. The results of the speciﬁc heat capacity

CPare subsequently ﬁtted to the given TCby the position of the peak in CP.

We start the discussion by investigating CP. Its computed temperature dependence is shown in

Fig. 3.29 in comparison with experimental data obtained under ambient pressure. We ﬁrst focus on

24We use the experimental values: 1044 K (Fe), 1388 K (Co), and 631 K (Ni).

3.4 Nearest-Neighbor Approach: Application to Fe, Co, and Ni 85

0500 1000 1500

vib+el+mag

vib+el

vib

Experiment

5,4

5,7

Temperature (K)

Specific Heat Capacity (kB)

bcc Fe

fcc Co

fcc Ni

rMC

Figure 3.29: The speciﬁc heat capacity CP(T)vs. temperature in comparison with experimental

data [107, 144–146]. The magnetic contribution is computed in a QMC calculation for N= 5832

spins.

the sum of the vibronic and electronic contributions (orange line). The fact that this sum reproduces

the experimental heat capacities very well up to ≈300 K nicely indicates that the systems are still

close to the magnetic saturation in this temperature regime. However, above room temperature the

importance of the magnetic contribution is clearly visible for all three magnets in agreement with

our previous ﬁndings for bcc iron. Please note that the ﬁrst peak in the case of Co at ≈700 K

occurs due to the hcp-fcc transition, which is not accounted for in the present approach (since both

phases are very similar [147] and the fcc phase is the dominant one at TC, we only consider fcc Co

here). It is worth mentioning that the electronic contribution increases from Fe to Co and is highest

for Ni. This is inherently related to the smallest band splitting in case of Ni (∼0.2−0.3eV),

which is in the order of the thermal energy at the considered temperatures. In this case, electrons

of both 3d sub-bands (spin-up and spin-down) can be thermally activated and contribute to the

total electronic entropy.

Within the QMC calculations the magnetic contribution to the heat capacity is directly com-

86 3.4 Nearest-Neighbor Approach: Application to Fe, Co, and Ni

-1200

-1000

-800

-600

-400

-200

-1200

-1000

-800

-600

-400

-200

500 1000 1500

-1200

-1000

-800

-600

-400

-200

0vib+el+mag

vib+el

vib

600 1200 1800

-50

-25

100

600 1200 1800

-50

-25

100

600 1200 1800

-50

-25

100

Temperature (K)

Free Energy (meV)

vibronic

electronic

magnetic

bcc Fe

fcc Co

fcc Ni

Figure 3.30: The free energy F(T)vs. temperature T(solid line) in comparison with calphad data

(symbols). The electronic and magnetic contributions are indicated by orange and blue shading,

respectively. The inset shows the diﬀerence between the computed and calphad values. For bcc

Fe, the results obtained by the rMC approach are shown.

puted employing the ﬂuctuation-dissipation theorem. All three ferromagnets are well described,

both below and above TC(see Fig. 3.29). Compared to Fe and Co, the magnetic contribution to CP

is smallest in the case of Ni due to the smallest spin quantum number (˜

S≈0.3), see Fig. 3.29. The

most signiﬁcant deviations from experiment are observed close to the critical temperature and are

most pronounced for bcc Fe. Increasing the system size in the QMC calculations yields a more pro-

nounced peak in CP(see inset in Fig. 3.29) but has only marginal inﬂuence on the overall shape.25

We note that from the numerical point of view it is challenging to estimate the maximum peak via

a ﬁnite size scaling ansatz (see e.g. Refs. [148–151]).26 At the same time the experimental scatter

is strongest around TCdue to the rapid change in CParound the critical point [144].

25The impact of the system size Non the free energy turns out to be negligible (i.e. well below 1meV) for the

considered Nand temperatures T.

26There is a long-standing debate whether the heat capacity of the 3D Heisenberg model shows a finite cusp or

diverges. This is inherently connected to the critical exponent αfor which different signs (cusp/divergent) have been

found. [148–151] However, even if the heat capacity shows a finite cusp, finite-size scaling approaches are very difficult

in practice [148].

3.4 Nearest-Neighbor Approach: Application to Fe, Co, and Ni 87

0,2

0,4

0,6

0,8

0,2

0,4

0,6

0,8

0,2

0,4

0,6

0,8

0 0,2 0,4 0,6 0,8 1

0,2

0,4

0,6

0,8

QMC,N=512

QMC,N=1728

QMC,N=5832

cMC,N=5832

Reduced Temperature

bcc Fe

fcc Co

fcc Ni

Reduced Magnetization

quantum effects

universal curve

taken from [111]

Figure 3.31: The reduced magnetization m(t) = M(t)/M0vs. the reduced temperature t=T/TC

for Fe, Co, and Ni in comparison with experimental data [152]. The numerical calculations are

carried out for three diﬀerent system sizes (N= 512,1728,5832).

We now turn to the discussion of the Helmholtz free energy Eq. (3.4). The results are shown in

Fig. 3.30 in comparison with calphad data. In order to compare the calphad with the ab initio

derived values we align them at 200 K as before. Similarly to the behavior of the heat capacity, the

electronic contribution to the total free energy increases from Fe to Co and reaches its maximum

for Ni (≈70 meV at 1800 K) due to the smallest band splitting in Ni as discussed above. The

remaining diﬀerences to the calphad free energy (lower orange line in insets of Fig. 3.30) clearly

reveal the necessity to account for the magnetic contribution Fmag.

The quantitative impact of Fmag on Fis reversed compared to the electronic contribution Fel

(see Fig. 3.30). The strongest magnetic contribution is obtained in the case of bcc Fe (≈100 meV

at 1800 K) due to the highest spin quantum number ( ˜

S≈1.1). Although Co has a larger spin

quantum number than Ni, the magnetic contributions in these systems is comparable. The reason

for this is that the lower TCrenders the Ni spin system close to the fully disordered state at the high

temperatures considered here, whereas in Co a noticeable amount of magnetic free energy remains in

the magnetic short-range order. Comparing the total free energies (including all contributions) with

88 3.4 Nearest-Neighbor Approach: Application to Fe, Co, and Ni

the calphad values we observe a very good agreement for all three ferromagnets. The strongest

deviations are observed for bcc Fe (≈30 meV at 1800 K) due to the poor description of CP

around TC(see Fig. 3.29) whereas a nearly perfect quantitative agreement is obtained for Co and

Ni (≈10 meV at 1800 K).

In Sec. 3.2.1 it has been shown that even sophisticated approaches like the RPA which give a

very good qualitative description of CPand F, can fail to accurately describe the magnetization

shape. We made similar observations for classical MC approaches [58]. It is therefore remarkable

that our eﬀective model, which reduces the long-ranged magnetic interactions responsible for the

magnetization to nearest-neighbor interactions only, also performs very well for m(t). The results

are shown in Fig. 3.31 for diﬀerent system sizes (N= 512,1728,5832) in comparison with exper-

imental data. The overall shape is in very good agreement for all three ferromagnets. For low

temperatures the exact solution of the quantum Heisenberg model ensures that the Bloch T3/2-law

for M(T)is fulﬁlled in agreement with experiment. For T→TCﬁnite size eﬀects are responsible

for the remaining tail in M(T)in the numerical results whereas the second order phase transition

(spontaneous symmetry breaking) is obtained in the thermodynamic limit of N→ ∞ only. Increas-

ing the system size Nreveals a steeper decrease of M(T→TC)(similar to the more pronounced

peak in CPin Fig. 3.29) but has only marginal eﬀect on the overall shape.

The fact that classical MC fails to accurately describe the magnetization shape indicates that

ﬁnite temperature magnetization is substantially aﬀected by spin quantization eﬀects. By perform-

ing classical MC calculations using a nearest-neighbor Hamiltonian, which reproduces the same TC,

we have a unique opportunity to quantify the quantum eﬀects. As can be seen in Fig. 3.31, the

magnetic contributions (light blue regions) are substantial even at high temperature (i.e. close to

TC) where quantum eﬀects are often assumed to be negligible.

Conclusion

The approach provides a straightforward mapping procedure that transforms the full magnetic

Hamilton operator of realistic systems (with long-ranged and oscillating interactions) onto an ef-

fective Hamiltonian with nearest-neighbor eﬀective interactions only. It allows to decouple the

estimation of TC- for which approximative methods (e.g., RPA or cMC) work very well - from

the determination of other thermodynamic properties. Particular important, the scheme provides a

high degree of transferability and performs equally well for Fe, Co, and Ni. This indicates that not

only the magnetization shape follows in a wide temperature window an universal behavior [111],

but also other thermodynamic properties such as, e.g., the speciﬁc heat capacity.

Whereas the considered systems are experimentally well studied, the proposed approach can also

be employed to study more complex magnetic systems for which experimental data is missing. In

the next section we apply our approach to such a system, i.e., cementite [153, 154], an industrially

important ferromagnetic Fe-C alloy.

3.5 Impact of Magnetism on the Phase Stability of Cementite 89

3.5 Impact of Magnetism on the Phase Stability of Cementite

In the previous section we employed a nearest-neighbor Heisenberg model to compute the mag-

netic entropy. In combination with the quasiharmonic approximation (vibronic entropy) and ﬁnite-

temperature DFT (electronic entropy), this approach provided accurate predictions of thermody-

namic data for (simple) ferromagnetic systems. In addition, the transferability of the method was

shown (all 3d ferromagnets were well described). The approach becomes therewith an ideal candi-

date to study even more complex ferromagnetic systems. In the following we employ this method

to gain a deeper insight into the phase diagram of Fe-C by investigating cementite (Fe3C).

Fe3Cis one of the key structures which appears in the Fe-C phase diagram, being the predom-

inant carbide in carbon steels [155]. Among others, it allows to control the mechanical properties

of steels by mixing the hard and brittle cementite phase with the soft and ductile phase of α−Fe

(ferrite). Fe3Chas also attracted signiﬁcant attention in geophysics due to a recent suggestion

that it can be contained in the Earth core [156, 157]. An experimental evaluation of the cementite

properties is, however, non-trivial since at ambient pressure it is a metastable phase up to ∼1000 K

[158], and experimental samples are typically contaminated with ferrite and bulk carbon phases

[159]. This leads to a sizable scatter of the measured data between diﬀerent experimental studies

when considering thermodynamic properties of cementite, such as the thermal volume expansion,

the heat capacity, and the heat of formation. Consequently, an unambiguous parametrization of

thermodynamic databases within calphad concept becomes cumbersome. In this situation, sup-

plementing experimental datasets, which cannot be easily characterized with existing experimental

setups, with results of the ﬁrst-principles calculations becomes particularly important. So far, how-

ever, published ﬁrst-principles studies are focused on the T=0 K properties only [160–164], and do

not take temperature eﬀects into account.

We employ in the following the nearest-neighbor approach to study the thermodynamic prop-

erties of cementite. Afterwards we compare our results with a recent thermodynamic assessment

of Fe3C, with an emphasize on the accurate description of the free energy. In addition we discuss

the thermodynamic stability of cementite versus the spinodal decomposition into bulk carbon and

iron phases by calculating the T-dependent heat of formation of Fe3Cwith a particular focus on

the magnetic contribution.

3.5.1 Thermodynamic Properties of Fe3C

To compute the free energy Fwe apply the adiabatic approximation as we did to compute the

thermodynamic properties of bcc and fcc Fe, Co, and Ni. The convergence parameters are chosen

such that they ensure a convergence of the ﬁnal contributions to Fto less than 1meV per atom.27

In contrast to experimental data for quantities such as the heat capacity, the critical temperature

of Fe3Cis experimentally well known [165]. These preconditions for Fe3Callow us to employ the

nearest-neighbor approach, i.e., it is not necessary to calculate the exchange coeﬃcients for the

real system. However, this requires additional approximations for the considered system: Since the

27The convergence of the calculated vibrational contribution is ensured by using a supercell consisting of 2x2x2

Fe3Cunit cells (128 atoms), a plane-wave energy cutoff Ecut = 400 eV, ∼13000 k-points per atom in homogeneous

Monkhorst-Pack sampling of the Brillouin zone. For the electronic contribution the k-point mesh with ∼40000

k-points per atom was used.

90 3.5 Impact of Magnetism on the Phase Stability of Cementite

Figure 3.32: Dependence of the simulated Fe3Cmagnetization on temperature as compared to the

experimental data from Ref. [159]. The experimental data are acquired at a magnetic ﬁeld of 5 kOe

for a mixture of Fe3Cnanoparticles of size 40±10 nm with residual carbon and about 2 wt% α-Fe.

approach is developed for systems containing single magnetic species only, we neglect the magnetic

contribution of C atoms. This is justiﬁed since their local magnetic moments are much smaller

(≤0.01µB) compared to that of Fe species (≈2µB). The eﬀective nearest-neighbor Heisenberg

model contains therefore Fe-type atoms only. Further, we assume that the Fe atoms are magnetically

indistinguishable (carry the same spin).28

To employ the nearest-neighbor approach three parameters are needed: the positions of the

eﬀective spins, the eﬀective nearest-neighbor exchange parameter J, and the eﬀective spin quantum

number. The cell geometry and the atomic positions are taken from the DFT calculations and

a model containing a 6x6x6 supercell (3456 atoms) is build up. The eﬀective nearest-neighbor

exchange parameter Jis determined by ensuring the correct experimental value of Texp

C= 483 K

[165]. From the DFT calculations the equilibrium local magnetic moments of Fe atoms are computed

as M0

Fe ≈1.9µBproviding an eﬀective spin ˜

S≈0.95 of the Fe3Csystem. According to (3.49) the

QMC simulations are performed using the spin quantum number S= 1 and S= 0.5and an

interpolation afterwards.

The theoretical temperature dependence of the magnetization in comparison with the corre-

sponding experimental data is shown in Fig. 3.32. Please note the diﬀerent scales on the right and

left axis. A quantitative comparison of the theoretical and the experimental data is hindered by

the presence of a mixture of Fe3Cand carbon in the experiments. Figure 3.32 shows, however, that

qualitatively the magnetization in Fe3Cis well reproduced consistent with the results for the pure

ferromagnets Fe, Co, and Ni (see also Fig. 3.31).

28Strictly speaking, there are two types of Fe-atoms in the unit cell which differ in the local magnetic moment by

less than 0.1µB.

3.5 Impact of Magnetism on the Phase Stability of Cementite 91

In Fig. 3.33 the calculated speciﬁc heat capacity, the corresponding available experimental data

and the results of diﬀerent empirical assessments are presented. The theoretical heat capacity is

further split into the vibronic, electronic and magnetic contributions.29 As for the other studied

systems (Fe, Co, Ni), we ﬁnd that the dominating contribution originates from the vibronic ex-

citations showing a typical Debye behavior at low temperatures (∼T3) and a linear increase at

high temperatures. We also observe a signiﬁcant contribution due to electronic excitations which

is due to the high density of electronic states close to the Fermi level in Fe3C[160]. The magnetic

contribution of the heat capacity is largest at TCdue to the break down of the long-range magnetic

order. Below 200 K (close to ferromagnetic saturation) as well as high temperatures (paramagnetic

state), it smoothly vanishes to zero. We can conclude that for temperatures larger than ≈1000 K,

the electronic and vibronic contributions are responsible for the linear increase in the theoretical

CP.

Comparing our results with the available measurements of the heat capacity at constant pressure

(Fig. 3.33) below the Curie temperature we ﬁnd that our predicted theoretical curve falls within

the scatter interval between diﬀerent measurements30. Above TC, our curve is located practically

on-top of the measurements of Naeser [166], and reproduces the presence of a local minimum at

T∼600 K and the increase at higher temperatures. These results are also consistent with the

empirical calphad assessments by Chipmann [158] and Darken et al. [167]. There is a strong

deviation between the above described behavior and the recent measurements by Umemoto and

coworkers where CPsteeply increases with the temperature above TC[165]. The latest empirical

assessment of CPby Hallstedt et al. [153], where the heat capacity above 600 K is assumed to be

constant and which is on-top of the experiments by Seltz et al. [168] below the Curie temperature,

is within the 10 J mol−1K−1from ﬁrst-principles results up to T=1500 K.

In Fig. 3.34 we compare the calculated free energy FFe3C(T)with the calphad data of Hallst-

edt et al. [153].31 Due to the diﬀerent description of CPfor the computed data and the empirical

assessment (Fig. 3.33), both datasets slightly and linearly deviate above ≈200 K. Considering the

magnetic free energy of Fe3C, substantial diﬀerences in the balance exist compared to our results

of pure bcc Fe. In the considered temperature window (up to 1500 K) the magnetic contribution

in pure iron is weaker than the one of Fe3C. Although both systems exhibit a similar maximum

magnetic entropy Smax

mag (due to the similar magnetic moments), Smax

mag in Fe3Cis reached at lower

temperatures due to the signiﬁcantly lower TCwhen compared to that of bcc iron. In general we

obtain an excellent agreement of the theoretical and the experimental free energies with deviations

below 15 meV/atom up to T=1500 K.

29Data for vibronic and electronic contributions for Fe3Cand diamond are taken from Alexey Dick.

30The large discrepancy at T > TCof the measured heat capacity in Ref. [165] is unclear and is likely to do with

the quality of the samples and/or their preparation procedure.

31We note, that due to the large scatter in the experimental data that underlies the Hallstedt et al. assessment,

these results which we use as benchmark are subject to an uncertainty that should be reduced in the future by

improved experimental setups.

92 3.5 Impact of Magnetism on the Phase Stability of Cementite

0500 1000 1500

Temperature (K)

100

150

Heat Capacity (J mol-1 K-1)

vib+el+mag

vib+el

vib

Chipmann 1972

Hallstedt 2010

Darken 1951

Gustafson 1985

Seltz 1940

Naeser 1934

Mazur 1969

Andes 1936

Umemoto 2001

Figure 3.33: Calculated heat capacity of cementite in comparison with available experimental data

(open black symbols) and thermodynamic assessments (ﬁlled green symbols). The diﬀerent the-

oretical contributions are indicated by shaded areas (vibronic grey, electronic orange, magnetic

blue).

200 400 600 800 1000 1200 1400

Temperature (K)

-800 -800

-600 -600

-400 -400

-200 -200

0 0

Free Energy (meV/atom)

400 800 1200

100

Fe3C

Figure 3.34: Calculated free energy for Fe3Cin comparison with the available thermodynamic

assessment from Ref. [153]. The ﬁrst-principles (blue line) and the calphad curves (open black

squares) are aligned at T= 200 K. The calculated electronic and magnetic contributions to the free

energy are indicated by shaded orange and blue areas correspondingly. Inset: Diﬀerence between

ab initio and calphad result.

3.5 Impact of Magnetism on the Phase Stability of Cementite 93

0500 1000 1500

Temperature (K)

-100 -100

0 0

100 100

200 200

∆Ff (meV/f.u.)

Fe3C

align

Figure 3.35: Calculated free energy of formation of cementite from α-Fe and diamond as a function

of temperature. For comparison the calphad assessment as given in [153] is shown. For an easier

comparison of the temperature dependence, the present theoretical curves aligned to the empirical

value at T= 0 K are also shown. The calculated electronic and magnetic contributions to the free

energy are shown as shaded orange and blue areas correspondingly.

3.5.2 Phase stability of Fe3C

In order to analyze the stability of Fe3Cversus bulk iron and carbon phases, and to compare it

with available calphad assessments, we consider the free energy of formation of Fe3C:

∆FFe3C

f(T) = FFe3C(T)−(3FFe(T) + FC(T)) .(3.50)

Here, FFe(T)is the free energy of bcc Fe taken from Sec. 3.4. FC(T)in the above equation de-

notes the free energy of bulk carbon. Typically, in experimental and empirical models only the

thermodynamically stable phase of graphite is considered [153]. It is well known, however, that

conventional exchange-correlation functionals used in DFT calculations fail to accurately describe

the weak correlation forces of van de Waals type present in graphite. Therefore, instead of graphite

only the diamond phase is considered here. The results obtained using graphite as reference reveal,

however, qualitatively the same ﬁndings. For a detailed discussion we refer to [154].

Employing (3.50) we analyze now the phase stability for Fe3C. A positive formation energy

indicates that Fe3Cis thermodynamically unstable versus spinodal decomposition into bulk Fe and

C phases. The resulting calculated heat of formation ∆FFe3C

f(T)is shown in Fig. 3.35. At T= 0 K

the predicted energy of formation is signiﬁcantly lower (by ≈250 meV per formula unit) than

the empirical assessment. The calculations predict Fe3Cto be thermodynamically stable at all

temperatures.32

In order to separate the accuracy of describing ground-state properties from the ﬁnite tempera-

ture eﬀects, we align the theoretical ∆Hfwith the empirical curves at T=0 K. The aligned curves

32In [154] it is shown that the results are not improved by taking all-electron methods into account. The failure

can be traced back to the used xc-functional which performs particular bad for the cohesive energies of Fe and Fe3C.

94 3.5 Impact of Magnetism on the Phase Stability of Cementite

are also shown in Fig. 3.35. We separate the magnetic and electronic contributions in order to

illustrate their impact on the stable-unstable transition temperature ∆TFe3C.

Considering the aligned ∆Hfand taking only vibronic contributions into account, we observe

a stabilization of Fe3Cat about 1300 K which is signiﬁcantly larger than the empirical predicted

∆Temp

Fe3C= 860 K. The electronic excitations are almost vanishing for ∆Hfand have practically

no inﬂuence on ∆TFe3C. This is due to the almost compensating contributions for α-Fe and Fe3C,

whereas the diamond phase has a negligible electronic contribution.33 In contrast to the electronic

excitations, a strong impact of the magnetic entropy contribution can be expected, since the critical

temperatures of α-Fe and Fe3Cstrongly diﬀer. Adding the magnetic excitations to the vibronic

and electronic heat of formation leads to a signiﬁcant decrease of the transition Tby more than

∼350 K and we obtain ultimately a theoretical predicted temperature ∆T theo

Fe3C= 950 K at which

the cementite becomes thermodynamically stable over the α-Fe and bulk carbon.

Conclusion

Using our developed combined ab initio approach employing the nearest-neighbor ansatz (Sec. 3.4)

it is possible to calculate the thermodynamic properties of Fe3Cwith a high degree of accuracy.

The quantitative deviations from the available experimental measurements are found to be smaller

than 0.1% and 10 mol−1K−1for temperatures up to 1500 K for the atomic volume and the heat

capacity, respectively. Noteworthy, that the magnitude of these deviations is comparable with the

one between diﬀerent experimental data themselves. Comparing the results for the free energy

with the recent empirical assessment shows that up to T=1500 K these energies diﬀer by less than

10 meV/atom over the whole considered temperature range.

To investigate the impact of magnetic contributions on the thermodynamic stability of Fe3C

against decomposition into bulk Fe and C phases the free energy of formation is calculated. Align-

ing the heat of formation at T=0 K to the corresponding experimental value we show that the

thermal excitations are accurately described within the present theory, providing a prediction of

the transition temperature ∆TFe3Conly ∼100 K above the corresponding empirical value. Such an

agreement is remarkable when taking into account the uncertainty which persists in the empirical

assessment.

The calculations further indicate that the largest deviation from the empirical assessment occur

for the zero-temperature heat of formation (in the order of 200 meV per formula unit). On the

other hand the T-dependence of the theoretical heat of formation is described with a high degree

of accuracy. A similar result was found in Ref. [80] by investigating a large set of non-magnetic

materials employing similar methods for vibronic and electronic excitations. It shows that the

accuracy of the existing exchange-correlation functionals for prediction of the heat of formation

at T=0 K should be improved in the ﬁrst place. For complex ferromagnetic systems, when the

temperature-dependent properties and trends are of importance, the existing approach can already

be used as complementary (or substituting) method to existing empirical and experimental schemes.

An incorporation into the calphad approach is planned.34

33Since diamond is an insulator, the electronic entropy contribution is much smaller than for α-Fe and Fe3C.

34Private discussions with Dr. B. Hallstedt, RWTH Aachen.

Chapter 4

Summary and Outlook

In the present work a major challenge in computationally guided materials design has been addressed

and to a large part solved: the ab initio prediction of finite-temperature properties of magnetic

materials. The accurate description of the magnetic contribution makes an interplay of various

theoretical treatments mandatory. For this purpose, a careful and profound analysis of diﬀerent

solution techniques of the Heisenberg model has been performed. The developed methods have been

evaluated for a well studied magnetic material, namely bcc iron, but have been also extended to

more complex magnetic materials like cementite. The presented techniques build up a ﬁrm basis to

compute the magnetic entropy contribution from ﬁrst principles and to address present questions

in materials design.

For iron-based materials, the Heisenberg model has proven to yield a numerically eﬃcient yet

suﬃciently accurate description of magnetic excitations. It provides an elegant way to couple ground

state ab initio calculations, which are used to determine the exchange integrals, with concepts of

many-body theory to describe the full temperature dependence of excited magnetic states and

to predict critical temperatures (Sec. 3.1). Beyond eﬃciency, the big advantage of the analytic

solutions is that they allow to decouple the diﬀerent physical mechanisms contributing to the Curie

temperature TC. Utilizing this, we revealed that the dependence of TCon pressure results from

an interplay between two competing pressure-dependent mechanisms: the increase of the magnetic

coupling between the localized magnetic moments and the simultaneous decrease of their magnitude

under the increase of the external pressure. While the ﬁrst mechanism tends to increase the Curie

temperature, the second mechanism acts against it. In the particular case of bcc iron they nearly

completely compensate each other. It is expected that such a nearly perfect compensation for

bcc iron is not an universal feature of all ferromagnetic materials, but is, likely, an exception.

Nevertheless, the proposed picture of two competing magnetic mechanisms can be used for an

intuitive interpretation of the physical reasoning underlying diﬀerent behaviors of the pressure-

dependent Curie temperature in other ferromagnetic materials.

In order to use the combination of DFT and the Heisenberg model for realistic materials, a

number of analytical and numerical approaches have been implemented and carefully checked with

respect to numerical eﬃciency and predictive power. In this context, the speciﬁc heat capacity

Cphas been identiﬁed as a critical quantity, which provides also a close link to experiments and

in addition to empirical calphad simulations. An analytical approach to the free energy of the

Heisenberg spin system (Sec. 3.2.1), which is based on the random phase approximation (RPA), in

combination with vibronic (quasiharmonic approximation) and electronic (ﬁnite temperature DFT)

contributions, has allowed us to obtain a remarkably accurate description of Cpbelow the critical

temperature. For a proper description of the region above the critical temperature, however, short

range correlations are critical and needed to be taken into account. Therefore, various Monte Carlo

(MC) schemes have been explored. Classical MC calculations that neglect quantum eﬀects, yield

signiﬁcant short comings below the critical temperature (Sec. 3.2.2). This implies an important

(and surprising) conclusion: Quantum eﬀects are important even at temperatures up to TC, which

is for bcc iron already about 1000 K!

Therefore, a rescaling scheme has been derived from an extensive set of quantum and classical

Heisenberg model systems, which provides the accuracy of quantum MC calculations, yet keeping

the full ﬂexibility of the classical MC approach. The proposed analytic rescaling function fulﬁlls

exact high-temperature and high-spin limits and transforms at low temperatures into an analytic

mean ﬁeld solution. Employing the rescaling function signiﬁcantly improves the deﬁciencies of the

original classical MC data at low temperatures even for realistic systems such as bcc iron. Due

to its straightforward applicability, it could be used whenever classical MC calculations are carried

out and allows a good estimation of the expected quantum eﬀects. In contrast to the analytic RPA

as well as the classical MC results, the rescaled MC approach has proven to allow a complete and

consistently good description of the magnetic entropy at temperatures below as well as above TC.

Applying the methods to antiferromagnetic fcc iron, classical MC calculations and the analytic

RPA method have provided similar predictions for the N´eel temperature. Combining the derived

magnetic entropy contribution with electronic and vibronic excitations, a free energy of fcc iron

is obtained similar to the results provided by the calphad approach. Due to the possibility of

deriving free energies of single phases (here, bcc and fcc iron), the investigation of phase transitions

has become feasible. The impact of diﬀerent entropy contributions on the phase transitions of iron

has been analyzed. It was thus possible to settle down a long standing debate regarding the impact

of vibronic contributions. They are (together with electronic contributions) equally important for

the stabilization of the fcc crystal phase as magnetic contributions. Based on the magnetic entropy

contributions, an intuitive explanation of the phase transitions could be provided: on the one

hand, the austenite is stabilized at intermediate temperatures due to the lower critical temperature

compared to ferrite. At high temperatures, on the other hand, the backward transition from the

austenite to the paramagnetic δ-iron is driven by the larger magnetic entropy of bcc iron due to the

larger spin moment compared to the fcc phase.

Analyzing the classical and quantum results for various model systems of the previous compar-

ison in more detail (Sec. 3.2.2), a universal behavior of the thermodynamic properties with respect

to the kind of the interaction (long- vs. short-ranged, ferro- vs. antiferromagnetic) and the lattice

structure has been discovered for the Heisenberg model. This universality has motivated the deﬁ-

nition of an eﬀective nearest-neighbor Heisenberg model with just one interaction parameter. The

latter is, e.g., related to TC, which can be either determined by ab initio calculations (Sec. 3.1) or

taken from experiment. Treating this eﬀective model with quantum MC has provided for bcc Fe

an excellent agreement with experimental Cpdata and free energies (Sec. 3.4). The application of

the same approach to the 3d ferromagnets cobalt and nickel has demonstrated its excellent trans-

ferability to other magnetic systems containing a single magnetic species. Therefore, the developed

methodology has been applied to more complex systems, such as, e.g., cementite (Fe3C), which is

a highly relevant phase for steel design (Sec. 3.5). Due to a large scatter in the experimental Cp(T)

an unambiguous thermodynamic assessment (e.g., calphad) of this system is challenging. As it

has been demonstrated in Sec. 3.5, the method has not only provided a clear prediction of Cp(T)

particularly in the high temperature regime, it also has provided a hitherto not achieved insight

into the impact of magnetic contributions to the phase stability of Fe3C.

In conclusion, the set of ab initio based methodologies developed in this work has allowed a

reliable and parameter-free prediction of the magnetic contribution to the free energy of materials, a

systematic separation from other physically relevant contributions to the thermodynamic properties,

and can be applied to a wide range of materials.

The developed methods in this study open up a huge number of potential investigations which

can now be tackled. Still, several aspects of the developed methodologies can be improved or

extended. In order to give a ﬂavor of the possible next steps we summarize in the following the

most important issues:

Multi-magnetic species systems:

In order to treat systems containing more than one kind of magnetic atom (such as FeMn-compounds),

the developed theories have to be extended. Concerning the analytic approaches, a corresponding

theory for obtaining TCbased on the RPA scheme already exists [169]. To extend the developed

theory to describe thermodynamic properties such as free energy and heat capacity, principally

a similar approach as in Sec. 3.2.1 could be used. Such an extension should be straightforward,

although likely requiring a considerably longer analytic derivation.

Regarding the rescaled Monte Carlo approach, the model calculations should be extended to the

case of multi-spin lattices in the ﬁrst place. After a detailed model comparison between the classical

and quantum results, a modiﬁed rescaling ansatz may be derived. In order to apply the nearest-

neighbor approach to a multi-magnetic-species system, most likely an additional eﬀective exchange

coeﬃcient is needed. This could be problematic regarding that the mapping of the long-ranged onto

the eﬀective nn-model is currently done by just one conserving parameter i.e. TC. As a ﬁrst order

approximation one could assume that the exchange between diﬀerent kind of magnetic atoms as

well as the diﬀerent spin moments can be described by one eﬀective (averaged) exchange parameter

and one eﬀective spin quantum number, respectively. Under this assumption, both methods can be

straightforwardly applied without further methodological extensions, and are equally well suited to

treat such systems.

Noteworthy, that regarding the DFT input, the computation of exchange-parameters for multi-

component-systems is already well established [169].

Inclusion of longitudinal spin fluctuations:

The impact of longitudinal spin ﬂuctuations, i.e., the change of the magnitude of the local magnetic

moments, is not clariﬁed. One should distinguish between the explicit and implicit temperature

dependence of the local moments and exchange coeﬃcients. The ﬁrst one is determined by explicit

temperature dependent excitation processes such as, e.g., single-spinﬂip processes of electrons or

ﬂip processes of electrons between diﬀerent bands. Their inclusion into the proposed Heisenberg

model ansatz is conceptually diﬃcult and requires most likely band-model based theories (Hubbard

model) [39, 170] and an additional mapping onto eﬀective exchange coeﬃcients. The extension

to multi-magnetic material systems is expected to be rather diﬃcult. A straightforward ﬁrst step

could be, to employ the already developed methods using diﬀerent electronic temperatures (Fermi

smearing) within the DFT formalism.1

The implicit dependence of the local moments and eﬀective exchange coeﬃcients, i.e., their

change due to volume expansion or magnetic disorder, can be simulated by DFT calculations di-

rectly. Diﬀerent methodologies are nowadays available (e.g., the so-called disordered local moment

method [171]) to simulate the paramagnetic state at T= 0 K. Based on these calculations, eﬀective

and implicitly temperature dependent (via the coupling of the local moment to the magnetic state)

exchange coeﬃcients can be computed. A recent incorporation of the latter into the classical Heisen-

berg model reveals promising results for Fe and Ni [48]. The inclusion of these implicit LSF could

be straightforwardly done in the rescaled MC approach (via the spin-dependent rescaling function).

Based on the ﬁndings in [48], it is expected that incorporating the LSF into the methodologies

might be potentially important for systems showing small band splitting like Ni, but should not

have signiﬁcant impact on results of strong-moment systems such as Fe.

Coupling terms:

In the present study, the diﬀerent types of excitations are treated in an adiabatic approximation.

In addition, the electronic and vibronic contributions were performed for the magnetic ground-state

conﬁguration, e.g., for the ferromagnetic state of bcc Fe. This implicitly ignores the dependence of

vibronic excitations on the magnetic state.

To account for this coupling (phonons in diﬀerent magnetic states etc.), a very promising route

was very recently proposed by Shang and Liu [172]. They propose to perform total energy and

phonon calculations for a number of diﬀerent magnetic states and to average afterwards following

a Boltzmann statistic. Physically, this corresponds to the situation that the system at a given tem-

perature is determined by the occupation of diﬀerent (magnetic) micro states. Naturally, coupling

terms between electronic-magnetic and vibronic-magnetic contributions are in such an approach

implicitly included. Whereas the recent approach ignores the explicit magnetic contributions [172],

a complete approach can be immediately employed using our developed methods.

Impurity dependent trends of the critical temperature:

Since TCis taking an important role in the developed methods, it would be interesting to investigate

how it is aﬀected by diﬀerent impurities such as, e.g., carbon. This would thus provide the next

step towards a complete theoretical prediction of the Fe-C phase diagram. For the study of the

thermodynamic properties of cementite for example, the ab initio determination of TC, which here

has been taken from experiment, is a natural next step to achieve a fully ﬁrst-principles Fe3C

prediction. This itself would be already an important issue since very little is known about the

dependence of TCof Fe3Con diﬀerent impurities, particularly, considering the fact that the diﬀerent

critical temperatures of ferrite and cementite play a crucial role for the stabilization of Fe3C. Since

the underlying methods have been developed, such a study could be immediately started.

1This could only simulate the different occupations within the spin bands. However, such an ansatz can not

account for temperature induced electronic intra- or inter-band flips.

Pressure-temperature phase diagrams:

Although an absolute quantitative agreement be-

0 2 4 6 8 10 12 14 16 18

Pressure (GPa)

400

800

1200

1600

2000

2400

2800

Temperature (K)

α-iron

γ-iron

ε-iron

tween theory and experiment for total free ener-

gies of pure Fe is not achieved, an investigation

of physical trends is possible such as, e.g., the

dependence of diﬀerent physical mechanisms on

pressure. Regarding the rather complex depen-

dencies for the stability of diﬀerent Fe crystal

phases on pressure, a simple mean-ﬁeld ansatz

could in principle already provide important in-

sights. Employing MF, the required ingredients

for the diﬀerent crystal phases λare the critical

temperatures Tλ

C(V)as well as the local magnetic

moments Mλ(V)(both may depend on volume).

Assuming as a ﬁrst order approximation Tλ

C(V)≈Tλ

C, the MF free energy Fmag

MF (T, V, λ), Sec. 3.2.1,

can be straightforwardly employed and combined with the ground state energies Eλ(V). Preliminary

results are already available and shown in the ﬁgure above. The dashed line indicates the theoretical

ferrite to austenite transition temperatures. Further development in this direction is planned, in

particular to elucidate the impact of the diﬀerent contributions with respect to pressure.

100

Appendix A

Appendix

A.1 Mean Field Approximation for Ferromagnets

The employed mean ﬁeld (MF) theory is derived as, e.g., found in Ref. [56].

A.1.1 The Mean Field Hamiltonian

The mean ﬁeld approximation is realized by the following treatment for operator products AB

AB =AhBi+hAiB−hAihBi+ (A−hAi)(B−hBi)(A.1)

→AhBi+hAiB−hAihBi.(A.2)

Within MF, ﬂuctuations of A and B are neglected, i.e.

(A−hAi)(B−hBi)→0.(A.3)

Thus, the Heisenberg Hamiltonian transforms into

H=−X

JijSiSj(A.4)

=−X

Jij(S+

iS−

j+Sz

iSz

j)(A.5)

−→ −X

Jij(hSz

iiSz

j+Sz

ihSz

ji) + X

JijhSz

iihSz

ji.(A.6)

Due to spin conservation at each lattice site, hS+

ii=hS−

ii ≡ 0, collective spin-wave excitations

(correlations) are supressed

hS+

iS−

jiMF

−→ 0.(A.7)

Assuming homogeneous magnetization,

hSz

ii ≡ hSzi,(A.8)

101

102 A.1 Mean Field Approximation for Ferromagnets

the Hamiltonian reads as

HMF := J0hSzi2(T)−BEX(T)X

i,(A.9)

with an eﬀective (mean ﬁeld) exchange ﬁeld,

BEX(T) := 2J0hSzi(T),(A.10)

which is proportional to the magnetization hSzi.

A.1.2 The Mean Field Magnetization

Within mean ﬁeld, the partition sum Zseparates into one-spin partition sums z1, i.e.

Z=Tre−βH=zN

1C, (A.11)

with

z1=

i=−S

eβBEX(T)Sz

i(A.12)

and the constant

i=−S

e−βJ0hSzi2(T)= (2S+ 1)e−βJ0hSzi2(T),(A.13)

which cancels when calculating observables (see derivation below).

The procedure is equivalent to Weiss theory of ferromagnetism (see e.g. [56]),

hSz

ii/S =TrSz′

ie−βHMF

STre−βHMF (A.14)

=P+S

−SSz′

ie−β(J0hSzi2(T)−BEX(T)PiSz′

SP+S

−Se−β(J0hSzi2(T)−BEX(T)PiSz′

i)(A.15)

=P+S

−SSz′

ieβBEX(T)PiSz′

SP+S

−SeβBEX(T)PiSz′

(A.16)

=BS(2S2J0

kBThSzi

S),(A.17)

where the Brillouin function BS(x)is given by

BS(x) := 2S+ 1

2Scoth(2S+ 1

2Sx)−1

2Scoth( x

2S).(A.18)

A.1.3 Mean Field Solution of f(t, S)for Low Temperatures

Within mean ﬁeld, the internal energy is proportional to the squared magnetization,

UMF(t, S)∼mMF(t, S)2,(A.19)

A.1 Mean Field Approximation for Ferromagnets 103

where the latter is given by

mMF(t, S)∼BS3S

S+ 1

m(t, S)

t,(A.20)

including the Brillouin function

BS(x) = 2S+ 1

2Scoth 2S+ 1

2Sx−1

2Scoth 1

2Sx.(A.21)

For S→ ∞, above function transforms into the Langevin function

BS→∞(x) = L(x) = coth(x)−1

x.(A.22)

For low temperature, i.e. x≫1, one obtains

BS(x≫1) = 1 −exp (−x/S)

S(A.23)

L(x≫1) = 1 −1/x. (A.24)

Combining Eqs. (A.19)-(A.24), one obtains straightforwardly the low temperature expansion for the

heat capacity,

CMF(t, S) = −3

2kB

S+ 1

dm2(t, S)

dt ,(A.25)

CMF(t→0, S)≈kB2tMF

S/t

exp(tMF

S/t)2

(A.26)

CMF(t→0, S → ∞)≈kB

2 1 + 1

p1−4/3t!,(A.27)

with the spin-dependent variable

1/tMF

S≡2/3(S+ 1).(A.28)

Therewith, the ratio of the quantum and classical (S→ ∞) heat capacity in the low temperature

limit t→0,

lim

t→0fMF(t, S) = lim

t→0CMF(t, S)/CMF(t, S → ∞),(A.29)

can be readily computed and we obtain

fMF(t→0, S) = 2tMF

S/t

exp(tMF

S/t)2

.(A.30)

104 A.2 Random Phase Approximation for Ferromagnets

A.2 Random Phase Approximation for Ferromagnets

The magnetization is derived as in Ref. [61] employing the random phase approximation (RPA) and

using a parametrized Green’s function.

We start with the following identity:

hS−

iS+

ii=~2S(S+ 1) + ~hSz

ii−h(Sz

i)2i.(A.31)

The operator product (Sz

i)2on the right site of the equation above requires a Green’s function which

provides the relation between the expectation values of the power set of (Sz

i)n. For this purpose,

we employ a parametrized Green’s function and follow a procedure proposed by Callen [61].

To determine the magnetization we introduce the following retarded Green’s function:

G(a)

ij =hhS+

i;B−(a)

jii =hhS+

i;eaSz

jS−

jii.(A.32)

Thus, we obtain for the equation of motion

EG(a)

ij (E) = ~h[S+

i;B−(a)

j]−i+hh[S+

i,H]−;B−(a)

jiiE.(A.33)

For the solution we need the inhomogeneity

h[S+

i;B−(a)

j]−i=hS+

i, eaSz

iS−

i−iδij =η(a)δij,(A.34)

as well as the commutator [S+

i,H]−. Employing the RPA yields

eiqRi[S+

i,H]−=S+

q,H−

RP A

−→ 2 (J0−Jq)hSzi=: Eq.(A.35)

With Eqs. (A.34) and (A.35) we can write down the equation of motion

~G(a)

q=η(a)

E−Eq

.(A.36)

The spectral density, S(a)

q=−1

πℑG(a)+−

q, can be written as a delta-function:

S(a)

q=~η(a)δ(E−Eq).(A.37)

Using the spectral theorem and Eq. (A.37) we can determine the expectation value

p(q, a) = 1

eiq(Ri−Rj)heaSz

jS−

jS+

ii(A.38)

=η(a)

eβE(q)−1.(A.39)

A.2 Random Phase Approximation for Ferromagnets 105

We sum over all wave vectors qand end up with:

p(a) = 1

p(q, a) = heaSzS−S+i=: η(a)ϕ, (A.40)

where we deﬁne the eﬀective Magnon occupation number

ϕ=1

eβE(q)−1.(A.41)

We obtain the desired Eigenvalue hSziusing a diﬀerential equation

Ω(a) = heaSzi(A.42)

hSzi=d

daΩ(a)a=0

.(A.43)

To derive the relation between p(a),η(a)and awe use the identity (proof by induction)

S+

i,(Sz

i)n−= ((Sz

i−~)n−(Sz

i)n)S+

i(A.44)

which yields

S+

i, eaSz

i−=X

n!anS+

i,(Sz

i)n−(A.45)

n!an((Sz

i−~)n−(Sz

i)n)S+

i(A.46)

=e−a~−1eaSz

iS+

i.(A.47)

We can now write down η(a)

η(a) = hS+

i, eaSz

iS−

i−i(A.48)

=heaSz

iS+

i, S−

i−+S+

i, eaSz

i−S−

ii(A.49)

= 2~heaSz

iSz

ii+e−a~−1heaSz

iS+

iS−

ii.(A.50)

Using the identity

~Sz

i=~2S(S+ 1) −(Sz

i)2−S−

iS+

i(A.51)

106 A.2 Random Phase Approximation for Ferromagnets

it follows

η(a) = e−a~−1heaSz

ii~2S(S+ 1) +

+heaSz

iSz

ii2~+~e−a~−~−

−e−a~−1heaSz

i(Sz

i)2i(A.52)

=~2S(S+ 1) e−a~−1Ω(a) + (A.53)

+~e−a~+ 1d

daΩ(a)−e−a~−1d2

da2Ω(a).(A.54)

With (A.51) its possible to write (A.40) as a function of Ω(a):

p(a) = heaSzS−S+i(A.55)

=~2S(S+ 1)Ω(a)−~d

daΩ(a)−d2

da2Ω(a) = η(a)ϕ(A.56)

Using (A.54) we can set up a diﬀerential equation for Ω(a):

h−e−a~−1ϕ+ 1id2

da2Ω(a) + h~e−a~+ 1ϕ+~id

daΩ(a)+

+h~2S(S+ 1) e−a~−1ϕ−~2S(S+ 1)iΩ(a) = 0 (A.57)

⇒d2

da2Ω(a)−~e−a~+ 1ϕ+ 1

(e−a~−1) ϕ−1−~2S(S+ 1)Ω(a) = 0.(A.58)

The solution of this diﬀerential equation of second order requires two limiting conditions. One

follows directly from the deﬁnition (A.42)

Ω(0) = 1.(A.59)

With Sz=−~S, . . . , ~Swe obtain the second condition

ms=−S

(Sz−~mS) = 0 ⇒

ms=−Sd

da −~mSΩ(a)a=0

= 0.(A.60)

Having the limiting conditions (A.59) and (A.60) we ﬁnally derive the solution (see e.g. [56]):

Ω(a) = ϕ2S+1e−a~S−(1 + ϕ)2S+1e~(S+1)a

(ϕ2S+1 −(1 + ϕ)2S+1)((1 + ϕ)e~a−ϕ).(A.61)

We now obtain the desired expectation value hSzi

hSzi=d

daΩ(a)a=0

=(1 + ϕ)2S+1(S−ϕ) + ϕ2S+1(S+ 1 + ϕ)

(1 + ϕ)2S+1 −ϕ2S+1 ,(A.62)

A.3 Random Phase Approximation for Noncollinear Magnets 107

with

ϕ=1

eβE(q)−1.(A.63)

A.3 Random Phase Approximation for Noncollinear Magnets

A.3.1 Magnetization and Néel temperature

In analogy to the ferromagnetic case we set up the equation of motion for the Green’s function

G+−

ij , but now in the locally rotated coordinate system Σ′. Therefore, the following commutators

have to be solved for the higher order Green’s functions:

S+′

m,H′−=−~X

j{Kjm(Sz′

mSx′

j+Sx′

jSz′

m) + iLjm(Sz′

mSy′

j+Sy′

jSz′

−Mjm(S+′

mSz′

j+Sz′

jS+′

m)−2X

Njm(−Sx′

jS+′

m+Sz′

mSz′

j)}.(A.64)

Applying the random phase approximation (RPA)

S+′

m,H′−

RPA

−→ −2~hSziX

jnKjmSx′

j+iLjmSy′

j−MjmS+′

mo(A.65)

=−~hSziX

jnKjm(S+′

j+S−′

j) + Ljm(S+′

j−S−′

j)−2MjmS+′

mo.(A.66)

Analogously one obtains:

S−′

m,H′−

RPA

−→ −~hSziX

jn−Kjm(S+′

j+S−′

j) + Ljm(S+′

j−S−′

j) + 2MjmS−′

mo.(A.67)

Equation of Motion: Writing down the equation of motion for Green’s functions G+−

ij and

G−+

ij , and employing a Fourier transformation, e.g. 1

NPieiqRiSi+′=S+′

q, one ends up with the

following coupled system of equations:

~ G+−

qG−−

G++

qG−+

q!= 2hSz′i0

0−2hSz′i!EI −Ω−1,(A.68)

where the excitation energies of the system are determined by the eigenvalues of the matrix

Ω(q) = hSz′i 2M−Kq−LqKq−Lq

−(Kq−Lq)−(2M−Kq−Lq)!= A B

−B−A!,(A.69)

and explicitly given by

Eq±=±2~hSz′iq(M−Kq) (M−Lq)≡E±.(A.70)

108 A.3 Random Phase Approximation for Noncollinear Magnets

Here, the exchange coeﬃcients are given by

Kq=Mq=1

2(Jq+Q+Jq−Q)(A.71)

and

Lq=Jq.(A.72)

Magnetization and critical temperature:

The magnon excitation energies and corresponding weights directly determine the eﬀective magnon

occupation number ϕ[56], given by

ϕ(T) = 1

χq+

eβEq+−1+χq−

eβEq−−1,(A.73)

where the ﬁrst term describes magnon creation and the second magnon annihilation, respectively.

The weights χq+are determined by the eigenvectors of (A.69) and given by

χq±=1

21±a

b,(A.74)

with

a:= M−1/2Kq−1/2Lq(A.75)

b:= q(M−Kq)(M−Lq).(A.76)

Eq. (A.73) can be rewritten as

ϕ(T) = 1

qa

bcoth [1/2βEq]−1,(A.77)

which turns out to be faster and more eﬃcient for numerical evaluation. The magnetization for the

general spin system (S≥1/2) is given by

hSz′i=(1 + ϕ)2S+1(S−ϕ) + ϕ2S+1(S+ 1 + ϕ)

(1 + ϕ)2S+1 −ϕ2S+1 ,(A.78)

with ϕas deﬁned in (A.73).

For S= 1/2above equation simpliﬁes to

hSz′i=S

1 + 2ϕ.(A.79)

A.3 Random Phase Approximation for Noncollinear Magnets 109

Critical temperature:

For T→Tcrit, the magnetization vanishes and ϕ→ ∞. Hence we can expand the right-hand site

of Eq. (A.78) in terms of O(ϕ−n)and obtain:

hSzi=1

3S(S+ 1)ϕ−1+O(ϕ−2)(A.80)

hSzi→0

−→ 1

3S(S+ 1) X

β4hSzi1

M−Kq

M−Lq!−1

(A.81)

⇒kBTcrit =4

3S(S+ 1) X

q1

M−Kq

M−Lq!−1

(A.82)

3 X

q"1

M−˜

Lq#!−1

,(A.83)

which corresponds to the Néel temperature TNemployed in Ref. [140, 141]. In the last step we have

absorbed the factor S(S+ 1) into the exchange constants.

A.3.2 Internal energy

To derive the internal energy we use, similar as for the ferromagnetic case, exact S= 1/2relations.

In the ﬁrst step we rewrite the Hamiltonian (2.70) in terms of correlation functions, which we access

by the Green’s function matrix Eq. (A.68).

From Eq. (A.68) we obtain:

~ G+−

qG−−

G++

qG−+

q! E−A−B

+B E +A!= 2hSz′i0

0−2hSz′i!(A.84)

Yielding:

G+−

q=E+A

E2−(A2−B2)2hSzi=E+A

(E−E+)(E+E+)2hSzi(A.85)

= 2hSzi1/2[1 + A/E+]

E−E+

+1/2[1 −A/E+]

E+E+(A.86)

From the matrix solution (e.g. ﬁrst row times second column) we obtain:

G+−

qB= (E+A)G−−

q(A.87)

and with (A.85):

G−−

q= 2hSziB

2E+−1

E−E+

E+E+(A.88)

and

G++

q=−2hSziB

2E+−1

E−E+

E+E+.(A.89)

110 A.3 Random Phase Approximation for Noncollinear Magnets

The internal energy

In order to rewrite the internal energy in terms of the correlation functions hSi±Sj±i, we replace all

operators Sik (k=x, y, z) in the Hamiltonian Eq. (3.29),

H′=−X

ij KijS′

ixS′

jx +LijS′

iyS′

jy +MijS′

izS′

jz+ 2 X

NijS′

ixS′

jz,(A.90)

by creation and annihilation operators:

H=−X

ij {1/4(Kij −Jij)S+′

iS+′

j+ 1/4(Kij −Jij)S−′

iS−′

j+ 1/2(Kij +Jij)S+′

iS−′

+KijSz′

iSz′

j−NijS+′

iSz′

j−NijS−′

iSz′

j}.(A.91)

Please be aware of the deﬁnitions of the N, K, J, i.e.

Jij =Jji Kij = cos(φi−φj)Jij =Kji Nij = sin(φi−φj)Jij =−Nji.(A.92)

For S= 1/2, the longitudinal correlation function can be expressed in terms of the transversal

correlation function:

S−′

iS+′

i,H−=−X

m{1/2(Kim −Jim)S−′

mS−′

i+ 1/2(Kim +Jim)S+′

mS−′

i−

2KimSz′

m(S−Sz′

i) + NmiS+′

m(S−Sz′

i) + NmiS−′

m(S−Sz′

i)−NimSz′

mS−′

i}(A.93)

and

−X

KimhSz′

mSz′

ii=1

ihS−′

iS+′

i,H−i+X

im {1/4(Kim −Jim)hS−′

mS−′

ii+

+ 1/4(Kim +Jim)hS+′

mS−′

ii−KimhSz′

miS−1/2NmihS+′

mSz′

ii−NmihS−′

mSz′

ii}.(A.94)

Inserting above relation into U=hHi, only the term 1/2Pim NmihS+′

mSz′

iiremains unknown. We

derive the missing correlation function using

S+′

iS+′

i,H−=−X

m{−1/2(Kim −Jim)S−′

mS+′

i−1/2(Kim +Jim)S+′

mS+′

i+NimSz′

mS+′

i}.(A.95)

and obtain (please be aware of Nij =−Nji)

1/2X

NmihS+′

mSz′

ii=−1/2X

ihS+′

iS+′

i,H−i+

im {1/4(Kim −Jim)hS−′

mS+′

ii+ 1/4(Kim +Jim)hS+′

mS+′

ii}.(A.96)

A.3 Random Phase Approximation for Noncollinear Magnets 111

Now we are able to express the internal energy purely in terms of the correlation functions hS±′S±′i

and hS±′S∓′i:

U=hHi =X

ij −1/4(Kij −Jij)hS+′

iS+′

ji−1/4(Kij −Jij)hS−′

iS−′

ji−

−1/2(Kij +Jij)hS+′

iS−′

ji−KijhSz′

iSz′

ji+NijhS+′

iSz′

ji+NijhS−′

iSz′

(A.94)

z}|{

ihS−′

iS+′

i,H−i+X

ij {−1/4(Kij +Jij)hS+′

jS−′

ii−KijhSz′

jiS+

ij −1/4(Kij −Jij)hS+′

iS+′

ji+ 1/2Nij hS+′

iSz′

(A.96)

z}|{

ihS−′

iS+′

i,H−i+X

ij {−1/4(Kij +Jij)hS+′

jS−′

ii−KijhSz′

jiS+

ij −1/4(Kij −Jij)hS+′

iS+′

ji−1/2X

ihS+′

iS+′

i,H−i+

ij {1/4(Kij −Jij)hS−′

jS+′

ii+ 1/4(Kij +Jij)hS+′

jS+′

ii}

ihS−′

iS+′

i,H−i− 1

ihS+′

iS+′

i,H−i+X

ij −KijhSz′

jiS+

+ 1/2JijhS+′

iS+′

ji−1/2JijhS+′

iS−′

ji.

With hSz′

ii=hS−S−′

iS+′

ii, and using S= 1/2, we ﬁnally end up with:

U=hHi =1

ihS−′

iS+′

i,H−i− 1

ihS+′

iS+′

i,H−i+1

JijhS+′

iS+′

ji+E0+

(K0δij −Jij)hS−′

jS+′

ii,(A.97)

where the ground-state energy E0is deﬁned as

E0=−K0S2.(A.98)

To compute terms such as 1

2PihS−′

iS+′

i,H−i, we employ that for the non time-dependent Hamil-

ton operator Hthe derivative S+′

i(t)can be written as:

i∂

∂tS+′

i(t) = S+′

i,H−(t).(A.99)

112 A.3 Random Phase Approximation for Noncollinear Magnets

Using above, the ﬁrst term on the right-hand site of Eq. (A.97) can be reformulated as (see e.g.

Ref. [28]):

ihS−′

iS+′

i,H−i=i

i ∂

∂t Z+∞

−∞

dE S+−

q(E)

eβE −1exp −iE(t−t′)!t=t′

(A.100)

qZ+∞

−∞

dE ES+−

q(E)

eβE −1.(A.101)

Analogously one obtains

−1

ihS+′

iS+′

i,H−i=−1

qZ+∞

−∞

dE ES++

q(E)

eβE −1(A.102)

and ﬁnally an expression for the internal energy of the noncollinear magnet as:

U=U0+1

qZ+∞

−∞

dE E+ 2S(K−Jq)

eβE −1S+−

q(E) + E−2SJq

eβE −1S++

q(E).(A.103)

The spectral densities are given as

S+−

q(E) = 2hSz′iχ+−

1δ(E−E+) + χ+−

2δ(E+E+),(A.104)

S++

q(E) = 2hSz′iχ++

1δ(E−E+) + χ++

2δ(E+E+),(A.105)

with

χ+−

1,2=1

2 1±K−1/2Kq−1/2Jq

p(K−Kq)(K−Jq)!,(A.106)

χ++

1,2=±1

1/2Kq−1/2Jq

p(K−Kq)(K−Jq),(A.107)

E+= 2hSz′iq(K−Kq)(K−Jq).(A.108)

For the ferromagnetic solution, i.e. Q= (0,0,0), Eq. (A.103) transforms into the known ferromag-

netic solution [28].

Bibliography

[1] Ironware piece unearthed from turkey found to be oldest steel, The Hindu Newspaper, March

26, 2009.

[2] U.S. Geological Survey, 2009, Iron and Steel statistics, in Kelly, T.D., and Matos, G.R., His-

torical statistics for mineral and material commodities in the United States: U.S. Geological

Survey Data Series 140, available online at http://pubs.usgs.gov/ds/2005/140/.

[3] Prof. Dr.-Ing. Dieter Ameling, Präsident Wirtschaftsvereinigung Stahl, Vorsitzender Stahlin-

stitut VDEh, Talk, 2007.

[4] T.Chandra, N.Wanderka, W.Reimers, and M.Ionescu, Materials Science Forum 638, 3134

(2010).

[5] B. Grabowski, PhD thesis, University of Paderborn, 2009.

[6] X. Sha and R. E. Cohen, Phys. Rev. B 73, 104303 (2006).

[7] L. Kaufman, E. V. Clougherty, and R. J. Weiss, Acta Metal. 11, 323 (1963).

[8] H. Hasegawa and D. G. Pettifor, Phys. Rev. Lett. 50, 130 (1983).

[9] U. von Barth and L. Hedin, Journal of Physics C: Solid State Physics 5, 1629 (1972).

[10] K. Capelle, Braz. J. Phys. 36, 1318 (2006).

[11] W. Heisenberg, Zeitschrift für Physik 38, 411 (1926).

[12] S. V. Halilov, A. Y. Perlov, P. M. Oppeneer, and H. Eschrig, Europhys. Lett. 39, 91 (1997).

[13] J. Rusz, L. Bergqvist, J. Kudrnovský, and I. Turek, Phys. Rev. B 73, 214412 (2006).

[14] G. Y. Gao, K. L. Yao, E. Şaşioğlu, L. M. Sandratskii, Z. L. Liu, and J. L. Jiang, Phys. Rev.

B75, 174442 (2007).

[15] M. Pajda, J. Kudrnovsky, I. Turek, V. Drchal, and P. Bruno, Phys. Rev. B 64, 174402 (2001).

[16] F. Körmann, A. Dick, B. Grabowski, B. Hallstedt, T. Hickel, and J. Neugebauer, Phys. Rev.

B78, 033102 (2008).

[17] F. Körmann, A. Dick, T. Hickel, and J. Neugebauer, Phys. Rev. B 79, 184406 (2009).

113

[18] N. Bogolyubov and S. Tyablikov, Soviet. Phys.-Doklady 4(1959).

[19] M. Marsman and J. Hafner, Phys. Rev. B 66, 224409 (2002).

[20] R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford Science

Publications, 1988.

[21] K. Capelle and G. Vignale, Phys. Rev. Lett. 86, 5546 (2001).

[22] H. Eschrig and W. E. Pickett, Solid State Communications 118, 123 (2001).

[23] N. Argaman and G. Makov, Phys. Rev. B 66, 052413 (2002).

[24] N. I. Gidopoulos, Phys. Rev. B 75, 134408 (2007).

[25] W. Kohn, A. Savin, and C. A. Ullrich, International Journal of Quantum Chemistry 101, 20

(2004).

[26] K. Capelle, C. A. Ullrich, and G. Vignale, Phys. Rev. A 76, 012508 (2007).

[27] T. Gál, P. W. Ayers, F. D. Proft, and P. Geerlings, J. Chem. Phys. 131, 154114 (2009).

[28] W. Nolting, Fundamentals of Many-body Physics, Springer, 2009.

[29] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).

[30] L. M. Sandratskii and P. G. Guletskii, Journal of Physics F: Metal Physics 16, L43 (1986).

[31] J. Kübler, Theory of Itinerant Electron Magnetism, Oxford University Press, 2009.

[32] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).

[33] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).

[34] C. S. Wang, B. M. Klein, and H. Krakauer, Phys. Rev. Lett. 54, 1852 (1985).

[35] M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Joannopoulos, Rev. Mod. Phys.

64, 1045 (1992).

[36] D. E. Eastman, F. J. Himpsel, and J. A. Knapp, Phys. Rev. Lett. 44, 95 (1980).

[37] S. Y. Savrasov, Phys. Rev. Lett. 81, 2570 (1998).

[38] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).

[39] E. Şaşıoğlu, A. Schindlmayr, C. Friedrich, F. Freimuth, and S. Blügel, Phys. Rev. B 81,

054434 (2010).

[40] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74, 601 (2002).

[41] W. Nolting, Grundkurs Theoretische Physik 7, Springer-Verlag, Berlin, 2005.

[42] L. Hedin, Phys. Rev. 139, A796 (1965).

114

[43] N. M. Rosengaard and B. Johansson, Phys. Rev. B 55, 14975 (1997).

[44] J. Kübler, J. Phys.: Condens. Matter 18, 9795 (2006).

[45] L. Landau and E. Lifshitz, Phys. Z. Sowietunion 8, 153 (1935).

[46] S. Shallcross, A. E. Kissavos, V. Meded, and A. V. Ruban, Phys. Rev. B 72, 104437 (2005).

[47] A. V. Ruban, S. Shallcross, S. I. Simak, and H. L. Skriver, Phys. Rev. B 70, 125115 (2004).

[48] A. V. Ruban, S. Khmelevskyi, P. Mohn, and B. Johansson, Phys. Rev. B 75, 054402 (2007).

[49] L. M. Sandratskii, Advances In Physics 47, 91 (1998).

[50] J. Hubbard, J. Hubbard, Proc. R. Soc. London, Ser. A 276, 238 (1963) Ser. A 276, 238

(1963).

[51] M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).

[52] J. Kanamori, Prog. Theor. Phys. 30, 275 (1963).

[53] C. Santos and W. Nolting, Phys. Rev. B 65, 144419 (2002).

[54] M. Ležaić, P. Mavropoulos, and S. Blügel, Appl. Phys. Lett. 90, 082504 (2007).

[55] N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).

[56] W. Nolting, Quantentheorie des Magnetismus, Teil 2, Modelle, Teubner Verlag, 1997.

[57] M. Exler, On Classical and Quantum Mechanical Energy Spectra of Finite Heisenberg Spin

Systems, PhD thesis, University Osnabrück, 2006, URN: urn:nbn:de:gbv:700-2006051610.

[58] F. Körmann, A. Dick, T. Hickel., and J. Neugebauer, Phys. Rev. B 81, 134425 (2010).

[59] J. Oitmaa and W. Zheng, Journal of Physics: Condensed Matter 16, 8653 (2004).

[60] J. Thoene, S. Chadov, G. Fecher, C. Felser, and J. Kübler, Journal of Physics D: Applied

Physics 42, 084013 (2009).

[61] H. Callen, Phys. Rev. 130, 890 (1963).

[62] S. H. Liu, Phys. Rev. 139, A1522 (1965).

[63] S. Tjablikov, Methods in the Quantum Theory of Magnetism,, Plenum Press, New York, 1967.

[64] J.-B. Fouet, A. Läuchli, S. Pilgram, R. M. Noack, and F. Mila, Phys. Rev. B 73, 014409

(2006).

[65] T. Fabritius, N. Laﬂorencie, and S. Wessel, Phys. Rev. B 82, 035402 (2010).

[66] N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn, Physics Letters A 238, 253 (1998).

[67] H. G. Evertz, Adv. Phys. 1, 1 (2003).

115

[68] O. F. Syljuåsen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002).

[69] A. W. Sandvik and J. Kurkijärvi, Phys. Rev. B 43, 5950 (1991).

[70] M. Troyer, S. Wessel, and F. Alet, Phys. Rev. Lett. 90, 120201 (2003).

[71] F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001).

[72] A. Albuquerque et al., J. Magn. Magn. Mater. 310, 1187 (2007).

[73] P. Henelius and A. W. Sandvik, Phys. Rev. B 62, 1102 (2000).

[74] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar,

Phys. Rev. B 41, 9301 (1990).

[75] K. Siemensmeyer, K. Habicht, T. Lonkai, S. Mat’as, S. Gabï£¡ni, N. Shitsevalova, E. Wulf,

and K. Flachbart, J. Low Temp. Phys. 146, 581 (2007).

[76] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94, 170201 (2005).

[77] F. D. Murnaghan, Proc. Natl. Acad. Sci. U.S.A. 30, 244 (1944).

[78] N. D. Mermin, Phys. Rev. 137, A1441 (1965).

[79] W. Nolting, Grundkurs Theoretische Physik 6, Springer-Verlag, Berlin, 2005.

[80] B. Grabowski, T. Hickel, and J. Neugebauer, Phys. Rev. B 76, 024309 (2007).

[81] B. Grabowski, L. Ismer, T. Hickel, and J. Neugebauer, Phys. Rev. B 79, 134106 (2009).

[82] J. Harris, R. O. Jones, and J. E. Mï£¡ller, J. Chem. Phys. 75, 3904 (1981).

[83] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).

[84] S. Boeck et al.,The s/phi/nx package, http://www.sphinxlib.de/.

[85] V. Heine, A. I. Liechtenstein, and O. N. Mryasov, EPL (Europhysics Letters) 12, 545 (1990).

[86] L. Kaufman and H. Bernstein, Computer Calculation of Phase Diagrams, Academic, New

York, 1970.

[87] N. Saunders and A. P. Miodownik, Calphad (Calculation of Phase Diagrams): A Comprehen-

sive Guide, Pergamon, Oxford, 1998.

[88] J. M. Leger, C. Loriers-Susse, and B. Vodar, Phys. Rev. B 6, 4250 (1972).

[89] S. Morán, C. Ederer, and M. Fähnle, Phys. Rev. B 67, 012407 (2003).

[90] E. Şaşıoğlu, L. M. Sandratskii, and P. Bruno, Phys. Rev. B 71, 214412 (2005).

[91] L. Bergqvist, B. Belhadji, S. Picozzi, and P. H. Dederichs, Phys. Rev. B 77, 014418 (2008).

[92] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).

116

[93] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).

[94] P. E. Blöchl, Phys. Rev. B 50, 17 953 (1994).

[95] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).

[96] J. Zhang and F. Guyot, Phys. Chem. Miner. 26, 206 (1999).

[97] O. Grotheer, C. Ederer, and M. Fähnle, Phys. Rev. B 63, 100401 (2001).

[98] C. Loong, J. Carpenter, J. Lynn, R. Robinson, and H. Mook, J. Appl. Phys. 55, 1895 (1984).

[99] J. Lynn, Phys. Rev. B 11, 2624 (1974).

[100] J. Crangle and G. M. Goodman, Proceedings of the Royal Society of London, Series A,

Mathematical and Physical Sciences 321, 477 (1971).

[101] F. M. Wang and R. Ingalls, Phys. Rev. B 57, 5647 (1998).

[102] D. C. Wallace, Thermodynamics of Crystals, Dover, New York, 1998.

[103] R. F. Sabiryanov and S. S. Jaswal, Phys. Rev. Lett. 83, 2062 (1999).

[104] Landolt-Börnstein, editor, Metals: Phonon States, Electron States and Fermi Surfaces, vol-

ume 13, Springer-Verlag, Berlin, 1981, New Series, Group III.

[105] A. F. Guillemet and P. Gustafson, High Temp. High Press. 16, 591 (1985).

[106] D. C. Wallace, P. H. Sidles, and G. C. Danielson, J. Appl. Phys. 31, 168 (1960).

[107] Y. S. Touloukian, R. K. Kirby, R. E. Taylor, and P. D. Desai, Thermophysical Properties of

Matter, volume 12, IFI/Plenum, New York, 1975.

[108] R. A. Tahir-Kheli, Phys. Rev. 159, 439 (1967).

[109] J. A. Copeland and H. A. Gersch, Phys. Rev. 143, 236 (1966).

[110] S. H. Liu and D. B. Siano, Phys. Rev. 164, 697 (1967).

[111] M. D. Kuz’min, Phys. Rev. Lett. 94, 107204 (2005).

[112] M. D. Kuz’min, M. Richter, and A. N. Yaresko, Phys. Rev. B 73, 100401 (2006).

[113] M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics, Oxford

University Press, USA, 1999.

[114] The critical temperature within the quantum case does not strictly follow TC(S)∼JS(S+1) =

Jas was, e.g., recently shown using high temperature series expansion [J. Phys.: Condens.

Matter 16, 8653 (2004)]. The complex dependence of TC(S)is, however, beyond the scope of

the present manuscript.

117

[115] The exclusions are TC(S)as discussed in the text and the integration constant UQMC(0) for

the internal energy, which cannot be determined solely from the heat capacity. To determine

UQMC(0) we use UcMC(Th)≡UQMC(Th), which holds true for high temperatures Th. For all

considered systems the choice Th= 3TCwas in this respect suﬃcient.

[116] The nearest neighbor antiferromagnetic fcc lattice is known to be frustrated, see e.g. [J. Low

Temp. Phys. 146, 581 (2007)], and is therefore not considered in this study.

[117] By taking 10 shells into account the average sign hσiTdecreases to 10−3for temperatures

below 2TC. To achieve statistical errors of less than 5%, the calculation time on a conventional

computer is in the order of years for a single temperature calculation.

[118] If the heat capacity shows a ﬁnite cusp with a magnitude CP(TC) = Cmax

P, a ﬁnite size ansatz

can be used to extrapolate the dependence Cmax

P(L)on system size L. Depending on the ﬁnite

size scaling ﬁt, the system size required to obtain, e.g., 50% of Cmax

Pmay result in system

sizes containing up to 107atoms (see e.g. Phys. Rev. B 42, 6087 (1990)), which renders the

problem to be computationally unfeasible.

[119] The calphad data has been obtained with the thermocalc program and the SGTE unary

database [105].

[120] R. J. Weiss and K. J. Tauer, Phys. Rev. 102, 1490 (1956).

[121] H. C. Herper, E. Hoﬀmann, and P. Entel, Phys. Rev. B 60, 3839 (1999).

[122] J. Zarestky and C. Stassis, Phys. Rev. B 35, 4500 (1987).

[123] C. Engin, L. Sandoval, and H. M. Urbassek, Modelling and Simulation in Materials Science

and Engineering 16, 035005 (2008).

[124] H. C. Hsueh, J. Crain, G. Y. Guo, H. Y. Chen, C. C. Lee, K. P. Chang, and H. L. Shih, Phys.

Rev. B 66, 052420 (2002).

[125] M. Y. Lavrentiev, D. Nguyen-Manh, and S. L. Dudarev, Phys. Rev. B 81, 184202 (2010).

[126] M. Lavrentiev, S. Dudarev, and D. Nguyen-Manh, Journal of Nuclear Materials 386, 22

(2009).

[127] D. J. Dever, J. Appl. Phys. 43, 3293 (1972).

[128] Y. Tsunoda and N. Kunitomi, Solid State Commun. 54, 547 (1985).

[129] Y. Tsunoda and N. Kunitomi, J. Phys. F: Met. Phys. 18, 1405 (1988).

[130] Y. Tsunoda, J. Phys.: Condens. Matter 3, 7231 (1991).

[131] B. Jonker, K. Walker, E. Kisker, G. Prinz, and C. Carbone, Phys. Rev. Lett. 57, 142 (1986).

[132] W. Macedo and W. Keune, Phys. Rev. Lett. 61, 475 (1988).

118

[133] H. Magnan, D. Chandesris, B. Villette, O. Heckmann, and J. Lecante, Phys. Rev. Lett. 67,

859 (1991).

[134] J. Thomassen, F. May, B. Feldmann, M. Wuttig, and H. Ibach, Phys. Rev. Lett. 69, 3831

(1992).

[135] S. Müller, P. Bayer, C. Reischl, K. Heinz, B. Feldmann, H. Zillgen, and M. Wuttig, Phys.

Rev. Lett. 74, 765 (1995).

[136] R. Ellerbrock, A. Fuest, A. Schatz, W. Keune, and R. Brand, Phys. Rev. Lett. 74, 3053

(1995).

[137] D. Keavney, D. Storm, J. Freeland, I. Grigorov, and J. Walker, Phys. Rev. Lett. 74, 4531

(1995).

[138] M. Straub, R. Vollmer, and J. Kirschner, Phys. Rev. Lett. 77, 743 (1996).

[139] Y. Tsunoda, N. Kunitomi, and R. M. Nicklow, Journal of Physics F: Metal Physics 17, 2447

(1987).

[140] I. Turek, J. Kudrnovský, M. Diviš, P. Franek, G. Bihlmayer, and S. Blügel, Phys. Rev. B 68,

224431 (2003).

[141] J. Kuneš and R. Laskowski, Phys. Rev. B 70, 174415 (2004).

[142] W. Huang, Calphad 13, 243 (1989).

[143] N. Singh, Phys. Status Solidi B 156, K33 (1989).

[144] P. D. Desai, J. Phys. Chem. Ref. Data 15, 967 (1986).

[145] Q. Chen and B. Sundman, J. Phase Equilib. 22, 631 (2001).

[146] A. F. Guillemet, Int. J. Thermophys. 8, 481 (1987).

[147] J. E. Saal, S. Shang, Y. Wang, and Z.-K. Liu, Journal of Physics: Condensed Matter 22,

096006 (2010).

[148] P. Peczak, A. M. Ferrenberg, and D. P. Landau, Phys. Rev. B 43, 6087 (1991).

[149] R. G. Brown and M. Ciftan, Phys. Rev. Lett. 78, 2266 (1997).

[150] R. G. Brown and M. Ciftan, Phys. Rev. Lett. 76, 1352 (1996).

[151] C. Holm and W. Janke, Phys. Rev. Lett. 78, 2265 (1997).

[152] Landolt-Börnstein, editor, Metals: Phonon States, Electron States and Fermi Surfaces, volume

19a, Springer-Verlag, Berlin, 1981, New Series, Group III.

[153] B. Hallstedt, D. Djurovic, J. von Appen, R. Dronskowski, A. Dick, F. Körmann, T. Hickel,

and J. Neugebauer, Calphad 34, 129 (2010).

119

[154] A. Dick, F. Körmann, T. Hickel., and J. Neugebauer, submitted, 2011.

[155] C. M. Fang, M. H. F. Sluiter, M. A. van Huis, C. K. Ande, and H. W. Zandbergen, Phys.

Rev. Lett. 105, 055503 (2010).

[156] H. P. Scott, Q. Williams, and E. Knittle, Geophys. Res. Lett. 28, 1875 (2001).

[157] B. J. Wood, Earth and Planetary Science Letters 117, 593 (1993).

[158] J. Chipman, Metall. Transactions 3, 55 (1972).

[159] E. Duman, M. Acet, T. Hülser, E. F. Wassermann, J. P. I. B. Rellinghaus, and P. Munsch, J.

Appl. Phys. 96, 5668 (2004).

[160] J. H. Jang, I. G. Kim, and H. Bhadeshia, Comp. Mat. Science 44, 1319 (2009).

[161] A. K. Arzhnikov, L. V. Dobysheva, and C. Demangeat, J. Phys.: Condens. Matter 19, 196214

(2007).

[162] G. Miyamoto, J. Oh, K. Hono, T. Furuhara, and T. Maki, Acta Mater. 55, 5027 (2007).

[163] H. I. Faraoun, Y. D. Zhang, C. Esling, and H. Aourag, J. Appl. Phys. 99, 093508 (2006).

[164] K. O. E. Henriksson, N. Sandberg, and J. Wallenius, Appl. Phys. Lett. 93, 191912 (2008).

[165] M. Umemoto, Z. Liu, K. Masuyama, and K. Tsuchiya, Scr. Mater. 45, 391 (2001).

[166] G. Naeser, Mitt. Kais.-Wilh.-Inst. Eisenforschg. 16, 207 (1934).

[167] L. S. Darken and R. W. Gurry, J. Metals 3, 1015 (1951).

[168] H. Seltz, H. McDonald, and C. Wells, Trans. AIME 140, 263 (1940).

[169] J. Rusz, I. Turek, and M. Diviš, Phys. Rev. B 71, 174408 (2005).

[170] W. Nolting, A. Vega, and T. Fauster, Z. Phys. B 96, 357 (1995).

[171] J. Staunton, B. Gyorﬀy, A. Pindor, G. Stocks, and H. Winter, Journal of Magnetism and

Magnetic Materials 45, 15 (1984).

[172] S.-L. Shang, Y. Wang, and Z.-K. Liu, Phys. Rev. B 82, 014425 (2010).

120

TEXT IN WHITE in order to keep this page empty

Publication List:

• A. Dick, F. Körmann, T. Hickel, and J. Neugebauer, “Ab initio based determination of ther-

modynamic properties of cementite including vibronic, magnetic and electronic excitations”,

submitted (2011).

• F. Körmann, A. Dick, T. Hickel, and J. Neugebauer, “Role of spin quantization in determining

the thermodynamic properties of magnetic transition metals”, Phys. Rev. B 83, 165114 (2011).

• M. Friák, T. Hickel, F. Körmann, A. Udyansky, A. Dick, J. von Pezold, D. Ma, O. Kim,

W. A. Counts, M. Šob, T. Gebhardt, D. Music, J. Schneider, D. Raabe, and J. Neugebauer,

”Determining the elasticity of materials employing quantum-mechanical approaches: From the

electronic ground state to the limits of materials stability“, Steel Res. Int. 82, 86 (2011).

• F. Körmann, A. Dick, T. Hickel, and J. Neugebauer, “Rescaled Monte Carlo approach for

magnetic systems: Ab initio thermodynamics of bcc iron“, Phys. Rev. B 81, 134425 (2010).

• B. Hallstedt, J. von Appen, A. Dick, F. Körmann, T. Hickel, and J. Neugebauer, “Thermody-

namic properties of cementite (Fe3C)“, Calphad 34, 129 (2010).

• J. Neugebauer, B. Grabowski, F. Körmann, A. Dick, J. von Pezold, M. Friak, and T. Hickel,

“Ab Initio Based Multiscale Modeling of Engineering Materials: From a Predictive Thermo-

dynamic Description to Tailored Mechanical Properties“, Asia Steel (Proceeding), (2009).

• F. Körmann, A. Dick, T. Hickel, and J. Neugebauer, “The pressure dependence of the Curie

temperature in bcc iron studied by ab initio simulations“, Phys. Rev. B 79, 184406 (2009).

• T. Hickel, A. Dick, B. Grabowski, F. Körmann, and J. Neugebauer, “Steel design from fully

parameter-free ab initio computer simulations“, Steel Res. Int. 80, 4 (2009).

• F. Körmann, J. Kienert, S. Schwieger, and W. Nolting, “Cu cap layer on Ni8/Cu(001): reori-

entation and TC-shift“, Eur. Phys. J. B 65, 499-504 (2008).

• F. Körmann, A. Dick, B. Grabowski, B. Hallstedt, T. Hickel, and J. Neugebauer, “Free energy

of bcc iron: Integrated ab initio derivation of vibrational, electronic, and magnetic contribu-

tions“, Phys. Rev. B “78, 033102 (2008).

• S. Henning, F. Körmann, S. Schwieger, J. Kienert, and W. Nolting, “Green function theory

versus quantum Monte Carlo calculations for thin magnetic films“, Phys. Rev. B 75, 214401

(2007).

• F. Körmann, S. Schwieger, J. Kienert, and W. Nolting, “A new type of temperature driven

reorientation transition in magnetic thin films“, Eur. Phys. J. B 53, 463-469 (2006).