LDA and GGA Investigations of Some Ground-State Properties and the Compton Profile of Copper with the All-Electron MAPW Method

Abstract

The electronic structure is self-consistently evaluated both in the local-density and the generalized-gradient approximations. The convergence with respect to the intrinsic parameters of the modified augmented plane-wave scheme is thoroughly investigated. Besides the Ar-core, the less-extended Ne-core and the most appropriate exchange-correlations functionals are considered. It is studied how some characteristic ground-state properties and the Compton profiles depend on these choices. Apart from the long wave limit of the Compton profiles, most of the other properties are found in good agreement with the experimental facts.

1. Introduction

The electronic structure of the noble metal Cu has been investigated extensively in the last four decades and reasonable agreement between theory and experiment was found. In most theoretical examinations the core state is chosen to be of Ar type which considerably simplifies the numerical work as the band structure scheme has only to be applied to the 4s and 3d electrons. However, the wave functions of the 3s and 3p electrons do not completely vanish at the surface of the sphere lying fully within the atomic cell with the consequence that these outer core electrons dynamically interact with the valence electrons. (Typically their wave functions differ by the factors 80 and 37 from the corresponding values of the 4s and 4p functions, resp.). It is interesting to examine to what extent the agreement is improved not only by including the outer core electrons in the band structure calculations but also by the use of more advanced theoretical schemes combined with the use of more powerful computers which allow a considerable increase of various parameters inherent in the computational schemes.

The outline of the paper is as follows. In Section 2 we start with a review of the modified augmented plane wave (MAPW) scheme [1, 2]. Then we will confine ourselves to some specific features as well as the proper choice of some intrinsic parameters of the MAPW scheme are described. Section 3 is devoted to the evaluation of the most important ground-state properties within the frame work of the local density approximation (LDA) and the generalized gradient approximation (GGA). In both approximation, the most appropriate exchange-correlation functionals [3–5] are applied. Section 4 resumes former investigations of the Compton profiles (CP) by the author [6] on a far more refined level nowadays possible. The MAPW scheme turns out to be best suited to meet with this problem without any truncations. At first we study the influence of the spherically averaged electron moment density (EMD) on the Compton profiles. Because it can be done quite efficiently it allows extensive investigations of different models and assumptions. In the following subsections, the directional Compton profiles are investigated. A critical comparison with previous Kohn-Korringa-Rostocker (KKR) studies by Sakurai et al. [7] concludes the work.

2. The Electronic Structure

2.1. Basic Concepts of the MAPW-Method

Since the MAPW method has received ample documentation elsewhere [1, 2], we will confine ourselves to reviewing only its most important features. As to the original APW-method [8], the trial functions outside the atomic spheres are approximated by a superposition of plane waves𝜒_⃗𝑘𝑟⃗=_⃗𝐾𝑣𝐾⃗𝑒^{⃗⃗𝑖(𝑘+𝐾)⋅⃗𝑟},(1) where the summation over the reciprocal fcc lattice vectors is restricted by the inequality𝑘⃗+𝐾⃗²≤QMXSQ2𝑎𝜋².(2)QMXSQ is a dimensionless magnitude which controls the number of plane waves 𝑁_pw. This choice of the reciprocal vectors guarantees the full point group symmetry of the Bloch energies 𝐸_⃗𝑘 at any 𝑘⃗ in the Brillouin zone. The optimal choice of QMXSQ is discussed below. To comprehend the form of the trial functions within the atomic spheres, the right-hand side of (1) is rewritten by use of the plane waves expansion in spherical coordinates with origin at the central atomic polyhedron_⃗𝐾𝑣𝐾⃗𝑒^{⃗⃗𝑖(𝑘+𝐾)⋅⃗𝑟}=4𝜋^∞_{𝑙=0𝑚=𝑙𝑚=−𝑙}×_⃗𝐾(−1)^𝑚𝑖^𝑙𝑌_𝑙,𝑚𝑟⁰⃗𝑣𝐾⃗𝑌_𝑙,−𝑚𝑘⃗+𝐾⃗⁰𝑗^𝑙||𝑘⃗+𝐾⃗||𝑟.(3)𝑌_𝑙,𝑚 are normalized spherical Harmonics and 𝑗_𝑙 are spherical Bessel functions of first kind. Unit vectors are denoted by the superscript 0. It is clear that the terms with small values of 𝑙 do not properly describe the influence of the crystal potential near the nuclei. Following Slater [8], we substitute the spherical Bessel functions for 𝑙≤𝐿 by a superposition of properly chosen linear combination of radial functions 𝑅_𝑠𝑙(𝑟). The first index 𝑠 counts the elements of the set. This augmentation idea leads finally inside the central atomic sphere to the Ansatz𝜒_⃗𝑘𝑟⃗=4𝜋^𝐿_{𝑙=0𝑚=𝑙𝑚=−𝑙}+4𝜋_∞(−1)^𝑚𝑖^𝑙𝑌_𝑙,𝑚𝑟⁰⃗_𝑠𝐴_𝑠𝑙𝑚𝑅_𝑠𝑙(𝑟)_{𝑙=𝐿+1𝑚=𝑙𝑚=−𝑙}×_⃗𝐾(−1)^𝑚𝑖^𝑙𝑌_𝑙,𝑚𝑟⁰⃗𝑣𝐾⃗𝑌_𝑙,−𝑚𝑘⃗+𝐾⃗⁰𝑗^𝑙||𝑘⃗+𝐾⃗||𝑟(4) which alternatively may be rewritten as𝜒_⃗𝑘𝑟⃗=_⃗𝐾𝑣𝐾⃗𝑒^{⃗⃗𝑖(𝑘+𝐾)⋅⃗𝑟}+4𝜋^𝐿_{𝑙=0𝑚=𝑙𝑚=−𝑙}×⎡⎢⎢⎣_𝑠(−1)^𝑚𝑖^𝑙𝑌_𝑙,𝑚𝑟⁰⃗𝐴_𝑠𝑙𝑚𝑅_𝑠𝑙(𝑟)−_⃗𝐾𝑣𝐾⃗𝑌_𝑙,−𝑚×𝑘⃗+𝐾⃗⁰𝑗^𝑙||𝑘⃗+𝐾⃗||𝑟(5) by use of the plane waves expansion.

