Wannier functions from Bloch Orbitals in Solids

Examensarbete 15 hp Juni 2013

Wannier functions from Bloch Orbitals in Solids

Anders Stangel

Kandidatprogrammet i Fysik

Department of Physics and Astronomy

Abstract

Wannierfunctions are a superposition of the Blochorbitals in a Brillouin zone

belonging to a manifold of energy bands. These Wannier functions have several uses regarding the analysis of the crystal on a local level. Since the Bloch orbital has a gauge indeterminacy and the Wannier functions therefore is strongly non-unique, the natural choice is the maximally localized Wannier funcition. These can be calculated from the standard Bloch orbital using unitary transformation by a steepest descent algorithm as proposed by Nicola Marzari and David Vanderbilt. Here the argument for this algorithm is discussed and explained.

Sammanfattning

I en perfekt kristall finns upprepade celler som ser likadana ut. I en sådan kristall kan elektronerna beskrivas med hjälp av Blochorbitaler. Dessa Blochorbitaler är dock inte begränsade i rummet utan har spridning i hela kristallen. Ett annat sätt att beskriva elektronerna på är då Wannierfunktionerna, som fås genom att lägga ihop flera

Blochorbitaler så att de får en begränsad utsträckning i rummet. Wannierfunktionerna är dock kraftigt icke-unika eftersom de beror på Blochorbitalerna och det naturliga är att välja de Wannierfunktioner som har så liten utsträckning i rummet som möjligt.

Dessa kan sedan användas för en analys av kristallen på lokal nivå. I detta arbete sammanfattas och förklaras de argument som föreligger för en metod att framställa dessa Wannierfunktioner.

Contents

Background 1

Definition of Wannier functions 1

Smooth Bloch orbitals through projection 6

Measuring the spread of the Wannier function 10

Maximally localized Wannier functions 22

The algorithm for calculating maximally localized Wannier Functions 27

Example: crystalline Si 30

Discussion 31

References 32

In this work a number of expressions are taken from other works, they are denoted by asterisks after the expression number as follows:

* This expression is taken from (Marzari, Vanderbilt (2012))

** This expression is taken from (Marzari, Vanderbilt (1997))

*** This expression is taken from (Blount (1962))

1

Background

When dealing with solids with periodic boundary conditions the electronic ground state can be expressed in terms of Bloch orbitals, labelled by their position in reciprocal space and their energy band or energy bands. These Bloch orbitals are localized in energy in the sense that they are the eigenstates to the Hamiltonian. The Bloch orbitals have an indeterminacy regarding both the overall phase and the choice of gauge, meaning that any arbitrary unitary transformation of the Bloch orbitals will yield the same physics (Marzari, Vanderbilt (2012) p.1420).

This can then be used to express the Wannier functions, a superposition of Bloch orbitals calculated through a series of unitary transformations introduced by Gregory Wannier in 1937. These Wannier functions were used as a starting point in calculations of, for example, semi-classical electron dynamics (Marzari, Vanderbilt (1997) p.12 847) but the actual Wannier functions were rarely required to actually be calculated. This changed when the Wannier functions was found to be useful in, for example, calculations of model Hamiltonians such as tight-binding Hamiltonians as well as providing an insight into the chemical bonding of the crystals (Marzari, Vanderbilt (2012) p.1421). These Wannier functions are not eigenstates of the Hamiltonian and as such trade localization in energy for localization in space (Marzari, Vanderbilt (2012) p.1420), and are labelled by their cell number and their energy band or energy bands.

Like in chemistry where similar orbitals, called localized molecular orbitals, are used to find a description that belongs to each individual molecule and are

approximately translational images of the localized molecular orbitals of other

molecules (Boys (1960) p.296), the Wannier functions can be found that belongs to a single cell in the solid, where the Wannier function of each cell is the translational image of the Wannier functions of the other cells in the crystal (Marzari, Vanderbilt (2012) p.1420).

These Wannier functions are strongly non-unique due to the indeterminacy of the Bloch orbitals and the natural choice is to choose the Wannier functions that are maximally localized. A steepest descent approach as an iterative method for

calculating the maximally localized Wannier functions was introduced by Nicola Marzari and David Vanderbilt in 1997.

Definition of Wannier Functions

For a perfect crystal, there is a periodicity in the reciprocal lattice so that in a point in the lattice, the electric potential should look the same as the corresponding point in an adjacent Brillouin zone, meaning that the electric potential also has the same periodicity as the lattice. If the translation operator TR is then defined as the operator that translates the Hamiltonian exactly one cell in some direction r then the

Hamiltonian with the translation operator acting on it will be the same as without the operator and also:

= , = − = 0 (1)^*

2 So the Hamiltonian and the translation operator commute. And the electric potential has periodicity R as shown in figure 1.

figure 1: an example of a periodic potential in a section of a one dimensional crystal where the red dots are the position of the ions and the lines indicate the electric potential.

The Bloch function is then defined as:

= ^∙ (2)^*

where has the same periodicity as the electric potential, r is the location in real space, k is the wave vector in the Brillouin zone, n is the band number and ^∙ is called a envelope function (Marzari, Vanderbilt (2012) p. 1421).

The states | is defined for every n so that | is the eigenstates of the Hamiltonian. The envelope function constructs a different wave function for every k and, due to the orthogonality of the sine and cosine functions and thus of the Bloch functions we should be able to construct other wave functions by superposition of two or more Bloch functions. One such possible wave function that should be expected is a localized wave packet. To construct these we start by superpositioning every Bloch function in the Brillouin zone.

