text:
# An ab-initio many-body method for electronic structure calculations of solids. I. Description of the method

## I Introduction

Most existing ab-initio (first-principles) methods for the numerical calculation of the electronic properties of solids are based on density functional theory (DFT)HK64 (1), which in principle is exact and properly takes into account many-body effects involving the Coulomb interaction between the electrons; for an overview on the present status of DFT we refer to the booksDrGross90 (2, 3). But, in general, the functional dependence of the kinetic energy and the exchange and correlation part of the Coulomb (interaction) energy on the electron density are not known explicitly, and hence additional approximations and assumptions are necessary. A well established additional approximation is the local density approximation (LDA)KS65 (4) (or local spin-density approximation, LSDA, for magnetic systems), which assumes that the electron density of the interacting system is the same as that of an effective non-interacting system and that the exchange-correlation potential depends only on the electronic density locally. Even then, the functional dependence of the exchange-correlation energy on the density is not known in general, and it is usually necessary to make an ansatz for the exchange-correlation functional, which is based on the homogeneous electron gas. The LDA goes beyond the simplest electron-gas approximation, the Hartree-Fock approximation (HFA), in that correlation energy (i.e. the part of the interaction energy beyond HFA) is explicitly taken into account. More recent generalized gradient approximation (GGA) calculations also make some improvements in the correlation energy that slightly improve on the homogeneous electron gas approximation. On the other hand, the exact HFA exchange potential is non-local, i.e., it depends on the electron wave functions and density at all other positions, an effect which the local LDA exchange potential misses. Also, LDA misses some cancellation between the exchange and correlation parts of the total energy (self-interaction corrections). However, in practice LDA-treatments are simpler than HFA-calculations (especially for metals), because local exchange is easier to treat than non-local exchange, and are usually in better agreement with experiment. Therefore, DFT-(LDA-like) treatments have been far more common than HFA during the past few decades, even in quantum chemistry (with a long tradition of methods based on HFA).

Ab-initio DFT calculations based on the L(S)DA have been very successful for many materials and ground-state properties such as crystal structure, ground state and ionization energy, lattice constant, bulk modulus, crystal anharmonicityRose84 (5), magnetic moments, and some photo emission spectra. However, there are also important limitations. For example, LDA predicts a band gap for semiconductors that is almost a factor of two too small. In addition, for many strongly correlated (narrow energy band) systems such as high-temperature superconductors, heavy fermion materials, transition-metal oxides, and $`3d`$ itinerant magnets, the LDA is usually not sufficient for an accurate description (predicting metallic rather than semiconducting behavior, failing to predict quasi-atomic-like sattelites, etc).

Therefore, it is justified and important to look for new, better ab-initio methods and improvements that go beyond L(S)DA. In fact, there have already been several attempts to improve LDA, such as gradient corrections, non-local density schemes, self-interaction corrections, and the GW approximation. Gradient correctionsDrGrossKap7 (6) approximately account for the fact that the electron density is not constant but $`𝐫`$-dependent in an inhomogeneous electron gas; estimates on the influence of perturbations of the homogeneity on the ground state energy of the homogeneous electron gas are then used to construct a new functional for the exchange-correlation potential containing $`n(𝐫)`$ terms. The non-local density schemes go beyond LDA by considering that the exact exchange-correlation potential $`V_{\text{xc}}(𝐫)`$ cannot depend only on the density $`n(𝐫)`$ at the same position $`𝐫`$ but should depend also on the electron density at all other positions $`n(𝐫^{})`$. Usually the new ansatz for the functional of the exchange-correlation energy contains the pair correlation function or the interaction of the electrons with the exchange-correlation hole.Gunnarsson80 (7, 6). The recently developed exact exchange (EXX) formalismVogl97 (8) cancels the spurious (unphysical) electronic self-interaction present in LDA and gradient corrected exchange functionals. A standard method for ab-initio calculations of excited states is the GW-approximation (GWA)Hedin (9, 10). Denoting the one-particle Green function by $`G`$ and the screened interaction by $`W`$, the GWA is an approximation for the electronic selfenergy $`\mathrm{\Sigma }GW`$, which is correct in linear order in $`W`$ and can diagramatically be represented by the lowest-order exchange (Fock) diagram. The one-particle Green function $`G`$ is usually obtained for the effective one-particle LDA Hamiltonian.

Recently there have been several attempts to combine ab-initio LDA calculations with many-body approximations.AZA91 (11, 12, 13, 14, 15, 16, 17, 18, 19, 20) All of these recent developments add local, screened Coulomb (Hubbard) energies U between localized orbitals to the one-particle part of the Hamiltonian obtained from an ab-initio LDA band-structure calculation, but differ in how they handle the correlation part. In the earliest attempts, the LDA+U methodAZA91 (11) used essentially a static mean-field-like (or Hubbard-I-like) approximation for the correlation. The simplest approximation beyond Hartree-Fock, second-order perturbation theory (SOPT) in U, was usedSASS91 (12, 13, 16, 19) to study the electronic properties of 3d-systems (like Fe and Ni) and heavy fermion systems (like UPt3). The LDA++ approachLK98 (15, 17, 18) has a similar strategy, but uses other many-body approximations to treat the correlation problem, namely, either the fluctuation-exchange approximation (FLEX) or the dynamical mean-field theoryGKKR96 (21) (DMFT); for a review on this recently very successful LDA+DMFTNHBPAV00 (20) approach see Ref. HeldVollhardtetal, 22. These approaches, including the LDA+U, have in common that they have to introduce a Hubbard U as an additional parameter and hence are not really first-principles (ab-initio) treatments. Although they use an LDA ab-initio method to obtain a realistic band structure, i.e., single-particle properties, Coulomb matrix elements for any particular material are not known, and the Hubbard U remains an adjustable parameter. In addition, since some correlations are included in LDA as well as by the Hubbard U, it is unclear how to separate the two effects and double-counting of correlation may be included in these approximations.

In this paper we suggest a different, alternative approach and propose a new type of ab-initio treatment. Because we want to avoid DFT-LDA due to its problems with double counting and self-interactions, we instead suggest starting from a realistic Hamiltonian in second quantization and directly applying many-body methods to this Hamiltonian. The underlying many-body theory in DFT-LDA involves the homogeneous electron gas, from which a functional for the exchange-correlation potential is derived. But we want to apply the many-body theory to a Hamiltonian tailored to the material of interest, i.e., study a problem of interacting electrons on a lattice and therefore an inhomogeneous electron gas. To get the realistic, microscopic Hamiltonian in second quantization the one-particle and Coulomb matrix elements have to be calculated. For that purpose one needs a suitable one-particle basis. For a solid it is reasonable that either Bloch or Wannier functions should form a suitable one-particle basis set. We suggest starting from a set of maximally localized first-principles Wannier functions, because the good localization of these wave function should make the one-particle (tight-binding) and two-particle (Coulomb) matrix elements important only for a finite number of near-neighbor shells (on-site, nearest, next nearest neighbors, etc.) such that only a few matrix elements have to be actually calculated. It is possible to determine these Wannier states for an effective (auxiliary) one-particle Hamiltonian. For example, one could use the true non-interacting Hamiltonian, the Hartree Hamiltonian, or even the LDA Hamiltonian. In principle, the final results should not depend on the choice of this auxiliary one-particle Hamiltonian, because any set of Wannier functions should span the whole one-particle Hilbert space. But, in practice, one typically works only within a restricted (finite dimensional) subspace. Therefore, the auxiliary one-particle Hamiltonian and the resulting one-particle basis need to be chosen in such a way that the physically most important subspace is spanned. In practice, we have found the Hartree Hamiltonian to be a suitable choice for the auxiliary, effective one-particle Hamiltonian.

Our suggested procedure is schematically outlined as follows: First, a traditional band-structure calculation is used for the auxiliary one-particle Hamiltonian to calculate the eigenstates in the form of Bloch functions. Wannier functions are closely related to these Bloch functions via a unitary transformation, and thus span the same one-particle space as the Bloch functions. However, since the phases of the Bloch functions are arbitrary, Wannier functions are not unique. Their non-uniqueness (gauge freedom) is used to construct “maximally localized Wannier functions” using a method proposed by Marzari and VanderbiltMV97 (23). A proper localization of the Wannier functions is important, because only then do the standard assumptions that are frequently used in model treatments hold, e.g., that both one-particle (tight-binding) and two-particle (Coulomb) matrix elements are important only on-site and for a few neighbor shells. Next, in order to describe the physical Hamiltonian in second quantization, one has to evaluate a restricted number of hopping and Coulomb matrix elements within the basis of the well localized Wannier functions. Making use of the special representation of the wave functions inherent to the specific band-structure method used (for instance in terms of spherical harmonics, plane waves or Gaussian orbitals) can simplify the actual computation of these matrix elements; for the Coulomb matrix elements we propose, in particular, an algorithm based on fast Fourier transformation, which works independent of the special band-structure method. The result of this procedure is a multi-band, second-quantized Hamiltonian in a Wannier representation (a kind of extended, multi-band Hubbard model) describing interacting electrons on a lattice. Unlike model studies of correlated lattice electrons, the matrix elements of this Hamiltonian are not free adjustable parameters, but are obtained from first principles for the given material. This Hamiltonian should be evaluated by an appropriate many-body methods. The simplest approximation is the mean-field Hartree-Fock approximation (HFA), which exactly treats exchange and excludes self-interaction effects. Unlike first-quantization treatments of HFA, which require the solution of a Schrödinger equation with a non-local potential, in second-quantization only the selfconsistent evaluation of expectation values is required (once an appropriate one-particle basis has been chosen). But, of course, HFA is certainly not sufficient for realistic materials properties; higher order correlations have to be taken into account, for example, perturbationally by resummations of Feynman diagrams (RPA-bubble diagrams leading to screening, ladder diagrams, etc.) or employing recently developed non-perturbational methods like dynamical mean-field theory (DMFT) (if the basic assumptions, for instance concerning the locality of the most important interaction terms, hold for the specific many-body Hamiltonian).