Summarizing, each trial function 𝜒 of the MAPW scheme is a linear combination of plane waves which are, for any angular quantum number 𝑙≤𝐿, augmented by the square bracket of (5) within the atomic spheres. The radial functions 𝑅_𝑠𝑙 are solutions of the spherical radial differential equation with properly spherically averaged crystal potential within the atomic spheres. For each value of 𝑙≤𝐿 a set of 𝑅_𝑠𝑙 counted by the first index 1≤𝑠≤𝑁𝑅(𝑙) is used. 𝑅_𝑠𝑙 is generated by requiring a given value of the logarithmic derivative (𝑑𝑅_𝑠𝑙/𝑑𝑟)/𝑅_𝑠𝑙 on the surface of the atomic sphere 𝑟_APW. It turned out to be mostly sufficient to choose these logarithmic derivatives to be either 1 or −1.

In contrast to other versions of APW each trial function as well as its gradient are made continuous everywhere in the atomic polyhedron by requiring that the corresponding square brackets of (5) vanish for 𝑟=𝑟_APW_𝑠𝐴_𝑠𝑙𝑚𝑅_𝑠𝑙(𝑟)=_⃗𝐾𝑣(𝐾)𝑌_𝑙,−𝑚∀(−𝑙≤𝑘⃗𝑚+𝐾≤⃗𝑙,⁰0𝑗≤^𝑙𝑙||𝑘≤⃗+𝐿𝐾)⃗||𝑟(6) and the analogous equation with the first derivative with respect to 𝑟. In the variational procedure which makes the expectation value of the crystal Hamiltonian 𝐻 extremal, these equations play the role of additional constraints with the consequence that a general Hermitian eigenvalue problem of rank 𝑁_PW+3∑^𝐿𝑙=0𝑁𝑅(𝑙) has to be solved. For a given value of 𝑘⃗ it gives Bloch states with energies 𝐸_𝑛⃗𝑘 and wave functions ⟨𝑛𝑘⃗∣𝑟⟩⃗ where 𝑛≤𝑁_PW+∑^𝐿𝑙=0(𝑁𝑅_𝑙−2). Approximately only a quarter of them is useful in further applications.

With respect to all-electron investigations the present schema has the advantage to allow the incorporation of the cores in a simple and rigid way. As core wave functions are completely localised within the atomic spheres, the left-hand side of (4) vanishes even when the coefficients 𝐴_𝑠𝑙𝑚 are nonzero. Consequently no plane waves exist outside the atomic spheres. By use of the radial differential equation, it can be shown that the core and valence electron wavefunctions are strictly orthogonal on each other.

A critical comparison in a recent paper [9] had the result that the MAPW scheme guarantees a higher accuracy than the APW schemes mostly used by other groups especially with regard to the EMD and CP.

2.2. Evaluation of the Electronic Structure of Cu

The valence electrons of Cu are of s-, p-, and d-type. Due to the dominant role of the d-electrons, we choose 𝑅_1,𝑙=2 to be almost atomlike, for example, both 𝑅_1,2 and 𝑑𝑅_1,2/𝑑𝑟 are made small on the surface of the atomic sphere by a properly chosen value of the energy in the radial differential equation. All other radial functions are generated with the recipe of Section 2.1. The electronic structure was self-consistently determined by use of the MAPW formalism yielding Bloch functions ⟨𝑟⃗∣𝑛𝑘⃗⟩ which due to their continuity within the atomic polyhedron are especially suited for the evaluation of the electron momentum density (EMD). For details we refer to previous investigations [10–13]. Complete self-consistent calculations for both the [Ne] and [Ar] cores are performed which treat 19 or 11 electrons, respectively, by the band structure scheme. In the present investigation it proves to be quite advantageous that, different from other schemes, the MAPW method yields Bloch states extending over an energy range of 8 Ry which are strictly orthogonal to each other.

As a peculiarity of the MAPW formalism in the angular momentum decomposition of the wave function within the APW sphere only 5 different properly chosen radial functions are sufficient for 𝑙≤2.

The density of the electrons 𝜌(𝑟)⃗ obtained from the occupied states is, as a consequence of the MAPW ansatz, a symmetric combination of plane waves everywhere and angular-dependent contributions within the atomic spheres. We have found it to be sufficiently accurate in the case of Cu to evaluate the latter contributions along 12 specially chosen directions [14–18] in the 1/48 of the unit sphere and to construct the crystal potential in the warped muffin-tin approximation by averaging over all angular-dependent contributions of these directions. In the sense of Fehlner et al. [15, 16] this set of directions allows the exact evaluation of the expansion coefficients of all Kubic harmonics (KHs) up to 𝑙=16. Outside the APW sphere the Coulomb potential is obtained without any further simplification, whereas the exchange and correlation as well as the Lam-Platzman(LP) [19] corrections could be treated in any accuracy wanted by use of a suitably chosen fine mesh of 𝑟⃗ points [12]. Consistent with the other intrinsic inaccuracies 916 points in the irreducible wedge are considered.

