Electronic structure and self energies of randomly substituted solids using density functional theory and model calculations
We describe procedures to obtain the electronic structure of disordered systems using either tight binding like models or quite directly from ab inito density functional band structure calculations. The band structure is calculated using super cells much larger than those containing a single minority component atom. We average over a large number of different super cell calculations containing different randomly positioned minority component atoms in the super cell as well as a varying total number of minority component atoms, weighted by the statistical probability. We develop an efficient and simple algorithm for unfolding of these bands, based on Bloch’s theorem. The unfolded band-structure obtained in this way exhibits momentum and energy broadened structures replacing the gaps observed in often used single super cell calculations. Using the super cell averaged band-structure one can introduce a self-energy, resulting from the scattering of randomly positioned alloy components. The self-energy is causal, and shows strong energy and some momentum dependence. The self-energy shows rather non-trivial behavior and is in general non-zero at the Fermi-energy, resulting in an ill or undefined Fermi surface. The real-part of the self-energy at the Fermi-energy relates to an apparent violation of Luttinger’s theorem. There is no simple relation between the apparent Fermi-surface volume and the electron count. Examples introducing these effects both for model and real binary alloy systems are presented.
pacs:71.15.-m, 71.23.-k, 79.60.Ht, 71.20.Gj, 71.20.Be, 71.20.Eh
The material specific properties of many systems depend crucially upon chemical substitution or doping, i.e. the intentionally introduction of impurities in an otherwise pure crystal structure. Substitutions can be chosen to introduce or change, the number of charge carriers, local magnetic moments, specific optical properties, or vortex pinning centers, just to mention a few commonly known applications. These substituted atoms also introduce a destruction of translational symmetry and therefore scattering into the band structure as an important deviation from the pure system. The effects of disorder resulting from random substitutions alter material properties in specific ways and the vast field of materials applications is strongly dependent on random substitutions.James52 (); Ling88 (); Bouhafs95 (); Orner97 (); Persson01 (); Maier02 (); Zorman05 (); Kakehashi10 (); Alam10a (); Alam10b (); Wadati10 (); Ku10 (); Berlijn11 (); Konbu11 (); Nakamura11 (); Marmodoro11 ()
Although we have powerful ab initio band structure methods to treat pure transitionally symmetric systems at least for systems in which electron correlation effects are less important these methods are generally not applicable directly to randomly substituted systems or disordered alloys. There are several useful approximation schemes available to treat randomly substituted systems or alloys, which can treat the disorder to a high level of accuracy.Faulkner82 (); Klauder61 (); Jarrell01 (); Janis01 (); Laad05 (); Batt06 (); Potthoff07 (); Kodderitzsch07 (); Rowlands09 (); Marodoro11 () It is desirable to create a method able to describe disordered systems with the same level of understanding and approximations as one has developed for periodic systems. Density Functional Theory (DFT)Hohenberg1964 (), within the Local Density Approximation (LDA) is an ab-initio method proved to be useful for a large class of materials. As long as the materials under consideration are not too correlated, the one electron wave functions in a Kohn and ShamKohn1965 () implementation can account for the material specific band-structure. First principle methods used for the description of crystals use the periodicity of the system in order to make calculations tractable. Two distinct problems arise when DFT is used to study the electronic structure of real materials: 1) how to account for disorder breaking the translational symmetry and 2) how to present the results in a meaningful way.
Often, ab-initio calculations including disorder are done assuming periodic boundary conditions within a so-call super cell approximation, i.e. using a cell which is a multiple of the primitive unit cell. It is also typical to choose the smallest cell size possible, containing one impurity atom, in order to reduce computational costs. Of course, the smallest cell with only one impurity or dopant, combined with periodic boundary conditions, can not account for a random impurity distribution in real materials. Such calculations basically result in the band structure of a fictitious new perfectly ordered material. Although in some cases such as the extreme dilute limit where the impurity impurity interference effects can be neglected this may be a reasonable approximation to the real problem it is important to develop theoretical methods based on the same approaches and approximations used for pure systems to check for the importance of randomness. For a detailed discussion of development in ab-initio theory of alloys we refer to a recent review paper by Ruban and Abrikosov.Ruban2008 () We too use the super cell approximation but allow for randomness in the impurity distribution across the super cell. We then propose a simple procedure to present the unfolded calculated band structures in a way as seen by angle resolved photo emission spectroscopy (ARPES) and demonstrate the effect of impurity scattering on the band structure and the self energy.
We first present the method used to unfold the band-structure in the first Brillouin zone of the super cell to the first Brillouin zone of the original cell. Which is basically based on Bloch’s theorem and symmetrization of the wave-functions. The method can be readily applied to all first-principle codes able to calculate partial densities of states and local symmetry or orbital projected band structures. We then demonstrate the basic concepts of the random impurity calculations using a simple tight binding like model and on site impurity potentials. In the last section before the conclusion we describe the use of these concepts in the DFT calculation of a real alloy system composed of Li substituted with either H or Mg.
I Unfolding and randomness of impurities
In order to recover the original band structure from a super cell calculation one needs to unfold the bands back to the first Brillouin zone of the original cell. Several methods are available, which either look at the momentum dependent Hamiltonian,Pavarini05 () or at the eigenfunctions.Ku10 () In case of some additional distortion, i.e. the super cell is not a perfect repetition of the original primitive-cell unfolding can become a rather non-trivial task. A possible method to define unfolded bands is to assign a certain weight to each eigenstate so that the eigenstates with a weight of 1 represent the original band structure while the eigenstates with a weight of zero identify the folded bands. In this section we first demonstrate how this can be accomplished with a rather simple procedure, after the example we present a more thorough derivation.
Consider a one-dimensional tight binding model with one s orbital per site and a nearest neighbor hopping integral . The solution to this problem is the wave function that is symmetric at the zone center and antisymmetric at the zone boundary:
As a super cell we take a doubling of the unit cell. This could model an infinite chain of dimers with intra-dimer hopping and inter-dimer hopping parameter . Assuming equal on-site energies and constant inter-atomic distances of size , the matrix representing the Hamiltonian is given by:
The Hamiltonian can be understood by looking at Fig. (1). There are two sites within the unit-cell, labeled and . (The unit cell, indicated by the black box has length .) Starting from the left site in the unit-cell at the atom labeled one can hop to the right with the strength or to the left with strength . Both processes represent a hopping from site to site , and thus enter as off-diagonal elements between and . In the Hamiltonian the hopping is multiplied by a momentum and directional dependent phase, whenever the hopping connects two species in a different unit-cell. For a hopping crossing the super cell boundary the additional phase is for a hopping to the right, or for a hopping to the left. With , the size of the super cell. For the eigenvalues and the eigen-vectors are:
The main band has the energy of (remember that ) at the zone center and at the zone boundary. The wave-function of the unfolded band within the super cell consists of the sum of and at the zone center and the difference at the zone boundary (). The second solution, i.e. the folded replica, shows the opposite behavior in energy and its wave-function at the zone center (boundary) is made by the difference (sum) of and . Both bands can be seen in panel (b) of Fig. (1).
In a true multi-component system one often uses a so-called ’fat’ band representation in order to separate various orbital contributions into the band structure at given energy and -point. There the weights (i.e. dependent band thickness) are determined by the squares of the eigen-vector components:
with the band index, and some orbital or site index, the crystal momentum (k) dependent eigen-vector component, and a sum over equivalent sites within the unit-cell. It is clear that weights calculated using the eigen-vectors and weights as defined above results in a constant weight throughout the whole Brillouin zone for both bands.
A simple change in the definition of the weighting-function allows one to include information about the relative phase in the wave-functions between different sites. If one defines the weight as:
one basically changed the sum of the norm-squared of the eigenstate pre-factors into the norm of the Fourier transform of the eigenstate pre-factors. Defining the positions of the atoms such that and and applying Eq. (5) to the wave-functions as found in the dimer model immediately shows that the band belonging to the original, unfolded system, with the wave-function gets the weight . The folded band, with the wave-function obtains the weight, .
As a proof of principle, we show in Fig. (2) the band structure of Li metal calculated with TB-LMTO density functional codeAndersen75 () using a super cell. The intensity is calculated according to Eq.(5) in both super cell and primitive cell calculations. Li is a metal with a body centered crystal structure and a lattice constant of 3.51 Å. The primitive-cell contains one Li atom, leading to 4 bands per k-point of different orbital character. At , for example, the Li derived band is at -4eV. The three fold degenerate bands are at about 14 eV above the Fermi energy. In other regions these bands cross and intermix. The variable labels the orbital character: . The index sums over the equivalent sites within the super cell. The band-index labels the different eigenstates found at each k point within the super cell. To highlight the complete band structure, we use the total intensity, traced over the orbital index : .
There is a complication with the definition for unfolding presented so-far. In principle one is free to choose the phase of each basis function separately. Multiplying the phase of the basis function in the dimer labeled as in Fig. (1) should not change any physical results. However it does change the unfolded weight as calculated with Eq. (5). The solution is simple, one should keep track of the phases of the basis functions and introduce a complex-conjugated of this phase in the definition of Eq. (5). Below we will derive the equation for the weight of an unfolded band-structure more rigorously. The here sketched phase-problem will then come out naturally as additional form-factors in the sum of Eq. (5).
i.1 Derivation of general unfolded weight of a band
The weight of a band can be calculated in similar fashion as one would calculate the momentum (inter-, and intra-Brillouin-zone) depended photo-emission spectra. Neglecting polarization effects, assuming a high photon energy such that the final states can be approximated by plane waves, and neglecting the fact that the action of the electric field of the light on the electrons is given by a momentum operator one can write the spectral weight () of a one particle eigenstate as:
whereby labels the momentum, k the corresponding reduced or crystal momentum, and the band-index. (A momentum dependent mixture of , , , or orbital character in the case of Li as shown in Fig. (2).)
For a super cell calculation of an undistorted structure it is straight-forward to show that the spectral weight function as defined in Eq. (6) will be zero for bands that do not belong to the unfolded bands of the primitive unit cell. Take the primitive cell to be defined by the matrix , consisting of three vectors , , and , which span the parallel piped defining the primitive unit cell. Using Bloch’s theorem for a periodic system we know:
Here is an integer vector () and is a periodic function with the period of the primitive unit cell, is periodic up-to a phase:
For a calculation done in the primitive unit cell one can use the block diagonal in k characteristic of the Hamiltonian. Solving for each k vector separately yields the band-structure. Within a super-cell calculation one does not obtain the momentum of a wave-function directly. If one takes a super-cell without perturbing the Hamiltonian the eigenstates are non-the-less given by Bloch’s theorem as given in Eq. (7). The Fourier transform, of such a function, which defines the spectral weight as defined in Eq. (6) is given as:
with k the crystal momentum of the wave-function in the primitive unit cell. The reduced momentum k, is an eigenstate property according to Bloch’s theorem, but undetermined from the super-cell calculation. The reduced momentum index within the super-cell has been suppressed as it is up to a super-cell reciprocal lattice vector implicitly labeled by k. labels an integral over all space and an integral over the original primitive unit cell, and the sum over all possible independent shifts by an integer lattice vectors () of the primitive unit cell. The first term of the unfolded spectral weight in the primitive unit cell () contains a sum over all unit-cells. This term leads to the selection rule:
with an integer vector and a matrix representing the reciprocal lattice: . This term is responsible for the fact that only those bands that belong to the principle unit cell will gain weight, i.e. it unfolds the folded band-structure. The second term of the unfolded spectral weight is a material and Brillouin zone dependent form factor. In photo-emission experiments it is this term (together with the here neglected polarization dependence) that is responsible for the Brillouin zone dependent spectral weight.
The unfolding of bands as described in Eq. (6) can be simplified further by approximating the Brillouin zone dependent form factor. We will first briefly describe model calculations and then extend it to a formalism that can readily, without much effort, be implemented in most ab initio density functional programs.
Unfolding of bands in tight-binding model calculations
For model calculations one often expresses the band-structure in terms of some form of tight-binding representation. If we take the variable as the index for the local orbital and or spin character of the tight-binding basis set (in the case of Li is , , , or ) then one can define the spectral weight matrix for a given one-particle eigenstate () as:
whereby the operator annihilates an electron at momentum and of orbital (spin) character . The momentum dependent annihilation operator can be written as a sum of position dependent operators as:
which describes the relation between the Fourier transform of a one-particle wave-function and the expectation value of the momentum dependent occupation number operator. If the orbital character of a certain band is of no importance one can write the total unfolded spectral weight of a certain band as the trace over the spectral weight matrix:
Unfolding of bands for ab initio density functional calculations.
In the case of an ab initio density-functional theory calculation one could naturally first create a tight binding model and then use Eq. (12) in order to calculate the unfolded spectral weight. Such a scenario has recently been proposed by Wei Ku et al. Ku10 () and successfully used on several systems.Berlijn11 (); Konbu11 () The unfolded band-structures found in this way are expected to be rather similar to the unfolded band-structure one obtains by defining the weight of a band via a Fourier transform (see Eq. (6)). For large super cells with several impurity atoms, defining a tight-binding Hamiltonian might become troublesome. An alternative method that allows one to directly obtain the unfolded band-structure from an ab initio calculation therefore is useful.
One can simplify the Brillouin zone dependent form factor of the unfolded spectral weight as it appears in Eq. (I.1) by looking at the one-particle wave-function character within a certain (muffin-tin) sphere defined around each atom. These projections are implemented in most ab initio DFT codes, where these projections are used to calculate the partial character of a given eigenstate needed to plot partial density of states, or bands of a certain orbital character. The projected wave-function can be written as:
with a band index of the one-particle eigenstates, an orbital and spin index, the sum over all sites that contain a muffin tin sphere onto which one projects, the position of the sphere labeled with the index , the functions projected onto (with orbital character and centered at position ), and the projected pre-factor. Similarly as done for model calculations it makes sense to define an orbital () specific unfolded weight, starting from:
Unfolding such a projected one-particle wave-function with the use of Eq. (6) yields:
defines a Brillouin zone dependent form factor, and although this is a physical quantity one can simplify the unfolding of bands substantially by neglecting this term. If the orbital basis set is defined such that the main outer part of is real and positive for all then a powerful approximation is to set the integral in Eq. (18) equal to 1. In this case the unfolded weight becomes:
In the case that the super cell is undistorted, i.e. a pure repetition of the system one can understand that this definition of the spectral weight indeed unfolds the band-structure with the argumentation based on Bloch’s theorem. This can be done by following the argumentation as given in the beginning of this section from Eq. (7) to Eq. (I.1). When the super cell is not a perfect repetition of primitive basis cells then the weight of the folded bands starts to deviate from zero. The weight in the folded bands is related to the mixing of the folded and original bands in the band-structure and therefore a measure of the distortion.
One should be slightly careful how the form-factor as presented in Eq. (18) is neglected. The phase of this integral depends on the phase of the basis function , which can be chosen freely. In many calculations the projections are done onto atomic like wave-functions. These functions are real, which fixes most of the phase, but still can either be positive or negative. The radial part of atomic like wave-functions has nodes. The important outer part of the radial wave-function is either positive or negative, depending on the number of nodes. The number of nodes for an atomic like wave-function is given as the phase of the basis function therefore is . Naturally different implementations might use different phases for the basis functions. (A complex phase of can for example be a useful alternative choice.) Using our phase definition, including the factor for the basis functions, the Brillouin zone dependent form factor will be approximated as:
whereby is the principle quantum number of the orbital with index and position and is the angular momentum of the orbital with index and position . The unfolded weight then becomes:
The partial character is the projected weight of the eigenstate with band index , orbital and spin character , at site in the super cell for the crystal momentum k. The unfolded weight of this band is calculated by taking a Fourier transform summed over (pseudo) equivalent sites () within the super-cell. This can either be the same atom at positions in the super cell slightly deviating from a perfect repetition of the primitive unit cell, or different atoms in the super cell in the case of substitution. The norm-squared of this Fourier transform defines the unfolded weight of a certain band and orbital character.
i.2 configuration average for the description of randomness
There is a crucial difference between a calculation for an ordered structure and a system with randomly substituted atoms. In a system with ordered impurities gaps will open in the band-structure, whereas for a random system this may not be the case and gaps can become smeared out structures, as will be shown by several examples later. In order to simulate the situation of an infinite crystal with a random impurity distribution one would like to create very large super cells. In order to do these calculations with the use of tractable cluster sizes it is important to average over different random distributions of impurities with equal impurity density and over different impurity densities within the cluster. Lets define the chance to have a substituted atom at site to be equal to . For a cluster of total sites the probability to have a number of substituted atoms within the cluster is given by the binomial distribution:
In order to illustrate the impurity distribution further we show in Fig. (3) a plot of this distribution for a doping of 5% () and a cluster size of 4, 8, 16, 32, 64, and 128 total atoms. One should realize that even for a reasonable cluster size there still is a finite chance to have no substituted atoms within the cluster. For each total number of substituted atoms in the cluster there are still many different ways to arrange these impurities. It is important to take this into account, since the scattering of nearest neighbor impurities is in general very different from the scattering of isolated impurities. In general one should determine the band structure in a super cell method with as large a cell as possible in real space and average over a varying number of impurities and a varying distribution of these impurities. Each of these calculations should be weighted with the probability to find this particular kind of configuration, including both the chance to have this specific density within the cluster and this specific order within the cluster. In practice it is not needed to calculate all possible impurity configurations (including variation of density and distribution) within the cluster. A random subset will already give a reasonable approximation to the real solution, provided that the cell is large enough and a large number of random configurations are taken.
We define a configuration to be a specific realization of a super cell with a random number of substituted atoms at random positions. The average density of substituted sites is given by .
Ii Random onsite energy impurity model
Here we exemplify several aspects of random impurity calculations, using a simple model system. The model we use is well studied which allows us to compare the super cell calculations to theoretical approximations used elsewhere.Klauder61 (); Faulkner82 (); Janis01 (); Jarrell01 (); Laad05 (); Batt06 (); Rowlands09 () We use a single band model with constant nearest neighbor hopping () and a onsite impurity potential at a fraction of all atoms randomly placed within the solid. The Hamiltonian is given as:
with the sum over all sites and nearest neighbor sites . is the hopping-strength, resulting in a band-width of . is the dimension of the system. is a potential with the value at a fraction of randomly placed sites and zero at other sites. is the density of substituted atoms.
As the Hamiltonian only includes one-particle interactions it is natural to Fourier transform it to momentum () space:
with and the position of atom . Within this representation the hopping part becomes diagonal in . The random potential allows for scattering between different -vectors. In the limit of large cluster sizes or for a potential averaged over many different cluster configurations one finds that and that for scales as . The averaged Hamiltonian (), is thus defined as . One might be tempted to state that for infinitely large clusters or for a cluster of finite size, averaged over a large number of configurations the Hamiltonian becomes diagonal in momentum and thus momentum is always a good quantum number. Care has to be taken as there is indeed only an infinitesimal small coupling to different momenta for infinitely large clusters, but there are an infinite number of different momenta possible to which one can scatter thus leading to a finite perturbation.
The one particle Green function is defined as:
The spectral function whose poles define the band-structure, is given as . The off diagonal elements of in , i.e. those elements for which scale as a function of cluster size, in the same way as the off diagonal elements of the Hamiltonian, namely as . For large clusters or for an average over many different random configurations the Green function becomes diagonal in momentum . Calculating the Green function, using the average Hamiltonian (which is diagonal in momentum) is however different from using the full Hamiltonian, even if one averages the Green function afterward. In the following whenever we write we will assume that this is the Green function averaged over several configurations ().
ii.1 Moment analyzes and series expansion
The random onsite energy impurity model has been extensively studied with the use of series expansions.Klauder61 (); Faulkner82 () The -th moment () of the momentum integrated spectral function (density of states) is given as:
with the total number of sites in the cluster. The Hamiltonian is given as with and . Analytical expressions can immediately be written down, if one neglects contributions to the moment which arise from products of and . The moments of are given as:
which are the moments of the original, unperturbed system. The non-zero lowest order moments are , , and . The moments of are:
Contributions due to products of and are important for the moments with .
The momentum dependent Green function can be expressed in terms of a series in with the use of the Dyson equations. For a Hamiltonian only containing one electron terms one can write:
For the random configuration averaged Green function, one finds to first order in the onsite impurity potential :
It is often useful to define the self energy for the averaged Green function as:
If one compares the definition of the self energy with Eq. (29) one finds a striking similarity. Multiplying Eq. (31) with and , one reproduces the form of Eq. (29), but now for the configuration averaged, single momentum dependent Green function, instead of the original Green matrix. Using Eq. (30) one finds that to first order in the self energy is given as . To linear order in the perturbing potential one can thus replace the random impurity system by a periodic system, which has an averaged potential on each site. This approximation is known as the virtual crystal approximation. It can correctly describe the averaged band-structure and total density of states to first order in . It can not describe changes to the one-particle wave-functions important for the band-character, the spatially varying charge density,Wadati10 () the broadening of bands, or the possible appearance of impurity bound states, which arise in quadratic and higher order of the perturbing potential.
ii.2 Cluster size and configuration averaging
In order to fully describe the effects of random impurities it is necessary to go beyond a series expansion in the perturbing potential. For such cases doing exact diagonalization of large clusters is a powerful method, especially with the current computational resources. Within this section we compare several calculations done for different cluster sizes. We furthermore show that the cluster sizes can be considerably reduced if one averages over several random configurations, i.e. both varying the impurity density as well as the impurity distribution within the cluster.
In the left most column of Fig. (4) we show calculations for the band-structure of a single random configuration, with fixed concentration. The different rows show the calculation for a different total cluster size. The calculations are done for a 25% () substitution and . The smallest cluster possible which allows for a substitution of 25%, namely with a total of sites and one impurity atom, shows the opening of large gaps in the band-structure. This is to be expected as for a cluster with substitution one has a perfectly ordered super-structure, which is very different from a random system. The larger the cluster becomes the smoother the band-structure becomes. The gaps become smaller and smaller and more and more distributed to random places in momentum space. One however needs very large clusters to find convergence if one does not average over different random configurations. Even for a single calculation of a super cell there is still some gaping visible in the band-structure (left column bottom row of Fig. (4)). This becomes much better if one averages over different random configurations (i.e. averaging over both a varying number of substituted atoms within the cluster and averaging over the possible different ways to substituted these atoms). The forth column shows the same band-structure calculations, but now averaged over 1000 random configurations. As one can see these look much more smooth. The gaps disappear when an average over different configurations is taken. Even for a super cell calculation the averaged band-structure looks quite reasonable, whereas the calculation of a single configuration is far from converged. What happens is that each different random configuration has gaps at different energies and vectors, leading to an averaged broadening of the bands and not to a gaping of the bands. A similar effect can be seen for the density of states. The calculations for a single configuration are rather noisy; the average over 1000 configurations however shows a nice continues behavior.
Averaging is crucial for the calculation of the Fermi-surface. The Fermi-surface is calculated by assuming a filling of 1 electron per atom in the unit cell after averaging all calculations done for different configurations. A single configuration shows only some points on the Fermi-surface, whereas an average over several configurations shows a very nice, though broadened Fermi-surface. Note that in the Fermi-surface plots for small super cell calculations one can see the periodicity of the super cell very clearly as features with the periodicity of the super cell Brillouin zone. The use of small super cell calculations can lead to spurious shadow bands.
ii.3 Comparison to other approximations
Next we compare the direct method of calculating the band-structure, Fermi-surface and density of states as presented here with other methods like the coherent potential approximation and its non-local version.Faulkner82 (); Jarrell01 () in Fig. (5) we present calculations as a function of the impurity potential for a 50% () substituted system. The model and parameter range used are the same as presented in Fig. (5) of the paper by Batt and Rowlands Batt06 (), which allows for a direct comparison between the different approximations. For a potential of one finds a clear Fermi-surface at a filling of one electron per site. The Fermi-surface is rather broad. These findings are basically equivalent to the results found with the coherent potential approximation and for the non-local coherent potential approximation calculated using a cluster. Increasing the impurity potential further will eventually lead to two separate bands separated by a gap. In principle one could state that the system becomes insulating at half filling. However, due to the large broadening of the bands, one still finds for a very large potential range a finite intensity at the Fermi-energy. It is in this range that one finds differences between the coherent potential approximation and their cluster extensions Batt06 (); Rowlands09 (). For one finds a new Fermi-surface topology. Around one finds an electron pocket which is related to the host band. Around one finds a hole pocket derived from the impurity band. It is note-worthy that the coherent potential approximation is not able to reproduce this effect. Batt06 (); Rowlands09 () The non-local coherent potential approximation, for large cluster sizes should naturally exactly reproduce the current calculations, under the condition that the cluster size taken into account is large enough. For a cluster size of one can see that some features of the gaping are indeed correctly reproduced. The exact potential at which the Fermi-surface topology change takes place, as well as some details of the Fermi-surface are however different. For accurate calculations it is important to take relatively large cluster sizes, within the non-local coherent potential approximation. The effect of the bath in the non-local coherent potential approximation would be relatively minor for such large clusters sizes. It is therefore relatively straight forward to extend the non-local coherent potential approximation to the current cluster calculations. For large cluster sizes both methods are equivalent.
ii.4 Self energy
A powerful concept for the understanding of random impurity calculations is the self energy . The self energy is a measure of how-much the band-structure of the system with impurities deviates from the system without substituted atoms. The full Green function is given in terms of the original Green function () and the self energy as:
The real part of defines an energy shift of , with respect to . The imaginary part of is related to an energy broadening of the Green-function. The Self-energy is casual, meaning that in our representation the imaginary part is strictly negative, and goes to zero as for large absolute values of . The real-part is related to the imaginary part by Kramers-Kronig. Note that for a system of impurities there are no specific relations for the self-energy at the Fermi energy. As seen in Fig. (4) and (5) there is a substantial broadening of the Fermi surface (non-zero imaginary part of the self energy) as well as a possible change of Fermi surface topology (energy and momentum dependent real and imaginary part of the self energy).
It is generally excepted that a system with random impurities has a broadened band-structure compared to a pure system. The imaginary part of the self energy must be integrable and goes to zero for large and small energies. By virtue of the Kramers-Kronig relations this means that one then also must have a downward energy-shift for energies lower then the main broadening and an upward energy shift for energies higher then the energy at which the broadening takes place. The broadening of bands will thus go hand in hand with shifts in the band-structure. This will become more clear with the use of the example given in Fig. (6).
In Fig. (6) we show the band-structure, density of states, Fermi-surface (for a filling of one electron per atom) and real and imaginary part of the self energy at as well as their momentum dependence for a cluster as a function of the impurity potential at . The real part of the self energy at (sixth column) has a constant part, independent of the energy, equal to , which one would expect in linear order in . In the forth column one can see the imaginary part of the self, energy, which peaks around . This leads by Kramers-Kronig to a real-part of the self-energy, which is related to a shift of the band-structure.
The shift of the band-structure due to the real-part of the self energy can be seen very clearly in the Fermi-surface plots. For a small potential of the order of or smaller the band-structure, density of states and Fermi-surface are, aside from some broadening, not changed much compared to the original host electronic structure. The Fermi-surface of the host is shown in the third column, by the small green pockets around . For the Fermi-surface of the system with impurities overlaps with the original Fermi-surface. For a larger impurity potential the band-broadening becomes more substantial, and a corresponding shift of the band-structure can be seen. For the Fermi-surface this results in apparently larger pockets around . The electron count does not change between the different calculations shown in the different rows in Fig. (6), The Fermi-energy is defined such that there are 1000 electrons in the cluster containing atoms. It appears that Luttinger’s theorem is violated. Calculating the number of electrons in the system by looking at the Fermi-surface as defined as the maximum intensity of at the Fermi energy gives an incorrect number of electrons. One can understand these findings by realizing that the calculations are not done for a periodic system, but averaged over many different randomly ordered periodic systems, within a super cell. The Green function is obtained by unfolded these super cell calculations back to the original Brillouin zone. The conditions for Luttinger’s theorem to be applicable are not fulfilled.
The change in Fermi-surface can be understood in a different way. For large negative impurity potentials the substituted atoms will build an impurity band, at a lower energy then the original host band. This impurity band, will be doubly occupied, i.e. full. This doubly occupation of the impurity band removes electrons from the Host band, leading to an apparent change in electron occupation. In order to test this simple picture we added the Fermi-surface one would expect for an unperturbed system with holes, i.e. on average only 0.7 electrons per unit cell. These are the large Green Fermi-surface contours shown in the third column of Fig. (6). Indeed it seems that for large potentials the apparent hole doping can be understood in terms of doubly occupied bound states around the substituted atoms, which do not contribute to the Luthinger’s count. There is no clear border between the two regimes where for small potentials the virtual crystal approximation is reasonably fulfilled to a regime for large potentials where the system can be understood in terms of an impurity band. There is always some deviation from the virtual crystal approximation, however only to quadratic order in the impurity potential.
Iii ab initio results
There is an important difference between model calculations and self-consistent first principle studies. Within a model the potential of the host and impurity are fixed, irrespectively of their electron count. Within a self-consistent calculation the potential depends on the local electron count. This later effect will lead to screening, and different impurity potentials when two impurities are nearest neighbors compared to impurities infinitely far separated from each other. The model presented before will go astray if one places a large patch of impurity atoms next to a large patch of host atoms. In this case there will be a charge donation of the host into the impurity, independent of the distance between them, thus even if the impurity atom is very far away from the host. In reality this will cost Coulomb energy, which scales with the distance and should prevent too much of a charge flow. It is thus important to do self consistent calculations. Density Functional Theory captures these effects reasonably well. The method, language and general features described in the previous sections are directly transferable to density-functional super cell calculations.
iii.1 Cluster shape dependence for LiMg
Mean-field approximated model calculations including several thousands of sites are possible and finite size effects therefore can be largely eliminated. Density functional theory calculations are more involved and although modern computers are very powerful it is important to try to optimize cluster sizes and shapes. Naturally Mean-field theories like most implementations of Density functional theory are still numerically easier to solve then true many-body calculations, from which some knowledge about optimal cluster sizes can be obtained.Betts99 () In Fig. (7) we show three calculations for the LDA band-structure of Li (middle panels) and LiMg (top panels), calculated with the use of LMTO.Andersen75 () The calculations all correspond to different ordered structures of the same chemical composition. The difference between the three rows is the shape of the super cell. In the middle row one can see in orange that the folded band-structure depends on the choice of the super cell used. The unfolded bands are independent of this choice (when no impurities are introduced). In the top row one can see the unfolded bands when one Li atom is replaced by Mg. If one compares the bands of LiMg with the folded bands of Li then a striking feature appears. Looking for example to the lowest derived band at -4 to 0 eV going from to in the unfolded Brillouin zone. Whenever the Li folded bands cross the Li main bands (see middle rows) a gaping of the bands at that place in space appears when one Li atom is replaced by Mg.
This is strongly related to the difference between the calculation of an ordered structure and the calculation of a randomly substituted material. For ordered superstructures gaps appear when the folded bands cross the original bands. The dependence of this folding on the shape of the super cell opens the possibility to not only average over different configurations, i.e. different impurity densities and different orderings of the impurities, but also to average over different cluster shapes.
The beauty of ab initio calculations over model calculations is that one for each calculation assuming a certain impurity configuration and super cell shape, not only obtains the electronic properties, but also a total energy. This total energy can be used as a weighting factor when the different cluster sizes are averaged. On the bottom of Fig. (7) we show the total energy of the 3 different super cell calculations per site, with respect to the energy of the lowest energy state. The material studied will be grown or annealed at a certain temperature, where one can assume the system to be in thermal equilibrium. Then one could rapidly cool the system. The distribution of the different clusters can be assumed to be governed by the Bolzman statistics one would have at the annealing temperature. Annealing the LiMg crystal at 300K would give a distribution of 0.29, 0.33 and 0.38 for the three different super cells as shown in Fig. (7). If the sample would be in thermal equilibrium at 30 K the distribution function would be 0.05, 0.23 and 0.72. Naturally the before sketched numbers are just exemplary. A real calculation needs much larger cluster sizes and the average over more then three different configurations or super cell shapes in order to give realistic results.
iii.2 Cluster size and averaging over the number of impurities in different super cells for LiMg
The model calculations as shown in the previous section show averages over configurations where both the possible number of substituted atoms as well as the random placement of these atoms is varied. One would expect that for large enough cluster sizes the variation of the number of impurities becomes less important. In Fig. (8) we show in panel (a) a calculation done for several configurations where for each configuration the number of substituted atoms is chosen to be 4, which for a super cell of 16 sites results in a density of . In panel (b) we show a calculation for the same super cell, but now we also vary the total number of substituted atoms. The calculation in panel (b) includes averaging over all configurations according to their distribution function as given in Eq. (22) with and . The average density is in both cases. One can see that including an average over different number of substituted atoms for small cluster sizes improves the calculation. One can compare these results to panel (c) and (d) of the same Fig. (8). Here the cluster size is increased from to . Again in panel (c) the total number of substituted sites is kept fixed at for each configuration which varies the random positions of these substituted sites. In panel (d) also the number of substituted sites is varied between the different configurations. For larger clusters (at these impurity levels) it becomes less important to have the exact random distribution as given in Eq. (22). One furthermore observes that for a cluster of a bcc structure containing sites the averaged band-structure already looks very reasonable.
iii.3 Concentration dependence in H substituted Li
To illustrate the effect of disorder on the band structure as obtained from density functional theory, we perform DFT calculation of Li metal substituted with Hydrogen using TB-LMTO code and super cell approximation. The size of the super cell is (128 sites due to bcc structure). A random number generator is used to generate 100 configurations for any given H concentration. The unfolding procedure is used to present the calculated band structures in the original BZ. All 100 band structures are plotted in the same plot. Fig. (9) shows such plots for 6.25, 12.5, 25 and 50% of Hydrogen. The impurity effect is twofold, first, the host band structure is broadened in both energy and momentum, and, second, there is an upward shift in energy with respect to the Fermi energy. The shifts can be understood in terms of on-site energy differences between Li and H states. Since the ionization energy of H is about 8.2eV higher than that of Li, the H -states are expected to have higher binding energy than that of Li. This means that the Hydrogen impurity is an acceptor. The shifts, however, are compensated somewhat by the changes in the band widths. At about 50 percent, the band structure of LiH alloy turns over to that of pure H which is also broadened and shows the shifts in the Fermi energy which is now due to Li impurities as can be seen in the right bottom panel of Fig. (9).
iii.4 Order versus disorder in the Li H alloy
The specific case of LiH is rather interesting since it allows us to study the effect of order on the electronic structure. The electronic structure of the ordered alloy is shown in the left columns of Fig. (10). It is calculated with the conventional bcc unit cell that corresponds to two interpenetrating simple cubic lattices. The material is an insulator with the gap of about 1eV which separates the occupied H 1s band with the conduction band formed by the states of mainly Li character. The random alloy is, however, metallic (see right columns of Fig. (10)). We note the very large smearing effect of the band structure especially in the region dominated by lithium based states above the Fermi energy. There is, of course, no background nor smearing of the bands in the ordered case.
In conclusion we have discussed a band structure based procedure which can be successfully used to describe the electronic structure of disordered alloys and randomly substituted systems. We have used a simple tight binding based model to describe the physics of the differences between minimum size super cell calculations and large super cells with random distributions of impurity atoms in the cell. In order to visualize the changes we present a simple unfolding procedure which is demonstrated to work very well indeed. We have shown that with small impurity potentials the band structure is qualitatively described by a rigid shift as in a virtual crystal approximation. However even here the bands become broadened in both energy and momentum space. This broadening becomes very large for larger impurity potentials. We also show that the Fermi surfaces are strongly smeared out for medium size impurity potentials and that the peak position of the broadened structure as a function of momentum would yield a Fermi surface determined by the crossing at an energy corresponding to the electron count i.e. the ”Fermi energy” which does not agree with the Luttenger theorem size of the Fermi surface.
In order to quantify the effects of the impurity scattering effects we present the real and imaginary parts of the self energy which clearly demonstrates the various effects impurity scattering and randomness has on the electronic structure. The smearing effects at the Fermi energy can cause a semiconducting density of states in the ordered case turn into a metallic density of states. This effect is quite spectacular in the comparison of a density functional based ab initio calculation of the electronic structure of ordered and disordered LiH. The results we show for Mg Li and Li H alloys using density functional methods demonstrate that ab initio methods can be used which opens the door to the study of a large number of real material systems.
We would like to thank O. K. Andersen and M. Berciu for stimulating discussion. We acknowledge financial support from the Candian funding agencies NSERC,CFI, and the CRC program. This work was done within the Max Planck/ UBC funded Center for Quantum Materials at UBC and was enabled in part by the use of computing resources provided by WestGrid and Compute/Calcul Canada.
- H. M. James, and A. S. Ginzbarg, J. Phys. Chem. 57, 840 (1953).
- M. F. Ling, and D. J. Miller, Phys. Rev. B 38, 6113 (1988).
- B. Bouhafs, F. Benkabou, M. Ferhat, B. Khelifa, J. P. Dufour, and H. Aourag, Infrared Phys. Tech. 36, 967 (1995).
- B. A. Orner, and J. Kolodzey, J. Appl. Phys. 81, 10 (1997).
- C. Persson, R. Ahuja, and B. Johansson, Phys.Rev. B 64, 033201 (2001).
- Th. A. Maier, and M. Jarrell, Phys. Rev. Lett. 89, 077001 (2002).
- B. Zorman, S. Krishnan, D. Vasileska, J. Xu, and M. van Schilfgaarde, J. Comp. Elec. 3, 351 (2004).
- Y. Kakehaski, M. Atiqur R. Patoary, and T. Tamashiro, Phys. Rev. B 81, 245133 (2010).
- A. Alam, T. Saha-Dasgupta, and A. Mookerjee, Phys. Rev. B 81, 054201 (2010).
- A. Alam, and A. Mookerjee, Phys. Rev. B 81, 184205 (2010).
- H. Wadati, I. Elfimov, and G. A. Sawatzky, Phys. Rev. Lett. 105, 157004 (2010).
- W. Ku, T. Berlijn, and C.-C. Lee, Phys. Rev. Lett. 104, 216401 (2010).
- T. Berlijn, D. Volja, and W. Ku, Phys. Rev. Lett. 106, 077005 (2011).
- S. Konbu, K. Nakamura, H. Ikeda, and R. Arita, arXiv 1108.0595 (2011).
- K. Nakamura, R. Arita, and H. Ikeda, Phys. Rev. B 83, 144512 (2011).
- A. Marmodoro, and J. B. Staunton, J. Phys.: Con. Series 286, 012033 (2011).
- J. S. Faulkner, Prog. Mat. Science 27, 1 (1982).
- J. R. Klauder, Ann. Physics 14, 43 (1961).
- M. Jarrell, and H. R. Krishnamurthy, Phys. Rev. B 63, 125102 (2001).
- V. Janis, Phys. Rev. B 64, 115115 (2001).
- M. S. Laad, and L. Craco, J. Phys.: Condens. Matter 17, 4765 (2005).
- G. M. Batt, and D. A. Rowlands, J. Phys.: Condens. Matter 18, 11031 (2006).
- M. Potthoff, and M. Balzer, Phys. Rev. B, 75, 125112 (2007).
- D. Ködderitzsch, H. Ebert, D. A. Rowlands, and A. Ernst, New J. Physics 9, 81 (2007).
- D. A. Rowlands, Rep. Prog. Phys 72, 086501 (2009).
- A. Marmodoro, and J. B. Staunton, J. Phys.: Conf. Ser. 286, 012033 (2011).
- P. Hohenberg, and W. Kohn, Physical Review 136, B864 (1964).
- W. Kohn, and Lu J. Sham, Physical Review 140, A1133 (1965).
- A. V. Ruban, I. A. Abrikosov, Rep. Prog. Phys. 71, 046501 (2008).
- E. Pavarini, A. Yamasaki, J. Nuss and O. K. Andersen, New J. Phys. 7, 188 (2005).
- O. K. Andersen, Phys. Rev. B 12, 3060 (1975).
- D. D. Betts, H. Q. Lin, and J. S. Flynn, Can. J. Phys. 77, 353 (1999).