= (3)

where is called a Wannier function and A is a normalization factor (Marzari, Vanderbilt (1997) p.12 848).

The Brillouin zone is periodic so that is the same for every Brillouin zone. It therefore makes sense to normalize the integral over the entire Brillouin zone. The volume of the Brillouin zone is:

3 = 2" ^#

$ =1

(4)

where G is the volume of the Brillouin zone and V is the volume of the real space primitive cell.

The Wannier function then becomes

= $

2" ^# (5)^*

More Wannier functions can be constructed by inserting a phase factor ^{& ∙} for some real vector R into the integral since the Wannier functions in different cells are translational images of each other. We can then get Wannier functions like

' , ₍ , …

= $

2" ^{# & ∙} (6)^*

and in Dirac bra-ket notation.

| = $

2" ^# | ^{& ∙} = | * (7)^**

Where | * = | for every value of n and R. This is a more convenient notation for the Wannier function in the Dirac bra-ket notation (Marzari, Vanderbilt (1997) p.12 848).

Note that the Bloch functions are normalized to one Brillouin zone and form an orthogonal set.

+ | _, ´- ≝ 011 ^∗∗ _, ´ 23456

= 2" ^#

$ 7 ^,7# ´

(8)^*

where 7 is the Dirac delta function.

4 The same calculations then for the Wannier functions yield:

+ *| ´8- = 9 $

2" ^# : | ^∙ ; 9 $

2" ^# | _, ^{& ∙ ´} ;

= < $

2" ^#=⁽ : | _, ^{∙ & ∙ ´}

= < $

2" ^#=⁽∗ 2" ^#

$ 7 ^{, ∙ & ∙ ´} = $

2" ^#7 _,∗ 2" ^#