The Brillouin zone (BZ) integrations over all occupied states yielding the electron density 𝜌(𝑟)⃗, the Fermi energy 𝐸_𝐹, the total energy 𝐸_tot, and the EMD are approximated by a sum over MG(MG+1)(MG+2)/12 points properly chosen in the irreducible wedge [20–23]. These points constitute a simple cubic lattice of length (1/2MG)(2𝜋/𝑎). The characteristic ground-state properties are calculated by use of the LDA functional [3] and of the GGA functionals [4, 5, 24, 25].

Extensive numerical investigations with QMXSQ extending up to 24.0 and MG up to 40 have been performed to analyze how these intrinsic parameters influence the ground-state properties. We have found that the asymptotic behavior, that is, MG→∞ and QMXSQ→∞, of the total energy minimized in the LDA scheme is well-described by the expression𝐸_tot=𝑃₁+𝑃²MG^𝜈+𝑃₃𝑒^−𝜇QMXSQ.(7) Because of the nonlinear parameters 𝜇 and 𝜈, the fitting procedure is divided into two steps: (1) For fixed values of 𝜇 and 𝜈 the parameters 𝑃₁,𝑃₂, and 𝑃₃ have been determined by the least-squares method, and (2) the minimum of the sum of the squared residuals over the (𝜇,𝜈)-plane gives the optimal values of the nonlinear parameters. In the case of the [Ne] core, LDA and 𝑎=6.8312171Bohr 𝑃₁=−3275.9028707Ry, 𝑃₂=−1.220±0.160Ry, 𝑃₃=13.002±0.443Ry, 𝜈=2.726±0.044, and 𝜇=0.58658±0.00211 are found. The positive sign of 𝑃₃ follows from the minimal property of 𝐸_tot, whereas the negative sign of 𝑃₂ is probably caused by the partitioning of the 𝑘⃗ space relative to the Fermi surface. According to this formula the 𝐸_tot is obtained with an accuracy of 2 mRy provided QMXSQ is chosen larger than 16.4. Then about 75 plane waves constitute a Bloch function. This number certainly is a lower limit for other plane wavelike schemes. For a proper rate of the magnitude of this energy a comparison with the zero-point energy of Cu is helpful which within the Debye approximation amounts to 2.66 mRy (Debye temperature Θ = 336 K and atomic volume 𝑉=79.70Bohr³). With regard to the evaluation of the Compton profiles we have chosen the relatively large value MG=40. The corresponding 𝑘⃗-grid contains 5740 points within the irreducible wedge and gives the total energy 𝐸_tot and the Fermi energy 𝐸_𝐹 with an accuracy of 50 𝜇Ry and 58 𝜇Ry, respectively.

This choice of the parameters together with the use of more elaborate exchange-correlation functionals guarantees results which can be considered as best obtainable approximation to the LDA at least in the warped-muffin-tin approximations. Remaining differences to the experimental results must be ascribed to many-body effects not properly approximated by LDA or GGA. Finally a useful remark: as the correction terms in (7) have opposite sign, their total contribution can be made zero by a proper choice of both MG and QMXSQ, for example, MG=20 requires QMXSQ=17.96. Then 𝑃₁ is a good approximation to the asymptotic value of 𝐸_tot.

3. Some Ground-State Properties

Self-consistent calculations both in LDA and GGA were performed for up to 7 different lattice constants in the range [6.48, 7.10] a.u. in the fcc-structure. We have found that the total energy as function of the volume can be approximated with an absolute error of 10𝜇Ry or less by the Murnagham expression [26]𝐸_tot=𝑝₁+𝑝₂𝑉+𝑝₃𝑉^1−𝛾(8) with the four adjustable parameters 𝑝₁,𝑝₂,𝑝₃, and 𝛾. It describes a wide variety of behaviors: (i) 𝛾<0 yields a more or less strong decrease at small values of V and (ii) 𝛾=3 gives a symmetric energy curve around a certain minimum. The volume 𝑉₀, the isothermal bulk modulus 𝐵₀, and the total energy 𝐸_tot,0 at the equilibrium configuration are determined by𝐵₀𝑉₀=(𝛾−1)𝑝₃/𝑝₂^1/𝛾,=𝛾𝑝₂,𝐸_tot,0=𝑝₁+𝑝₂𝑉₀+𝑝₃𝑉^01−𝛾.(9) The nonlinear parameter 𝛾 in the Murnagham expression numerically coincides with the pressure derivative of the bulk modulus 𝑑𝐵/𝑑𝑝. As the latter is the third derivative of the total energy, a proper fit to the Murnagham expression requires energies 𝐸_tot with high intrinsic accuracy, say 1𝜇Ry, in the case of Cu. As displayed in the case of Al [12] the energy-volume dependence is almost rigidly shifted upwards in accord with (7).

In Table 1 the result of the fit for the different LDA and GGA functionals and for both core models are listed. The corresponding experimental values are given in the last row. The lattices constant at 𝑇=0K was derived by use of the linear expansivity measured by Himmler et al. [27] using its value at 18°C [28]. The bulk modulus is taken from ultrasonic experiments at 𝑇=4.2 K [29]. 𝛾=𝑑𝐵/𝑑𝑝 is a mean value of different measurements according to [30].

