Theoretical understanding and calculation of the Edelstein effect

(1)

TVE-F 17 017 juni

Examensarbete 15 hp Juni 2017

Theoretical understanding and calculation of the Edelstein effect

Gustav Eriksson

Hampus Nyström

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Theoretical understanding and calculation of the Edelstein effect

Gustav Eriksson, Hampus Nyström

The main topic of this project is the so called Edelstein effect.

This recently discovered effect consists in the possibility of converting an electric field (a current) into a magnetization in materials that fulfill specific characteristics, more specifically materials where an effective Rashba spin-orbit coupling is present.

The Edelstein effect is appealing to the scientific community from the fundamental physics point of view as well as from the technological point of view. In fact the possibility of efficiently converting an electric signal into a magnetic signal could revolutionize the current information storage technology.

In this project, after a study of basic concepts of solid state physics: crystal structure, Bloch's theorem, spin-orbit coupling; we addressed the study of the basics of a powerful numerical tool, called density functional theory (DFT), for predicting the electronic properties of solids. This tool provides us with all the needed quantities for numerically calculating any kind of linear response, which we show that the Edelstein effect is a specific form of. Using a specific implementation of DFT, called augmented spherical wave (ASW), we calculate the Edelstein effect in iron and copper (where no effect is expected) and manganese silicide (where the effect is

expected to appear). We also perform a systematic study on how the Edelstein effect depends on the symmetry of the material and the magnitude of the spin-orbit coupling. The calculations showed

promising results from which we concluded that the numerical methods used could clearly distinguish between the presence of the Edelstein effect or not in mentioned materials.

Ämnesgranskare: Martin Sjödin

Handledare: Marco Berritta, Peter M. Oppeneer

(3)

Populärvetenskaplig sammanfattning

Användningsområdena för datorer har ökat konstant med ett kontinuerligt ökande krav på komponenterna både när det kommer till prestanda men även platseffektivitet. I dagens datorer används elektroners laddning som som informationsbärare, dessa kan represen- teras av ett binärt talsystem, ettor och nollor, och har länge fungerat med tillförlitliga tekniska implementationer. I framtiden kommer dock nya tekniker behövas vilket för med sig ett krav på nya teoretiska modeller. Ett förslag innebär att använda elektroner- nas spinnmagnetisering som informationsbärare. Om spinnmagnetiseringen i ett material skulle användas i komponenter, som ersättning till dagens RAM minnen till exempel, så skulle mer information kunna lagras på mindre yta. Tekniker för att manipulera mag- netiseringen i dessa komponenter tyder också på att informationen kan ändras snabbare, vilket möjliggör mer effektiv hantering av information i datorer. En utav dessa teorier är den så kallade Edelsteineffekten.

Edelsteineffekten omvandlar ett elektriskt fält, eller en strömtäthet, till en spinnpolar- isation i material med vissa specifika egenskaper. I det här projektet har vi teoretiskt undersökt Edelsteineffekten, vad som krävs för att den ska uppstå, och med numeriska metoder försökt beräkna den för tre olika material. Målet var att få en teoretisk förståelse för effekten och dess härkomst men också att verifiera om programmet, som de numeriska metoderna var implementerade i, kan användas i fortsatta studier. Alltså undersöka om de resultat programmet ger oss stämmer överens med de förväntningar teorin ger.

Beräkningarna genomfördes på UPPMAX superdatorer i Uppsala och visade sig ge goda resultat. Vi kunde verifiera att programmet identifierade Edelsteineffekten i material där den var förväntad och gav indikationer på obefintlig effekt i material där den inte förväntades. Med hjälp av resultaten kunde vi estimera det numeriska felet och dra slut- satser gällande dessa metoders användningsområden när det kommer till experimentella mätningar av Edelsteineffekten. Vi noterade att programmet, och därmed de numeriska metoderna, fungerar väl och kan, med relativt få modifikationer, utökas till att undersöka fler material och därmed användas för att i framtiden identifiera material lämpliga för implementationer i tekniska komponenter.

(4)

1 Introduction

1.1 Background

Nowadays the random-access memories (RAM) are based on devices where information is stored in the form of electric currents and electric charges. An alternative to the current technology is based on the so called spintronic devices, where the information units are encoded through the magnetic states of a magnetic material. Implementing such a device requires a way to read and write (or manipulate) the magnetic information units. Meth- ods have already been developed for reading the information on such devices and can, in the case of antiferromagnetic (AFM) states, be based on different magnetoresistance effects.¹ These methods basically measure effects that describe correlations between magnetic fields and the resistance of a material in those fields.

The problem that is still being heavily investigated is instead how to efficiently manipulate the magnetization. This can be done by applying an external magnetic field but an electrical method is more desirable, of which there are two approaches of current interest. One is to utilize the spin-transfer torque which requires a spin-polarized current, in essence transferring spin-angular momentum from the spin-polarized current to the material. The other way is using spin-orbit torque which instead is caused by spin-orbit coupling in the material. This requires a material with a broken inversion symmetry but the current does not need to be spin-polarized. The Edelstein effect is one of the effects that lead to such a spin-orbit torque and is therefore of great interest. With an electric field along a certain direction, this effect induces in certain materials a magnetization perpendicular to the electric field itself.

A way to quantify the Edelstein effect in different materials is introduced by the Kubo linear response theory, which can be used to determine the relationship between the induced spin polarization, M, and an applied electrical field, E, such that M = χE. Where χ is a material dependent tensor. In previous studies a DFT (density functional theory) based method, which accounts for the specifics of the material under consideration, has been implemented in a program, called ’MOKE’, to calculate the optical conductivity and the magneto-optical Kerr effect (MOKE) in different materials.² Similar methods can be used for calculating the mentioned χ tensor caused by the Edelstein effect. In this project we will use a modified version of the MOKE program, called ’EDEFF’, to calculate certain elements of the χ tensor. In particular we will run the EDEFF program for iron, copper and manganese silicide and study the results with different expectations due to the symmetry properties of the materials. Comparing the results of the calculations with these expectations will then provide us with a general idea whether the program works and if so how it can be used. Additionally, we will conduct a theoretical study of the required basics of solid state physics to obtain an understanding of the origins of the Edelstein effect. This also includes the basics of the needed tools to formulate and numerically calculate the χ tensor.

