Electronic structure from equivalent differential equations of Hartree–Fock equations
Lin Hai
State Key Laboratory of High Field Laser Physics, Shanghai Institute of Optics and Fine Mechanics, Shanghai 201800, China

 

† Corresponding author. E-mail: linhai@siom.ac.cn

Abstract

A strict universal method of calculating the electronic structure of condensed matter from the Hartree–Fock equation is proposed. It is based on a partial differential equation (PDE) strictly equivalent to the Hartree–Fock equation, which is an integral–differential equation of fermion single-body wavefunctions. Although the maximum order of the differential operator in the Hartree–Fock equation is 2, the mathematical property of its integral kernel function can warrant the equation to be strictly equivalent to a 4th-order nonlinear partial differential equation of fermion single-body wavefunctions. This allows the electronic structure calculation to eliminate empirical and random choices of the starting trial wavefunction (which is inevitable for achieving rapid convergence with respect to iterative times, in the iterative method of studying integral–differential equations), and strictly relates the electronic structure to the space boundary conditions of the single-body wavefunction.

1. Introduction

Exactly solving the wavefunction of a quantum many-body system is the kernel task of several main branches of modern physics, such as that of nuclear structures in nuclear physics,[1,2] and that of the electronic structures of multi-electron atoms, molecules[3] and solids[4] in condensed matter physics. The difficulty in solving true/exact wavefunctions has promoted the progress of various approximation methods.[537] The Hartree–Fock (HF) equation is among the most widely used of these methods. Its kernel approximation is that the wavefunction of a many-fermion system can be expressed by the Slater determinant of the single-fermion wavefunction[57] as well as its various linear combinations.[8,9] Consequently, solving a multi-body wavefunction is simplified into solving the single-body wavefunction in a self-consistent field.

As an integral–differential equation, the HF equation is naturally attacked by the iterative method, the most popular method of solving integral–differential equations. Although this is a universal method, its efficiency is limited because achieving fast convergence of trial solutions with respect to the iterative times requires empirically and skillfully choosing the starting trial solution. If the convergence requirement is expressed as for any x, where is the initial trial function of variable x and represents the derived function after the i-th iteration, one can find this requirement to be the uniform convergence requirement which is the most severe in various convergence requirements. Clearly, whether convergence can be reached is closely related to the starting trial solution which needs to be chosen empirically and randomly in order to rapidly converge with respect to the iterative times.

This difficulty promotes solution of the integral–differential equation through its corresponding partial differential equation (PDE). The merit of doing so is to avoid having to empirically and randomly choose the starting trial solution and instead relate a strict solution to the conditions for determining a solution. But such a corresponding PDE is often based on approximations. For example, in density functional theory (DFT),[1016] the kernel approximation is to assume the integral term in the HF equation equals a summation of the products of single-body wavefunctions in finite number and hence transforms an integral–differential equation into a PDE with the same order (still a 2nd-order PDE). As shown below, if the single-body wavefunction is described with a PDE, the PDE should be at least 3rd-order, rather than 2nd-order.

Compared with the iterative method, the stricter method of solving the integral–differential equation solves its strictly equivalent PDE. For an integral–differential equation whose partial differential operatorʼs maximum order is N, it should be mathematically equivalent to at least an (N+1)-th-order PDE. Therefore, the HF equation, the maximum order of whose differential operator is 2, should be mathematically equivalent to at least a 3rd-order PDE of the single-body wavefunction. In the following section, a strict theory indicates that the HF equation is equivalent to a 4th-order PDE.

2. Theory and method
2.1. Main idea

The single-electron wavefunction in solids is identified by wavevector k. The HF equation of reads as

where the two real-valued functions UL and Uc satisfy the Poisson equations (Ql and Rl are the charge and the position of the l-th ion) and the complex-valued function reads as
where w is defined as . Here, because the summations over q and over k + q are different expressions of the same object, w is k-independent. Besides the HF equation, each also satisfies
where is the volume of a cell, a is the lattice constant, VFS is the volume of the Fermi sphere, and is the total ionic charge per cell. Moreover, the are orthogonal to each other.