Table 1 

Ground-state magnitudes for different exchange correlation functionals and both choices of the core. 5740 𝑘⃗ points within the irreducible wedge, QMXSQ=15.5 a.u. vwn: LDA functionals according Vosko et al. [3], pbe: GGA functionals according Perdew et al. [5], pw91: GGA functionals according Perdew et al. [4], and KKR: KKR results [31].

As in previous investigations [11], the GGA overestimates the equilibrium lattice constant 𝑎₀, whereas LDA underestimates it. Band structure calculations with the smaller core [Ne] shift 𝑎₀ closer to the experimental values. Compared to other theoretical investigations [4, 32], which mostly underestimate 𝑎₀ by 2%, a considerable progress is achieved. The specific choice of the core has only a slight influence on the value of the bulk modulus 𝐵. Again LDA and GGA approach the experimental value from different sides where in LDA the difference is slightly smaller than in GGA, apart from the [Ar]pbe case which seems to be exceptional. In accord with the extremal principle 𝐸_tot assumes in the [Ne] configuration a smaller value than in the corresponding [Ar] configuration; in LDA this difference is rather small. With regard to the volume dependence of the 𝐵 described by the Murnagham parameter 𝛾 there is no unique behavior which is worth to be further discussed due the experimental uncertainty. In general the ground-state properties listed in Table 1 do not allow to favour one of the GGA functionals.

The Fermi surface of Cu consists of a distorted sphere at Γ and is multiply connected at L by necks. Of special interest are the extremal orbits: the belly orbits B[100], B[111], the dog’s bone DB[110], and the rosette R [100]. For further details see the review article [33]. In Tables 2 and 3 we present extremal areas A and their dilatation dependence 𝑑ln𝐴/𝑑ln𝑉 evaluated for different models and 𝑎=3.6030 Å = 6.80868 Bohr, the experimental lattice constant at zero temperature. It is remarkable how close the results are found with different exchange-correlation functionals and with both choices of the core. For comparison experimental values of 𝐴 and of 𝑑ln𝐴/𝑑ln𝑉 obtained by high precision de Haas van Alphen experiments [34, 35] are given in the last column. Apart from the neck orbit they differ from corresponding theoretical values by less than 1% that means that the topology of the Fermi surface is quite insensitive to details of the potential. The relatively large error in the neck result is due to the fact that the extension of this orbit sensitively depends on the locations of the state 𝐿₂′ relative to the Fermi level. A shift of approximately 8 mRy would completely reduce this deficiency. Interestingly, a recent study [36] showed that a similar order of magnitude of the shift needed to be applied to calculations of the Fermi surface of Ag. Compared to the FLAPW investigations [37] which give the extended orbits with errors up to 8% both in LDA and SIC a considerable improvement is achieved. As the authors of this investigation concede the extremely small area of the neck orbit is a distinct disadvantage of the SIC approximation.

Table 2 

Extremal areas of cross-sections of the Fermi surface in units of the free electron cross section evaluated for 𝑎=3.6030 Å.

Table 3 

Dilatation dependence of the cross-section 𝑑ln𝐴/𝑑ln𝑉 evaluated for 𝑎=3.6030 Å.

The cyclotron masses also turn out to be rather insensitive to the exchange-correlation functionals and the choice of the core. Typical values are 1.347, 1.404, 0.3931, −1.210, and −1.278, for the ⟨100⟩belly,⟨111⟩belly,⟨111⟩neck,⟨110⟩dogsbone, and the⟨100⟩rosette, respectively.

Using the Fermi surface radii along 48 suitably chosen directions on the unit sphere the coefficients of the expansion in symmetrized combinations of plane waves as defined in [39] have been evaluated by use of a singular value decomposition. The result listed in Table 4 confirms the previous results that the shape of the Fermi surface only weakly depends on the choice of the core and of the exchange correlation functionals. It is not surprising that there does not exist the faintest similarity with the values derived by Coleridge and Templeton [34, 40] from their high-precision data as their system of linear equations involved in the least-squares fitting procedure is not well conditioned. The preceding results are further supported by the fact that the energetic locations of some special states only weakly depend on the various approaches apart from the widths of the 3s and 3p core states which are 0.0017 and 0.0033 Ryd, respectively, in the case of the [Ne]-core.

Table 4 

Coefficients for seven term fits of the Fermi surface for different exchange-correlation functionals and both choices of the core. 5740 𝑘⃗ points within the irreducible wedge, QMXSQ=15.5 a.u. Notation following Halse [38]. Last column: root mean of sum of squared energy residua.

Besides the total energy the charge density is the fundamental quantity of the density-functional formalism. Fortunately, it can directly be measured by the form factors of X-ray scattering or by 𝛾-ray diffractometry. In Table 5 theoretical values obtained for different models with the lattice constant 3.614924 Å which is close to that of the room temperature value [28] are compared with measurements by Schneider et al. [41] and Temkin et al. [42]. Again we find that this quantity is quite insensitive to the special model which is certainly due to the fact that the charge density is globally controlled by the requirement of charge neutrality. It is within the experimental margin that the data by Schneider slightly favour the model with the [Ne] core, whereas those of Temkin favour more the other core choice.

Table 5 

Form factors for different exchange- correlation functionals and both choices of the core. 5740 𝑘⃗ points within the irreducible wedge and QMXSQ=15.5 a.u. (a) experimental values at 𝑇=50 K (b); experimental values at room temperature.

4. The Compton Profile

4.1. Basic Considerations

