Correlated Band Structure of a Transition Metal Oxide ZnO Obtained from a Many-Body Wave Function Theory

TL;DR

This paper demonstrates the first application of the BiTC many-body wave function method to compute the band structure of ZnO, showing improved accuracy over conventional methods for correlated materials.

Contribution

It introduces the BiTC method for solid-state calculations with d electrons, providing a parameter-free, systematically improvable approach for correlated materials.

Findings

BiTC yields band gaps closer to experimental values

BiTC outperforms conventional methods in accuracy

First successful application of BiTC to a transition metal oxide

Abstract

Obtaining accurate band structures of correlated solids has been one of the most important and challenging problems in first-principles electronic structure calculation. There have been promising recent active developments of wave function theory for condensed matter, but its application to band-structure calculation remains computationally expensive. In this Letter, we report the first application of the biorthogonal transcorrelated (BiTC) method: self-consistent, free from adjustable parameters, and systematically improvable many-body wave function theory, to solid-state calculations with d electrons: wurtzite ZnO. We find that the BiTC band structure better reproduces the experimental values of the gaps between the bands with different characters than several other conventional methods. This study paves the way for reliable first-principles calculations of the properties of strongly correlated materials.

Tables

Table 1. Table 1: Some characteristic values in the band structures of ZnO, as calculated by several methods  Tab1_note . The bottoms of the Zn- 3 ​ d 3 𝑑 3d and O- 2 ​ p 2 𝑝 2p bands are evaluated at the Γ Γ \Gamma point where the valence-band top is set to zero for each method. The errors between the calculated and experimental values shown in Ref.  ZnOexpt_old are presented in parentheses. For LDA, the O- 2 ​ p 2 𝑝 2p bottom and Zn- 3 ​ d 3 𝑑 3d -averaged levels are not presented here because of the overlap among the O- 2 ​ p 2 𝑝 2p and Zn- 3 ​ d 3 𝑑 3d bands. For G 0 ​ W 0 superscript 𝐺 0 superscript 𝑊 0 G^{0}W^{0} (HSE03), the Zn- 3 ​ d 3 𝑑 3d bottom position was read from the density of states presented in Ref.  ZnOHSEGW . All values are in eV. a Ref.  ZnOHSEGW . b Ref.  GWZnO . c Ref.  ZnOARPES . d Ref.  ZnOQSGW . e Ref.  ZnOscGW . f Ref.  ZnOAFQMC . g Ref.  ZnO_DMC . h Ref.  ZnOexpt_old . i Ref.  ZnO3d . j Ref.  ZnO3d_3 . k Ref.  ZnO3d_2 .

		Band gap [Error]		O-$2p$ bottom [Error]		Zn-$3d$ averaged	Zn-$3d$ bottom
DFT	LDA	0.7	[$-$2.7]	$-$	$-$	$-$	$-$5.8
	HSE03a	2.1	[$-$1.3]	$-$4.9	[$+$0.4]	$-$	$-$6.5
$GW$	$G^0{W}^0$(LDA)b	2.4	[$-$1.0]	$-$5.2	[$+$0.1]	$-$	$-$6.5
	$G^0{W}^0$(HSE03)	3.2a, 3.46c	[$-$0.2, $+$0.06]	$-$	$-$	$-$6.21c	$-$7.2a
	$G^0{W}^0$ (GGA+$U$)c	2.94	[$-$0.46]	$-$5.6	[$-$0.3]	$-$6.33	$-$7.1
	($U-J=$ 6 eV)
	$G^0{W}^0$+$V_d$ (GGA+$U$)c	3.30	[$-$0.1]	$-$5.5	[$-$0.2]	$-$7.45	$-$8.0
	($U-J=$ 6 eV, $V_d=$ 1.5 eV)
	QS$GW$d	3.87	[$+$0.47]	$-$5.3	[$\pm$ 0]	$-$	$-$7.2
	sc$GW$ (RPA) e	3.8	[$+$0.4]	$-$	$-$	$-$6.4	$-$
	sc$GW$ ($e$-$h$)e	3.2	[$-$0.2]	$-$	$-$	$-$6.7	$-$