This paper is organized as follows: In Section II we shortly repeat some basic notations, give the Hamiltonian in first and second quantization, describe DFT and LDA, Wannier- and Bloch-states as possible basis states spanning the one-particle Hilbert space and the necessary truncation. Section III describes how a suitable restricted basis set is obtained by solving the Schrödinger equation for an auxiliary one-particle Hamiltonian, which should be chosen appropriate for the material of interest and incorporating information about the lattice structure, lattice constant, total electron number, etc. In Section IV we discuss how another equivalent basis set of maximally localized Wannier functions can be obtained from this one-particle basis. Within this Wannier basis ab-initio one-particle (tight-binding) and two-particle (Coulomb) matrix elements can be calculated, as outlined in section V. The resulting many-body Hamiltonian in second quantization with ab-initio parameters is then ready to be studied by suitable many-body approximations, which is the subject of Section VI. We then discuss and compare the advantages and disadvantages of our proposed ab-initio method relative to other ab-initio methods in Section VII. The following paper (II) contains an initial simple application (an unscreened HFA) of this approach to the 3d transition metals Fe, Co, Ni and Cu.

## II Hamilton Operator in first and second quantization

A system of $`N_e`$ interacting (non-relativistic) electrons can be described by the Hamiltonian

$$H=T+V+W=\underset{i=1}{\overset{N_e}{}}\frac{𝐩_i^2}{2m}+\underset{i=1}{\overset{N_e}{}}V(𝐫_i)+\underset{i>j}{}\frac{e^2}{|𝐫_i𝐫_j|}.$$
(1)
The first part $`T`$ is the kinetic energy of the electrons. The $`V(𝐫)`$ describes the external one-particle potential, which for molecules and solids is solely caused by the positive nuclei or ions and–within the adiabatic or Born-Oppenheimer approximation–can be regarded as a static potential. For solids, i.e., for electrons in a crystal, $`V(𝐫)`$ is a periodic potential. The third part $`W`$ is the Coulomb interaction between the electrons.

For the interacting system it is impossible to exactly solve the full Schrödinger equation and to calculate the many-body eigenenergies (or even the ground state energy) and the totally antisymmetric many-particle wave functions, which can be represented as a linear combination of Slater determinants (a configuration-interaction or CI approach to correlation). This is practical only for a very small number of electrons. The basic idea of DFT is to avoid determining the antisymmetric (ground state) wave function and instead to focus on the (ground state) density $`n(𝐫)`$. This requires knowing the functional dependence of the kinetic energy $`T`$, the one-particle potential $`V`$, and the interaction energy $`W`$ on the density $`n`$. From this the ground state density can be determined from a variational principle. But in general this DFT-concept can only be approximately applied, since the functional dependence of the kinetic energy and the exchange-correlation energy on the density is unknown. Therefore additional assumptions and approximations are necessary, such as the LDA where the density $`n(𝐫)`$ is the same as a non-interacting electron system in an effective one-particle potential and the exchange-correlation energy depends only locally on the electron density. From a fundamental point of view these additional assumptions are hard to justify and uncontrolled, and an estimate of the error or systematic corrections are very difficult to obtain.

In this paper we choose yet another approach that employs the (very elegant) formalism of “second quantization”, which automatically accounts for the antisymmetry through the fermion anticommutation relations. In second quantization the full many-body Hamiltonian can be written as:

$$H=\underset{i,j,\sigma }{}t_{ij}c_{i\sigma }^{}c_{j\sigma }+\frac{1}{2}\underset{i,j,k,l,\sigma ,\sigma ^{}}{}W_{ij,kl}c_{i\sigma }^{}c_{j\sigma ^{}}^{}c_{k\sigma ^{}}c_{l\sigma }$$
(2)
Here $`i,j,k,l`$ denote a complete set of one-particle orbital quantum numbers and $`\sigma ,\sigma ^{}`$ are the spin quantum numbers. The states $`|i`$ and the corresponding wave functions $`\phi _i(𝐫)=𝐫|i`$ form a basis of the one-particle Hilbert space. The matrix elements in Eq. 2 are defined by:

$`t_{ij}`$ $`=`$ $`i|{\displaystyle \frac{𝐩^2}{2m}}+V(𝐫)|j`$
$`=`$ $`{\displaystyle d^3r\phi _i^{}(𝐫)\left(\frac{\mathrm{}^2}{2m}^2+V(𝐫)\right)\phi _j(𝐫)}`$
$`W_{ij,kl}`$ $`=`$ $`i|j|{\displaystyle \frac{e^2}{|𝐫𝐫^{}|}}|k|l`$
$`=`$ $`{\displaystyle d^3rd^3r^{}\phi _i^{}(𝐫)\phi _j^{}(𝐫^{})\frac{e^2}{|𝐫𝐫^{}|}\phi _k(𝐫^{})\phi _l(𝐫)}`$
Of course these matrix elements depend on the one-particle basis $`\{|i\}`$ that is chosen. But because of the completeness relation the physical results should, in principle, not depend on the choice of the one-particle basis. Because of the lattice periodicity an obvious choice for a one-particle basis is a Bloch basis $`\{|n𝐤\}`$; then the orbital one-particle quantum numbers $`n,𝐤`$ are the band index $`n`$ and the wave number $`𝐤`$ (within the first Brillouin zone). The Bloch wave functions can be written as a product of a plane wave with a lattice periodic Bloch factor

$`\psi _{n𝐤}(𝐫)`$ $`=`$ $`𝐫|n𝐤=e^{i\mathrm{𝐤𝐫}}u_{n𝐤}(𝐫)`$ (4)
$`u_{n𝐤}(𝐫)`$ $`=`$ $`u_{n𝐤}(𝐫+𝐑)`$
where $`𝐑`$ denotes a lattice vector. In practice one can work only on a truncated, finite-diemnsional one-particle Hilbert space. Here the truncation consists of including only a finite number of bands and a set of $`𝐤`$-values from a discrete mesh in $`𝐤`$-space. A Bloch basis can be obtained by applying a traditional band-structure method to solve the Schrödinger equation for an effective one-particle Hamiltonian with a periodic potential. But, because the Bloch states are delocalized, a very large number of Coulomb matrix elements (depending on four quantum numbers) between all possible $`𝐤`$-states would have to be calculated. Therefore, it seems that a more appropriate basis would be to use well localized wave functions. Although the one-particle part is not diagonal within such a localized basis, it is expected that a short-range tight-binding assumption will hold, i.e., that the on-site and the intersite matrix elements for only a few neighbor shells are sufficient. It is well known from elementary solid state theory that Wannier states provide an alternative orthonormal basis set, which spans the same space as the Bloch states. The Wannier states are related to the Bloch states by the unitary transformations:

$`w_{𝐑n}(𝐫)`$ $`=`$ $`𝐫|𝐑n={\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}e^{i\mathrm{𝐤𝐑}}\psi _{n𝐤}(𝐫)`$ (5)
$`|\psi _{n𝐤}`$ $`=`$ $`{\displaystyle \underset{𝐑}{}}e^{i\mathrm{𝐤𝐑}}|𝐑n`$
Now our strategy is the following:

* Perform a traditional band-structure calculation for an effective one-particle Hamiltonian $`H_{\text{eff}}`$ with lattice periodicity to obtain a Bloch basis of the Hilbert space. Only a finite number of band indices will be considered and the calculations will be done for a discretized, finite mesh in $`𝐤`$-space, i.e., we will work only on a reduced, truncated Hilbert space.
* Determine well-localized Wannier functions spanning the same (truncated) Hilbert space as the Bloch basis from the canonical transformation (5) described above.