In the so-called impulse approximation [43], the Compton profile is defined by the double integral in the momentum space of the ground state momentum density 𝜌(𝑝)⃗𝐽_⃗𝑛(𝑞)=^∞−∞𝜌𝑝_𝑥,𝑝_𝑦,𝑝_𝑧=𝑞𝑑𝑝_𝑥𝑑𝑝_𝑦.(10) For simplicity we use a Cartesian coordinate system with the 𝑝_𝑧 axis along the scattering direction 𝑛⃗. 𝑞 is the momentum transfer. 𝜌(𝑝)⃗ is the ground state momentum density defined by the sum 𝜌𝑝⃗=2_𝑛⃗𝑘|||𝑝⃗∣𝑛𝑘⃗|||²𝑓_𝑛⃗𝑘(11) extended over the first Brillouin zone and over all bands. The preceding factor is due to the spin degeneracy, 𝑓_𝑛⃗𝑘 denotes the Fermi-Dirac function.⟨𝑝⃗∣𝑛𝑘⃗⟩ is the Fourier transform of the Bloch function |𝑛𝑘⃗>|||𝑛𝑘⃗=1(2𝜋)^3/2_𝑉𝑐𝑒^{𝑖⃗𝑝⋅⃗𝑟}𝑟⃗∣𝑛𝑘⃗𝑑³𝑟,(12) where 𝑉_𝑐 denotes the volume of the elementary cell. As usual, the Bloch functions are normalized within the elementary cell. ⟨𝑝⃗∣𝑛𝑘⃗⟩ is only nonzero if 𝑝⃗ is equal to 𝑘⃗ up to a reciprocal lattice vector. The definition according equation (12) implies that the integral of EMD over the momentum space gives the number of electrons within the unit cell (EMD sum rule) _∞𝜌𝑝⃗𝑑³𝑝=2_𝑛⃗𝑘𝑓_𝑛⃗𝑘.(13) Correspondingly the CP sum rule reads ^∞0𝐽_⃗𝑛(𝑞)𝑑𝑞=2_𝑛⃗𝑘𝑓_𝑛⃗𝑘(14) for any direction 𝑛⃗. Because the MAPW functions as well as their first derivatives are strictly continuous everywhere their Fourier transform asymptotically decreases at least with 1/|𝑝|⃗⁴ as the exact Bloch functions. The APW or KKR schemes yield ⟨𝑝⃗∣𝑛𝑘⃗⟩ decreasing at best with 1/|𝑝|⃗³ [44]. Nevertheless, the evaluation of the incomplete integral equation (10) needs some care even in the present case. We proceed in two steps. In Section 4.2 the influence of the spherically averaged EMD on the Compton profiles is investigated which turns out to be the dominant part. Section 4.3 deals with the directional CP’s.

4.2. Contribution from the Spherically Averaged EMD

In momentum space the mean value of 𝜌(𝑝)⃗ over the sphere of radius 𝑝=|𝑝|⃗ gives the spherically averaged EMD 𝜌_av(𝑝)=14𝜋_{sphere|⃗𝑝|=𝑝}𝜌𝑝⃗𝑑𝜔_⃗𝑝.(15)𝑑𝜔_⃗𝑝: unit surface element. It likewise satisfies the EMD sum rule 4𝜋^∞0𝜌_av(𝑝)𝑝²𝑑𝑝=2_𝑛⃗𝑘𝑓_𝑛⃗𝑘.(16) Its contribution to the Compton profile is given by the one-dimensional integral 𝐉_av(𝑞)=2𝜋^∞𝑞𝜌_av(𝑝)𝑝𝑑𝑝,(17) which is evidently spherically symmetric, too. To avoid confusion with the spherically averaged Compton profile 𝐽_⃗𝑛(𝑞) we have used a bold font. As 𝐉_av(𝑞) also satisfies the CP sum rule ^∞0𝐉_av(𝑞)𝑑𝑞=2_𝑛⃗𝑘𝑓_𝑛⃗𝑘(18) we expect that it is a good approximation for the spherically averaged value of CP, 𝐽_⃗𝑛(𝑞).

Using the concept of special directions described in a previous paper [18] 𝜌_av(𝑝) is directly evaluated from the EMD along 48 directions in the irreducible wedge which are defined by requiring that 111 Kubic harmonics 𝐾_𝐿(𝑝⁰) with 𝑙≤68 have nonvanishing residua in the sense as required by Bansil [14–17]. For any direction 𝑝⁰ the contributions from the valence electrons are evaluated on a fine mesh up to 129.0 a.u. by use of MAPW results where for each value of 𝑝⃗ a Bloch vector 𝑘⃗ is found by subtracting a suitably chosen reciprocal vector 𝐾⃗, for example, 𝑝⃗=𝑘⃗+𝐾⃗. With the proper weights of the 48 directions formula the spherical average 𝜌_av and by use of (17) 𝐉_av(𝑞) are obtained. In contrast to the free electron metals Li [13] and Al [12], 𝜌_av has no discontinuities due to possible Fermi breaks. Figure 1 shows 𝐉_av(𝑞) augmented by the core profile on a semilogarithmic scale up to 𝑞= 129.0 a.u. In accordance with the general considerations on the continuity of the MAPW functions [23] the contribution from the valence electrons decays by 9 orders of magnitude in the momentum range up to 129 a.u. It is expected that valence profiles obtained by KKR or FLAPW show a weaker decay. A reliable test of the numerical accuracy is the CP sum rule which gives the value 11.002 in the momentum range considered and for the [Ar] core. The corresponding value in the case of the [Ne] core is 19.002. In Figure 1 the dashed-dotted curve is the Compton profile in LDA simply obtained from the SCF-MAPW density by use of 𝐽_LDA(𝑞)=2_𝑉𝐶𝜌𝑟⃗𝑗_𝑖𝜌𝑟⃗;𝑞𝑑³𝑟.(19)𝑗_𝑖(𝜌₀;𝑞) is the CP per electron of an interacting electron gas of density 𝜌₀. As for case Li [13] and Al [12] it roughly coincides on the semilogarithmic scale with the all electron CP in the case of Cu up to 𝑞=12 a.u. Again we find confirmed that the LP corrections which are based on the same functional 𝑗_𝑖 apply to the all-electron profiles.