In DFT,[1016] the exchange–correlation energy is assumed to be in a simple functional of n and , where n is the total electronic density and then each HF equation, which is an integral–differential equation of , is approximated to be a differential one with the same order, , a 2nd-order PDE. As previously mentioned, if the single-body wavefunction is described with a PDE, the PDE should be at least 3rd-order, rather than 2nd-order.

Actually, it is unnecessary to depend on those assumptions for . There exists a strict and straightforward method of solving the HF equations. By denoting as , we can re-write each HF as

where Ek is space-independent and also time-independent. All terms in this form belong to one of two classes: one is k-dependent and the other k-independent. The condition for an equation in such a form having solutions at any k-value is that all k-independent terms in Eq. (4) form a closed relation and the k-dependent terms are governed by a few k-independent terms,
where E0 is the corresponding atomic level of those extended-states and ** is a constant . Equation (5) is a 2nd-order linear PDE of w, and its solution is determined by a physical boundary condition which means the density of valence electrons at the cation position is 0. Moreover, the linearity of Eq. (5) can warrant to be satisfied easily by multiplying the solution of Eq. (5) by a constant.

Due to the (1/r)-dependent part of the operator in Eq. (5), exists if . But the mathematical property of the isotropic part of UL, or the fact that is a power series of , determines the expansion coefficients wi of as 0 for all odd i in the case of . This might be unable to warrant the physical boundary condition to be fulfilled when . In contrast, if , it is possible for to appear, which is favorable for the fulfillment of the physical boundary condition.

After knowing the profile of the total particle density from Eq. (5), each , or each “relative probability amplitude” , can be solved from Eq. (6). The real part and the imaginary part of Eq. (6) yield

where
Having obtained the expressions of VS and VC in terms of sk and from Eqs. (7) and (8), we can substitute them into Eqs. (9) and (10) and have
Comparing the -dependent terms, as well as the -dependent terms, in Eqs. (14) and (15), we obtain 4 equations respectively for their -dependent terms and -dependent terms. The four equations are in fact 2 pairs and in each pair the two equations are the same because one is merely the other multiplied by a constant. Thus, we obtain two equations:

Here, we utilize the formulas

Note that equations (16) and (17) are a 4th-order linear PDE set of sk. Therefore, equation (3) can be easily satisfied by multiplying the solution of Eqs. (16) and (17) by a constant.

Equations (16) and (17) are of the same general form: and , where and are space–time functions and is a space–time-independent constant satisfying . In principle, equations (16) and (17) and can finally lead to a self-consistent higher-order PDE of sk, whose solution can lead to a self-consistent value of . In the following subsection, we explain in detail how to derive self-consistently an electron structure from Eqs. (16) and (17) and .

2.2. Detailed mathematics

Solutions of Eqs. (16) and (17) are determined by the boundary conditions . For the simplicity of symbols, we first present a detailed formula for the isotropic case where all quantities are -independent. By introducing the power series expressions and and comparing terms proportional to in Eqs. (16) and (17) order-by-order, we can find that equations (16) and (17) lead to a recurrence formula of (see later Eqs. (30) and (44)). Because the isotropic parts of and of read as

Equations (16) and (17) are of the same general form , where are analytic functions of r. Note that equation (16) does not contain any terms that are proportional to and hence exists. In contrast, the term in Eq. (17) will have a contribution proportional to , , through in . But another term has an opposite contribution proportional to , through in . Therefore, equation (17) yields

Because of Eqs. (20) and (21), there exist
and

Using the power series of , , , and , we find that equations (16) and (17) lead to , or . Likewise, the fact that all terms that are proportional to in Eqs. (16) and (17) should be 0, leads to and .

and read as

and and read as
Moreover, from Eq. (17), we can express as
where analytic functions C0 and C1 are defined as

Thus, can be written as

or