$ 7^{# ´}

= 7 _,7# ´ (9)

and the Wannier functions also form an orthogonal set.

| * has the form of a Fourier transform, and the inverse can thus be constructed.

| = ^∙ | * (10)

This is the general formula for continuous R. But R are the cells in the crystal and as such should be seen as a discrete mesh and the Bloch functions then becomes:

| = > ^∙ | *

(11)^*

Note here that there is a gauge freedom in how to choose the Bloch functions

|~ = ^@| |~ ^∙ = | ^{@ ∙} |~ = | ^@

(12)

where A is a real function that has the same periodicity as the Hamiltonian.

This gauge freedom does not change the physics of the system. By choosing a convenient gauge the Bloch function can be made to be smooth so that it is partially differentiable in every point, that is, ∇_C| ~ is defined for every k. This will change the Wannier functions into,

5

| * = $

2" ^# | ^{@ & ∙} (13)^*

which will not be the same function even though it is still a valid Wannier function.

The Wannier function above is given for a specific band n on the Bloch

functions. If the Wannier function is instead taken to be a superposition of J Wannier functions from several bands.

| * = > $

2" ^# D₁ | ₁ ^{@ & ∙}

E 1F'

(14)^*

where D_, is a square matrix of size J with the periodicity of the Hamiltonian since the physics of the system stay the same under a unitary transformation of the Bloch orbitals.

However, this should not change the normalization given earlier.

G~ H~_,

´I = > : ₁ |D₁ ^∗D_1,| ₁

E 1F'

= 2" ^#

$ 7 ^,7# ´ > D₁ ^∗D_1,: ₁ | ₁

E

1F'

= >: ₁ | ₁

E 1F'

D₁ ^∗D_1, = J ∗ 7_, (15)

where I is the identity matrix of size J, and thus the D₁ matrix is unitary. Note that for |~ n is no longer an isolated band, but rather a manifold of J bands.

The k-space of the Brillouin zone can now be seen not as a continuum but as a discrete mesh consisting of a number of uniformly distributed points. The Wannier functions of equation (6) can then be written as,

| = > ^∙ | * K | * = 1

L > ^{& ∙} | (16)^*

where N is the number of discrete points in the Brillouin zone.

6 Using this convention means that N is treated like the volume of the Brillouin zone, since the mesh points are uniformly distributed and a larger volume therefore means more mesh points. It is also possible to use the convention,

| = 1

√L> ^∙ | * K | * = 1

√L > ^{& ∙} | (17)^*

This gives,

+ | _, ´- = 1

L >: *| ^{& ∙} > ^´∙ | 8 =L

L + *| 8-7^, 7# ´ = 7_, 7# ´

(18)

But | is still normalized to the Brillouin zone with N discrete points.

+ | - = : | ∗ | = L (19)

where A is a normalization constant.

( = L | = √L| ^∙ (20)

Smooth Bloch orbitals through projection

In order to find the maximally localized Wannier function, the Bloch orbitals should be smooth. One way of creating smooth Bloch orbitals is via projection.

In the previous section the Wannier functions were denoted by the energy bands * where * is a manifold of energy bands. In order to project the Wannier functions onto these bands the projection operator is defined.

N = | : | (21)

where A is a normalization constant (Cloizeaux (1964) p.A685).

7 But is still normalized over the entire Brillouin zone and N should normalize to 1 (Cloizeaux (1964) p.A685).

G~ H~

, ´I = ∗ 2" ^#

$ 7 ^,7^# _´= 7 _,7^# _´ = $ 2" ^# N = $

2" ^# | : | = $

2" ^# P> ^∙ | * > : *| ^{& ∙} Q

(22)

but | * forms an orthogonal set.

N = $

2" ^# P> ^∙ | * : *| ^{& ∙} Q = $

2" ^# P> | * : *|Q

= $

2" ^# > | * : *| = $

2" ^# > | * : *| 2" ^#

$

= > | * : *| (23)

This projection operator is invariant to the multiple band choice in the sense that,

N = $

2" ^# |~ :~ | = $

2" ^# D₁ | :~ |D₁ ^∗

= $

2" ^# J ∗ | : | = $

2" ^# | : | (24)

And treating the Brillouin zone like a uniform mesh the projection operator becomes

N = $

2" ^#> | : | (25)^*

Due to the gauge freedom of the Bloch orbitals, there exists a gauge | ^@ so that ∇_C|~ is defined for every and the Bloch orbitals become smooth (Marzari, Vanderbilt (2012) p.1422). In order to find this gauge a trial function R = |R is created. This function consists of a smooth set of J localized orbitals in the location where the Wannier function is expected to be centred (Marzari, Vanderbilt (2012) p.1424) where the manifold of J bands are isolated. Using the projection operator on this trial function yields:

8 N|R = $

2" ^#>| + |R - (26)^*

But n is a manifold of J bands.

N|R = $

2" ^# > >| _, + _, |R -

E ,F'

(27)

And by this a new Bloch orbital can be defined as

>|S = 2" ^#

$ N|R = > >| ^, + _, |R -

E ,F'

(28)

Which can be defined in a single point on the grid.

|S = >| _, + _, |R -

E ,F'

(29)^*

The scalar product +S |S - can then be calculated.

+S, |S - = >TR,U _V W: _V |

E VF'

>| ₁ + ₁ |R -

E

1F'

= > >TR_,U _V W

E 1F'

T _V U ₁ W+ ₁ |R -

E

VF' (30)

But ₁ is orthogonal as stated before.

+S, |S - = > >TR,U _V W

E 1F'

T _V U ₁ W+ ₁ |R -

E

VF'

= > >TR_,U _V W

E 1F'

2" ^#

$ 7^V1+ ₁ |R -

E VF'

= 2" ^#

$ >TR^,U _V WT _V UR W

E

VF' (31)

9 This scalar product is over all space. The factor ^(XY

Z is the volume of one primitive cell in the reciprocal space and it can be removed by instead defining the scalar product over one cell.

+S, |S -Z ≝ S_, ^∗ = >TR_,U _V WT _V UR W

E VF'

(32)

But TR U _V W is the scalar product for two manifolds n and j that can be described in a * × \ matrix and through singular value decomposition can be decomposed as

T _V UR W = TR U _V W^† = ^ _V_ _V`^† _V (33)^*

where Z and W are unitary matrices and D is real and diagonal (Klema, Laub (1980) p.166) giving

_ = _^†

(34)

Then the scalar product of +S_, |S -_Z is given as +S, |S -Z = >TR_,U _V WT _V UR W

E VF'

= `<sub>,V</sub>_<sub>,V</sub><sup>†</sup> ^<sub>,V</sub><sup>†</sup> ^ <sub>V</sub>_ <sub>V</sub>`^† _V (35)

But ^ is unitary

^_,V^† ^ _V = J ∗ 7_,

(36)

And the scalar product becomes

+S_, |S -_Z = `<sub>,V</sub>_<sub>,V</sub><sup>†</sup> ^<sub>,V</sub><sup>†</sup> ^ <sub>V</sub>_ <sub>V</sub>`^† _V7 _,= ` <sub>V</sub>_<sup>†</sup><sub>V</sub>_ <sub>V</sub>`^† _V = ` <sub>V</sub>_<sup>(</sup><sub>V</sub>`^† _V

(37)

10 and since ` is unitary it is given for any given value of n and j equation (35)

becomes

+S |S -Z = _⁽_V

(38)

Equation (29) can then be expressed in terms of these matrix indices.

|S = >| _V T _V UR W

E VF'

= >| _V

E VF'

^ _V_ _V`^† _V (39)

It can then be clearly seen that

|S ∗ +S_, |S -_Z ^&'(= >| _V

E VF'

^ _V_ _V`^a _V∗ _^&'_V = >| _V

E VF'

^ _V`^† _V (40)

But Z and W are unitary, so this is a unitary transformation of | _V and we can define a new Bloch orbital.

|~ = >|S_, ∗ +S_, |S -_Z^&'(

E ,F'

= > >| _V T _V UR_,W

E ,F'

∗ +S_, |S -_Z^&'(

E VF'

(41)^*

which has a smooth gauge in since it consists of smooth orbitals projected onto a manifold of smooth energy bands and is given by a unitary transformation from | _V .

|~ is uniquely defined by the choice of the trial function |R but since the transformation is unitary they should still yield valid Wannier functions.

Measuring the spread of the Wannier function

In order to find the Wannier function that is maximally localized in the set of all

possible Wannier functions using the calculation above it must first be defined what it means to be maximally localized. To do this, a localization functional can be defined.

Ω = > +c*|d⁽|c*- − +c*| |c*-⁽ (42)^**

11 Where Ω is the localization functional (Marzari, Vanderbilt (1997) p.12 849) and is the location in real space.

This is a measurement of the spread of the Wannier function centralized

in = c. The maximally localized Wannier function is therefore the Wannier function where Ω has a minimum. But the Wannier function is not invariant to the gauge of the Bloch orbital. One way of treating the localization functional is therefore to transform it into a function of Bloch orbitals rather than as a function of Wannier functions. That way the Wannier function doesn’t have to be calculated until after the Bloch orbitals associated with the maximally localized Wannier function has been found. To accomplish this, the localization functional can also be written as:

Ω = > P+c*|d⁽|c*- − +c*| |c*-⁽ − > |+ 8| |c*-|⁽

ec

+ > |+ *| |c*-|⁽

ec

Q (43)

If this localization functional can be shown to be positive definite then there exists a global minimum associated with the maximally localized Wannier function.

The second term of expression (43) can be written as:

+c*| |c*-⁽ = g> < $

2" ^#=⁽ : ₁ |D₁ ^{∗ & @} ∗ ∗ D₁ | ₁ ^@

E 1F'

h

(

(44)

and is therefore positive definite. The localization functional can also be written as:

Ω = > P+c*|d⁽|c*- − >|+ *| |c*-|⁽+ >|+ *| |c*-|⁽

ec

Q (45)

The middle term of this can now be achieved by using the projection operator.

>+c*|d_iNd_i|c*-

i

= > jc*kd_i> | * : *|d_ikc*l

i

= > >+c*|d_i| * : *|d_i|c*-

i

= > >+ *|d_i|c*-^a+ *|d_i|c*-

i

= > >|+ *|d_i|c*-|⁽

i

= >|+ *| |c*-|⁽ (46)

where m is the dimensions for d and N is the projection operator.

If the first term is then expressed as

12

>+c*|d⁽|c*- = >+c*|d_i∗ d_i|c*-

i (47)

the first two terms can be defined as Ω_n = >+c*|d_i 1 − N d_i|c*-

i

= >+c*|d⁽|c*- − >|+ *| |c*-|⁽ (48)^**

Which can be expanded into Ω_n = >+c*|d_i 1 − N d_i|c*-

i

= >+c*|d_i⁽|c*-

i

− >+c*|d_iNd_i|c*-

i (49)

that has terms equal to the diagonals on the matrices:

>+c*|d_i⁽|c*-

i

= > d₅ Nd_i⁽

i (50)

and

>+c*|d_iNd_i|c*-

i

= > d₅ Nd_iNd_i

i (51)

where d₅ is the trace per unit cell of the matrices, defined as the sum of the diagonal elements per cell of the matrices.

This yields

Ω_n = >+c*|d_i⁽|c*-

i

− >+c*|d_iNd_i|c*-

i

= > d₅ Nd_i⁽− Nd_iNd_i

i

= > d₅ Nd_i 1 − N d_i

i (52)^*

But these matrices can also be written as:

Nd_i 1 − N d_i = Nd_i⁽− Nd_iNd_i= Nd_i− Nd_iNd_i+ Nd_iNd_i− Nd_iNd_i

= Nd_i⁽− 2Nd_iNd_i+ Nd_iNd_i (53)

13 But the projection operator is idempotent (Marzari, Vanderbilt (2012) p.1425) so that

NNd_i= $

2" ^#> | : | > | : |d_i

= $

2" ^#> | + | - : |d_i == $

2" ^#> | : |d_i

= Nd_i (54)

which is also intuitively true considering projection onto a manifold of energy bands of a function that is already projected onto the same manifold should yield the same projection. The same therefore holds true for 1 − N since it acts as a projection operator on every band that N does not project on. The projection operator is also a self-adjoint operator so that:

N = > | * : *| = N^† (55)

Therefore:

Nd_i 1 − N d_i= Nd_i⁽− 2Nd_iNd_i+ Nd_iNNd_i = Nd_i 1 − N 1 − N d_i

= Nd_i 1 − N 1 − N d_i N = Nd_i 1 − N oNd_i 1 − N p^†

= qNd_i 1 − N q (56)^**

But N is gauge invariant and qNdi 1 − N q is positive definite. Therefore Ω_n is also gauge invariant and positive definite.

Then ~

Ω can be defined so that:

Ω = Ω_n+ ~Ω ~

Ω = >|+ *| |c*-|⁽

ec (57)^**

This must then be the gauge dependent term that can further be divided into a diagonal term ~

Ω_r and an off-diagonal term ~ Ω_sr.

~Ω = ~

Ω_r+ ~Ωsr (58)^**

14 where

~Ωr = > >|+ *| |c*-|⁽

ec (59)^**

and

~Ω_sr = > >|+ 8| |c*-|⁽

e, (60)^**

Due to the absolute value, both ~

Ωr and ~

Ωsr are positive definite and thus the localization functional is positive definite.

In order to find the Bloch orbitals associated with the maximally localized Wannierfunction, the localazation functional should then be rewritten in terms of the Bloch orbitals.

Due to the definition of the Wannier function in equation (5), d can be expressed as an operator on |c* using the convention of

| = √L| ^∙ t |c* = > u $

2" ^# | (61)

for some manifold of bands n.

|c* = u $

2" ^# >| = u $

2" ^# >| (62)^***

Because is not a function of (Blount (1962) p.357). But the gradient of the Bloch orbital with respect to is

∇ | = ∇ | ^∙ = v | + ^∙ ∇ |

(63)

Solving this for | yields

| =∇ | − ^∙ ∇ |

v = v ^∙ ∇ | − v∇ | (64)

15 And equation (62) becomes

|c* = u $

2" ^# >|

= u $

2" ^#P v ^∙ ∇ >| − v∇ >| Q (65)^***

But the second term corresponds to a plane of discontinuity (Blount (1962) p.358) and it vanishes since | is continuous everywhere.

|c* = u $

2" ^# v ^∙ ∇ >| (66)

This can then be redone for the d⁽ operator.

d⁽|c* = d⁽u $

2" ^# >| = u $

2" ^# d⁽>| (67)

This can, with the same argument as for |c* be calculated using the second gradient of | .

∇⁽| = ∇ ∙ ∇ | = ∇ ∙ ov | + ^∙ ∇ | p =

= ∇ ∙ ov | ^∙ + ^∙ ∇ | p

= v ov | + ^∙ ∇ | p + v ^∙ ∇ | + ^∙ ∇⁽|

= 2v ^∙ ∇ | + ^∙ ∇⁽| − d⁽| (68)

Solving this for d⁽| yields

d⁽| = 2v ^∙ ∇ | + ^∙ ∇⁽| − ∇⁽|

(69)

which yields an expression for d⁽|c* .

16 d⁽|c* = u $

2" ^# d⁽>|

= u $

2" ^#P 2v ^∙ ∇ > | + ^∙ ∇⁽> |

− ∇⁽ | Q (70)

But the last term corresponds to a plane of discontinuity in the same way as in equation (65) (Blount (1962) p.358), and vanishes because | is continuous everywhere.

d⁽|c* = u $

2" ^#P2 v ^∙ ∇ > | + ^∙ ∇⁽> | Q (71)

Another expression for d⁽|c* can be gotten from the expression for |c* . d⁽|c* = ∙ |c* = ∙ u $

2" ^# v ^∙ ∇ >|

= u $

2" ^# v ^∙ ∇ >|

= u $

2" ^#P2 v ^∙ ∇ > | + ^∙ ∇⁽> | Q (72)

v ^∙ ∇ | = 2v ^∙ ∇ | + ^∙ ∇⁽|

v ^∙ ∇ | = − ^∙ ∇⁽| (73)

17 d⁽|c* = u $

2" ^#P2 v ^∙ ∇ > | + ^∙ ∇⁽> | Q =

= u $

2" ^#P ^∙ ∇⁽> | − 2 ^∙ ∇⁽> | Q =

= −u $

2" ^{# ∙} ∇⁽>| (74)

An expression for + *| |c8- and + *|d⁽|c8- can then be gotten by multiplying these expression with : *| (Blount (1962) p.326).

+ *| |c8- = u $

2" ^#∗ u $

2" ^{# & ∙} : |v ^∙ ∇ >| _,

,

(75)

+ *|d⁽|c8- = u $

2" ^# ∗ u $

2" ^{# & ∙} : | ^∙ ∇⁽>| (76)

And the expression for the matrix elements 8 and * then becomes:

+ *| |c8- = v $

2" ^{# ∙} + |∇ | _, - (77)^***

+ *|d⁽|c8- = − $

2" ^{# ∙} T U∇ ⁽U _, W (78)^***

Since | = | ^∙ is a smooth function, then | is also a smooth function and both ∇ | and ∇⁽| is defined for all .

In order to calculate the gradient the Brillouin zone can be discretized into a

Monkhorst-Pack mesh, defined as a uniform mesh of points with the Bloch orbitals determined at those points (Monkhorst, Pack (1976) p.5188). For each point stars can defined as the neighbours to the points in the previous star (Monkhorst, Pack (1976) p.5188) so that the first star consists of the points first neighbours, the second star consists of the points second neighbours and so on. In another convention,

18 these stars are also called shells (Chadi, Cohen (1973) p.5748) but here they shall be called stars.

If the following condition:

> _wx_ix_y = 7_iy

z (79)^**

can be fulfilled for some number of stars, where z_i is the vector connecting the point to a point m, and _w is a weight that normalizes the condition (Marzari, Vanderbilt (1997) p.12 862), then the gradient of a smooth function { can be approximated as the function in a point going from the point to another point + z, and the

difference in value for the function going between these points so that:

∇ { ≈ > _wzo{ + z − { p

z (80)^*

Notice the similarities to this to the definition of the derivative.

This also yields:

|∇ { |⁽ = k> _wzo{ + z − { p

z

k

(

= k> _wzo{ + z − { p

z

∗ > _wzo{ + z − { p

z

k

= }> g> _w⁽x_i⁽o{ + z − { p

i

( z

+ > > _w⁽

yei

x_ix_yo{ + z − { p^~

i ;} (81)

But due to the condition above the second term is zero and the first term yields:

|∇ { |⁽ = k> P> _w⁽x_i⁽o{ + z − { p

i

(Q

z

k = > _wo{ + z − { p

z

(

(82)^*

Now using the convention from equation (17) the second term of Ω_n can be gotten.

19 + *| |c8- = v

L > ^∙ : |

,z

wz | _{, •z} − | _,

= v

L >_,z ^{∙ w}z + | _{, •z}- − 1 (83)

which gives

+c*| |c8- = v

L > ^wz + | _•z- − 1

,z (84)^**

The other part of Ω_n can also be gotten.

+ *|d⁽|c8- = −1

L > ^∙ : | _w 2 ∗ € | _{, •z} − | _, ⁽

,z

= −1

L >_,z ^{∙ w}€ 2+ | _{, •z}- − 2

= 1

L >_,z ^{∙ w} 2 − 2 ∗ € + | _{, •z}- (85)

Which gives:

+c*|d⁽|c8- = 1

L > ^w 2 − 2 ∗ € + | _•z-

,z (86)^**

The overlap + | _{, •z}- between the Bloch orbitals in neighbouring mesh points can be defined as.

•_,^,z = + | _{, •z}-

(87)^**

But

+ _•z| _•z- = 1 = + | -

(88)

And we can define • ^,z to have the form of a Taylor expansion as:

20

• ^,z = 1 + v‚x +1

2 ƒx⁽+ „ x^# (89)^**

For some ‚, ƒ (Marzari, Vanderbilt (1997) p.12 851).

And since + _•z| _•z- = 1 and the form of the gradients to a function above, both

‚, ƒ ∈ ℝ. This yields.

• ^,z − 1 = v‚x + „ x⁽

(90)^**

which can be simplified through a power series of ln • ^,z .

ln • ^,z = lno1 + v‚x + „ x⁽ p ≈ 1 + v‚x + „ x⁽ J8‰ln • ^,z Š ≈ ‚x + „ x⁽ (91)

for small b. But b is the distance between points in the Monkhorst-Pack mesh that is an approximation of a continuum and as such are small.

A result isolating

2 − 2 ∗ € ‰• ^,z Š = 2 − 2 ∗ 91 +1

2 ƒx⁽+ „ x^# ; = −ƒx⁽+ „ x^# (92)

which yields

+c*| |c8- = v

L > ^wz + | _•z- − 1

,z

= v

L > ^wz‰• ^,z − 1Š

,z

≈ v

L > ^wz v‚x

,z

= −1

L > ^wz ‚x

,z

≈ −1

L > ^wz ∗ J8 ‰ln‰• ^,z ŠŠ

,z

(93)

and

+c*|d⁽|c8- = 1

L > ^w<2 − 2 ∗ € ‰• ^,z Š=

,z

≈ 1

L > ^w −ƒx⁽

,z (94)

21 But

H• ^,z H = 91 +1

2 ƒx⁽+ „ x^# ;⁽+ o‚x + „ x^# p⁽

= 1 + ƒx⁽+1

4 ƒ⁽x^Œ+ ‚x ⁽+ „ x^# = 1 + ƒx⁽+ ‚x ⁽+ „ x^# (95)

Solving for – ƒx⁽ yields

−ƒx⁽ = 1 − H• ^,z H + ‚⁽x⁽+ „ x^# ≈ 1 − H• ^,z H + <J8 ‰ln‰• ^,z ŠŠ=⁽ (96)

which gives +c*|d⁽|c8- = 1

L > ^w −ƒx⁽

,z

≈ 1

L > ^w91 − H• ^,z H + <J8 ‰ln‰• ^,z ŠŠ=⁽;

,z (97)^**

The localization functional can then be expressed in terms of • ^,z .

Ω_n = > P+c*|d⁽|c*- − >|+ *| |c*-|⁽Q

= >1

L P> ^w91 − H• ^,z H + <J8 ‰ln‰• ^,z ŠŠ=⁽;

,z

+ > _wz ∗ J8 ‰ln‰• ^,z ŠŠ

,z ; (98)^*

~Ω_r = > >|+ *| |c*-|⁽

ec

= 1

L > ^w> g−J8 ‰ln‰• ^,z ŠŠ − z

,z

∙ P1

L > ^wz ∗ J8 ‰ln‰• ^,z ŠŠ

,z Qh

(

(99)^*

22

~Ωsr = > >|+ 8| |c*-|⁽

e,

= −1

L > ^w > H•_,^,z H⁽

,e

,z (100)^*

Maximally localized Wannier functions

In order to find the maximally localized Wannier function, the Wannier function in which the localization functional reaches its global minimum is sought.

This can be found by first considering an infinitesimal unitary transformation operator D_, = 7_, + `_,

(101)^**

where W is anti-Hermitian so that

`<sub>,</sub> <sup>†</sup> = − `_,

(102)^**

and D_, is unitary

D_, | = | + > `_, | _,

, (103)^**

The global minima can then be sought by determining the gradient of the localization functional.

9 Ω

`; _,= Ω

`_, (104)^*

But Ω is positive definite and a condition for a minimum is that this gradient vanishes.

Then for some operator B that can be expressed as

Ž = d `Ž

` (105)^**

the superoperators can be defined

23 Ž =oŽ − Ž^†p

2 = € o d `Ž p

` (106)^**

• Ž =oŽ + Ž^†p

2v = J8o d `Ž p

` (107)^**

(Marzari, Vanderbilt (1997) p.12 852) and then the operators can be defined

• ^,z = J8 ‰ln‰• ^,z ŠŠ − z ∙ P1

L > ^wz ∗ J8 ‰ln‰• ^,z ŠŠ

,z

Q (108)^**

€_,^,z = •_,^,z • ^{,z ∗}

(109)^*

~€,

,z = •_,^,z

• ^,z (110)^**

,,z

= ~€,

,z • ^,z =•_,^,z

• ^,z • ^,z (111)^*

• ^,z can then from equation (103) be found to be (Marzari, Vanderbilt (1997) p.12 852)

• ^,z = + | _•z- = −‰ ` • <sup>,z</sup> Š + ‰• <sup>,z</sup> ` Š

(112)^**

but it follows from equation (87) that the hermitian conjugate of • ^,z can be given as

• ^,z = + | _•z- = + _•z| -^†= • •z, &z& †

= • ^{•z,&z †} (113)

24

• ^,z = −‰ ` • <sup>,z</sup> Š + ‰• <sup>,z</sup> ` Š

= −‰ ` • <sup>,z</sup> Š − ‰ ` • ^•z,&z Š^∗ (114)

but

Ω_n+ ~Ω_sr = 1

L > ^w> 1 − H• ^,z H⁽

,z

= 1

L > ^w> 1 −

,z

• ^,z • ^{,z ∗}

= 1

L > ^w> 1 −

,z

€_,^,z = 1

L > ^w> 1 −

,z

€_,^,z (115)

‰Ω_n+ ~Ω_srŠ = P1

L > ^w> 1 − €_,^,z

,z

Q = 1

L > ^w> − €_,^,z

,z

= 2

L > ^w> − `_, €_,^,z

,z

= 4

L > ^w€ < d‰ `€_,^,z Š=

,z

(116)

For the diagonal part ~

Ω_rfirst • ^•z,&z can be calculated as

• ^•z,&z = J8 ‰ln‰• ^{,z ∗}ŠŠ − −z ∙ P1

L > ^w −z ∗ J8 ‰ln‰• ^{,z ∗}ŠŠ

,z

Q

= −J8 ‰ln‰• ^,z ŠŠ − z ∙ P1

L > ^wz ∗ <−J8 ‰ln‰• ^,z ŠŠ=

,z

Q

= −• ^,z (117)

and from this ~ Ω_r is

~ Ω_r = 1

L > ^w> • ^{,z (}

,z (118)

25

~

Ω_r = P1

L > ^w> • ^{,z (}

,z

Q

= 2

L > ^w> • ^,z J8 P− `• ^,z

• ^,z + `• ^{,z ∗}

• ^{,z ∗}Q

,z

= 2

L > ^w> • ^,z J8 <− `~€ <sup>,z</sup> + `~€ ^•z,&z =

,z

= 2

L > ^w> J8 <− `~€ <sup>,z</sup> • <sup>,z</sup> + `~€ ^•z,&z • ^,z =

,z

= 2

L > ^w> J8 <− `~€ <sup>,z</sup> • <sup>,z</sup> − `~€ ^•z,&z • ^•z,&z =

,z

= −2

L > ^w> J8‰ ` <sup>,z</sup> + ` ^•z,&z Š

,z

=

= −4

L > ^w J8 < d‰ ` ^,z Š=

,z

(119)

The gradient of the localization functional then becomes

Ω

`_, = 4

L > ^w

€ < d‰ `€_,^,z Š=

,z `

− 4

L > ^w

J8 < d‰ ` _,^,z Š=

,z `

Ω

`_, = 4 > _w‰ o€ ^,z p − •o ^,z pŠ

z

(120)^*

There is a possibility that this gradient being zero does not correspond to a global minimum but instead corresponds to a false local minimum. This can be

experienced when dealing with isolated bands in fine k-point meshes. In a false local minimum the Wannier functions are typically complex, while at a global minimum they are, apart from an overall phase, always real. At a good global minimum the values of

• ^,z are also such that

H• ^,z H ≪ "

(121)^**

while at a false local minimum some H• ^,z H approaches " (Marzari, Vanderbilt (1997) p.12 855).

Using the projection operator and a trial function to create a set of smooth orthonormal Bloch orbitals, this tend not to be a problem for actual calculations. If the

26 problem occurs, then the trial function was either not smooth enough or placed too far away from the centre of the Wannier functions, and a different choice of trial function will solve the problem.

In order to find the point for which the gradient of the localization functional is zero ` can be defined as an infinitesimal step in one direction of the gradient so that

` = ’ Ω

`_, (122)^**

where ’ is a positive infinitesimal.

This then means that a step Ω in the localization functional can be written as:

Ω= > Ω

`<sub>,</sub> ` (123)^**

but ` is defined to be anti-Hermitian so that

` = ’ Ω

`_, = −’ Ω

`_,

†

(124)

which gives for the step in the localization functional

Ω= > d P Ω

`<sub>,</sub> ` Q = −’ > d P Ω

`_,

Ω

`_,

†Q

= −’ > “ Ω

`_, “

(

(125)

Since both ’ and ”_•–^•Ω_—˜ ”⁽ is positive then Ω is always negative and will decrease the localization functional.

In practical computation it is not possible to treat ’ as infinitesimal. Instead it can be treated as a fixed finite step.

’ = m

4 ∑_{z z} (126)^**

27 where m is a constant (Marzari, Vanderbilt (1997) p.12 854).

Changing the infinitesimal ` into

` t ∆` = m

∑_{z z}> _w‰ o€ ^,z p − •o ^,z pŠ

z (127)^**

But from equation (102)

∆– ∗ ^{∆– †} = ^∆–∗ ^&∆– = 1

(128)

and ^∆– is unitary.

To find a value for m the factor _{∑ ›}ⁱ

œ

z should be viewed as an approximation of an infinitesimal. It can then be expected that the smaller m becomes, the more accurate the calculation will be. For the study of perfect cubic and hexagonal crystals

structures, m = 1 is enough (Marzari, Vanderbilt (1997) p.12 854).

The Bloch orbital can then be updated by ^∆–. The overlap then becomes updated with

•_,^,z = + | _•z- t ^&∆– + | _•z- ^{∆– •z} = ^&∆– •_,^{,z ∆– •z}

(129)^**

The algorithm for calculating maximally localized Wannier Functions

The steepest descent algorithm proposed by Nicola Marzari and David Vanderbilt then becomes:

First the Bloch orbital is defined

| = | ^∙

(130)

A smooth trial function |R is defined (for example as a Gaussian times a spherical harmonic).and a new Bloch orbital |S is defined.

|S = >| _, : _, |R

E ,F'

(131)

28 This can then be normalized into a smooth, orthonormal Bloch orbital | ~

| ~

= >|S_, ∗ +S_, |S -_Z^&'(

E ,F'

(132)

From this, the new | can be extracted.

| = | ~ _{& ∙}

(133)

From this the overlap •_,^,z can be calculated.

•_,^,z = + | _•z-

(134)

The operators €_,^,z and _,^,z is then defined.

€_,^,z = •_,^,z •_,^{,z ∗}

(135)

,,z = •_,^,z

• ^,z ∗ gJ8 ‰ln‰• ^,z ŠŠ − z ∙ P> _wz ∗ J8 ‰ln‰• ^,z ŠŠ

z

Qh (136)

The crystal structure is then divided into stars at each point in k-space and the vector z_i is defined as the vector going between the central k-point and a neighbour in star m.

A number of stars are then chosen so that the condition

> _wx_ix_y= 7_iy

z (137)

can be fulfilled for some weight _w . If the crystal is a cubic lattice, then only one star needs to be used and _w = ^#_wž where Z is the number of neighbours in the star.

A starting point for the unitary transformation steps is then defined.

29 D_, = 7_,

(138)

A step for the unitary transformation is then calculated.

∆` = m

∑_{z z}> _w‰ o€ ^,z p − •o ^,z pŠ

z (139)

For some constant m.

The unitary transformation is then updated using this step.

D_, t D_, ^∆–

(140)

The overlap is then updated using the unitary transformation.

•_,^,z t D_, ^a•_,^,z D_,

(141)

A new step can now be calculated using the new overlap. This is iterated until

convergence has been obtained. The overlap does then correspond to the maximally localized Wannier function and the Wannier function can be calculated using

|c* = $

2" ^# | ^C∙ (142)

30

Example: crystalline Si

Crystalline Si is a semiconductor and as such has a gap between the s-like bands and the p-like bands as shown in figure x (Marzari, Vanderbilt (2012) p. 1440).

Figure x: the band structure of Si. the blue circled lines represent a manifold of eight energy bands that are not entangled with any energy bands outside the manifold.

The zero energy is the Fermi energy.

Using the manifold of energy band for Crystalline Si as shown in figure x, a maximally localized Wannier function can be calculated as shown in figure y.

Figure y: The maximally localized Wannier functions of the lowest eight energy bands of crystalline Si, Placed on top of a representation of a section of the crystal.

31

Discussion

In order for this algorithm to be applicable, the Bloch orbitals need to be defined in such a way that there are no degeneracies in the energy band manifold chosen. This is obviously not possible in all cases. Some electric conductors have degeneracies at least somewhere in the crystal for all energy bands in such a way that choosing a degeneracy free manifold may prove unsuccessful. In such cases the energy bands need to be deentangled by some process. One such process was suggested by Nicola Marzari, David Vanderbilt and Ivo Souza in 2001 (Marzari, Vanderbilt, Souza (2001) p.035109-3). After deentanglement The Wannier functions are calculated using similar arguments as given here.

32

References

Blount E I (1962) Formalisms of band theory, In: Seitz F & Turnbull D eds.

Solid state physics vol 13 Pp. 305-373 1th ed. New York: Academic press inc.

Boys S F (1960) Construction of some molecular orbitals to be approximately invariant for changes from one molecule to another, Reviews of modern physics vol 32 no 2:296-299

Chadi D J, Cohen M L (1973) Special points in the Brillouin zone, Physical review B vol 8 no 12:5747-5753

Cloizeaux J des (1964) Energy bands and projection operators in a crystal: analytic and asymptotic properties, Physical review vol 133 no 3A:A685-A697

Klema V C, Laub A J (1980) The singular value decomposition: its computation and some applications, IEEE Transactions on automatic control vol AC25 no 2:164-176 Marzari N, Vanderbilt D (1997) Maximally localized generalized Wannier functions for composite energy bands, Physical Review B vol 56 no 20:12847-12865

Marzari N, Vanderbilt D et. al. (2012) Maximally localized Wannier functions: Thory and applications, Reviews of modern physics vol 84:1419-1475

Monkhorst H J, Pack J D (1976) Special points for Brillouin-zone integrations, Physical review B vol 13 no 12:5188-5192

Souza I, Marzari N, Vanderbilt D (2001) Maximally localized Wannier functions for entangled energy bands, Physical review B vol 65:035109:1-13

In this work a number of expressions are taken from other works, they are denoted by asterisks after the expression number as follows:

* This expression is taken from (Marzari, Vanderbilt (2012))

** This expression is taken from (Marzari, Vanderbilt (1997))

*** This expression is taken from (Blount (1962))