Figure 1 

𝐉av(𝑞) of Cu on a semilogarithmic scale. Full curve: contribution from the valence electrons only, and dashed curve: all electron CP, dashed-dotted curve: CP in LDA using (19). In units of electrons/Bohr.

The Compton profiles of the valence electrons along the three principal directions have been obtained by Sakurai et al. [7] by subtracting the KKR core contributions from the measured values. Using the weights quoted in Section 2 of the previous work [12] the spherical average of the valence profiles has been derived up to 𝑞= 8 a.u. It is remarkable that KKR and the present MAPW core contributions do not differ significantly from the results of Biggs et al. [45] based on Hartree-Fock calculations. In Figure 13 of the previous investigations [18], theoretical profiles defined by (12) are displayed after having been folded with a Gaussian of FWHM=0.12 a.u. to reflect the experimental resolution. The overall shapes of both profiles are quite similar. To make this more clear in Figure 2 the differences between the theoretical, either corrected according Lam Platzman or not, and the measured profiles are plotted. The uncorrected results, displayed by the broken curve behave quite similarly as the profiles derived with the KKR scheme [7]: at low momenta the difference amounts almost up to 0.28, after a crossover around 1.0 a.u. the deviations reach −0.07 and approach the value 0.02 near 8.0 a.u, all in electrons/Bohr. The weights just mentioned were also used to evaluate another spherical average from the CPs along the principal directions considered in details in Section 4.3. In the low momentum region it runs up to 0.033 electrons/Bohr below the broken curve. A distinct step in the direction of the experimental results is due to the LP corrections using the momentum density function of the interacting free electron gas either as proposed by Farid et al. [46] or by Ortiz and Ballone [47]. At small values of 𝑞 we observe a reduction by roughly 0.05 and near 𝑞=1.0 a.u. an increase by 0.02 electrons/Bohr. Beyond 2.0 a.u. these corrections are almost negligible. Different from the cases of Li [13] and Al [12], they are smaller but not negligible as assumed in [7] and both parameterizations produce similar corrections. The results displayed in Figure 2 are obtained within the LDA and the [Ne]-core but the other cases behave similarly, within the plotting accuracy. More quantitative insights give Table 6 which lists the deviations 𝛿𝐽 at 𝑞=0 and the square root of the variance defined by 𝜎=^07.85𝐉_av(𝑞)−𝐽_av,exp(𝑞)²𝑑𝑞(20) for each model and for both LP parameterizations. As for the form factors the Compton profiles turn out to be only slightly influenced by details of the model, for example, LDA, GGA, or the choice of the core, certainly as consequence of the CP sum rule. The [Ne]-vwn model best competes with all others in the long-wave limit of 𝐽(𝑞), whereas in the overall agreement characterized by the variance the [Ne]-pw91 model is the favorite. In all cases the momentum density proposed by Farid et al. [46] give a better agreement. Once more we want to emphasize that the preceding results have been obtained without any truncation, nevertheless they are considerably less CPU-time expensive than the investigations in the following section. But they fully reflect the essential trends how the CPs depend on different models and assumptions.

Table 6 

Deviations between theoretical and experimental CPs, 𝛿𝐽, at 𝑞=0 and the square root of the variance evaluated in the range 0≤𝑞≤7.85 a.u. These values have to be compared with the profile integrated over the same range amounting 10.896 a.u. Lower index 𝐹 and 𝑂 parameterizations according to Farid et al. [46] and Ortiz and Ballone [47], respectively. All in units electrons/Bohr.

Figure 2 

Differences between the calculated and the experimental spherically averaged valence electron Compton profiles. Full curve: LP corrections according to [46] included, dotted curve: LP corrections according to [47] included, broken curve: LP corrections excluded, dashed-dotted curve: spherical average of the CPs along the principal directions, and long-broken curve: APW result [48]. CPs in electrons/Bohr.

4.3. Directional Compton Profiles

As in a previous investigation [11] we make use of the fact that outside a circle with sufficiently large radius 𝑃_max the momentum density 𝜌(𝑝)⃗ is in a good approximation spherical with the consequence that the contribution from outside the circle is defined by Δ𝐽_⃗𝑛(𝑞)=2𝜋^∞_{𝑞2+𝑃2max}𝜌(𝑡)𝑑𝑡=_𝐉av𝑞²+𝑃^2max(21) with the 𝐉_av shown in Figure 1. This procedure only produces a truncation error from the nonspherical contribution of the EMD which is roughly by the factor 10 smaller than the mean value. This consideration also allows a rough estimate of the accuracy of the KKRcalculations [7] which are based on the EMD evaluated on a mesh generated by 2421 reciprocal lattice vectors extending to about 13.0 a.u., an upper limit markedly different from 20.0 quoted by these authors. (In the fcc reciprocal lattice the last shell includes 2421 reciprocal lattice vectors of the type (10,6,6) yielding a 𝑝⃗ vector of maximal length √(172.0)(2𝜋/𝑎)). According to Figure 12 of the previous paper [18] the neglected contribution Δ𝐽(𝑞) amounts to 0.08 and 0.06 electrons/Bohr for 𝑞=0 and 𝑞=8.0 a.u., respectively. From this estimate follows that Sakurai et al. [7] have considerably overestimated their accuracy by claiming that the final Compton profiles are accurate to a few parts in 10⁴.