Equation (35) requires its -related part to equal 0:

This requirement leads to

2.3. Main theoretical results

Equations (27)–(29) and (37) can lead to a close relation to and the boundary conditions

Because of the formulas
we finally obtain a quadratic equation of as follows:
where

Equation (30) is the lowest-order version, or the i = 0 version, of the recurrence formula

2.4. Feasibility of extension to more realistic anisotropic case

The same procedure can be extended to more realistic anisotropic cases. UL is -dependent,

where
with being dependent on ** . Note that ** and the spherical harmonic function satisfy a common equation
Therefore, corresponding to the expression in Eq. (45), w and sk have similar expressions . Namely, when UL is -dependent, equation (4) will correspond to a set of equations, each of which is for all terms associated with . The radial coefficient function in each can be solved from this set, where . Note that those -components should warrant
Thus, for -dependent , we can solve, in a procedure similar to that presented previously, and from the -related terms of Eqs. (5) and (6).

3. Numerical results and discussion

Because an HF equation is a differential–integral equation, its Green function (GF) will contain infinite terms related to the unperturbed GF which is associated with a differential equation. Although many terms are small enough to be negligible, their number is so large that their summation is not negligible. Therefore, the GF of an HF equation is difficult to exactly solve. If the HF equation is approximated as a differential equation by making assumptions for , such as DFT,[1021] it will be a troublesome task to verify the self-consistency of these assumptions because of a fundamental discrepancy: the differential equation that the single-body wavefunction satisfies should be at least 3rd-order, rather than the 2nd-order that the approximation assumes.

In contrast, our method does not require any assumption for . It is based on the inherent mathematical properties of the coupled HF equations under the same constraint condition and hence self-consistency can be warranted fully.

We use this method to calculate the 3d-band of a transition metal, i.e., the dimensionless parameter is chosen to be . Two elements of the body-centered cubic (bcc) crystal type and two elements of the face-centered cubic (fcc) type are chosen, and their lattice constants can be found in some typical textbooks.[38,39] Their 3d-bands all contain multiple electrons and hence ought to be studied by many-body methods. For the bcc type, the origin of the space coordinate frame is chosen to be located at the center of a face of the cell. For the fcc type, the origin is chosen to be located at the center of the cell. Obviously, applying LCAO to the 3d-band cannot warrant the HF equation to be satisfied, because the Hamiltonian of the molecule/solid is not a linear combination of the atomic Hamiltonian.

According to the solutions of Eqs. (41)–(43) in the case of , is , where G and are both constant. The curves of the density of states (DOS) are plotted in Fig. 1 and the Fermi level EF is shown in Fig. 2. Because what is calculated in the present research is for approximate cases where UL is taken as isotropic and non-magnetic, it is unrealistic to expect the calculated results to accord complet ely with experimental measurements. But Figure 2 indicates that the density of states at the Fermi surface is nearly inversely proportional to the number of the 3d-valence. The when the 3d-shell is half-filled is higher than the density of states when the 3d-shell is nearly completely filled. This is a qualitatively reasonable result for a non-magnetic Hamiltonian.

Fig. 1. DOS curves.
Fig. 2. Values of DOS at Fermi level.
4. Conclusions

Solving an integral–differential equation through its strictly equivalent PDE is stricter and more efficient than directly applying an iterative method to it. For the HF equation, that its Coulomb-type integral kernel function obeys the Poisson equation determines its equivalence to a 4th-order nonlinear PDE. According to the standard theory of the PDE,the solution of such a 4th-order PDE depends on its boundary conditions. Having converted the differential–integral equation into its equivalent higher-order differential counterpart according to the standard strict procedure described in the text, we can find that the equivalent differential equation contains “potentials” of space singularity and . As shown in Eqs. (16) and (17), this space singularity will affect the allowed values of the space boundary conditions of an analytic wavefunction. By analyzing the allowed values of the space boundary conditions from Eqs. (22)–(26), we obtain a detailed formula reflecting the dependence of the wavefunction energy on the space boundary conditions of the wavefunction modulus, and finally obtain the dispersion relation of modulated plane-waves, which are the allowed solutions of the HF equation. This corresponds to a stricter method of calculating electronic structures within the framework of the single-electron problem.