Of course, the important energy bands (and corresponding band indices) are those that determine the electronic properties of the system, i.e., the bands near to the Fermi level. Because the Hilbert space is truncated, we do no longer work with a complete basis set. Hence, it is important to start from Bloch functions obtained from a band-structure calculation for a well chosen effective one-particle Hamiltonian. Possible and suitable choices for the effective (auxiliary) one-particle Hamiltonian will be discussed in the next Section III. Furthermore, the Wannier functions obtained according to (5) are not unique, because the phases of the Bloch functions can be chosen arbitrarily. In fact, for a given set of randomly chosen phase factors the Wannier functions may not be localized at all. On the other hand, one may use this gauge freedom to construct optimally localized Wannier functionsMV97 (23), as will be described in Section IV.

## III Band structure calculation for the auxiliary one-particle problem

Our Bloch state basis should be suitably chosen for the material of interest, and should be obtained from a standard band-structure calculation of an effective one-particle Hamiltonian with a lattice periodic potential:

$$V_{\text{eff}}(𝐫)=V_{\text{eff}}(𝐫+𝐑)$$
(6)
for any lattice vector $`𝐑`$. The corresponding one-particle Schrödinger equation is

$$\left(\frac{𝐩^2}{2m}+V_{\text{eff}}(𝐫)\right)\psi _{n𝐤}(𝐫)=\epsilon _n(𝐤)\psi _{n𝐤}(𝐫).$$
(7)
Different choices for the effective periodic potential are possible. The simplest choice would be the bare one-particle potential $`V(𝐫)`$. However, since the Coulomb interaction is strong, in general, the resulting eigenenergies of this non-interacting Hamiltonian will be much lower than the relevant energies of the interacting system. For example, in the 3d-systems like Ni or Cu the 3d- and 4s-states form bands close to the Fermi energy. However, without any Coulomb repulsion the 3d-states become very strongly bound atomic-like (core) states, which would be pushed well below the Fermi energy, and therefore the corresponding Bloch eigenfunctions are not a good starting point to describe the electronic bands close to the Fermi level. Because the Hilbert space is truncated, it is extremely important to start from a band Hamiltonian $`T+V_{\text{eff}}`$ that gives eigenfunctions as close as possible to those which are expected to form the relevant many-body states of the interacting system. The bare one-particle potential is consequently a bad choice.

Another possible choice for the effective periodic one-particle Hamiltonian is the Hartree Hamiltonian. Then the Bloch basis is obtained by solving the one-particle Schrödinger equation

$$\left(\frac{𝐩^2}{2m}+V(𝐫)+V_\text{H}(𝐫)\right)\psi _{n𝐤}(𝐫)=\epsilon _n(𝐤)\psi _{n𝐤}(𝐫)$$
(8)
where the Hartree potential is given by

$$V_\text{H}(𝐫)=d^3r^{}\frac{e^2\rho (𝐫^{})}{|𝐫𝐫^{}|}$$
(9)
and the particle density

$$\rho (𝐫)=\underset{l𝐤\sigma }{}f(\epsilon _l(𝐤))|\psi _{l𝐤}(𝐫)|^2$$

$$\text{where }f(E)=\frac{1}{e^{\beta (E\mu )}+1}$$

is the Fermi function (step function for the ground state $`T=0`$ or $`\beta =\mathrm{}`$). Since the occupied states $`|\psi _{l𝐤}`$ determine the particle density and Hartree potential, Eq. 8 has to be solved selfconsistently, which is usually done by iteration procedure. The advantage of including the Hartree potential into the effective potential, i.e. $`V_{\text{eff}}(𝐫)=V(𝐫)+V_\text{H}(𝐫)`$, is that the approximate effects of the Coulomb interaction are included (in the mean-field approximation). Therefore, the eigenenergies (energy bands) will be about the right magnitude and the resulting basis functions can be expected to be more suitable in the energy regime around the Fermi level.

Since the only purpose in solving the effective one-particle Schrödinger equation (7) is the construction of a suitable basis set of Bloch functions, we will not make use of the eigenenergies $`\epsilon _n(𝐤)`$ obtained in (7) or give these solutions any physical interpretation. For this reason one could also use a variety of different artificial effective potentials. For example, we have studied a “weighted Hartree potential”

$$V_{\text{eff}}(𝐫)=V(𝐫)+xV_\text{H}(𝐫)$$
(10)
with an additional parameter $`x`$, which can be varied so that the choice of the basis functions minimizes the resulting Hartree- or Hartree-Fock ground state energy. Another possible choice for the effective (auxiliary) one-particle potential would be the LDA-potential. Because, in this context, LDA would only be used to construct basis functions, there would be no problem with double counting interaction terms.

After choosing the appropriate effective one-particle potential one has to employ a numerical method to solve the Schrödinger equation (7) for the periodic potential. For this established and well known band-structure methods can be used, such as the “linearized muffin-tin orbital” (LMTO) methodOKAndersen (24, 25), the (full potential) “linearized augmented plane wave” (LAPW) method (using the WIEN2k computer codeWIEN2K (27)), or the “full-potential local-orbital” (FPLO)Eschrig99 (28) minimum-basis method.

## IV Determination of maximally localized Wannier functions

The band-structure calculation diagonalizes the auxiliary one-particle Hamiltonian to determine the Bloch states $`\psi _{n𝐤}(𝐫)=e^{i\mathrm{𝐤𝐫}}u_{n𝐤}(𝐫)`$ labeled by band indices $`\{n\}`$. After truncating this one-particle Hilbert space by considering only a finite number $`J`$ of band indices, unitary transformations are possible for each fixed $`𝐤`$ that lead to a new basis set

$$\stackrel{~}{\psi }_{n𝐤}(𝐫)=\underset{m}{}U_{mn}^𝐤\psi _{m𝐤}(𝐫),$$
(11)
which are still Bloch functions but with other band indices (within which the auxiliary one-particle Hamiltonian is, in general, no longer diagonal). One can use this gauge freedom to construct a Bloch basis (still of the same restricted one-particle Hilbert space) for which the corresponding Wannier functions determined according to (5) are maximally localized.

A suitable measure for the localization of the Wannier function is the spread functional

$$\mathrm{\Omega }=\underset{n}{}\left[r^2_n𝐫_n^2\right],$$
(12)
where the notation $`A_n=\mathrm{𝟎}n|A|\mathrm{𝟎}n`$ for any operator $`A`$ has been used. Among all the possible gauges, we will use that gauge which minimizes this spread functional $`\mathrm{\Omega }`$.

A method for minimizing Eq. (12) has been developed by Marzari and VanderbiltMV97 (23). This method has already found widespread applications recentlyMV2002 (29, 30). It starts from a decomposition into invariant, diagonal, and off-diagonal contributions:

$`\mathrm{\Omega }`$ $`=`$ $`\mathrm{\Omega }_\text{I}+\stackrel{~}{\mathrm{\Omega }}`$
$`=`$ $`{\displaystyle \underset{n}{}}\left[r^2_n{\displaystyle \underset{𝐑m}{}}\left|𝐑m|𝐫|\mathrm{𝟎}n\right|^2\right]`$
$`+`$ $`{\displaystyle \underset{n}{}}{\displaystyle \underset{𝐑m\mathrm{𝟎}n}{}}\left|𝐑m|𝐫|\mathrm{𝟎}n\right|^2`$
The first term is “gauge-invariant”, i.e., independent of the choice of unitary transformations among the bands. The second term can be decomposed into band-off-diagonal and band-diagonal components:

$$\stackrel{~}{\mathrm{\Omega }}=\underset{\mathrm{\Omega }_{\text{OD}}}{\underset{}{\underset{mn}{}\underset{𝐑}{}\left|𝐑m|𝐫|\mathrm{𝟎}n\right|^2}}+\underset{\mathrm{\Omega }_\text{D}}{\underset{}{\underset{n}{}\underset{𝐑\mathrm{𝟎}}{}\left|𝐑n|𝐫|\mathrm{𝟎}n\right|^2}}$$
(14)
In practice one needs the unitary matrices $`U_{mn}^𝐤`$ on a discretized $`𝐤`$-mesh, which can be chosen to be simple cubic. If $`𝐛`$ denotes the vectors connecting each $`𝐤`$-point to its nearest neighbors, one can define

$$M_{mn}^{𝐤,𝐛}=u_{m𝐤}|u_{n,𝐤+𝐛},$$
(15)
where $`u_{n𝐤}`$ are the Bloch factors appearing in (4). Using (15) one can express the expectation values of $`𝐫`$ and $`r^2`$ asMV97 (23):

$`𝐫_n`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤,𝐛}{}}w_b𝐛\text{Im}\mathrm{ln}M_{nn}^{𝐤,𝐛}w_b={\displaystyle \frac{1}{2b^2}}`$ (16)
$`r^2_n`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤,𝐛}{}}w_b\left\{[1|M_{nn}^{𝐤,𝐛}|^2]+[\text{Im}\mathrm{ln}M_{nn}^{𝐤,𝐛}]^2\right\}`$
These expressions are not unique. This non-uniqueness arises from the fact that it is only required that these expectation values become exact in the limit of a dense mesh. The present choice guarantees that under the gauge transformation $`|u_{n𝐤}e^{i\mathrm{𝐤𝐑}}|u_{n𝐤}`$ (corresponding to a translation by a lattice vector) the following physical properties are fulfilled:

$`𝐫_n`$ $``$ $`𝐫_n+𝐑`$ (17)
$`r^2_n`$ $``$ $`r^2_n+2𝐫_n𝐑+R^2,`$
Then one finds for $`\mathrm{\Omega }=\mathrm{\Omega }_\text{I}+\mathrm{\Omega }_{\text{OD}}+\mathrm{\Omega }_\text{D}`$:

$`\mathrm{\Omega }_\text{I}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤,𝐛}{}}w_b\left(J{\displaystyle \underset{mn}{}}|M_{mn}^{𝐤,𝐛}|^2\right)`$ (18)
$`\mathrm{\Omega }_{\text{OD}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤,𝐛}{}}w_b{\displaystyle \underset{mn}{}}|M_{mn}^{𝐤,𝐛}|^2`$
$`\mathrm{\Omega }_\text{D}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤,𝐛}{}}w_b{\displaystyle \underset{n}{}}\left(\text{Im}\mathrm{ln}M_{nn}^{𝐤,𝐛}𝐛𝐫_n\right)^2`$
The Marzari-Vanderbilt algorithmMV97 (23) now consists in an iterated application of a (in principle infinitesimal, in practice numerically discrete) canonical transformation chosen so that the gradient of $`\mathrm{\Omega }`$ is negative, i.e., the spread functional and thus the delocalization of the Wannier function decreases. In each iteration cycle $`\mathrm{\Omega }_\text{I}`$ should not change whereas $`\mathrm{\Omega }_{\text{OD}}+\mathrm{\Omega }_\text{D}`$ should decrease.

In practice the iteration cycle can schematically be written down asMV97 (23):

$$U^{(N)𝐤}=U^{(N1)𝐤}\mathrm{exp}\left[\mathrm{\Delta }W^𝐤(M^{(N1)𝐤,𝐛})\right]$$
(19)
$$M^{(N)𝐤,𝐛}=U^{(N)𝐤}M^{(N1)𝐤,𝐛}U^{(N)𝐤+𝐛}$$
(20)
with the initialization:

$$U_{mn}^{(0)𝐤}=\delta _{mn}$$
(21)
$$M_{mn}^{(0)𝐤,𝐛}=u_{m𝐤}|u_{n,𝐤+𝐛}$$
(22)
Here $`\mathrm{\Delta }W^𝐤(M^{𝐤,𝐛})`$ is an anti-hermitian matrix (so that $`\mathrm{exp}\left[\mathrm{\Delta }W^𝐤\right]`$ is unitary) which according to Ref. MV97, 23 is explicitly given by the following equations:

$$\mathrm{\Delta }W^𝐤=\frac{\alpha }{w}\underset{𝐛}{}w_b\left(𝒜(R^{\mathrm{𝐤𝐛}})𝒮(T^{\mathrm{𝐤𝐛}})\right)$$
(23)
where $`𝒜(B)=(BB^{})/2,𝒮(B)=(B+B^{})/2i`$ are the symmetrizing and antisymmetrizing operations for the operators or matrices $`B`$, $`w=_𝐛w_b`$, $`\alpha [0,1]`$ is a numerical parameter determining the discrete, finite step size, and the $`J\times J`$-matrices $`R^{\mathrm{𝐤𝐛}}`$ and $`T^{\mathrm{𝐤𝐛}}`$ are explicitly given by:

$`R_{mn}^{\mathrm{𝐤𝐛}}`$ $`=`$ $`M_{mn}^{\mathrm{𝐤𝐛}}M_{nn}^{\mathrm{𝐤𝐛}}`$ (24)
$`T_{mn}^{\mathrm{𝐤𝐛}}`$ $`=`$ $`{\displaystyle \frac{M_{mn}^{\mathrm{𝐤𝐛}}}{M_{nn}^{\mathrm{𝐤𝐛}}}}\left(\text{Im}\mathrm{ln}M_{nn}^{\mathrm{𝐤𝐛}}+𝐛𝐫_n\right)`$
The whole algorithm (19) can be considered as an iteration scheme to construct the unitary transformation $`U^𝐤`$. Once this iteration has converged one can perform the transformation to the new Bloch functions according to (11) and then calculate the optimally localized Wannier functions according to the definition (5).

In practice it is useful to prepare the Bloch orbitals to make the starting Wannier functions somewhat localized. This has two advantages: (i) the minimization procedure converges faster, and (ii) this helps to avoid getting trapped in local minima. For that purpose we make a rather trivial “gauge transformation” to each Bloch function obtained by the band structure calculation, namely a multiplication with a phase factor according to

$$\psi _{n𝐤}(𝐫)\mathrm{exp}\left[i\text{Im}\mathrm{ln}\psi _{n𝐤}(𝐫_n)\right]\psi _{n𝐤}(𝐫).$$
(25)
This gauge transformation has the property that $`\text{Im}\mathrm{ln}\mathrm{\Psi }_{n𝐤}(𝐫_n)`$ transforms to zero. So at the point $`𝐫_n`$, all the Bloch functions will have the same phase and $`𝐫_n|\mathrm{𝟎}n`$ will take a large value. To make the method work well, one should choose $`𝐫_n`$ where the Wannier functions are expected to be reasonably large and one can choose $`𝐫_n`$ individually for each band $`n`$.

## V One-particle and Coulomb matrix elements

After the maximally localized Wannier functions $`w_{𝐑n}(𝐫)`$ have been determined, the next task is to calculate the one-particle and Coulomb matrix elements of the Hamiltonian. The one-particle matrix elements in Wannier representation are explicitly given by:

$$t_{12}=d^3rw_1^{}(𝐫)\left(\frac{\mathrm{}^2}{2m}^2+V(𝐫)\right)w_2(𝐫)$$
(26)
where $`V(𝐫)`$ is the lattice periodic one-particle potential for the non-interacting electrons. Here, and in the following, we use the abbreviated notation 1 to mean $`𝐑_1n_1`$ and 2 to mean for $`𝐑_2n_2`$, etc. Note that the one-particle part of the Hamiltonian cannot be expected to be band diagonal because of the unitary transformation to another set of the band indices within which the Wannier functions are maximally localized. But by construction and definition the Wannier functions obey the following translational invariance property:

$$w_{𝐑n}(𝐫)=w_{\mathrm{𝟎}n}(𝐫𝐑)$$
(27)
where $`\mathrm{𝟎}`$ is the arbitrarily chosen zero lattice vector. Thus, the one-particle matrix elements depend only on the relative distance $`𝐑=𝐑_1𝐑_2`$:

$$t_{𝐑n,m}=𝐑n|\left(\frac{\mathrm{}^2}{2m}^2+V(𝐫)\right)|\mathrm{𝟎}m$$
(28)
Because of the strong localization of the Wannier functions one can safely expect that these one-particle (hopping) matrix elements decrease with increasing distance $`|𝐑|`$ and have to be evaluated explicitly only for a limited number of $`𝐑`$, in particular on-site, i.e. for $`\mathrm{𝟎}`$, and for a few neighbor shells $`𝐑`$. In practice it turned out that for properly localized Wannier functions (optimized according to the Marzari-Vanderbilt algorithm) the explicit evaluation of one-particle matrix elements up to the 5th neighbor shell is sufficient, whereas for non-optimally-localized Wannier functions (made only “somewhat localized” according to the prescription (25), for instance) up to 30 neighbors have to be taken into account.Schnell2002 (31) Depending on the specific band-structure method used for the effective (auxiliary) one-particle problem, the eigenfunctions (Bloch functions) are usually represented as linear combinations of certain elementary functions, for instance plane waves, Gaussians, or spherical harmonics. Thus, also the evaluation of the three-dimensional integral in (26) can usually be traced back to known integrals over these elementary functions, and usually at most a one-dimensional integral remains to be calculated explicitly numerically.

The Coulomb matrix elements are given by:

$$W_{12,34}=d^3𝐫d^3𝐫^{}w_1^{}(𝐫)w_2^{}(𝐫^{})\frac{e^2}{|𝐫𝐫^{}|}w_3(𝐫^{})w_4(𝐫)$$
(29)
Let us first take a brief look at general properties of the matrix elements which are useful for minimizing computing time and memory storage. From (29) and (27), it follows that

$$W_{12,34}=W_{(𝐑_1𝐑_4n_1)(𝐑_2𝐑_4n_2),(𝐑_3𝐑_4n_3)(\mathrm{𝟎}n_4)}$$
(30)
That is, we may always translate the lattice site indices in a way that $`𝐑_4\mathrm{𝟎}`$. Moreover, since $`𝐫`$ and $`𝐫^{}`$ in (29) can be interchanged, we have $`W_{12,34}=W_{21,43}`$ and since $`W`$ is Hermitian we have $`W_{12,34}=W_{43,21}^{}`$. For the practical evaluation of these 6-fold integrals a fast Fourier transformation (FFT) algorithm turned out to be very efficient and independent of the specific representation of the Wannier functions. Using