WFT	AFQMCf	3.26(16)	[$-$0.14]	$-$	$-$	$-$	$-$
	DMCg	3.8(2)	[$+$0.4]	$-$	$-$	$-$	$-$
	BiTC	3.1	[$-$0.3]	$-$5.1	[$+$0.2]	$-$9.1	$-$9.7
	HF	11.4	[$+$8.0]	$-$5.7	[$-$0.4]	$-$9.3	$-$9.9
Expt.		3.4h		$-$5.3h, $-$5.2(3)i		$-$7.5c, $-$7.5(2)j, $-$8.5(4)k,	$-$
						$-$8.6(2)k, $-$8.81(15)i

Equations

$$\mathcal{H}\Psi =E\Psi \iff {\mathcal{H}}_{\mathrm{TC}}\Phi =E\Phi ({\mathcal{H}}_{\mathrm{TC}}=F^{-1}\mathcal{H}F),$$

$$\left(-\frac{1}{2}{\nabla }_1^2+v_{\mathrm{ext}}(x_1)\right){\varphi }_i(x_1)$$

$$\displaystyle+\sum_{j=1}^{N}\int\mathrm{d}x_{2}\ \phi_{j}^{*}(x_{2})v_{\mathrm{2body}}(x_{1},x_{2})\mathrm{det}\left[\begin{array}[]{rrr}\phi_{i}(x_{1})&\phi_{i}(x_{2})\\ \phi_{j}(x_{1})&\phi_{j}(x_{2})\\ \end{array}\right]$$

$$-\frac{1}{2}\sum _{j=1}^N\sum _{k=1}^N\int \mathrm{d}{x}_2\mathrm{d}{x}_3{\varphi }_j^{\ast }(x_2){\varphi }_k^{\ast }(x_3)v_{3body}(x_1,x_2,x_3)$$

$$\displaystyle\times\mathrm{det}\left[\begin{array}[]{rrr}\phi_{i}(x_{1})&\phi_{i}(x_{2})&\phi_{i}(x_{3})\\ \phi_{j}(x_{1})&\phi_{j}(x_{2})&\phi_{j}(x_{3})\\ \phi_{k}(x_{1})&\phi_{k}(x_{2})&\phi_{k}(x_{3})\end{array}\right]=\sum_{j=1}^{N}\epsilon_{ij}\phi_{j}(x_{1}),$$

$$v_{2body}(x_1,x_2)\equiv \frac{1}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}+\frac{1}{2}\sum _{i=1}^2[{\nabla }_i^{2}u(x_1,x_2)$$

$$-({\nabla }_{i}u(x_1,x_2){)}^2+2{\nabla }_{i}u(x_1,x_2)\cdot {\nabla }_i],$$

$$v_{3body}(x_1,x_2,x_3)\equiv {\nabla }_{1}u(x_1,x_2)\cdot {\nabla }_{1}u(x_1,x_3)$$

$$+{\nabla }_{2}u(x_2,x_1)\cdot {\nabla }_{2}u(x_2,x_3)+{\nabla }_{3}u(x_3,x_1)\cdot {\nabla }_{3}u(x_3,x_2),$$

$$u(x,x^{\prime })=\frac{A}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}\left(1-\exp \left(-|\mathbf{r}-{\mathbf{r}}^{\prime }|/C_{\sigma ,{\sigma }^{\prime }}\right)\right),$$

$$E_{gap}(N_k)\simeq E_{gap}(N_k=\infty )+C/N_k,$$

Full text

To reveal fertile and nontrivial physics in condensed matter, first-principles electronic-structure calculation has established itself as an indispensable tool in recent studies. For this purpose, density functional theory (DFT) HK ; KS has played a leading role and has been applied to various materials; however, the limitations of this theoretical framework have come to light. One of the major problems is an inaccurate description of strong electron correlations, e.g., in transition metal oxides. The $GW$ method GW1 ; GW2 ; GW3 is a promising way to ameliorate the inaccuracy of the band structures and has been applied to several solids, including $d$-electron systems. However, because the $GW$ method is often applied without satisfying self-consistency, a nontrivial dependence on the initial DFT calculations is introduced. It has also been reported that the $GW$ method sometimes exhibits severe difficulty in obtaining converged results GWconvergence , owing to its perturbative nature. Another possible choice is to construct effective models from DFT that include correlation terms such as Hubbard $U$ and to solve the models using elaborated methodologies Models such as dynamical mean-field theory DMFTreview ; DMFTreview2 ; DMFT2 . However, the exact correspondence between the effective models and the first-principles Hamiltonian is a nontrivial problem.