Since iron and copper have centrosymmetric crystal structures we expect the Edelstein effect to be non-present, i.e. we expect the χ tensor calculated from the program to be

(6)

zero. According to the theory the effect is expected in materials without inversion symmetry, which manganese silicide is an example of, so consequently we expect a non-zero χ tensor from those calculations. The general aim of this project is to validate the program, to study the results of the calculations done with different settings and attempt to draw conclusions whether the program yields credible results. Furthermore, we want to gain a deeper understanding of the Edelstein effect and the physics behind it.

1.2 Theory

This section introduces the basic theory needed to understand the Edelstein effect and the calculation of its linear response (the χ tensor).

1.2.1 Crystal Structure

The structure of a crystal, how the individual atoms are ordered in a material, is a major factor that determines the material properties. This is true for how it conducts electricity and heat but also for its optical absorption and many other properties.³ The Edelstein effect is not different from these and therefore a basic understanding of crystal structures is needed to comprehend it.

What defines a crystal is a unit cell (a portion of space which contains a set of atoms) which periodically repeats itself in the whole space. It can be proven that there exist no more than 14 different ways in which a 3D structure can periodically repeat itself (5 in 2D), these are called Bravais lattices. We can then classify all the possible crystals in 3D according to the so called Bravais lattices. On top of these "translational symmetries" it is possible to have other transformations transforming the crystal in itself. For instance if the lattice that is formed from infinitely many cubes like the one in figure 1, where all the vertices are equivalent (for example when all vertices are formed by one atom), a rotation by 90° around the middle of any face would be a symmetry operation.

Figure 1: Simple cubic unit cell.

The composition of the translational symmetries of the Bravais lattices with the so called crystallographic point groups (of transformations which can transform the crystal in itself, for instance rotations, reflections on a specific plane etc.) allows us to define a total of 230 3D crystal space groups which all crystals adhere to. These 230 crystal structures always has an underlying Bravais lattice from which their unit cell is derived.

(7)

What is seen in figure 1 is known as a unit cell of the simple cubic Bravais lattice.

This is one way of defining a unit cell and it is good for visualizing how it looks but it is far from the only one. The only requirement of a unit cell is that it should be a space that when stacked produces the entire lattice and because of this definition it can be created in many different ways, one of which is the Wigner-Seitz cell. Unlike the one in figure 1 it only contains one lattice point and is defined as the space that is closer to that lattice point than any other. For the simple cubic Bravais lattice this just results in another cubic cell with one lattice point in the middle instead but for other Bravais lattices the changes are greater. For instance the body-centered cubic (BCC) Bravais lattice (figure 2) has the Wigner-Seitz cell seen in figure 3. Although the simple unit cell seen in figure 2 gives a clear image of how the lattice looks the Wigner-Seitz cell is more useful in many calculations.

Figure 2: Body-centered cubic unit cell.

Figure 3: Body-centered cubic Wigner-Seitz cell.

There are three lattices that are of particular interest in this paper. The first is the body centered cubic (BCC) lattice which has already been visualized in figure 2 and 3. The other two Bravais lattices are body centered tetragonal (BCT) and orthorhombic Bravais lattices. In the case of copper it is worth noting that the BCT lattice, seen in figure 4, is equal to the FCC (face centered cubic) lattice due to the relation c =√

2 a, where c the height of the cell and a is the side length of the square base. In all other ways the complete BCT lattice is formed in the same way as for the BCC lattice.

(8)

Figure 4: Body-centered tetragonal unit cell.

The orthorhombic unit cell in turn looks the same as the two previous ones but without a point in the middle where all the sides are of different length. What also can be stated is that iron has a BCC structure and copper a BCT (or FCC) structure which both have inversion symmetry, meaning that mirroring at a point in the lattice is a symmetry operator that brings the lattice into itself. Manganese silicide on the other hand only has the translational symmetry of an orthorhombic unit cell but consequently has an atom structure that breaks the full symmetry, in particular the inversion symmetry of the lattice, which is or great relevance to the Edelstein effect.

1.2.2 Reciprocal space

The reciprocal space, or k-space,³ is the Fourier transform of the Bravais lattice. In other words it is a representation of the periodicity of the crystal structure. The formal definition is given as follows:

The set of all wave vectors K that yield plane waves with the same periodicity as a given Bravais lattice, is known as it’s reciprocal space.

This give rise to the following condition; K belongs to the reciprocal space of a Bra- vais lattice if

e^iK·(r+R) = e^iK·r (1)

holds for any arbitrary vector r and all points R in the Bravais lattice.

There are several important qualities of the reciprocal space which are not easily re- alized at first. The first of these is the correlation between planes in the Bravais lattice and vectors in the reciprocal space. In any Bravais lattice it is possible to identify lattice planes which are planes that contain at least three lattice points. Each of these belong to a family of lattice planes which are parallel and spaced an equal distance apart which when added together contain all the points of the Bravais lattice. One of the interesting properties of these planes are that for every vector in the reciprocal lattice there is a plane in direct space that is orthogonal to the vector. This is helpful due to the relevance that planes have in the diffraction of photons. In other words the reciprocal space simplifies the calculations but they can still be made in the direct space. This effect transfers to many more calculations where the wave properties of particles are of interest.

(9)

As in the direct lattice, we can also define a unit cell in the reciprocal space. In the reciprocal space the Wigner-Seitz cell is instead called the first Brillouin zone but is oth- erwise defined the exact same way as for the direct lattice. Due to this it obviously has the same benefits and these benefits are intimately related to what is known as Bloch’s Theorem.

1.2.3 Bloch states