$$d^3𝐪\frac{e^{i\mathrm{𝐪𝐫}}}{q^2}=\frac{2\pi ^2}{|𝐫|}$$
(31)
one obtains:

$$W_{12,34}=\frac{e^2}{2\pi ^2}d^3𝐪\frac{1}{q^2}f_1(𝐪)f_2(𝐪)$$
(32)
where

$$f_1(𝐪)=d^3𝐫e^{i\mathrm{𝐪𝐫}}w_{𝐑_1n_1}^{}(𝐫)w_{𝐑_4n_4}(𝐫)$$
(33)
$$f_2(𝐪)=d^3𝐫e^{i\mathrm{𝐪𝐫}}w_{𝐑_2n_2}^{}(𝐫)w_{𝐑_3n_3}(𝐫)$$
(34)
These functions are just the Fourier transforms of a product of Wannier functions. They can be calculated very efficiently by evaluating the Wannier functions on a cubic mesh in $`𝐫`$-space with some finite spacing $`\mathrm{\Delta }x`$ and then applying a standard (numerical) FFT-algorithm. The result of this Fourier transformation is $`f_{1,2}(𝐪)`$ on a cubic mesh in $`𝐪`$-space with some $`\mathrm{\Delta }q`$. As $`F(𝐪)=f_1(𝐪)f_2(𝐪)`$ is a smooth function at $`𝐪=\mathrm{𝟎}`$, a divergence arises in the integrand from $`q^2`$. In order to treat this divergence, we split the integral by subtracting and adding $`F(\mathrm{𝟎})`$:

$$d^3𝐪\frac{F(𝐪)}{q^2}=d^3𝐪\frac{F(𝐪)F(\mathrm{𝟎})}{q^2}+F(\mathrm{𝟎})d^3𝐪\frac{1}{q^2}$$
(35)
The first non-vanishing term of a polynomial expansion of the numerator $`F(𝐪)F(\mathrm{𝟎})`$ is of order $`q^2`$. Hence, the polynomial expansion of the integrand starts with a constant and the divergence is avoided.

All integrals are over a cube with length $`2p=N\mathrm{\Delta }q`$. The first integral in (35) is evaluated by transforming the integral into a sum over little cubes with volume $`(\mathrm{\Delta }q)^3`$. At $`𝐪=\mathrm{𝟎}`$, the value of the integrand is calculated via the second derivative of $`F(𝐪)`$ at $`𝐪=\mathrm{𝟎}`$ numerically.

The remaining integral in (35) is simply a constant given by:

$`{\displaystyle _p^{+p}}𝑑q_x{\displaystyle _p^{+p}}𝑑q_y{\displaystyle _p^{+p}}𝑑q_z{\displaystyle \frac{1}{q^2}}`$ $`=`$
$`p{\displaystyle _1^{+1}}𝑑q_x{\displaystyle _1^{+1}}𝑑q_y{\displaystyle _1^{+1}}𝑑q_z{\displaystyle \frac{1}{q^2}}`$ $`=`$ $`pC`$ (36)
with $`C=15.34825`$, which we have evaluated numerically.

As argued already for the one-particle matrix elements, the Coulomb matrix elements need explicitly be evaluated only on-site, i.e. for $`𝐑_1=𝐑_2=𝐑_3=\mathrm{𝟎}`$, and for at most only a few neighbor shells, probably only nearest neighbors, because of the good localization of the Wannier functions used. Here again the application of the Marzari-Vanderbilt algorithm turns out to be important in order to reduce the number of explicitly necessary evaluations of Coulomb matrix elements, which are (computational) time consuming.

## VI Many-body treatment of second quantized Hamiltonian with ab-initio parameters

After the evaluation of the one-particle and two-particle matrix elements we have the Hamiltonian in second quantization, i.e., in the form (2) in a Wannier representation with all parameters determined from first principles. The only approximations made so far are the restriction to a finite number $`J`$ of band indices for a specific basis set and the approximations inherent to the band-structure calculation method used (for example, muffin-tin potentials, linearization approximations, discretization in $`𝐤`$\- and real $`𝐫`$-space, etc.). The Hamiltonian is in the form of a multi-($`J`$-)band extended Hubbard model. A one-band description is usually not sufficient even for the simplest solids. Furthermore, the hopping matrix elements are in general not restricted to nearest neighbors (but taken into account up to the 5th nearest neighbors) and direct and exchange Coulomb matrix elements on-site, inter-band and intra-band and, if necessary, also inter-site (but again only up to a few neighbor shells) are taken into account.

Now, in principle, one can apply any many-body method that has been developed for interacting electron systems on a lattice. Many of the standard methods rely on systematic perturbation theory with respect to the Coulomb interaction, and the terms can be represented by Feynman diagramsMahan (32). The lowest-order selfconsistent approximation within this diagrammatic approach is the Hartree-Fock approximation (HFA), which can also be derived from a variational principle. Within the HFA the selfenergy is explicitly given by:

$`\mathrm{\Sigma }_{12,\sigma }^{\text{HF}}`$ $`=`$ $`\mathrm{\Sigma }_{12}^{\text{Hart}}+\mathrm{\Sigma }_{12,\sigma }^{\text{Fock}}`$
$`=`$ $`{\displaystyle \underset{34\sigma ^{}}{}}\left[W_{13,42}\delta _{\sigma \sigma ^{}}W_{31,42}\right]c_{3\sigma ^{}}^{}c_{4\sigma ^{}}`$
Here the expectation values (or the density matrix)

$$A_{12}^\sigma =c_{1\sigma }^{}c_{2\sigma }$$
(38)
have to be determined selfconsistently for the Hartree-Fock Hamiltonian:

$$H_{\text{HF}}=\underset{12\sigma }{}\left(t_{12}+\mathrm{\Sigma }_{12,\sigma }^{\text{HF}}\right)c_{1\sigma }^{}c_{2\sigma }$$
(39)
These expectation values are easily calculated by going back to the Bloch states within which the effective one-particle Hamiltonian (39) is diagonal with eigenenergies $`E_m^\sigma (𝐤)`$. In this Bloch representation the expectation values for zero temperature (i.e., within the ground state) are given by

$$c_{n𝐤\sigma }^{}c_{n^{}𝐤^{}\sigma ^{}}=\delta _{nn^{}}\delta _{\mathrm{𝐤𝐤}^{}}\delta _{\sigma \sigma ^{}}\theta (E_FE_n^\sigma (𝐤)).$$
(40)
One finds

$$A_{12}^\sigma =\frac{1}{N}\underset{𝐤}{}e^{i𝐤(𝐑_1𝐑_2)}\underset{m}{}U_{mn_1}^𝐤U_{mn_2}^𝐤\theta (E_FE_m^\sigma (𝐤))$$
(41)
where (similar as in (11)) the $`U_{mn}^𝐤`$ denote the unitary transformation between the basis with band indices leading to maximally localized Wannier functions and the basis within which the effective one-particle Hamiltonian (39) is diagonal. For the total energy in the Hartree-Fock approximation one finds:

$$E_{\text{HF}}=H=\underset{12\sigma }{}\left(t_{12}+\frac{1}{2}\mathrm{\Sigma }_{12}^{\text{Hart}}+\frac{1}{2}\mathrm{\Sigma }_{12,\sigma }^{\text{Fock}}\right)A_{12}^\sigma $$
(42)
This follows formally from the original many-body Hamiltonian (2) using the decoupling

$$c_{1\sigma }^{}c_{2\sigma ^{}}^{}c_{3\sigma ^{}}c_{4\sigma }c_{1\sigma }^{}c_{4\sigma }c_{2\sigma ^{}}^{}c_{3\sigma ^{}}c_{1\sigma }^{}c_{3\sigma ^{}}c_{2\sigma ^{}}^{}c_{4\sigma }$$
(43)
which is valid according to Wick’s theorem, if the expectation values are calculated with respect to an effective one-particle Hamiltonian, which the Hartree-Fock-Hamiltonian (39) is.

Of course the HFA is only the simplest and basic many-body approximation which can be applied (mean-field). It is correct only up to linear order in the Coulomb interaction. The main advantage of the diagrammatic approach to the second-quantized many-body Hamiltonian (2) is that one can formulate systematic improvements and corrections to the basic approximations. One possibility to go beyond HFA is to take into account all skeleton diagrams up to second order in the Coulomb interactions. Probably better and more reliable many-body approximations are obtained by resummations of an infinite series of diagrams of a certain class. For instance the resummation of all bubble diagrams corresponds essentially to the random phase approximation (RPA) and means that in the Fock diagram the bare interaction has to be replaced by an effective screened interaction.