In our investigation we have chosen 𝑃_max=22.5 a.u., 26.0 a.u. and 21.8 a.u., for the direction 𝑛⃗ along [100],[110] and [111], respectively and have determined Δ𝐽 by spline interpolation of 𝐉_av. The evaluation of the contribution to the integral equation (21) is quite tedious since for each value of 𝑝⃗ the corresponding Bloch function |𝑛𝑘⃗> with 𝑘⃗+𝐾⃗ is needed. A considerable reduction of the numerical work is possible by use of a two-dimensional grid of points originally proposed to perform Brillouin zone integrals and which have been already used in Section 2. As shown in [6] this partition of the 𝑝⃗ space implies discrete values of 𝑞 defined by (1/2MG)(2𝜋/𝑎)𝑞_[𝑖]𝐼 with the integer 𝐼=0,1,2 and 𝑞_[100]=1, 𝑞_[100]=1/√2, and 𝑞_[111]=1/√3. To comply with the experiments [7] and the aim to reliably identify fine structures in the profiles we have chosen MG=40. This gives 5740 𝑝⃗ points in the irreducible wedge of the BZ and more than 870 million points lying within the sphere of radius 𝑃_max in contrast to the 48×1183×2421 points or to the 48×505×3143 points considered in [7] or in [37], respectively. Finally the integration over the two-dimensional grid is approximated by a scheme [23] similar to that proposed by Gilat and Raubenheimer [49, 50] which allows to take into account the exact position of the Fermi break within the corresponding parallelepipeds. To avoid spurious contributions, the extensions of these parallelepipeds in 𝑝_𝑧 are chosen small compared to (1/MG)(2𝜋/𝑎).

The Compton profiles of the valence electrons along the principal directions are displayed in Figure 3. LP corrections have been included in the full curve. Both theoretical curves are convoluted with a Gaussian with full width at half maximum (FWHM) of 0.12 a.u. according to Sakurai et al. [7]. For comparison the corresponding measured CPs are plotted which have been renormalized in [111]-direction by 0.438% according to the results [Ne]vwn with LP corrections included. Obviously the LP corrections slightly improve the agreement with the experimental results. But they still remain higher than the measured ones up to 𝑞=0.75 a.u. and lie slightly below them in the momentum range [1.32–5.0] a.u. Both theoretical profiles show a weak wiggly structure in the [110] profile.

(a)

(b)

(c)

Figure 3 

The valence electron Compton profiles along [100,110] and [111] obtained by subtracting all core contributions. Full curves: theoretical result with LP corrections included, dashed curves without LP corrections; both convoluted with a Gaussian with FWHM of 0.12 a.u. Crosses connected by thin curves: experimental results [7]. Dotted curves: APW results [48].

Figure 4 showing the differences between the theoretical and measured profiles along the principal directions makes the deviations more clear. The uncorrected results shown in panel (a) behave quite similarly to those derived with the KKR-scheme by Sakurai et al. [7]. As far as comparison is possible the MAPW calculations yield a slightly larger deviation near 𝑞=0. According to the estimate at the beginning of this section this is a consequence of the smaller value of the cut-off 𝑃_max used in the KKR investigations [7]. The LP-corrections based on the momentum density functions [46, 47] reduce the deviations in the momentum range up to 𝑞=0.8 a.u., especially near 𝑞=0 by 0.05 a.u. but have only slight influence beyond 1.5 a.u. As already found in Section 4.2 the momentum density of the interacting electron gas proposed by Ortiz and Ballone [47], Figure 4(b), and by Farid et al. [46], Figure 4(c), yield similar curves.

(a)

(b)

(c)

Figure 4 

Differences between the calculated and the experimental valence-electron Compton profiles. (a) LDA results without LP corrections, (b) LP corrections [47] included, (c) LP corrections [46] included. [100]: full curve, [110]: dashed curve, [111]: dotted curve. All in units of electrons/Bohr.