Bloch’s theorem (theorem 1) is a powerful tool for studying quantum effects in periodic potentials of any kind.³

Theorem 1 The eigenstates ψ of the one-electron Hamiltonian H = ^−¯_2m^h²∇² + U (r), where U (r + R) = U (r) for all vectors R in the Bravais lattice, can be chosen to have the form of a plane wave times a function with the periodicity of the Bravais lattice:

ψ_nk(r) = e^ik·ru_nk(r) (2)

where

u_nk(r + R) = u_nk(r) (3)

for all R in the Bravais lattice.

This can also be written in an alternative form:

ψnk(r + R) = e^ik·Rψnk(r) (4)

The wave vector k which is introduced in Bloch’s theorem is of particular interest due to the fact that we never have to consider any effects outside of the first Brillouin zone (the Wigner-Seitz cell of the reciprocal space). The reason for this is that any k not in the first Brillouin zone can be written as k = kb + K where K is a vector in the reciprocal lattice and kb lies in the first Brillouin zone. Due to equation 1 the relation e^iK·R = 1 holds for all vectors in the reciprocal lattice, which proves that if equation 4 holds for kb

it also holds for k.

This wave vector is also interesting in the regard that it plays a similar role as the wave vector in the free electron model. This is not to say that this wave vector is a measure of the momentum; this can be proven to be false as they do not commute. Instead, it is to say it has many of the properties for which momentum is so useful and it describes the translational symmetry of a periodic potential the same way that momentum is a characteristic of the translational symmetry of free space. Therefore, it is a useful way of labeling the eigenstates of the electronic wave function in a solid. This is also one of the reasons why calculations in the reciprocal space are so useful.

The index n in equations 2 to 4 is the result of the periodicity of the system and the complexity of the unit cell. For every such n we can vary k continuously which instead of constant energy levels gives us continuous energy functions εn(k + K) = ε_n(k) with the same periodicity as the reciprocal lattice. These functions are what defines the band structure of the material. One of the reasons for this name is quite apparent. Due to the

(10)

fact that all εn(k) are both periodic and continuous they create an upper and a lower bound which all values of εn(k) lies within, creating a band. These bands are a direct result of the Pauli exclusion principle and also gives rise to what is known as the Fermi surface.

1.2.4 Fermi surface

The ground state of a number of Bloch atoms³ is not as straightforward as for the free electron model. In the case of the free electron model the ground state is given by filling all the energy levels ε(k) = ^¯^h_2m²^k² with energies less than the Fermi energy εF. The Fermi energy in turn is defined by requiring that the number of energy levels less than εF is the same as the total number of electrons in the system. We can then say that εF characterizes the ground state of a system with a specific number of electrons. The trouble when we are observing Bloch electrons is that ε is no longer dependent only on k but on both the crystal and the band index n. This gives rise to two different scenarios when looking at the ground state of a number of Bloch electrons. The first is that there are a number of full energy bands and all others are empty. This is when we can define what is known as a band gap in the material. The band gap is the difference in energy between the highest occupied energy band and the lowest unoccupied one and in this band gap we find the Fermi energy. The value of this band gap is very important for how a solid behaves, especially regarding the electrical conduction of the material. The other scenario is that there are a number of partially filled bands. In this case we can define a surface in the k-space that separates the levels that are filled from the ones that are not. The set of all these surfaces is what is known as the Fermi surface and is important for understanding how the electrons in a material behaves, especially for the electrical conductivity.

Figure 5: Visualization of the band structure of an insulating material.

What can be seen in figure 5 is a visualization of the band structure of an unknown material. From this diagram we can draw the conclusion that it is not a metal, as the Fermi energy is not in the middle of one or more energy bands. Instead, it must be either an insulator or a semiconductor, depending on the magnitude of the band gap. This is a simple theoretical case, as the energy bands rarely are evenly spaced but it gives a good visualization of what the Fermi energy and Fermi surface says about a material. When an

(11)

electron is excited to the conduction band this allows charge to flow in the material, both through the electron in the conduction band but also through the hole it left behind in the valence band. What decides how easily a material can conduct electricity is therefore how easily an electron can pass from the valence band to the conduction band. As was previously stated the Fermi energy cuts right through one or more bands when it comes to metal, which can be seen in figure 6, and hence this band contain both occupied and unoccupied states. This is the reason why they conduct charges so easily; there is no energy difference between the highest occupied and the lowest unoccupied state.

Figure 6: Visualization of the band structure of a metal.

For us, studying the Edelstein effect, we assume that there will always be a Fermi surface as we are attempting to manipulate conducting magnetic materials. The fact that we are attempting to manipulate the magnetization with a current in conducting materials allows us rewrite the expression for the Edelstein effect. Using the relation ji = σ_ijE_j, where σij is the conductivity of the material and jy the current density, gives us Mx = χ_xyE_y = χ_xy(_σ^j^y

yy) where Mx is the spin-polarization and χxy is the material dependent linear response of the Edelstein effect. Thus, the spin-polarization can be expressed in terms of the current density as Mx = χ_xyρ_yyj_y, where ρyy is the resistivity of the material.

1.2.5 Spin-orbit coupling

It has been proven that the Edelstein effect is intimately related to the so called spin- orbit coupling,⁴ which can be pictured as the interaction between the magnetic moment of the spin of an electron and the magnetic moment due to its orbital angular momentum.

It can be shown that for an electron in a centrosymmetric potential we can model the spin-orbit coupling by adding a term in the Hamiltonian of the electron of the form

H_SOC = 2µB

¯ hm_eec²

1 r

∂U (r)

∂r

L · ˆˆ S (5)

where µB is the Bohr magneton, me the mass of the electron, e the electrons charge, c the speed of light, U(r) is the spherically symmetric potential, ˆL is the orbital angular momentum of the electrons and ˆS is its spin.

(12)

Depending on the symmetry of the system the spin-orbit coupling can assume different forms, in particular the Edelstein effect is intimately related to the Rashba spin-orbit coupling which has the form