Of course the application of many-body theory is not restricted to the standard perturbational methods in terms of Feynman diagrams. During the last decade much progress in many-body theory has been made due to the development and successful application of non-perturbative methods, such as the DMFTGKKR96 (21). In its original version the DMFT requires a local ($`𝐤`$-independent) selfenergy, which is usually justified for 3-dimensional systems. But even if this assumption should not be justified, the application of non-local (cluster) extensions of DMFT is possible. In any case, for realistic materials a multi-band system has to be studied; therefore a mapping on a multi-level single-impurity problem and new types of DMFT-selfconsistency cycles are probably required. Another possible non-perturbative many-body approach to our Hamiltonian (2) with ab-initio parameters is the application of variational methods, such as the generalized Gutzwiller ansatz, which has recently been applied to the 3d-ferromagnet nickelGebhardWeber (33).

## VII Comparison with other ab-initio methods

In this section we compare our ab-initio approach with existing and established first-principles methods for electronic-structure calculations.

DFT-based approaches encounter the problem that the density functional of the exchange-correlation energy is not known, and hence additional assumptions and ansatz are necessary. Within LDA the assumption is made that the particle density to be determined is the same as that of an effective one-particle system, and the ansatz for the (local) density dependence of the exchange-correlation energy is chosen so that in the case of the homogeneous electron gas the most accurately known results are reproduced. So results of many-body theory for the homogeneous electron gas enter, but many-body theory is not directly applied for the lattice system (an inhomogeneous electron gas). In LDA the (local) Hartree potential for the inhomogeneous electron system is explicitly taken into account when solving the Kohn-Sham equations, but the exchange (Fock) contribution is only approximately considered (in a form which is correct only for the homogeneous electron system). Therefore, whereas self-interaction terms exactly cancel, if the Hartree and Fock contributions are treated on the same level, they usually do not cancel in LDA. Our approach, on the other hand, is free from the problem of self-interaction, as can immediately be seen from (43), because the Hartree- and exchange terms are treated on the same level. Furthermore, many-body theory is directly applied for the lattice electron system instead of only for the homogeneous electron gas. From the very beginning we never assume a constant electron density or a dependence on the local density. The effects which the “gradient correction” and nonlocal density schemes corrections to LDA aim at are automatically considered in our approach.

The standard Hartree-Fock approximation in first quantization requires the solution of the equations

$`\left({\displaystyle \frac{\mathrm{}^2}{2m}}^2+V(𝐫)\right)\phi _i(𝐫)+{\displaystyle \underset{j\sigma ^{}}{}}{\displaystyle d^3r^{}\frac{e^2|\phi _j(𝐫^{})|^2}{|𝐫𝐫^{}|}\phi _i(𝐫)}`$
$`{\displaystyle \underset{j}{}}{\displaystyle d^3r^{}\frac{e^2\phi _j^{}(𝐫^{})\phi _i(𝐫^{})}{|𝐫𝐫^{}|}\phi _j(𝐫)}=\epsilon _i\phi _i(𝐫).`$ (44)
For equal spin and $`i=j`$ the second (Hartree) and the third (Fock or exchange) contributions exactly cancel each other so that no self-interaction occurs. The aim of the recent self-interaction corrected (exact exchange) schemesVogl97 (8) as LDA improvements is, therefore, to consider these exchange terms rigorously. But the effective exchange (Fock) potential is a non-local potential, and the one-particle wave-function $`\phi _i`$ to be determined enters not only at the position $`𝐫`$ but also at all other positions $`𝐫^{}`$. For that reason HFA treatments in first quantization and exact exchange schemes are computationally more complicated than LDA methods, especially for a periodic solid. Within our approach HFA-calculations are as simple as a Hartree calculation, since only expectation values have to be evaluated. But in our method self-interaction terms also exactly cancel. Therefore, our approach should be as good as the exact-exchange approach but easier in practical applications.

When using an RPA (screened HFA) as the many-body approximation, our approach should be similar to the GWA. But in GWA one usually uses as the Green function (G-line) the LDA (effective one-particle) result whereas we suggest using a selfconsistently calculated Green function (within the many-body approximation used). Also, when using the LDA Green function, only a part of the exchange is taken into account unlike the exact exchange of our approximation.

All the recent attempts to combine ab-initio and many-body methods (on different levels concerning the many-body treatment), namely, the LDA+U, the LDA++, and the LDA+DMFT-schemes, have in common that they start from an LDA ab-initio calculation and obtain their one-particle bands and density of states from this LDA treatment, and then add an interaction (correlation) term with an on-site Hubbard-U term. Then, on one hand, some of the correlations and interactions are already included in the effective one-particle energies, namely those on the LDA (homogeneous electron gas) level. On the other hand, since other correlations are treated in many-body theory for the Hubbard U, there may be a double counting of interaction (correlation) contributions, which is hard to justify and the magnitude of which is difficult to estimate. Furthermore, in these theories the Hubbard U usually is an additional free parameter, so that these are not really “ab-initio” (first-principles) approaches. In the our approach we do not start from the LDA but from an effective band structure calculation on the Hartree level. Therefore no correlation terms are implicitly included within the one-particle band structures and the problem of double counting of interaction terms does not occur. Furthermore, all Coulomb matrix elements are calculated from “first principles”.

## VIII Summary and Outlook

In this paper we have suggested a new approach for combining ab-initio and many-body methods for the calculation of the electronic properties of solids. The starting point is a traditional band-structure calculation for an effective (auxiliary) one-particle Hamiltonian, which can be the Hartree-Hamiltonian. This yields, in particular, the eigenfunctions in the form of Bloch functions. Keeping only a finite number of $`J`$ band indices restricts and truncates the Hilbert space for further calculations. We use the Marzari-Vanderbilt algorithm to construct maximally localized Wannier functions (within the truncated one-particle Hilbert space). Then all the one-particle (tight-binding) and two-particle (Coulomb) matrix elements between these Wannier functions can be calculated. The strong localization guarantees that only on-site matrix elements and near-neighbor inter-site matrix elements have to be calculated. We are left with a many-body Hamiltonian in second quantization but with parameters determined from first principles for any given material. This should be solved by a suitable many-body approximation. The simplest approximation is the (unscreened) HFA, but improved methods (e.g., RPA or non-perturbative many-body approximations like DMFT) are also feasible. Our suggested approach is free from the problems of double counting of correlation effects and self-interaction and considers exchange contributions exactly. It does not rely on assumptions based on the homogeneous electron gas or a dependence on the local electron density density. An inhomogeneous (lattice) electron system is considered right from the beginning.

An application of this scheme to the 3d transition metals Fe, Co, Ni and Cu using the unscreened HFA as the many-body method is presented in the following paper.Schnell2003-2 (34)

###### Acknowledgements.

This work has been supported by a grant from the Deutsche Forschungsgemeinschaft No. Cz/31-12-1. It was also partially supported by the Department of Energy under contract W-7405-ENG-36. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.


format: text

text:
# Detection of type-2 quasars in the radio galaxies B3 0731+438 and 3C 257

## 1 Introduction

An important issue in AGN research is whether heavily obscured powerful AGNs (type-2 quasars) are common. The unification of AGNs (Antonucci 1993), UV spectral shapes (Maiolino et al. 2001) and the modeling of AGNs infrared spectral energy distributions (Pier & Krolik 1993) all suggest the presence of numerous type-2 quasars in the Universe. These obscured AGN are thought to outnumber unobscured objects by a factor of $``$ 3:2 at redshifts out to z $``$ 3.5 (Willott et al. 2000; Fiore et al. 2001) and it is important to understand whether this ratio applies at higher redshifts and luminosities (Barcons et al. 1998; Almaini, Lawerence & Boyle 1999; Stern et al. 2002). Obscured AGN are also believed to be responsible for a significant fraction of the hard X-ray background (Comastri et al. 1995). Almost all of the soft X-ray background ($`<`$ 2 keV) has been resolved into discrete sources, the majority of which are broad line (BL) type-1 QSO’s. At higher energies, typical unobscured QSO spectra are too steep to explain the XRB and account for only $``$ 20 $`\%`$ of the sources (Hasinger et al. 1998). Obscured AGNs provide a natural explanation, since photoelectric absorption allows only the hard X-rays to penetrate. This implies the majority of the energy density generated by accretion in the Universe takes place in obscured AGN. These results have been confirmed by data from two deep Chandra observations (Brandt et al. 2001; Rosati et al. 2002) and recent analysis from an XMM-Newton observation of the Lockman Hole has also revealed similar results; Mainieri et al. (2002) found $``$ 27 $`\%`$ of the Lockman Hole sample are Extremely Red Objects, the majority of which contain an obscured AGN.