The orientation dependence of the Compton profile is quite small as shown in Figure 5 where the profile differences along the principal directions obtained by MAPW are compared with the experimental results. By this display not only the core contributions and the many-body aspects treated by the LP procedure but also the asymptotic tail of the momentum density in the integrand of (21) are completely absent in the theoretical curves since all these contributions only depend on the absolute magnitude of the momentum 𝑞⃗. Our results are quite similar to those of the previous investigations of the author [6] and the KKR work [7] and show that the anisotropy of the CPs is only slightly influenced by the band structure schemes and details of the potential. The agreement with the experimental results is satisfactory. In the both differences Δ𝐽_[111]-[110] and Δ𝐽_[110]-[100] not only the periodicity of the oscillations up to 𝑝=4.0 a.u. but also the locations of some fine structures are properly described. The heights of the oscillations are distinctly overestimated. The broken curve in Figure 5 is obtained by use of the Fourier-Hankel method [51] based on the expansion of the EMD in Kubic harmonics 𝐾_𝑙 in the present case up to 𝑙=28. With the notation of [13] Section 4 we get, for example, Δ𝐽_[110]-[100](𝑞)=2𝜋_𝐿𝐹_[110],𝐿−𝐹_[100],𝐿𝑔_𝐿(𝑞).(22) Because both (𝐹_[110],𝐿 and 𝐹_[100],𝐿)𝑔_𝐿 are defined by azimuthal integrals over Kubic Harmonics, see (11) of [13], their differences can be obtained without loss of accuracy. The functions 𝑔_𝐿(𝑞) are found by folding the EMD expansion in KHs with the Legendre polynomial of the same degree. From the close agreement of both theoretical curve we learn that the Fourier-Hankel method, which is considerably less CPU demanding than the two-dimensional integration according to equation (21), is also well suited for the evaluation of the directional profiles 𝐽_⃗𝑛 provided the momentum density 𝜌(𝑝)⃗ is known along a sufficiently dense grid of orientations 𝑝⁰.

(a)

(b)

Figure 5 

Directional Compton differences. Upper panel: [111]-[110] lower panel: [110]-[100]. Full curve: theoretical results, broken curve: Fourier-Hankel reconstruction, long-broken curve APW results [48] string of crosses: experimental data. All in units of electrons/Bohr.

The second derivatives of the theoretical and experimental profiles 𝑑²𝐽(𝑞)/𝑑𝑞² derived by spline interpolation and convoluted with a Gaussian of FWHM=0.12 a.u. show good overall agreement as to both the positions and the heights of the peaks. However, the distinct maxima occurring in[100] and [110] directions cannot be assigned to any Fermi diameter. In contrast to Li and Al, measurements of the Compton profiles are not suited to map out the Fermi surface in the case of Cu.

4.4. Another Explanation of the Disagreement Near 𝑞=0

Measured profiles are mostly corrected by computer simulation of the multiple Compton scattering events which need specific assumptions for the scattering process. Different from Sakurai et al. [7], Bauer and Schneider [52] used a semiempirical approach which is based on the theoretical valence electron CP derived by using a basis of Gaussian orbitals [53]. It is justified by the fact that the corrected profiles for samples of different thicknesses coincide within a margin of 0.03 a.u. Figure 6 shows the profiles for the [110] direction after being corrected compared with the theoretical profile. Both profiles mostly differ by less than 0.02 a.u. over the whole momentum range up to 5.0 a.u. Near 𝑞=0 the discrepancy is completely reduced making other explanations obsolete.

Figure 6 

Difference between experimental [52] and theoretical all-electron profiles along the [110] directions. All in units of electrons/Bohr, profile difference enlarged by a factor of 10.

4.5. A Comparison with the KKR Investigations [7]

Apart from the choice of some intrinsic parameters as already mentioned the KKR investigations [7] are essentially different from the present one in various points. (i) They are based on the muffin-tin approximations with the consequence that the wave functions outside the APW spheres are not variationally given and, therefore, are less accurate than those obtained by MAPW. (ii) Less sophisticated exchange and correlation functionals [54] are used in the SCF scheme and in the evaluation of the Lam-Platzman corrections. (iii) No details are given how the latter are determined, especially outside the APW-spheres. (iv) The Compton profiles are evaluated by an improper integral which has been truncated at a rather small value of |𝑝|⃗. (v) In contrast to a generalizations of the Gilat-Raubenheimer scheme which saves considerably computer time [55] the tetrahedron method has been used for the evaluation of the momentum density. Despite these differences with respect to methodology and computational accuracy the overall similarity of the profiles obtained by completely different schemes is remarkable. In the case of the noble metal Cu investigations of the Compton profiles are no proper mean for a decision in favour of one theoretical scheme.

5. Conclusions

The preceding considerations had the result that there exists no unique approach which gives the best agreement with the experimental facts. LDA with the [Ne]-core describes the energetics and the form factors considerably well whereas the Compton profiles are slightly favoured by the GGA functionals proposed by Perdew et al. [4] and the [Ne] core. In a consistent treatment the former approach requiring a larger set of variational functions is to be preferred. This dilemma is certainly due to many-body effects not properly taken into account. Undoubtedly, those models and approximations should be favoured which explain simultaneously different phenomena sufficiently well. The full-potential linearized APW investigations by Kubo et al. [37] do not meet with these criteria as they show a quite good agreement with the measured profiles [7] when self-interaction corrections are taken into account, as in the case of free atoms [56], but on the other hand, these corrections change the band structure, especially near the Fermi level, in a rather unconventional way.

Our hope that one of the approaches will essentially reduce the disagreement of the profiles near 𝑞=0 that has not been fulfilled. The present investigations rule out that the treatment of the core contributions is the origin of the disagreement: as mentioned the core contributions evaluated by MAPW differ from Biggs atomic profiles by less than 1%. Nevertheless, the experimentalists are recommended to report in future also the core CP which they use to shift their data. Thus it remains an open question where the origin of the discrepancy lies. It may either be caused by the handling of the measured data or by many-body effects which are not taken into account by the Lam Platzman corrections.

Acknowledgments

The author would like to thank Professor Y. Sakurai for providing results of their high-resolution Compton-scattering study. Dr. Herbert Stöhr has accompanied these and the preceding investigations with great engagement. For this support and the critical reading of the paper many thanks are due. The help of Dr. Reinhold Bader as well as the generous hospitality of Professor Jan von Delft and Professor U. Schollwöck are gratefully acknowledged.