H_R−SOC = α(ˆp × ˆS) · ˆz (6)

where α is a coupling constant, ˆp is the momentum of the electrons, ˆS their spin and ˆz a direction determined by the symmetry of the system. This additional term in the Hamiltonian is what causes the Edelstein effect, therefore we expect a larger effect if this term grows and a smaller effect, or non-present, if this term decreases towards zero. What can also be shown is that in crystals where there exists inversion symmetry the net effect of this coupling is not present and is therefore the reason for why the non-zero Edelstein effect does not arise in such materials.

1.2.6 Density functional theory

Even though the Born-Oppenheimer approximation simplifies the many-body Schrödinger equation, fixing the nuclei to a static external potential, the problem of solving it is still too complicated to be done conveniently.⁵ The interaction between electrons in the system prevents the equation to be divided into multiple simpler ones. So, in order to numerically solve the many-body Schrödinger equation efficiently, an approximation has to be used. This is where density functional theory (DFT) is useful.⁶ The DFT allows us to identify the ground-state properties of a system without having to deal with the many-electron state. Using the concept of a spatially dependent electron density the theory maps the many-body problem, including the electron interaction, onto a one-body problem. The Born-Oppenheimer Hamiltonian (in Gaussian units) is defined by

HBO = −X

i

¯ h²

2m∇²_i −X

i,A

e²Z_A r_iA +X

i>j

e² 1

r_ij + X

B>A

e²Z_AZ_B

R_AB (7)

where m is the electron mass, e the electron charge, ZA the atomic numbers of the atoms, ri the coordinates of the electrons and RAB are the coordinates of the nuclei. In comparison, the Kohn-Sham Hamiltonian is defined as

H_KS = −h¯²

2m∇²+ V_ext+ Z

dr⁰ ρ(r⁰)

|r − r⁰|+ V_xc(r) (8) and can be seen as a single electron Hamiltonian. In fact while ˆH_BOact on state functions of the form ψ(r1, r₂, r₃...), where r1, r₂, r₃... are the coordinates of the various electrons, H_KS depends only on r. This simplifies the calculations significantly. In the Kohn-Sham Hamiltonian Vxc is called exchange potential and is a crucial quantity in DFT, we will not discuss its details here however since it is out of the scope of this project. It is worth noting here however that the Kohn-Sham Hamiltonian must have the same symmetry of the crystal of the material it should model. Despite it can be proven that HKS is exact only for calculating the ground state of a material it provides a powerful tool to reliably calculate the band structure of a wide range of materials. In practice we now have a tool for calculating HKS|nki = ε_nk|nki where εnk are the bands and |nki the

(13)

eigenstates of HKS. We can then calculate, at least in principle, the matrix elements of any observable, in particular a spin-polarization. Using this we essentially convert the interaction potentials of the electrons into an outer density potential, providing us with a numerically solvable one-body Hamiltonian which in turn provides us with the ground-state properties needed to calculate the Edelstein effect.

1.2.7 Kubo linear response theory

The Kubo linear response theory is an approximation to treat how a time dependent perturbation affects the expectation value of an observable.⁷ For the derivation of the Kubo formula we start from the density matrix formalism which does not restrict us to the treatment of pure quantum states.⁸ For a given generic observable ˆO we can, once the density matrix ˆρ is defined, calculate its expectation value as

h ˆOi = T r{ρ(t) ˆO} =X

j

p_jjhψ_j|O_j|ψ_ji, T r{a} =X

i

a_ii (9)

ρ(t) = e^{−β ˆ}^H

Z , Z = T r{e^{−β ˆ}^H} and β = 1 kβ

T (10)

where kβ is the Boltzmann constant and T the temperature. When describing the re- action of a time independent system to a small perturbation it is helpful to use the interaction picture instead of the Schrödinger picture. In the Schrödinger picture the time dependence is strictly "attached" to the states while the observables are always time independent. There is also the Heisenberg picture in which the operators carry all the time dependence and the states are time independent. The interaction picture is a mix of the two. We consider a Hamiltonian on the form ˆH = ˆH₀+ ˆH₁(t)where ˆH₀ is the Kohn-Sham Hamiltonian and ˆH₁(t) a time dependent perturbation.

The time dependence of the density operator in the interaction picture is given by d ˆρⁱ(t)

dt = 1

i¯h[ ˆH₁ⁱ(t), ˆρⁱ(t)] (11) where the index i stands for "interaction picture". After integrating we get

ˆ

ρⁱ(t) − ˆρⁱ(−∞) = 1 i¯h

Z t

−∞

dt⁰[ ˆH₁ⁱ(t⁰), ˆρⁱ(t⁰)], ρˆⁱ(−∞) = ˆρ₀. (12) Combining equation 9 with equation 12 we obtain

h Ôi = T r{ ˆρ(t) Ô} = T r{ ˆρⁱ(t) Ôⁱ(t)} = T r{ ˆρ₀O} +ˆ 1 i¯h

Z t

−∞

dt⁰T r{[ ˆH₁ⁱ(t⁰), ˆρⁱ(t⁰)] ˆOⁱ} (13)

⇒ h ˆOi = T r{ ˆρ₀Oˆⁱ} + i

¯ h

Z t

−∞

dt⁰T r{ ˆρⁱ(t⁰)[ ˆH₁ⁱ(t⁰), ˆOⁱ]}. (14) Only considering the first order of the density operator gives the first order approximation of the expectation value

h ˆOi ≈ h ˆOi₀+ i h

Z t

−∞

dt⁰h[ ˆH₁ⁱ(t⁰), ˆOⁱ(t)]i₀. (15)

(14)

The change in the expectation value from the time independent solution can be written as

δh ˆOi ≈ i

¯ h

Z t

−∞