Recently, wave function theory (WFT), which had been mainly applied to molecular systems and established itself as the gold standard in theoretical chemistry Szabo , has become a promising alternative to DFT for accurate descriptions of electron correlation in solids FCIQMC . Among the most powerful frameworks in WFT are first-principles quantum Monte Carlo (QMC) methods QMCreview ; QMCreview2 , such as the variational Monte Carlo (VMC), diffusion Monte Carlo (DMC), auxiliary-field quantum Monte Carlo (AFQMC), and full-configuration-interaction (FCI) QMC methods. Other kinds of WFT, called post-Hartree-Fock (post-HF) methods, have also been applied to condensed matter in recent years MP2_Kresse_band ; MP2_Kresse_totE ; MP2_Shepherd ; CCD_Shepherd ; CCD_Shepherd2 . However, their targets are in most cases limited to solids with small unit cells, owing to their expensive computational cost. In addition, the correlated band structure, which is quite useful in various kinds of theoretical analyses, is not easily obtained in many WFTs. For example, calculation of the band structure in the framework of VMC or DMC requires a large number of single-point calculations of the excited states, which is a clear difference from a mean-field-like approach such as DFT, whereby the whole band structure is obtained at once.

From this viewpoint, the transcorrelated (TC) method BoysHandy ; Handy ; Ten-no ; Umezawa is a fascinating WFT that can be applied to solids with reliable accuracy and moderate computational cost Sakuma ; TCaccel ; TCjfo ; TCCIS ; TCMP2 ; TCPW . The TC method adopts the so-called Jastrow ansatz, which is based on a promising strategy often adopted in several WFTs such as QMC methods to describe strong electron correlations; i.e., the electron-electron distance is included into many-body wave functions. Explicitly correlated electronic-structure theory R12 in quantum chemistry also adopts this strategy. In fact, the Gutzwiller- and Jastrow-correlation factors have often been used to describe strong electron correlation, including the Mott physics in systems such as the Hubbard model VMC ; MottVMC1 ; MottVMC2 ; MottVMC3 . It is also important to note that, unlike several WFTs, the whole band structure is obtained at once by solving a one-body self-consistent-field (SCF) equation in the TC method, as described later in this Letter. Moreover, the TC method is deterministic, i.e., free from the statistical error, unlike the QMC methods. Accurate calculations for the Hubbard model TCHubbard ; LieAlgebra and molecular systems LuoTC ; LuoVTC ; CanonicalTC were also reported using the TC method or other theories that have a close relationship with the TC method. However, insofar as solid-state calculations are concerned, the TC method has so far been applied only to weakly correlated systems.

In this Letter, we present the first application of the TC method to the band-structure calculation of a $d$-electron system: wurtzite ZnO. $3d$ transition metal oxides have posed theoretical challenges for first-principles band-structure calculations, as it is well known that popular approximations such as the local-density approximation (LDA) fail to provide their band structures accurately, as we shall see later. We find that the TC method with the bi-orthogonal formulation (the BiTC method) BiTC ; TCMP2 successfully reproduces the experimental band structure of ZnO. We also clarify how the Jastrow factor improves the first-principles description of the correlated electronic states through comparison of the band structures and electron densities among the BiTC and other methods.

The central concept of the TC method is to make use of the similarity transformation of the many-body Hamiltonian with the Jastrow factor $F=\mathrm{exp}(-\sum_{i,j}u(x_{i},x_{j}))$,

[TABLE]