Reference
[1] Ring P Schuck P 1980 The nuclear many-body problem Springer Verlag
[2] Hodgson P 1971 Nuclear Reaction and Nuclear Structure Oxford University Press
[3] Levine I 1999 Quantum Chemistry Prentice Hall
[4] Ziman J M 1972 Principle of the Theory of Solids Cambridge University Press
[5] Hartree D R 1928 Proc. Cambridge Phil. Soc. 24 89
[6] Hartree D R 1928 Proc. Cambridge Phil. Soc. 24 111
[7] Fock V 1930 Z. Physik. 61 126
[8] Sherrill C D Schaefer H F 1999 The Configuration Interaction Method: Advance in Highly Correlated Approaches. Advances in Quantum Chemistry 34 Lowdin P O San Diego Academic Press 143
[9] Cramer C J 2002 Essentials of Computational Chemistry Chichester John Wiley & Sons 191
[10] Abrikosov A A Gorkov L P Dzyaloshinskii Y E 1965 Quantum Field Theoretical Methods in Statistical Physics New York Pergamon
[11] Georges A Kotliar G Krauth W Rozenberg M J 1996 Rev. Mod. Phys. 68 13 and references therein
[12] Kotliar G 2006 Rev. Mod. Phys. 78 865 and references therein
[13] Freericks J K Zlatic V 2003 Rev. Mod. Phys. 75 1333 and references therein
[14] Bickers N E Scalapino D J 1989 Ann. Phys. 193 206
[15] Aryasetiawan F Gunnarson O 1998 Rep. Prog. Phys. 61 237
[16] Hedin L Lundqvist S 1969 Solid State Phys. 23 New York Academic 1
[17] Hohenberg P Kohn W 1964 Phys. Rev. 136 B864
[18] Kohn W Sham L J 1965 Phys. Rev. 140 A1133
[19] Sham L J Kohn W 1966 Phys. Rev. 145 561
[20] Dreizler R M Gross E K U 1990 Density Functional Theory New York Springer
[21] Alonso J A Girifalco L A 1978 Phys. Rev. B 17 3735
[22] Janak J F Moruzzi V L Williams A R 1975 Phys. Rev. B 12 1257
[23] Perdew J P 1986 Phys. Rev. B 33 8822
[24] Perdew J P Wang Y 1986 Phys. Rev. B 33 8800
[25] Langreth D C Perdew J P 1980 Phys. Rev. B 21 5469
[26] Perdew J P Burke K Wang Y 1996 Phys. Rev. B 54 16533
[27] Perdew J P Levy M 1983 Phys. Rev. Lett. 51 1884
[28] Berke A D 1986 J. Chem. Phys. 85 7184
[29] Berke A D 1988 J. Chem. Phys. 88 1053
[30] Berke A D 1993 J. Chem. Phys. 98 1372
[31] Berke A D 1988 Phys. Rev. A. 38 3098
[32] Martin R M 2008 Electronic Structure Cambridge University Press
[33] Harrison W A 1966 Pseudopotentials in the theory of metals New York Benjamin
[34] Yu Z M Liu Y L 2014 Chin. Phys. Lett. 31 017103
[35] Peng P Jiang Y Q 2018 Acta Phys. Sin. 67 132101 in Chinese
[36] Fang Y Z Liu J H Kong X J 2018 Acta Phys. Sin. 67 117101 in Chinese
[37] Pan F C Chen H M Lin X L 2018 Acta Phys. Sin. 67 107701 in Chinese
[38] Ashcroft N W Mermin N D 1976 Solid State Physics Harcourt Inc 2
[39] Kittel C 1953 Introduction Solid State Physics New York John Wiley Sons 233 Table 12.1