dt⁰h[ ˆH₁ⁱ(t), ˆOⁱ(t)]i₀. (16) If we then consider ˆH₁ⁱ as two parts ˆH₁ⁱ = ˆP (t)F (t), where F is an external field and ˆP is the coupling between that field and an observable of the electronic system, we get

δh ˆOi ≈ i

¯ h

Z t

−∞

dt⁰h[ ˆPⁱ(t⁰), ˆOⁱ(t)]i₀F (t⁰) (17) If we define a response as, assuming a perturbation that is switched on at t’,

χ = T r{i

¯

hρˆ₀[ ˆPⁱ, ˆOⁱ]}Θ(t − t⁰) (18) where Θ(t) is the Heaviside function, we get

δh ˆOi ≈ Z t

0

dt⁰χ(t, t⁰)F (t⁰) (19)

χ = T r{i

¯

hρˆ₀[ ˆPⁱ, ˆOⁱ]}Θ(t − t⁰) = i

¯ h

X

n

ρ_nnhn|[ ˆPⁱ, ˆOⁱ]|niΘ(t − t⁰) (20)

= i

¯ h

X

n

ρ_nnhn|( ˆPⁱOˆⁱ − ˆOⁱPˆⁱ)|niΘ(t − t⁰) (21)

= i

¯ h

X

n

ρ_nnhn|(e^iH0t

0

¯

h P eˆ ^iH0(t−t

0)

¯

h Oeˆ ^−iH0t^¯^h − e^iH0t^¯^h Oeˆ ^iH0(t^¯^h^0−t)P eˆ ^−iH0t

0

¯

h )|niΘ(t − t⁰) (22)

= i

¯ h

X

nm

ρnn(e^iEnt

0

¯

h hn| ˆP |mie^iEm(t−t

0)

¯

h hm| Ô|nie^−iEnt^¯^h − eîEnt^¯^h hn| Ô|mieîEm(t0−t)^¯^h hm| ˆP |nie^−iEnt

0

¯

h )Θ(t − t⁰)

(23)

= i

¯ h

X

nm

ρ_nn(e^−iEn(t−t

0)

¯

h e^iEm(t−t

0)

¯

h P_nmO_mn− e^iEn(t−t

0)

¯

h e^iEm(t^¯^h^0−t)O_nmP_mn)Θ(t − t⁰) (24)

= i

¯ h

X

nm

(ρ_nnei(Em−En)(t−t0 )

¯

h P_nmO_mn− ρ_mmei(Em−En)(t−t0

)

¯

h P_nmO_mn)Θ(t − t⁰) (25)

= i

¯ h

X

nm

(ei(Em−En)(t−t0)

¯

h P_nmO_mn(ρ_nn− ρ_mm))Θ(t − t⁰) (26) Where ρnn is the Fermi-Dirac function. A Fourier transform with respect to δt = t − t⁰ and adding limδt→0e^−ηδt in order for the integral to converge gives us

χ(ω) = lim

η→0⁺

i

¯ h

X

nm

P_nmO_mn(ρ_nn− ρ_mm) Z ∞

0

e^i(ω^mn^+ω+iη)δtdδt, ω_mn= E_m− E_n

¯

h (27)

→ χ(ω) = lim

η→0⁺

−1

¯ h

X

nm

P_nmO_mn(ρ_nn− ρ_mm)

ωmn+ ω + iη . (28)

Additionally, it can be proven that a finite value of η can account for the so called line broadening (lifetime of a state).

(15)

1.2.8 Edelstein effect

To measure the Edelstein effect we want to determine the off-diagonal term of the tensor χand in order to do that we have to define the operator ˆP, the field F and the observable Oˆ that are involved. The Edelstein effect is an effect where the electrical field affects the spin polarization of a material. This means that F is the electric field, ˆP is the coupling between the electric field and the electron density and the observable ˆO is the spin of the material. This gives rise to the following equations that will be substituted in equation 28.

O_nm = S_nm^x , P_mn = er_mn^y , r^y_mn= −i¯h m_e

p^y_mn

E_m− E_n (29)

where Snm^x is the spin in the x-direction, e is the electron charge and rmn^y is the position of the electron. The third relation in equation 29 is due to the commutation relation ˆ

p_i = ^m_i¯_h^e[ ˆH, ˆr_i], where me is the mass of the electron, p^y_mn is the momentum and Em and E_n are the energies. This leads to the following expression:

χ(ω) = ie

¯ hm_e

X

nm

(ρ_nn− ρ_mm) ω_nm

p^y_mnS_nm^x

ω − ω_mn+ iη (30)

= ie

¯ hm_e

( X

n>m

(ρ_nn − ρ_mm) ω_nm

p^y_mnS_nm^x

ω − ω_mn+ iη +X

m>n

(ρ_nn− ρ_mm) ω_mn

p^y_mnS_nm^x ω − ω_mn+ iη

)

(31)

= ie

¯ hm_e

( X

n>m

(ρ_nn− ρ_mm) ω_nm

p^y_mnS_nm^x

ω − ω_mn+ iη +(ρ_mm− ρ_nn) ω_mn

p^y_nmS_mn^x ω + ω_mn+ iη

)

(32)

⇒ χ(ω) = ie

¯ hme

X

n>m

(ρ_nn− ρ_mm) ωnm

p^y_mnS_nm^x

ω − ωmn+ iη + p^y_nmS_mn^x ω + ωmn+ iη

(33)

In all of this derivation a summation of k (crystal momentum) has been considered implicitly, i.e. the Em term correspond to the previously mentioned Bloch energies εnk. In order to actually calculate the response tensor we need to consider a summation over k in the first Brillouin zone. Note also that all quantities in the summation carry the k dependence of |nki. In equation 34 the k dependence is incorporated in the relation:

χ(ω) = ie

¯ hm_e

X

k,n>m

(ρ_knn− ρ_kmm) ω_knm

p^y_kmnS_knm^x

ω − ω_kmn+ iη + p^y_knmS_kmn^x ω + ω_kmn+ iη

. (34)

Density functional theory is then used to derive the solid state Hamiltonian of the system.