where the correlated wave function is represented as $\Psi=F\Phi$ and $\mathcal{H}_{\mathrm{TC}}$ is called the TC Hamiltonian. Here, $x$ denotes a pair of space and spin coordinates: $x=(\mathbf{r},\sigma)$. By adopting the so-called Slater-Jastrow ansatz, $\Phi$ becomes a Slater determinant consisting of one-electron orbitals, $\phi(x)$: $\Phi=\mathrm{det}[\phi_{i}(x_{j})]$, and Eq. (1) yields an SCF equation for one-electron orbitals that experience the effective interaction described with the TC Hamiltonian:

[TABLE]

where $v_{\mathrm{ext}}(x_{1})$, $v_{\mathrm{2body}}(x_{1},x_{2})$, and $v_{\mathrm{3body}}(x_{1},x_{2},x_{3})$ are the external potential including the nucleus-electron interaction note_Vext and the two- and three-body effective interactions in the TC Hamiltonian, defined as

[TABLE]

and

[TABLE]

respectively. As is evident, the HF method can be regarded as the TC method with $u=0$. Owing to this effective one-body picture, it is possible to treat the many-body correlation with moderate computational cost. In addition, one can obtain the band structure of the quasiparticles by using the real part of the eigenvalues of the $\epsilon$ matrix on the right-hand side of Eq. (8). One of the authors proved in prior work Umezawa that such a use of the $\epsilon$ matrix as quasiparticle energies is consistent with Koopmans’ theorem. We note that one can systematically improve the accuracy of the TC method by utilizing quantum chemical methodologies such as the coupled-cluster and configuration interaction methods to go beyond a single Slater determinant (e.g. Refs. Ten-no ; BiTC ; TCMP2 ).

Here, the Jastrow function, $u(x,x^{\prime})$, is set to the following simple form without adjustable parameters: Sakuma ; QMCreview ; Ceperley ; CeperleyAlder

[TABLE]

where $A=\sqrt{V/(4\pi N)}$ ($N$ is the number of valence electrons in the simulation cell, $V$ is the volume of the simulation cell) and $C_{\sigma,\sigma^{\prime}}=\sqrt{2A}$ (spin parallel:$\ \sigma=\sigma^{\prime}$), $\sqrt{A}$ (spin antiparallel:$\ \sigma\neq\sigma^{\prime}$). The long-range asymptotic form of this function describes the screening effect of the electron-electron Coulomb interaction BohmPines . The short-range behavior of the exact Jastrow function should obey the cusp condition cusp ; cusp2 ; cusp_note . The Jastrow ansatz adopted here works well for state-of-the-art QMC methods QMCreview ; QMCreview2 . Although our choice of the Jastrow function is rather simple, we shall see that, nevertheless, it works well not only for weakly correlated systems Sakuma ; TCaccel but also for the 3$d$-electron system. Of course, it is possible to improve a quality of the many-body wave function by using a complicated Jastrow factor, but we adopted this simple trial wave function to realize moderate computational cost.

In this study calc_note , we adopted the BiTC method BiTC ; TCMP2 , in which the left one-electron orbitals, $\chi(x)$ in the left Slater determinant X $=\mathrm{det}[\chi_{i}(x_{j})]$, replace the bra orbitals in the SCF equation (8), while the ket orbitals remain $\phi(x)$. Because the TC Hamiltonian is non-Hermitian, $\phi(x)$ and $\chi(x)$ become different. We do not show the band structure calculated using the TC method without the bi-orthogonal extension here, because of the large imaginary part of the eigenvalues Sppl .

Figure 1 presents the band structures of ZnO calculated with the LDA, all-electron $G^{0}W^{0}$ starting from LDA using the LAPW LAPW method, BiTC, and HF methods. The characteristic energy values in these band structures, as evaluated at the $\Gamma$ point, are listed in Table 1. Table 1 also lists calculated values with various other methods. By focusing on the gaps between the bands with different characters as listed in Table 1, the BiTC band structure exhibits better accuracy than many other conventional methods, including the $G^{0}W^{0}$ method starting from LDA, and accuracy comparable to that for the most accurate varieties of the $GW$ scheme, i.e., the $G^{0}W^{0}$ method starting from HSE03 and the self-consistent $GW$ methods such as the QS$GW$ method QSGW . It is noteworthy that both the $GW$ and BiTC methods yield such successful results despite being based on conceptually different formulations. It is important that the BiTC method is based on the self-consistent formulation and thus is independent of DFT calculations, whereas the $G^{0}W^{0}$ method strongly depends upon the unperturbed DFT calculations, as seen in Table 1. We should also note that, while some methods employ parameters $U$ and $V_{d}$, which are difficult to determine in the ab initio way, the BiTC method does not use such parameters.