Recently a small number of type-2 quasars candidates have been identified, but these results are usually based on limited data (Fabian et al. 2000; Willott et al. 2001; Fabian et al. 2001). The best published type-2 quasars, all obtained using the Chandra X-ray Observatory, are IRAS 09104+4109 at z $`=`$ 0.44 (Iwasawa et al. 2001), originally observed with BeppoSAX (Franceschini et al. 2000), CDF-S 202 at z $`=`$ 3.7 (Norman et al. 2001), A18 at z $`=`$ 1.467 (Crawford et al. 2001), CX052, z $`=`$ 3.288 (Stern et al. 2002) and B2 0902+343, z $`=`$ 3.395 (Fabian, Crawford & Iwasawa 2002). Correcting for absorption, these objects have X-ray luminosities of the order 1045 erg s-1 and intrinsic absorption NH $``$ 1024 cm-2.

To identify likely type-2 quasars which are bright enough to provide reasonable X-ray spectra, we have observed two Narrow line radio galaxies (NLRGs) with XMM-Newton. The correlation between the mass of the central SMBHs and the galaxy spheroidal component (Magorrian et al. 1998) indicates that SMBHs are most likely to reside in massive spheroids. Radio galaxies almost always host (forming) massive spheroids (van Bruegel et al. 1998), and NLRGs lack broad components to their optical emission lines, so these galaxies are promising candidates for hosting obscured SMBHs (Barthel 1989). To identify such objects several parameters need to be measured; in particular the absorbing column of the obscuring material and the intrinsic X-ray luminosity. In the context of this paper Type-2 quasars are defined to have an intrinsic, hard X-ray (2 – 10 keV) luminosity $`>`$ 1044 erg s-1 (we assume cosmological values Ho $`=`$ 75 km s-1 Mpc-1 and qo $`=`$ 0.5 throughout); Less luminous objects are called Seyfert galaxies. The X-ray absorbing (gas and dust) column density of Seyfert 2s has NH $``$ 1022 \- 1025 cm-2 (Risaliti, Maiolino & Salvati 1999), and type-2 quasars are expected to have similar values.

We present XMM-Newton observations and spectral fitting of two powerful FRII radio galaxies, B3 0731+438 and 3C 257. Both sources are high redshift NLRGs with predicted X- ray luminosities well within the range of that expected for type-2 quasars. 3C 257 is the highest redshift (z $`=`$ 2.474) radio galaxy observed in the 3C catalog with a fairly smooth K band emission and morphology resembling an elliptical galaxy (van Breugel 1998). B3 0731+438 (z $`=`$ 2.429) is a typical FRII radio galaxy with double-lobed, radio emitting hot spots and a central core (Carilli et al. 1997; Motohara et al. 2000).

This paper is organized as follows; in Section 2 we discuss sample selection and the use of optical emission line luminosities to determine the predicted X-ray luminosity of both sources. This is followed in Sections 3 and 4 by the XMM-Newton observations and data analysis. In Section 5 we discuss the spectral fitting and in the final section we conclude our findings.

## 2 Target Selection

There is an observed correlation between L\[OIII\] ($`\lambda `$5007) and the X-ray luminosity LX (2 – 10 keV) in radio loud AGNs (Mulchaey et al. 1994; Sambruna, Eracleous & Mushotzky 1999). This relation indicates NLRGs with L\[OIII\] $`>`$$`\times `$ 1043 erg s-1 will have an X-ray luminosity LX (2 – 10 keV) $`>`$ 5 $`\times `$ 1044 erg s-1 at the 90 $`\%`$ probability level. The relationship was used to select target NLRGs which are likely to host luminous AGNs and the two most luminous in \[OIII\] were chosen for observation. From the \[OIII\] luminosities (Evans 1998) the predicted X-ray luminosities for B3 0731+438 (L\[OIII\] $`=`$ 12.4 $`\times `$ 1043 erg s-1) and 3C 257 (L\[OIII\] $`=`$ 14.0 $`\times `$ 1043 erg s-1) were estimated to be greater than 12 and 14 $`\times `$ 1044 erg s-1 respectively at the 90 $`\%`$ probability level. The width of the H$`\alpha `$ line in the K-band spectra was determined for both sources (Evans 1998; Jackson and Rawlings 1997); B3 0731+438 shows no broad H$`\alpha `$ component (FWHM$`_{\mathrm{H}\alpha +[\mathrm{NII}]}`$800 km s-1) and 3C 257 has a somewhat broader line (FWHM$`_{\mathrm{H}\alpha +[\mathrm{NII}]}`$1800 km s-1). Thus, for both objects their spectra indicate prominent emission from a narrow-line region while any nuclear emission is highly obscured.

For B3 0731+438 the diagnostic emission-line ratios are consistent with the presence of a Seyfert-2 nucleus (Evans 1998). A narrow-band optical image shows a large Ly$`\alpha `$ emission nebula (McCarthy 1991), suggestive of a hidden nucleus ionising the extended gas. This possibility is supported by a narrow-band 2.25 $`\mu `$m image obtained using the Subaru telescope showing biconical clouds extending out to 37 kpc, radiating H$`\alpha `$+\[NII\] emission lines with a total H$`\alpha `$+\[NII\] flux of 3.5 $`\times `$ 10-15 erg s-1 cm-2 (Motohara et al. 2000), corresponding to a luminosity of 6.7 $`\times `$ 1043 erg s-1, (correcting to a cosmological scale of Ho $`=`$ 75 km s-1 Mpc-1 and qo $`=`$ 0.5). The H$`\alpha `$+\[NII\] ionization cone, aligned with the axis of the radio hot spots, implies the existence of a hidden, more luminous central source not observed in the optical, consistent with Keck spectropolarimetry data which suggest the nucleus is extinguished and we are only seeing a young stellar population in direct optical light (Vernet et al. 2001). The combination of these data suggest B3 0731+438 is a powerful, radio-loud AGN with a large intrinsic absorption. Little is known about 3C 257.

## 3 Observations and Data Reduction

B3 0731+438 and 3C 257 were observed during revolutions 343 (22nd October 2001) and 370 (15th December 2001) of XMM-Newton respectively. The observations of the EPIC MOS and PN detectors for both sources were in full window mode using a medium filter. The data were screened with the XMM Science Analysis Software (SAS) and pipeline processed. X-ray events corresponding to patterns 0 - 12 for the two MOS cameras and patterns 0 - 4, corresponding to single and double events for the PN were selected.

We then extracted the spectra and lightcurves for both the source and background in each observation, separately for each EPIC detector. Circular source regions of 1 arcmin diameter for B3 0731+438 and 0.6 arcmin for 3C 257 were selected. The background regions were placed at positions close to the sources and any significant flaring in the source and background regions were then removed using time selections in XMMSELECT. We then coadded the MOS1 and MOS2 data into a single spectral file to maximize the signal-to-noise ratio and grouped the data into a minimum of 20 counts per bin.

The total effective observing time for B3 0731+438 is 22.75 ksec for PN and 27.40 ksec for both MOS1 and MOS2 with a source count rate of 0.01 cts/s for PN and 0.0033 cts/s for the coadded MOS data. The effective observing time for 3C 257 is 25.73 ksec for PN and 30.15 ksec for both MOS1 and MOS2, with a source count rate 0.0069 cts/s for PN and 0.0028 cts/s for the coadded MOS data. Finally, the background subtracted spectra were fitted using XSPEC version 11.1.0ae with the latest response matrices.

## 4 Spectral Analysis

The spectral analysis was carried out in the conventional way. The data were compared with a range of parametric models until the best fit (best reduced $`\chi `$2) was obtained. Seperate PN and MOS fits were carried out as a standard check for the fit and it was found that the data produced similar results to that of the coadded PN and MOS data. Hereafter, we present fits only to the coadded data with one sigma errors quoted throughout.

### 4.1 B3 0731+438

Initially we fitted the coadded PN and MOS spectra (0.3 – 10 keV) with a single powerlaw and neutral absorption (Galactic line-of-sight absorption) model; Model 1 in Table 1. Table 1 also includes values for the powerlaw photon index ($`\mathrm{\Gamma }`$), the associated intrinsic column density (NH), where appropriate, and values for $`\chi `$2 and number of degrees of freedom in the fit. Using a galactic column density NH = 6.2 $`\times `$ 1020 cm-2, Model 1 gives $`\mathrm{\Gamma }`$ $`=`$ 1.4 $`\pm `$ 0.15 ($`\chi `$2/d.o.f $`=`$ 29.15/32), suggesting some intrinsic absorption. The unabsorbed flux and luminosity (2 – 10 keV) were found from the fit to be 3.2 $`\times `$ 10-14 erg s-1 cm-2 and 2.8 $`\times `$ 1044 erg s-1 respectively.