In this Hamiltonian the structure of the material comes into play and the symmetries are of particular importance. If the crystal has complete inversion symmetry there also arises a symmetry in the Hamiltonian which completely cancels out the Edelstein effect. Instead, in order for the Edelstein effect to be present we need a material with either locally or globally broken inversion symmetries. In this computation of the Hamiltonian, Bloch’s theorem is also of great importance as it simplifies the computation. The eigenstates that are produced by the Hamiltonian are then used to compute p, S, ρ and ωkmn. This of course means that the indexes in equation 34 are the same as the band index in the Bloch states. Additionally, note that we do not account here for the intraband contribution. There are however some hints that these contributions might be zero,⁹ but further numerical and/or symmetry based investigations are needed in this direction.

(16)

1.3 Program

This section describes each part of the EDEFF program calculating the χ tensor (it can be divided into four parts) and describes in a general way the role each part played in the calculation. The calculations are done in Rydberg atomic units (R.a.u.) where the induced spin polarization is measured in units of the Bohr magneton, µB, where one Bohr magneton corresponds to a complete switch of the polarization, and the electric field is measured in units of ERy/a₀e, where a0 is the Bohr radius, e the electron charge and ERy

the energy in Rydberg atomic units. Consequently, the χ tensor is measured in units of µBa₀e/E_Ry and calculated versus the photon energy of the electric field measured in electron volts.

1.3.1 Reciprocal space discretization

As previously mentioned, to calculate the linear response (χ tensor) an integral over only the first Brillouin zone is required. This can be done numerically by creating a mesh in this volume in k-space and approximating the integrals with a weighted sum on these points. Since for symmetry reasons it is possible to restrict the Brillouin zone to a smaller portion of the k-space, that contains exactly the same information, we use a routine that allows us to choose a specific part of the Brillouin zone where to perform the integration.

After which we can use mentioned symmetries to enlarge the solution to the whole Bril- louin zone. In addition to choosing the integration volume the method includes a choice of number of k-points to be used. Using the tetrahedron method,¹⁰ the chosen part of the Brillouin zone is divided into smaller volumes, tetrahedrons, which correspond linearly to the number of k-points. The method includes a choice to divide the first tetrahedrons in thirds instead of the standard division by two. For example, a structure with 4-1 subdivisions would represent a total of four divisions where the first division is done in thirds and the remaining three in twos. Whereas 4-0 subdivisions corresponds to four divisions in twos.

By varying the number of divisions, and thus the number of k-points, a convergence study was done for a given choice of integration volume. These choices for iron, copper and manganese silicide respectively was called BCC (1/48th of the Brillouin zone), BCT (1/16th of the Brillouin zone) and ORTH2 (1/4th of the Brillouin zone). The options for iron and copper was complemented with BCCT and BCTT where the last T stands for total, i.e. integration over the whole Brillouin zone. These use a larger part of the Brillouin zone than the standard BCC and BCT options for the calculations and thus will require more computing power but might result in more accurate results.

1.3.2 Electron density and potential calculation

To approximate the electron density and potential of the system, which is used to formulate the Hamiltonian, the augmented spherical wave (ASW) method is used. The method is based on DFT, described in section 1.2.6, and provides the basis set ψnk that is used by the linear response program.

(17)

1.3.3 Edelstein effect calculation

The state properties of the system are calculated by solving the single-body Hamiltonian, given by the ASW method, on the k-points created by the discretization routine. These properties are in turn implemented in equation 34 and the off diagonal elements of the χ tensor, the interband contribution, versus the photon energy of the electric field is calculated. This routine includes an option to rescale the effect of the spin-orbit coupling in the calculations. How the Edelstein effect depends on the magnitude of the spin-orbit coupling might suggest whether the effect is contributing to the results.

1.3.4 Refining results

A method to rescale and interpolate the resulting calculations is used. The rescaling takes into account the integration volume used in the discretization routine, enlarging the result to the whole Brillouin zone, so that it is comparable for different structures and materials. The interpolation of the data points allows for better handling of the results, especially plotting.

1.4 Method

A modified version of the MOKE program,² called EDEFF (the main difference is in the implementation of the Kubo formalism), was ran on the Rackham cluster on UPPMAX.

The settings in the program chosen by us included the integration volume, the discretization in this volume, the magnitude of the of spin-orbit coupling and the rescaling of the result depending on the integration volume used. Initially, we ran the discretization routine where the integration volume and number of subdivisions where chosen. The electron density and potential calculations was already done by our supervisors and included in the material specific code given to us. The main code was run, solving the Hamiltonian and calculating the specific elements of the χ tensor. When the calculations finished the refining method was ran (with a rescaling factor modified corresponding to the integration volume used).

To determine how the calculated χ values depend on the discretization a convergence study was done comparing the results of the same material and integration volume calculated with multiple different discretizations. For most calculations a maximum of 5-1 subdivisions were used except for iron with the BCCT options were 6-0 and 6-1 subdivisions were used as well. This was done because the convergence of these calculations showed to be of particular interest, which will be discussed later. It was also noted that the calculations with 5-0 and 5-1 subdivisions for manganese silicide exceeded the time limit on the Rackham cluster, there were indications however that using 4-1 subdivisions was sufficient. The study of spin-orbit coupling dependency was done for iron, copper and manganese silicide with the respective options BCCT, BCTT and ORTH2. For these calculations 5-1 subdivisions was used for iron and copper and for manganese silicide 4-1 subdivisions were used. The χ tensor was compared for each material with the spin-orbit coupling scaled to 0, 1 and 10 times the normal magnitude (non present, normal and greatly enhanced spin-orbit coupling). This is of interest since the Edelstein effect is a consequence of the spin-orbit coupling and expected to depend linearly on its magnitude.

(18)

2 Results

In this section the results of the calculations are presented and compared for different materials, integration volumes, discretizations and spin-orbit coupling scaling. The graphs show the real part of the interband contribution values of the χ tensor plotted versus the photon energy of the electric field in electron volts.