Consistency of the calculated band gaps among WFTs that use similar trial wave functions (AFQMC, DMC, and BiTC) is also remarkable. We again stress that the whole band structure is not easily obtained in QMC simulations and requires many single-point excited-state calculations. We can see that the band structures calculated with the BiTC and HF methods are very similar, except for the band gap. Therefore, the main role of the Jastrow factor used in this study on the band structure seems to be improvement of the size of the band gap through the screening effect of the electron–electron interaction, which is described with the long-ranged asymptotic behavior of the Jastrow factor. To obtain a more accurate band structure, e.g., with respect to the depth of the Zn-3$d$ bands, more elaborated Jastrow factors, as used in QMC studies QMCreview ; QMCreview2 , will be necessary. Because the Zn-$3d$ bands are almost flat, a key point might be accurate description of the atomic states, which is an important issue for future investigation.

For the $GW$ method, it was pointed out that the shallow Zn-$d$ bands can strengthen the $p$-$d$ hybridization, and thus can result in underestimation of the band gap ZnOARPES ; ZnOscGW . A similar situation might also be realized in the BiTC band structure, whereas the Zn-$d$ bands are rather deep note_GW and so the band gap can be overestimated. However, the BiTC band gap is also affected by the $A$ parameter in the Jastrow factor, as mentioned in the above comparison between the BiTC and HF band structures. Because the $A$ parameter used in this study was determined by RPA analysis of the uniform electron gas, it can cause overscreening in the insulator TCjfo , thereby decreasing the band gap. One possibility is that these two factors are canceled here, but more detailed investigation on other materials is also an important future issue note_GW2 .

Figure 2 presents the electron densities calculated with the LDA, HF, and BiTC methods on the line shown in the crystal structure. The electron density obtained by the BiTC method is defined as $n(\mathbf{r})=$ Re$[\sum_{i=1}^{N}\chi_{i}^{*}(\mathbf{r})\phi_{i}(\mathbf{r})]$, where the condition $\int\mathrm{d}\mathbf{r}\ n(\mathbf{r})=N$ is satisfied due to the bi-orthonormalization condition $\langle\chi_{i}|\phi_{j}\rangle=\delta_{ij}$. We can see that the electron densities of these methods are almost the same. However, a slight increase of the electron density at the atomic sites is observed for the BiTC method compared with the others. Such a tendency is consistent with the fact that the strong divergence of the electron-electron Coulomb repulsion is alleviated for the effective two-body interaction in the similarity-transformed TC Hamiltonian, because the Jastrow factor satisfies the cusp condition. More concretely, the $\nabla_{1}^{2}u(x_{1},x_{2})$ and $\nabla_{2}^{2}u(x_{1},x_{2})$ terms in Eq. (9) yield $1/|\mathbf{r}_{1}-\mathbf{r}_{2}|$ divergence with a different sign than the electron-electron Coulomb repulsion. As can be seen in the proof of the cusp condition cusp ; cusp2 , the true many-body wave function should exhibit deformation described with the two-body degrees of freedom near the electron-electron coalescence point, which cannot be represented solely with one-body degrees of freedom. This is a characteristic advantage of the Slater-Jastrow-type wave function for the description of localized electronic states, such as in strongly correlated systems. It is noteworthy that the atomic calculations of the TC method also exhibit a similar tendency for localization Umezawa .