The fit is improved (93 $`\%`$ significance) when intrinsic absorption is included in the spectral model; Model 2. This produces a steeper underlying powerlaw, $`\mathrm{\Gamma }`$ $``$ 1.8 and NH $``$ 3.3 $`\times `$ 1022 cm-2 ($`\chi `$2/d.o.f $`=`$ 26.11/31). This result suggests the source has a possible intrinisically scattered component, typical of a type-2 AGN. Therefore a partial covering model was fitted; Model 3, consisting of intrinsically absorbed and unabsorbed (scattered) powerlaw components. Model 3 provides an improved spectral fit for B3 0731+438 at greater than 99.9 $`\%`$ significance, with $`\mathrm{\Gamma }`$ $`=`$ 2.6 $`\pm `$ 0.4 and NH $`=`$ 1.9 $`\times `$ 1023 cm-2 $`\pm `$ 0.7 ($`\chi `$2/d.o.f $`=`$ 16.89/30). The unabsorbed X-ray luminosity (2 – 10 keV) is 9.1 $`\times `$ 1044 erg s-1. Figure 1 displays the best fit (Model 3) for B3 0731+438; the partial covering model. Figure 2 illustrates a contour plot of Model 3, showing the column density and X-ray luminosity distribution. The scattering efficiency was found to be 11 $`\%`$ $`\pm `$$`\%`$. Constraining the scattering efficiency to 10 $`\%`$ we find the results are consistent with the unconstrained fit. These results are consistent with the low scattering efficiencies for Seyfert 2s (Lumsden et al. 2001).

The spectra of B3 0731+438 may be consistent with a reflection dominated component, but fitting the data with the reflection model PEXRAV (Magdziarz & Zdziarski 1995) requires a very steep photon index of $`\mathrm{\Gamma }`$ $``$ 3.3 $`\pm `$ 0.2 and no direct view of the incident continuum. If we assume the reflected component is $``$ 5 $`\%`$, the implied upper limit on the intrinsic X-ray luminosity (2 – 10 keV), is 5.6 $`\times `$ 1045 erg s-1 if the planar reflection subtends 2$`\pi `$ sr at the illuminating source. Fitting the data with a broad (1 keV) FeK$`\alpha `$ emission line, we find a 90 $`\%`$ confidence upper limit on the intensity to be 2.7 $`\times `$ 10-6 photon s-1 cm-2, corresponding to a rest-frame equivalent width of 263 eV. This equivalent width is low, particularly for a type-2 object (Nandra et al. 1997; Turner et al. 1997).

### 4.2 3C 257

The analysis of the PN data for 3C 257 was complicated by a residual spectral feature at $``$ 1.5 keV due to incomplete subtraction of an internal detector background aluminium K-alpha line (Lumb et al. 2002). This line was removed from the spectral fitting by ignoring energies in the range 1.4 – 1.6 keV for the 3C 257 PN data during the analysis.

The models used for B3 0731+438 were again used for 3C 257. Model 1, the single powerlaw and neutral absorption model, using a Galactic value of NH $`=`$ 4.4 $`\times `$ 1020 cm-2, produced a very hard photon index $`\mathrm{\Gamma }`$ $`=`$ 1.0 $`\pm `$ 0.1 ($`\chi `$2/d.o.f $`=`$ 7.15/12). The flux and unabsorbed luminosity (2 – 10 keV) obtained from the fit were 2.6 $`\times `$ 10-14 erg s-1 cm-2 and 1.5 $`\times `$ 1044 erg s-1 respectively. As in the case of B3 0731+438 the fit for 3C 257 is improved, with a 90 $`\%`$ significance, when intrinsic absorption is included in the spectral model; Model 2. This again produced a steeper underlying powerlaw with $`\mathrm{\Gamma }`$ $``$ 1.5 and NH $``$ 5.0 $`\times `$ 1022 cm-2 ($`\chi `$2/d.o.f $`=`$ 4.04/11).

We then used the partial covering model, Model 3, but it is not possible to distinguish between this model and Model 2. Since it was difficult to constrain the fit, we fixed the scattering efficiency to 10 $`\%`$, consistent with that of B3 0731+438. We found $`\mathrm{\Gamma }`$ $``$ 1.4 and NH $``$ 5.5 $`\times `$ 1022 cm-2 with a corrected, unabsorbed X-ray luminosity (2 – 10 keV) $`=`$ 2.6 $`\times `$ 1044 erg s-1. Figure 3 displays the spectrum of source 3C 257, fitted using Model 3. Both models 2 and 3 suggest there is significant intrinsic absorption in 3C 257.

Fitting the data with the reflection model PEXRAV requires a steep photon index $`\mathrm{\Gamma }`$ $``$ 2.7 with no direct view of the incident continuum. Making the same assumptions as for B3 0731+438 , the implied upper limit on the intrinsic X-ray luminosity (2 – 10 keV), is 3.0 $`\times `$ 1045 erg s-1.

## 5 Discussion

The best-fit partial-scattering model for B3 0731+438 gives an intrinsic absorbing column of NH $``$ 1.9 $`\times `$ 1023 cm-2 which is comparable to the column density suggested by Risaliti, Maiolino & Salvati, (1999) for type-2 Seyferts of NH $`>`$ 1023 cm-2. This model also gave a photon index of $`\mathrm{\Gamma }`$ $``$ 2.6. This value is steep, but is comparable to the typical quasar value of $`\mathrm{\Gamma }`$ $`=`$ 1.9, at 90 $`\%`$ significance. Some luminous quasars have steep photon indices (e.g. $`\mathrm{\Gamma }`$ $``$ 2 for the quasar PDS 456 (Reeves et al. 1999)). For 3C 257 it was difficult to constrain the fit as a result of the limited number of photons. However, the spectral fits are consistent with this source also having a large X-ray luminosity and intrinsic absorption.

Both sources appear less luminous than predicted based on their \[OIII\] optical fluxes, but are still within the quasar range. The optical to X-ray spectral index $`\alpha `$ox was determined for B3 0731+438 using IR spectroscopy from Eales & Rawlings (1993) and Evans (1998). Assuming an optical extinction factor of Av $`=`$ 1.9 (Motohara et al. 2000) and using an absorption corrected X-ray flux, we calculate
$`\alpha `$ox $`=`$ \- 1.87 $`\pm `$ 0.2, consistent with sources studied by Bechtold et al. (2002).

For B3 0731+438 we can compare the ionizing photon number derived by Motohara et al. (2000) with the observed X-ray continuum flux. Integrating a powerlaw of photon index 2.5, consistent with that used by Motohara et al. (2000), between 13.67 eV and infinity, normalised to the unabsorbed X-ray continuum, we derive a value of 1.7 $`\times `$ 1056 photons s-1. Repeating this calculation with our derived value of $`\alpha `$ox, extrapolating to 13.67 eV using the SED of Zheng et al. (1997), we derive a value of 1.0 $`\times `$ 1057 photons s-1. These photon numbers are consistent with the result obtained by Motohara et al. (2000), of $`>`$ 2.1 $`\times `$ 1056 photons s-1 based on the H$`\alpha `$ luminosity, assuming case B and a unity covering factor in the H$`\alpha `$ cone (corrected to a cosmological scale of Ho $`=`$ 75 km s-1 Mpc-1 and qo $`=`$ 0.5). We can also compare the X-ray and H$`\alpha `$ luminosities assuming the relation that the intrinsic X-ray luminosity (2 – 10 keV) is approximately ten times that of the H$`\alpha `$ luminosity (Ward et al. 1988; Ho et al. 2001). Again correcting to a cosmological scale of Ho $`=`$ 75 km s-1 Mpc-1 and qo $`=`$ 0.5, Motohara et al. (2000) found a H$`\alpha `$ luminosity for B3 0731+438 of 6.7 $`\times `$ 1043 erg s-1, giving a predicted X-ray luminosity of 6.7 $`\times `$ 1044 erg s-1. This is in good agreement with our value obtained for B3 0731+438 of 9.1 $`\times `$ 1044 erg s-1. These comparisons suggest the X-ray emission is fairly isotropic.

## 6 Conclusions

In this paper we have presented XMM-Newton spectral analysis for the two NLRGs B3 0731+438 and 3C 257. The analysis of B3 0731+438 from fitting the spectra with a partial covering model indicates a spectral index of $`\mathrm{\Gamma }`$ $`=`$ 2.6 and a large column density NH $``$ 2 $`\times `$ 1023 cm-2. As a result of few counts, it was difficult to constrain the fit for 3C 257. However, the results suggest a large column density with a hard photon index, in agreement with an obscured source. Both objects have large absorption corrected X-ray luminosities (2 – 10 keV of 9 $`\times `$ 1044 erg s-1 for B3 0731+438 and 2.6 $`\times `$ 1044 erg s-1 for 3C 257.

In summary, the analysis has shown that both objects have intrinsic column densities and large X-ray luminosities well within the type-2 quasar range. We therefore conclude B3 0731+438 and 3C 257 are examples of type-2 quasar and currently provide the best X-ray spectra obtained of such objects with XMM-Newton.

## 7 ACKNOWLEDGMENTS

The work in this paper is based on observations with XMM-Newton, an ESA science mission, with instruments and contributions directly funded by ESA and NASA. The authors would like to thank the EPIC Consortium for all their work during the calibration phase, and the SOC and SSC teams for making the observation and analysis possible. This research has made use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space Flight Center. PMD would also like to acknowledge support from a PPARC studentship.