2.1 Iron

The results for iron using the BCC option in the discretization routine with 5-1 subdivisions is shown in figure 7. This configuration corresponds to a total of 442368 k-points in the integration volume. We see a response of magnitude around 10^-3 (R.a.u.) except for a region between 2 and 3 eV where the value grows. The graph for iron with the BCCT options with 6-1 subdivisions (84934656 k-points) has a similar shape but with two or- ders of magnitude lower value, as seen in figure 8. The assumption that the Edelstein effect is zero for iron and this result suggests that using symmetries in the Brillouin zone to calculate the effect is not synonymous with using the whole zone, as the calculations taking the whole Brillouin zone into account perform better.

0 1 2 3 4 5 6 7

Photon energy [eV]

-14 -12 -10 -8 -6 -4 -2 0 2 4

[Ba0e/ERy]

10^-3

Figure 7: The calculated element of the χ tensor plotted versus photon energy for iron with the BCC option with 5-1 subdivisions.

(19)

0 1 2 3 4 5 6 7 Photon energy [eV]

-8 -6 -4 -2 0 2 4

[Ba0e/ERy]

10^-5

Figure 8: The calculated element of the χ tensor plotted versus photon energy for iron with the BCCT option with 6-1 subdivisions.

2.2 Copper

Similar calculations as the ones for iron were done for copper. In figure 9 the results for the BCT option with 5-1 subdivisions corresponding to 1327104 k-points is seen. Com- paring this with the results for the calculations done in the whole Brillouin zone, which can be seen in figure 10, further enhances the notion that the use of symmetries in the Brillouin zone results in a worse solution. For copper with a BCTT option with 5-1 subdivisions (21233664 k-points) we see an order of magnitude of 10^-19 (R.a.u.) while for regular BCT the magnitude is around 10^-13 (R.a.u.). Both are however very small and considered to be zero.

0 1 2 3 4 5 6 7

Photon energy [eV]

-4 -3 -2 -1 0 1 2

[Ba0e/ERy]

10^-13

Figure 9: The calculated element of the χ tensor plotted versus photon energy for copper with the BCT option with 5-1 subdivisions.

(20)

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5

[Ba0e/ERy]

10^-19

Figure 10: The calculated element of the χ tensor plotted versus photon energy for copper with the BCTT option with 5-1 subdivisions.

2.3 Manganese silicide

In figure 11 we see the results for manganese silicide with the ORTH2 option with 4-1 subdivisions (663552 k-points). For low photon energies we see a larger order of magnitude, around 10^-3 (R.a.u.), compared to both the corresponding values for iron and copper with the BCCT and BCTT options respectively, as expected if the program works as intended. In the result of the calculations of the individual atoms for manganese silicide, figure 12, we note that the values of each atom varies a lot but are of similar magnitude as the results calculated with the whole unit cell. This indicates that the Edelstein effect is present in all atoms and that they contribute with a similar amount in the unit cell.

0 2 4 6 8 10 12

Photon energy [eV]

-8 -6 -4 -2 0 2 4 6 8 10

[Ba0e/ERy]

10^-4

Figure 11: The calculated element of the χ tensor plotted versus photon energy for manganese silicide with the ORTH2 option with 4-1 subdivisions.

(21)

0 2 4 6 8 10 12 Photon energy [eV]

-1.5 -1 -0.5 0 0.5 1 1.5

[ Ba 0e/E Ry]

10^-3

Atom 1 Atom 2 Atom 3 Atom 4 Atom 5 Atom 6 Atom 7 Atom 8

Figure 12: The calculated element of the χ tensor plotted versus photon energy for each atom in the unit cell of manganese silicide.

2.4 Convergence

In figures 13 and 14 the results of the convergence study are shown (the four smallest discretizations are plotted here) for iron with the BCC and BCCT options respectively.

For BCC it is clear that the shape of the graph is preserved, the χ tensor converges towards non-zero values, for increased discretization. This is not the case for BCCT however, we note that the magnitude for these calculations is decreasing with more k- points, indicating that the correct solution is approached by increasing the discretization.

0 1 2 3 4 5 6 7

Photon energy [eV]

-14 -12 -10 -8 -6 -4 -2 0 2 4

[Ba0e/ERy]

10^-3

16384 k-points 55296 k-points 131072 k-points 442368 k-points

Figure 13: The calculated element of the χ tensor plotted versus photon energy for iron with the BCC option with varying subdivisions.

(22)

-1.5 -1 -0.5 0 0.5 1 1.5 2

[Ba0e/ERy]

10^-4

Figure 14: The calculated element of the χ tensor plotted versus photon energy for iron with the BCCT option with varying subdivisions.

Once again similar results as the ones for iron are achieved for copper, which can be seen in figures 15 and 16. The results with the BCT option quickly converges towards non-zero values while the results with BCTT decreases with increased discretizations.

Although, as mentioned earlier, the order of magnitude for copper is much smaller than for iron. The convergence study for both iron and copper indicates once again that the use of symmetries in the Brillouin zone causes the solution to not converge to zero as it should. While using the whole Brillouin zone in the calculations yields a more realistic solution that converges towards zero.

0 1 2 3 4 5 6 7

Photon energy [eV]

-4 -3 -2 -1 0 1 2

[Ba0e/ERy]

10^-13

Figure 15: The calculated element of the χ tensor plotted versus photon energy for copper with the BCT option with varying subdivisions.

(23)

-1 -0.5 0 0.5 1 1.5 2

[Ba0e/ERy]

10^-18

Figure 16: The calculated element of the χ tensor plotted versus photon energy for copper with the BCTT option with varying subdivisions.