Finally, we mention the computational effort required for the BiTC calculation. Computation takes place on time scales given by $\mathcal{O}(N_{k}^{2}N_{b}^{2}N_{pw}\log N_{pw})$, where $N_{k}$, $N_{b}$, and $N_{pw}$ are the numbers of $\mathbf{k}$-points, occupied bands, and plane waves, respectively note_computation . This is the same order as that for the HF or hybrid DFT calculations with a prefactor about 20 to 40 note_cmptime . The BiTC calculation involves neither the frequency index nor the convergence with respect to the number of conduction bands, unlike some perturbative methodologies such as the $GW$ method. As can be seen from the fact that the hybrid DFT calculations have now been applied to various periodic systems, the computational cost of the BiTC method is reasonable for solid-state calculations. One remaining obstacle for wide application of the BiTC method is that we use the norm-conserving pseudopotential with a very high cutoff energy to handle the semicore states at present. However, this problem is not inherent to the BiTC method and can be overcome in principle by the development of a pseudopotential formalism such as the PAW method PAW1 ; PAW2 adapted to the TC method, which is an important future issue note_PAW . We also note that the deformation of the electron density near atomic sites shown in Fig. 2 also implies the importance of careful treatment for the core states. This can also be a common problem for QMC calculations using the Jastrow correlation factor.

To conclude, we apply the bi-orthogonal version of the TC method to wurtzite ZnO and find that it well reproduces the experimental band structure. Our study encourages further investigation of other strongly correlated materials using the BiTC method.

I Supplemental Material

II Convergence with respect to the number of $\mathbf{k}$-points

We have performed band-structure calculations using a $2\times 2\times 2$ $\mathbf{k}$-mesh (32 atoms in the unit cell) for the HF and BiTC methods to confirm the convergence with respect to the number of $\mathbf{k}$-points. For the BiTC calculations, the band gap and the relative energy of the Zn-$d$-band bottom to the valence-band top are 2.92 and $-9.73$ eV, respectively, for a $2\times 2\times 2$ $\mathbf{k}$-mesh, whereas they are 3.09 and $-9.68$ eV, respectively, for a $3\times 3\times 3$ $\mathbf{k}$-mesh (108 atoms in the unit cell). For the HF calculations, the corresponding values are 11.53 and 10.07 eV, respectively, for a $2\times 2\times 2$ $\mathbf{k}$-mesh, whereas they are 11.42 and $-9.90$ eV, respectively, for a $3\times 3\times 3$ $\mathbf{k}$-mesh. The difference between the values calculated using the two $\mathbf{k}$-meshes is less than 0.2 eV for all of these values.

In our previous study on the TC method TCaccel , we saw that the calculated band gap $E_{gap}(N_{k})$ with $N_{k}$ $\mathbf{k}$-points obeys an approximate relation:

[TABLE]

where $C$ is a constant. Using this approximate relation for the band energies, the residual finite-size error in our calculations using a $3\times 3\times 3$ $\mathbf{k}$-mesh is estimated to be less than 0.1 eV. Therefore, we can say that the convergence with respect to the number of $\mathbf{k}$-points is sufficiently achieved. Although it is possible that the approximate relation that is used above does not hold for this case ZnO_DMC , we can expect that the convergence will be sufficient for our discussion because the finite-size error is expected to be significantly reduced in a $3\times 3\times 3$ $\mathbf{k}$-mesh from its value for a $2\times 2\times 2$ one.

III Imaginary part of the TC eigenvalues without the bi-orthogonal formulation

We do not show the band structure calculated using the TC method without the bi-orthogonal extension in the main text. This is because we found that the imaginary part of the eigenvalues of the Zn-$3d$ bands obtained by solving the one-body SCF equation became as large as 0.1 eV in the TC method, which is more than 1,000 times larger than that in the BiTC method. Because the non-Hermiticity in the (Bi)TC formulation originates from the similarity-transformation, the presumed equivalency between the TC solution and the true many-body eigenstate is suggested to break when the imaginary parts of the eigenvalues become too large. The difference between the band structures obtained with the TC and BiTC methods was small for weakly correlated systems except LiF that is rather localized-electron system TCMP2 .

IV Electron-density differences among the LDA, HF, and BiTC methods

The differences in the electron densities calculated with the LDA, HF, and BiTC methods are shown in Fig. S1. In Fig. S1, we can verify the tendencies explained in the main text, i.e., that the BiTC electron density tends to be localized near the atomic sites. In addition, we can see an increase of the BiTC electron density along the Zn-O bonds near the oxygen atoms.