For manganese silicide the results of the convergence study are shown in figure 17. Unlike for the BCCT and BCTT studies for iron and copper the results for manganese silicide with the ORTH2 option converges towards a non-zero shape. Since these calculations are done in a relatively large part of the Brillouin zone (1/4th), compared to the BCC and BCT options (1/48th and 1/16th), the error caused by an extensive use of symmetries in the Brillouin zone is likely not present, to the same degree, in this result. Providing a sign that the program is able to recognize the Edelstein effect in manganese silicide.

0 2 4 6 8 10 12

Photon energy [eV]

-8 -6 -4 -2 0 2 4 6 8 10 12

[Ba0e/ERy]

10^-4

Figure 17: The calculated element of the χ tensor plotted versus photon energy for manganese silicide with the ORTH2 option with varying subdivisions.

(24)

2.5 Spin-orbit coupling dependency

As previously mentioned the spin-orbit coupling dependency was determined with the BCCT, BCTT and ORTH2 options. For iron and copper (BCCT and BCTT) 5-1 subdivisions was used and for manganese silicide 4-1 subdivisions was used. The results for iron and copper can be seen in figures 18 and 19, where we note that the Edelstein effect with 10x scaling does not seem to follow the shape of the graph with normal (1x) spin- orbit coupling. The calculated values do not depend on the magnitude of the spin-orbit coupling in any obvious way. To be compared to the corresponding results for manganese silicide, seen in figure 20, where we note that the effect with good approximation scales linearly with the increased spin-orbit coupling magnitude. This reinforces that the program identifies the non-zero Edelstein effect in the calculations for manganese silicide but not for iron or copper. In addition to this we see that the program treats the Edelstein effect as a direct consequence of the spin-orbit coupling, as intended, since the χ tensor values are zero when the magnitude of the spin-orbit coupling is set to zero, for all materials.

0 1 2 3 4 5 6 7

Photon energy [eV]

-14 -12 -10 -8 -6 -4 -2 0 2 4 6

[Ba0e/ERy]

10^-5

0x SOC scaling 1x SOC scaling 10x SOC scaling

Figure 18: The calculated element of the χ tensor plotted versus photon energy for iron with the BCCT option with varying influence of spin-orbit coupling.

(25)

-3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5

[Ba0e/ERy]

10^-18

Figure 19: The calculated element of the χ tensor plotted versus photon energy for copper with the BCTT option with varying influence of spin-orbit coupling.

0 2 4 6 8 10 12

Photon energy [eV]

-3 -2 -1 0 1 2 3 4 5 6

[Ba0e/ERy]

10^-3

Figure 20: The calculated element of the χ tensor plotted versus photon energy for manganese silicide with the ORTH2 option with varying influence of spin-orbit coupling.

(26)

3 Discussion

As mentioned in section 1.2.8 we do not expect the Edelstein effect to appear for iron or copper. Nonetheless, the χ tensor for those materials assumes non-zero values regardless of which integration volume is used. For the iron BCCT and copper BCTT calculations this appears to be due to a numerical error that decreases with increased discretization.

This is clear from observing figures 14 and 16 as they both seem to converge towards zero. The magnitudes of the values for copper are better in these convergence studies but only by a constant factor as their dependence on the discretization seems to be the same. The results of the convergence studies for iron BCC and copper BCT in figures 13 and 15 are significantly worse. It is hard from these graphs to determine how they depend on the discretization but what we can say is that they do not converge towards zero as expected. This is due to the portion of the Brillouin zone taken into account not containing enough information for the expression we are calculating. The reason for this would be that the spin-orbit coupling breaks some symmetry, so when we apply these symmetries and enlarge the results to the entire zone we get another value that the program converges to. This argument is further reinforced by the fact that the BCCT and BCTT options do converge to the correct value. An obvious solution to this is to use as large part of the Brillouin zone as possible (preferable the whole zone) when doing the Edelstein effect calculations as they appear to give the most credible results.

In any regard, using the best performing settings (the largest possible integration volume with smallest discretization) we see a clear tendency towards low values of the Edelstein effect for iron and copper. In figure 11 the result for manganese silicide with 4-1 subdivisions show that for low photon energies the Edelstein effect is significantly stronger than for iron and copper with χ values peaking at 10^-3 (R.a.u.). This is also the case for the individual atoms in manganese silicide as can be seen in figure 12. The magnitude of the response tensor in the individual atoms is of a similar magnitude as the unit cell as a whole, especially for low photon energies. Confirming that all atoms in the unit cell contributes with a similar amount. Worth noting concerning the calculations for manganese silicide is also how close the results with varying discretizations are to each other, as can be seen in figure 17. Unlike for iron or copper the calculations for manganese silicide quickly converge towards specific non-zero values. Indicating that the Edelstein effect is noticed by the program and contributes to the result for manganese silicide.

As mentioned in section 1.2.5, spin-orbit coupling plays a crucial part of what causes the Edelstein effect. So consequently an interesting feature in the program is how the calculated response tensor depends on the magnitude of the spin-orbit coupling, or the absence of it. An obvious confirmation that the program treats the existence of the spin-orbit coupling as intended is seen in figures 18, 19 and 20 where the χ tensor for iron, copper and manganese silicide respectively is zero when the spin-orbit coupling is set to zero. The more interesting comparison is when its influence in the calculations is increased, in our cased scaled to ten times the normal effect. For iron this causes the χ tensor to oscillate in general, but not increase in magnitude. Thus, the shape of the graphs is not preserved with increased spin-orbit coupling as would be expected if the values were due to the Edelstein effect. The same is true for copper, as seen in figure

Theoretical understanding and calculation of the Edelstein effect

Examensarbete 15 hp Juni 2017

Theoretical understanding and calculation of the Edelstein effect

Gustav Eriksson

Hampus Nyström

Abstract

Theoretical understanding and calculation of the Edelstein effect

Populärvetenskaplig sammanfattning

Contents

1 Introduction

1.1 Background

1.2 Theory

1.3 Program

1.4 Method

2 Results

2.1 Iron

2.2 Copper

2.3 Manganese silicide

2.4 Convergence

2.5 Spin-orbit coupling dependency

3 Discussion