DFT+DMFT calculations of the complex band and tunneling behavior for the transition metal monoxides MnO, FeO, CoO and NiO
Long Zhang, Peter Staar, Anton Kozhevnikov, Yun-Peng Wang, Jonathan Trinastic, Thomas Schulthess, Hai-Ping Cheng
DDFT + DMFT calculations of the complex band and tunneling behavior for the transition metalmonoxides MnO, FeO, CoO and NiO
Long Zhang , Peter Staar , ∗ Anton Kozhevnikov , Yun-Peng Wang ,Jonathan Trinastic , Thomas Schulthess , , and Hai-Ping Cheng † Department of Physics and The Quantum Theory Project , University of Florida , Gainesville FL 32611 , USA Institute for Theoretical Physics , ETH Zurich , , Switzerland and Swiss National Supercomputing Center , ETH Zurich , , Switzerland
We report complex band structure (CBS) calculations for the four late transition metal monoxides, MnO, FeO,CoO and NiO, in their paramagnetic phase. The CBS is obtained from density functional theory plus dynamicalmean field theory (DMFT) calculations to take into account correlation e ff ects. The so-called β parameters,governing the exponential decay of the transmission probability in the non-resonant tunneling regime of theseoxides, are extracted from the CBS. Di ff erent model constructions are examined in the DMFT part of the cal-culation. The calculated β parameters provide theoretical estimation for the decay length in the evanescentchannel, which would be useful for tunnel junction applications of these materials. I. INTRODUCTION
Motivated by the application of transition metal oxides(TMO) in modern electronics, the charge transport throughTMO nano-junctions has been extensively investigated in thepast 20 years, theoretically and experimentally. A largeamount of literature focuses on the non-resonant tunnelingexperiments in which the tunneling current decays exponen-tially, I = I · exp ( − β L ), as the length of the tunnel junction( L ) increases. Although the β parameter depends on the inter-facial properties between the junction and the metallic elec-trodes, it is mainly determined by the electronic properties ofthe junction material itself. Since many TMOs are stronglycorrelated electronic systems, calculation of β from first prin-ciples with the inclusion of electron correlation would be nec-essary and important for understanding TMO nano-junctions.In this work we report ab initio calculations of β , based onDFT and single-site DMFT, for the late TMO monoxides.Existing studies have shown the β parameter is related to amaterial’s band gap, the hopping parameter t of the insulatingmaterial, and the alignment of the Fermi level in the metalelectrodes with the band gap of the insulating junction .One way to calculate β from first principle is to evaluate thecomplex band structure (CBS) rather than the ordinary realband structure (RBS). Complex band structure is the energyeigenvalues defined for complex values of (cid:126) k . The wavefunc-tion of a crystal structure has the well-known Bloch form of ψ = ψ e i (cid:126) k · (cid:126) r , where (cid:126) k is real. In fact, only the solutions ofSchr¨odinger’s equation with real wave vector are considered,and wavefunctions having complex wave vectors are often notconsidered because they would grow exponentially in somedirection which is physically unreasonable for periodic sys-tems. However, when dealing with surfaces and interfaces ofsemiconductors, i.e. finite or non-periodic solids, solutionswith complex wave vectors have physical meaning for ener-gies within the band gap. They represent the states exponen-tially decaying into the semiconductor, also called evanescentinterface-induced gap states. If (cid:126) k becomes a complex vari-able (cid:126) k = (cid:126) k Re + i (cid:126) k Im , then the wavefunction can be written as ψ = ψ e i (cid:126) k · (cid:126) r = ( ψ e − (cid:126) k Im · (cid:126) r ) e i (cid:126) k Re · (cid:126) r , yielding an exponential de- cay factor to the amplitude of the wavefunction. The decayis for the direction (cid:126) k Re , and corresponds to the electron non-resonant tunneling current decay I = I · exp ( − β L ) in thatdirection. The β parameter is associated with the imaginarypart, (cid:126) k Im , via the relation β = | (cid:126) k Im | . The energy bands aregeneralized to be defined on contours in the complex plane of (cid:126) k . A very useful feature of CBS is that one can directly read (cid:126) k Im thus β from the band structure without additional calcu-lations. By picking an arbitrary energy of the wavefunctionwithin the gap, one can trace sideways to the nearest complexband at that energy level and trace down to the corresponding (cid:126) k Im . In many cases, it is su ffi cient to apply the CBS approachwith the standard Kohn-Sham (KS) density functional theory(DFT). However, it could yield wrong results for materials inwhich electron correlation plays an important role. The self-energy due to correlation must be considered in such cases.There have been CBS studies based on beyond-DFT calcula-tions. For example, the GW approximation and hybrid den-sity functionals had been used to calculate β of simple organicmolecules and yielded better agreement with experimentallyknown values . In this study, we analyze the CBS and calcu-late the β of Mott insulating materials, using NiO, CoO, FeOand MnO as examples. The correlation e ff ect is taken into ac-count by carrying out DFT plus DMFT calculation and β isevaluated from the DMFT-corrected band structure.The series of 3 d transition metal monoxides with rock saltstructure present very rich physical properties. Early in the se-ries, TiO and VO are metallic materials, whereas later mem-bers, e . g . NiO, CoO, FeO and MnO, show clear insulatingproperties and antiferromagnetic (AFM) ordering below theNeel temperature T N . At room temperature, NiO is in AFMphase and the other three are in paramagnetic (PM) phase( T N = K , K , K , K for NiO, CoO, FeO andMnO, respectively) . Because most of the tunnel junctionapplications are operated at room temperature or even lowertemperatures, carrying out the calculation of wavefunction de-cay rate in the PM phase of these materials needs to be justi-fied, at least for NiO. This is supported by several experimentsas well as developments in the calculation side. It is wellknown that DFT in the local density approximation (LDA) a r X i v : . [ c ond - m a t . s t r- e l ] F e b fails to provide a band gap for these materials. Better match-ing between the experimental data and the calculated bandstructures was reached first by using the spin-resolved ver-sion of LDA, the local spin density approximation (LSDA) .Independent angle resolved photo-emission experiments in the 1990s studied the valence band property of bulk NiObelow 525K. The experiments demonstrated that the valenceband structure from a LSDA calculation with AFM orderingagreed with the experimental data better than LDA. Howeverthe band structure close to Fermi energy was still very di ff er-ent from the measured data and the calculated band gap wastoo small ( < > → PM transition of NiO ,which contradicts theories that the band gap of NiO is mainlydue to AFM ordering. On the calculation side, the develop-ment of the LDA + U method provided a much improved bandgap of 3 . , and the LSDA + U calculationsalso had success in describing the electronic structure of 3 d metal monoxides. These facts suggest that the band gap ofNiO is mainly due to electronic correlation rather than AFMordering. The decay rates of evanescent channel calculatedin the PM phase should not be significantly di ff erent than inthe magnetic ordered phase, because the β is mainly related tothe materials’ band gap and the band gap is not significantlya ff ected from the AFM → PM transition. For CoO, FeO andMnO, we do not find similar experiments studying whetherthe photo-emission spectra features change in AFM → PMtransition. All three are in PM phase at room temperature.The calculations presented here are carried out in a straight-forward way. The four materials’ ground state band struc-tures are first calculated using the full potential linearized aug-mented planewave (FP-LAPW) method. The obtained bandstructures are then used to construct e ff ective Hamiltonian inWannier orbital basis, and also used to compute Coulomb in-teraction matrices using the cRPA method . With the Hamil-tonian and the U matrices, we perform DMFT calculationsto get the k -resolved spectral functions and analyze the bandgaps. Using the DMFT self-energy, we construct the fullGreen’s function and use it to calculate the complex bandstructures (CBS) and extract the decay rate. Outline : The remainder of the article is organized as fol-lows. Section II introduces the calculation methods, includingthe DFT plus DMFT scheme and the way we obtained the β parameter from CBS. The essential step of computing theCoulomb interaction U matrices are grouped in Appendix A.The resulting spectral function, as well as the k -dependenceof the β parameter, are described in Section III. Section IVprovides the conclusion. II. METHODS
The CBS can be calculated using either wavefunctions orGreen’s functions. We used the Green’s function approachbecause it’s consistent with the DMFT formalism. The self-energy from DFT + DMFT is used to construct the full Green’sfunction, which is then used to evaluate the CBS and β . We will first describe our DFT and DMFT calculation, then ex-plain how we calculate CBS from Green’s function, and howwe find β from CBS.The four monoxides, especially NiO, have been extensivelystudied in the DMFT community and used as benchmark ma-terial for novel computational methods . The existingDFT + DMFT calculations of NiO were not done in the ex-actly same way. One di ff erence in our calculation is the use ofcRPA method to calculate the U matrices in the same Wannierorbital basis used for the Hamiltonian construction. Thus thehopping and interaction parameters of the e ff ective Hubbardmodel are consistently built from the same DFT ground state. II.A. DFT calculation
Our DFT calculation is done using the FP-LAPW method,as implemented in a modified version of the ELK code . Theground state is calculated within the generalized gradient ap-proximation (GGA) using the PBE functional. The mu ffi n tinsphere radii are, for example, 2.02 a and 1.72 a , for Niand O, respectively. The experimental values of lattice con-stants are used . A dense k-point grid of 16x16x16 wasused to perform Brillouin zone integration. Figure 1 displaysthe ground state band structure of NiO and orbital characters(the amount of overlapping between Bloch states and atomicorbital states). -8 -6-4-2 0 2 4 6 8K G X W L G E ne r g y ( e V ) DFT bands p-weight d-weight s-weight
FIG. 1: Non-magnetic NiO ground state band structure (solidlines), and band characters (open circles) calculated by pro-jecting Bloch states onto atomic orbital states. The radius ofthe open circles are proportional to the weight of the atomicstates. Fermi level is at zero.It is clearly seen that there are five d -like bands aroundthe Fermi level, representing the partially filled d states ofthe transition metal atom and giving the material a metal-lic state. Below them in the [ − . , − .
0] eV range are threebands showing p -orbital character. Above the d -like bands, inthe [ + . , + .
0] eV range, there is a single band of transitionmetal (TM) 4 s character. It is a common feature of the fourmaterials that the d -like bands and p -like bands are separatedby a small gap. The group of p -like and d -like bands are iso-lated from lower bands, but are very close to the s -like bandat the Γ point. Within the GGA-PBE calculation, as shown inFigure 2, we find that the s -like band is slightly gapped fromthe d -like bands in the cases of NiO and CoO but is over-lapping in energy with the d -like bands for FeO and MnO.Though not shown, we also find that, when using the LDAfunctional, the s -like band has overall more overlap with the d -like bands for these four materials. (a) (b) (c) (d) FIG. 2: Non-magnetic ground state band structures (solidcurves) of NiO, CoO, FeO and MnO. And the TM 4 s orbitalcharacter only (open circles). Horizontal solid lines are placedat the maximum of d -like bands, to make the separation oroverlap of s -like and d -like bands clearly seen. Fermi level isalways at zero. -8-6-4-2 0 2 4 6 8K G X W L G E ne r g y ( e V ) DFT bands Wannier orbital bands
FIG. 3: DFT band structure of NiO (solid line), and re-constructed bands (dots) in symmetry-preserving Wannier or-bital basis, which are identical to DFT bands by construction.Fermi level is at zero. The orbitals’ characters in Fig.1 and Fig.2 display a clear d - p mixing in these materials, which motivates our model con-struction explained in next section. The s orbital weight iswell located in the singe band in the [ + . , + .
0] eV range.We do not observe significant mixing between s and the groupof p and d . We will keep the s -like band in the analysisbecause it had been demonstrated in existing DFT + DMFTstudies that the TM s -like band has significant contribu-tion to the photo emission spectrum of these monoxides. Inaddition, as we will see in later sections, the extension of the s -like band in the complex domain goes across the gap region;thus it should be included for a correct determination of Fermilevel pinning position. So, in order to construct localized or-bital basis for all later calculations, we downfold the Blochbands in the energy window [ − . , + .
0] eV to the symmetry-preserving Wannier orbital basis that includes the TM d -likebands, the TM s -like band and the Oxygen p -like bands. Thereconstructed bands of NiO are shown in Fig. 3. II.B. DMFT calculation
Dynamical mean field theory is one successful way to moreaccurately capture the electronic correlation e ff ect and rem-edy the failure of DFT. Application of DMFT to TMOs orig-inated from the work of Peierls and Mott . Usually anappropriate correlation subspace was identified as those elec-tron states in the partially-filled, transition-metal d shell, andwas associated with interactions including the on-site intra- d and inter- d interactions. During the past two decades, theDMFT method has been developed for the low energy e ff ec-tive Hubbard model constructed for real materials’ d -like or f -like bands. The widely-adopted numerical scheme involvesselecting DFT bands near the Fermi energy as the correlationsubspace and fitting them to a tight-binding model using thedownfolding technique applied to localized orbitals, such asWannier orbitals . For each (cid:126) k point, the Bloch Hamilto-nian is downfolded to the Wannier orbital basis , H Wann ( (cid:126) k ).Through the Fourier transformation, H Wann ( (cid:126) k ) serves as sin-gle particle hopping t i j in the first term of Eq.(1) below. ThisHamiltonian contains contributions from the e ff ective poten-tial of the DFT calculation that also creates a double-countingissue, which is explicitly accounted for by a correction withinDMFT.The multi-orbital Hubbard model Hamiltonian with on-site Coulomb interaction can be expressed within the secondquantization framework as ˆ H Hubbard = ˆ H Kinetic + ˆ H Coulomb = (cid:88) i , j ,µ,ν t dpi j ,µν ˆ c + i ˆ c j + (cid:88) i ,α,β,γ,δ U di ,αβγδ ˆ c + i ,α ˆ c + i ,β ˆ c i ,γ ˆ c i ,δ (1)Here the indices i , j are site indices, and µ , ν are orbital indicesincluding spin for all orbitals within the correlation subspace.The indices { α, β, γ, δ } are a subset of { µ , ν } to indicate thoseorbitals associated with on-site Coulomb interactions. For thefour monoxides in this study, due to the d - p mixing men-tioned in Sec.II.A, the correlation subspace includes both d and p orbitals and the subset { α, β, γ, δ } of interacting orbitalsis limited to d only. The Coulomb interaction tensor U di ,αβγδ can be computed from first-principle or sometimes used asempirical parameters of the model. With the hopping andCoulomb interaction parameters at hand, the DMFT methoditeratively solves the model by mapping it to an e ff ective An-derson single-impurity model (ASIM). The impurity Green’sfunction is often expressed in the path integral formulation,with integration over Grassmann fields of second quantizationcreation and annihilation operators, ˆ c + and ˆ c , G imp ( i , τ ; i , τ ) = − (cid:82) D [ˆ c + ] D [ˆ c ] e − S [ˆ c + , ˆ c ] { ˆ c ( τ )ˆ c + ( τ ) } (cid:82) D [ˆ c + ] D [ˆ c ] e − S [ˆ c + , ˆ c ] (2)In Eq.(2), D [ . ] is the standard integration measure. S [ˆ c + , ˆ c ] isthe e ff ective action as defined in Eq.(3) below for the impurity. S [ˆ c + , ˆ c ] = − (cid:90) β d τ (cid:90) β d τ (cid:48) (cid:88) i , j ˆ c + i ( τ ) G − , i j ( τ − τ (cid:48) ))ˆ c j ( τ (cid:48) ) + (cid:90) β d τ ˆ H Coulomb (ˆ c + , ˆ c ) (3)In Equations (2) and (3), i and τ are the site index and imagi-nary time. G , i j is the bare propagator, which is also called the bath Green’s function. It plays a similar role as the Weiss fieldin classical mean-field theory. Specifically, it describes an ef-fective field coupled to the impurity that contains all non-localinformation of the underlying lattice, and the lattice is consid-ered as a reservoir of non-interacting electrons. The di ff erencefrom the classical Weiss field arises in its time dependence,which accounts for local dynamics. We refer readers to Ref. for explicit definitions of D [ . ] and S [ˆ c + , ˆ c ].Given the e ff ective action, there exists several well estab-lished numerical methods to solve for the impurity’s Green’sfunction. The family of quantum Monte Carlo (QMC) solversare widely accepted and numerically exact if the simulationtime is su ffi ciently long. We refer readers to Ref. for tech-nical details about general QMC impurity solver. In thiswork, we are using the continuous time hybridization ex-pansion (CT-HYB) QMC solver implemented in the DCA ++ code . The solver adopts the segment picture to take intoaccount density-density interactions. The coupling to the bath is diagonal only in orbital space. Our calculations are per-formed at inverse temperature 1 / kT =
20. The number ofMonte Carlo sweeps in the QMC calculation is 10 in eachsolver run. The continuous-pole-expansion method is usedfor obtaining the self-energy and impurity Green’s function inreal frequency domain.For material-specific calculations, the lattice Green’s func-tion is constructed, within the correlation subspace, from thedownfolded Hamiltonian: G ( i ω n ) = N k (cid:88) (cid:126) k i ω n + µ ) − H Wann ( (cid:126) k ) − Σ ( i ω n ) (4)where N k is the number of (cid:126) k points, ω n is the Matsubara fre-quency and µ the chemical potential. The self-energy Σ ( i ω n ) is supplied with an initial guess, then updated in each of theDMFT iteration. One uses the Dyson’s equation in each iter-ation to derive the bath Green’s function and the e ff ective im-purity problem numerically, i . e . G − ( i ω n ) = G ( i ω n ) − +Σ ( i ω n ),and solve the impurity problem. For the late TM monox-ides with strong d - p mixing, we include five d and three p orbitals in the correlation window, while limiting interac-tions to d orbitals. Thus the Hamiltonian H Wann ( (cid:126) k ) and lat-tice Green’s function G ( i ω n ) in Eq. (4) are eight-dimensional. Σ ( i ω n ) is always five-dimensional, and is added to the d -blockof H Wann ( (cid:126) k ). When constructing the bath G − ( i ω n ) for inter-acting orbitals, we use the d -block of G ( i ω n ) together with Σ ( i ω n ) in Dyson’s equation.The value of interaction parameters in ˆ H Coulomb (ˆ c + , ˆ c ) arecalculated using ab initio methods from the materials’ DFTground states. The constrained Random Phase Approxima-tion (cRPA) method, as explained in details in Appendix A, isadopted for this step. The important step in cRPA is to choosea screening window, within which the particle-hole polariza-tion are excluded. The d - p mixing gives some arbitrarinesshere because one cannot find a window of bands that includesexactly all d -weight and exclude all p -weight. There are nat-urally two choices: excluding both d - like and p - like bands in[-8.0, + d p model; exclud-ing only the five d - like bands which is called the d - d p model.We have calculated the on-site Coulomb interactions of thefive Wannier d orbitals for the two models. The results arediscussed in Appendix A. The later complex band analysis isbuild on the d - d p model only. We present the cRPA resultsof both models in this work for the partial purpose of bench-marking the current cRPA implementation.The DFT + DMFT calculation scheme and its variants havebeen widely used in the past two decades to study TMOs thathave a pronounced correlation e ff ect. It is worth briefly re-viewing the existing studies and pointing out the di ff erencesand limitations in the current work. One of the earliest appli-cations using DFT + DMFT for real materials was a study ofNiO , where a realistic gap and the near-gap spectra were ob-tained for a correlation subspace that contains d orbitals only.Shortly after, the oxygen p states of NiO were included ,in a way similar to that described in this section, to providefuller description of the valence-band spectrum. It was foundin these studies that doping holes leads to the filling of the cor-relation gap and a significant transfer of the d spectral weight.In these studies the low-energy Zhang-Rice bands were alsoobtained for the first time. Besides the paramagnetic state,magnetic state properties of NiO were also investigated withinthe framework of DFT + DMFT , where the Iterated Per-turbation Theory (IPT) solver and the numerical exact diago-nalization (ED) solver were used to solve the impurity prob-lem. NiO has also been actively used as a benchmark ma-terial for DFT + DMFT method development, e.g. new meth-ods related to the double counting correction and an-alytical continuation . The other members of the late TMmonoxides, CoO, FeO and MnO, together with NiO havebeen studied within the DFT + DMFT scheme in a system-atic investigation of the band gaps of these materials , in-vestigation of the fundamental quantum entanglement of in-distinguishable particles , and as benchmarking materials inmethod developments for determining the Coulomb correla-tion strength .Besides studies under ambient conditions, there are a sig-nificant number of first-principle calculations focusing onthe beyond-equilibrium properties of the late-TM monox-ides, particularly the changes of electronic structure relatedto high pressure and lattice distortions. It was first reported that MnO experiences a simultaneous moment collapse, vol-ume collapse, and metallization transition under a pressureof about 100 GPa. Upon compression of 60 −
70 GPa, theB1 structure of FeO has a spin-state transition accompaniedby an orbital-selective Mott metal-insulator transition , ingood agreement with the experimental result . A pressure-driven orbital selective insulator-to-metal transition is also ob-served in CoO . Similar to what is seen in FeO, the t g orbitals of Co become metallic first at about 60 GPa, and the e g orbitals remain insulating until the much higher pressure ofabout 170 GPa. It is found that the transition to fully metallicstate is driven by a high-spin to low-spin transition of the Co + ions. A systematic study of all four TM monoxides under highpressure reveals a remarkably high pressure of 430 GPa forthe insulator-metal transition in NiO, which is well out of therange 170 −
40 GPa of MnO–CoO.The full charge-density self consistent (CSC) DFT + DMFTscheme is often used in TM oxide calculations under highpressure because the charge density is subject to change withlattice distortion. In the CSC scheme , the DMFT iterationdescribed in this section is nested in an outer iteration of thecharge density. The many-body e ff ect within the correlationsubspace is self-consistently included in the entire system.It has been demonstrated that the CSC scheme is necessaryin studying the metal-insulator transition of V O , where astrong enhancement of the a g − e π g crystal-field splitting causesa substantial re-distribution of charge density and thereby in-fluences the lattice structure due to electron-lattice coupling.In a recent study of pressure-induced insulator-metal transi-tion in Fe O , a site-selective redistribution of the Fe 3 d charges between the t g and e g orbitals associated with spinstate transition was captured within a CSC DFT + DMFT cal-culation. The CSC scheme might be important for the cur-rent study mainly because the 4 s -like band enters the cor-relation window at the Γ -point. If the 4 s -like band is sig-nificantly shifted in a CSC DFT + DMFT calculation then itscomplex extension would be shifted too and a ff ect the com-plex band structure within the Mott gap. Indeed, in exist-ing CSC DFT + DMFT studies of late TM monoxides underpressure , the 4 s -like band is significantly lowered (byabout 3 eV ) to become much closer to the Fermi energy. How-ever, under ambient conditions, the 4 s -like band is not signifi-cantly moved in CSC DFT + DMFT calculations, which makessense because of the minimum hybridization between the 4 s -like band and the group of p -like and d -like bands (Fig. 1).Given the fact that a non-CSC scheme was successfully ap-plied in many studies of the late TM monoxides , wecarry out the calculations in the non-CSC scheme even whenthe 4 s -like bands of FeO and MnO enter the correlation win-dow (Fig.2). Though the 4 s -like band does not take part in the DMFT iteration, it is used in constructing the final latticeGreen’s function for complex band analysis. II.C. Evaluation of CBS and Decay Rate
As mentioned in the Introduction, the wavefunction decayrate in direction (cid:126) k can be estimated by supplying an imagi-nary part to it, i.e. (cid:126) k = (cid:126) k Re + i (cid:126) k Im , and studying the complexband structure (CBS). Here, (cid:126) k Re is the decay direction, and (cid:126) k Im yields β . CBS is always defined on the complex plane of (cid:126) k where (cid:126) k Re lies on the real axis. A grid or path of real (cid:126) k points inthe sense of traditional Brillouin zone sampling in DFT calcu-lations is not defined here. Contrary to the normal procedureof solving for eigenenergies after specifying real (cid:126) k points, oneneeds to first specify the value of the eigenenergy, then searchthe complex plane of (cid:126) k at that eigenenergy for poles of thefull Green’s function to locate the (cid:126) k (s) for that eigenenergy. Inpractice, it’s convenient to express any (cid:126) k in the 1st Brillouinzone as: (cid:126) k = C · ˆ k ⊥ + ( C · ˆ k + C · ˆ k ) ≡ (cid:126) k ⊥ + (cid:126) k (cid:107) (5)where the real unit vector ˆ k ⊥ is always the decay direction,and C is the complex coe ffi cient that defines the complexplane of searching for poles, i . e . ( C · ˆ k ⊥ ) ≡ (cid:126) k ⊥ . Here, thequantity ( Re [ C ] · ˆ k ⊥ ) is same as the (cid:126) k Re used at the beginningof this section. And, ˆ k and ˆ k are user-defined real unit vec-tors in the plane perpendicular to ˆ k ⊥ . C and C are both real . C · ˆ k + C · ˆ k ≡ (cid:126) k (cid:107) defines the parallel component of (cid:126) k .The DMFT self-energy in real frequency domain is usedto construct the lattice Green’s function of the subspace con-taining s , p and d orbitals (9-dimensional). The full Green’sfunction is expressed in the usual way: G ( (cid:126) k , ω ) = ω + µ − ˜ H Wann ( (cid:126) k ) − Σ DMFT ( ω ) (6)In Eq.(6), µ is the Fermi energy, Σ ( ω ) is the convergedDMFT self-energy after analytical continuation. The original H Wann ( (cid:126) k ) is Fourier transformed to real space hopping t i j ( (cid:126) R ),then Fourier transformed back to ˜ H Wann ( (cid:126) k ) at any (cid:126) k , real orcomplex, as defined in Eq.(5). For given ˆ k ⊥ and (cid:126) k (cid:107) , ˜ H Wann ( (cid:126) k )can be considered as a function of C only. Thus the Green’sfunction is a function of C and ω , i . e . G ( C , ω ). Fig.4 showsan example case of the complex plane defined by G ( C , ω ).In the example ˆ k ⊥ is chosen to be the unit reciprocal latticevector (cid:126) b / | (cid:126) b | and (cid:126) k (cid:107) =
0. The poles of G ( C , ω ) are resolvedby finding roots of the equation: det | G − ( C , ω ) | =
0. The setof roots for di ff erent values of ω gives the complex bands.In the following work, we study the complex band and decayrate in three directions: (a) ˆ k ⊥ = ( k x , k y , k z ) = (1 , , k ⊥ = (1 , , / √
2; (c) ˆ k ⊥ = (1 , , / √
3, where k x , k y , k z areCartesian coordinates in k -space.FIG. 4: Complex plane of C for searching poles of G ( C , ω ),for the direction ˆ k ⊥ = (cid:126) b / | (cid:126) b | . Shaded area is the searchingarea, which includes the real axis and boundaries of the 1stBrillouin zone. III. RESULTS AND DISCUSSION
The main purpose is to study the features of the CBS whenthe correlation e ff ect is included. In this section, we first re-port the band gaps and spectral functions from DFT + DMFTand compare it to experimental data and existing calculations.Then the CBS and β parameters are studied at the so-calledcharge neutrality level within the band gap, which is furtherrelated to the Fermi level pining position and the height ofSchottky barrier in tunnel junction applications. III.A. Band Gaps from DFT + DMFT
The resulting density of states (DOS) from the DMFT cal-culations are presented in Fig.5 for the d p and d - d p mod-els. The band gaps are measured between the widths at halfheight of the conduction band and valence band in DOS. Ta-ble I has the measured band gaps. Here we briefly compare tothe experiments. NiO and MnO have been extensively stud-ied in experiments. The band gap of NiO was determinedto be 3 . . . The experimental values of MnO gapis in the range of 3 . . . The present calculation ofthe d - d p model are in agreement with experiments for thesetwo materials. Experimental values of the band gap of CoOhas diverse values. Some experiments reported a band gapof 2 . . . Some other studies found higher values, e.g.5 . and indirectgap of 2 . . The value ofour calculated band gap based on the d - d p model falls in therange and is close to the conductivity experiment . Quasi-particle calculations using the hybridfunctional and the G W method yields similar value (3 . ofCoO. It seems the absorption experiments underestimates thegap compared to DMFT results. Unfortunately, there are verylimited experiments reporting the measured band gaps of FeOin PM phase. This is because the preparation of a pure FeO sample is di ffi cult due to Fe segregation . Thus, compari-son of theoretical calculation with experimental spectra is veryrare. The only reported experimental estimate that we foundis 2 . . We arenot sure if this is an underestimated value.FIG. 5: Total(black), t g (blue), e g (green) and p (red) spectralfunctions, A ( ω ), of the two models of NiO, CoO, FeO andMnO. Band Gap (eV) NiO CoO FeO MnODFT + DMFT, dp model + DMFT, d - dp model / a 3.8Exp(XAS-XES) 4.0 2.6 n / a 4.1Exp(PES-BIS) 4.5 2.5 n / a 3.9Exp(absorption) 4.0 2.8 2.4 3.6-3.8 TABLE I: The band gaps measured from DFT + DMFT den-sity of states. Sources of the experimental gaps are in text ofSec.III.A.We noticed our band gap result of the d - d p model of NiOis in general agreement with existing DFT + DMFT calcula-tions. Though the existing calculations are not exactly same,most of them are performed in the PM state and are not inCSC scheme. In Ref. a similar model construction wasconsidered and U = . J involved) calculated fromconstrained-LDA was used. A band gap of about 4 eV wasobtained. In Ref. , only the Ni- d orbitals were taken into ac-count (which should be called d - d model following the nam-ing convention used here) and U = . J = . . had used NiO to test a new doublecounting method, where they used the same U = . J = . d and p orbitals. The calcula-tion was carried at high T of 2300K. They found a band gapof about 4 . p orbital peak below Fermienergy depend on the double counting potential. (a) NiO (b) CoO(c) FeO (d) MnO FIG. 6: (cid:126) k -resolved spectral functions of the d - d p model, ofNiO, CoO, FeO and MnO.In our calculation, the correct representation of the U ma-trices is important to yield accurate band gaps for the fourmaterials. We have found results from the d - d p model arein better agreement with existing experiments than the d p model, which seems to overestimate the band gap. However,the preference of the d - d p model over the d p model shouldbe more carefully supported by taking into consideration ofother important factors. For example, it has been proposedthat the current cRPA method could be improved by includ-ing the Pauli exclusion principle in the formalism, and overallthe e ff ect would be a reduction of the interaction strength .There are also other factors within the DMFT, e.g. di ff erentdouble counting methods, that a ff ect the resulting band gapand spectral function. For example the double counting meth-ods introduced in Ref. and Ref. . Those factors are worthyof dedicated studies and are outside the scope of the current work. Thus we limit our discussion to be within the original d - d p cRPA scheme and with only fully localized limit (FLL)double counting in the DMFT part. The d p cRPA schemeis not used in later decay rate analysis, but exists in the Ap-pendix for purpose of benchmarking the cRPA implementa-tion. We would like to emphasize that, for consistency, oneset of projected Wannier orbital basis functions is used in boththe downfolded Hamiltonian H ( k ) and the cRPA-calculatedCoulomb U matrices, which is di ff erent from existing calcu-lations of these materials.In order to include the e ff ect of the TM 4 s band in the latercomplex band analysis, we prepare the full Green’s function,Eq. (6), in the Wannier orbital basis containing d , p , and s or-bitals. The DMFT self energy is associated with the d orbitals.The k -resolved spectral function A ( (cid:126) k , ω ) = ( − /π ) (cid:61) [ G ( (cid:126) k , ω )]of the d - d p model is shown in Fig. 6. The rest of the calcula-tions are based on this Green’s function. III.B. Complex band structure including DMFT self-energy
The complex band structure (CBS) and real band structure(RBS) are obtained from resolving poles of the Green’s func-tion, as explained in Section II. In this section we analyze theobtained complex band structure and argue about the pinningposition of the Fermi level for tunnel junction applications andcalculating the β parameter at that energy level.In a general tunnel junction setup, where the insulating ma-terial is connected to metal leads on both sides, the Fermi lev-els of the two materials are brought into coincidence. At theinterface, there are metal-induced gap states (MIGS) in theinsulator gap region decaying exponentially into the material,which are Bloch states of the insulator with complex wavevector. For semiconductors, this forms a continuum of statesaround the Fermi energy ( E F ) within the gap. The gap statescontinuously change (as function of energy) from valence- toconduction-band character, appearing as an arch-shaped com-plex band connecting the real valence band with the real con-duction band. The idea of local charge neutrality , proposedby Terso ff , says: when the density of MIGS is reasonablylarge, it is necessary to occupy those MIGS which are pri-marily valence-band character, while leaving those of mainlyconduction-band character empty. Therefore, E F should bepinned at or near the energy where the gap states cross overfrom valence- to conduction-band character. If the complexband within the gap is a smooth curve, the crossing over isnaturally found at dE / dk → ∞ , which is called the branchpoint , and the corresponding energy is called the charge neu-trality level E B .Fig.7 displays the CBS and RBS for ˆ k ⊥ = (cid:126) b / | (cid:126) b | and (cid:126) k (cid:107) =
0, where (cid:126) b is a reciprocal lattice vector. The CBSpart, ( Re [ C ] , Im [ C ]) = (0 . , − . → (0 . , .
0) in the lefthalf of Fig.7, contains important information for determin-ing the charge neutrality level E B and wavefunction decayrate β . One major observation is that the gap states transi-tioning from below to above the Mott gap are not continu-ous. The real conduction bands and the top of valence bandsare both primarily d -orbitals, as seen in the DOS. The exten-sions of these bands into the complex sector cross each otherrather than connect smoothly, which is di ff erent from typi-cal semiconductors. For NiO in Fig.7, the crossing happens atabout 3 . + DMFT Fermi level ( E DMFTF ) and C = (0 . , − . E B should be pinned at or near the crossing point,where dE / dk has a finite jump rather than → ∞ . This is moreclearly shown in Fig.8. We call this crossing point the Mottbranch point, and the corresponding energy level E Mott B . Theconduction band obtained from DFT + DMFT is very flat. Itsextension in the complex domain is almost a straight line ofsmall slope, while the complex band originated from valencebands goes up steeply from below the gap. This feature leadsto the E Mott B being very close to the conduction band minimum E cond.min (which is found at Γ in our case). The small di ff erencebetween E Mott B and E cond.min gives a small Schottky barrier height Φ B of the tunnel junction, i.e. Φ B = E cond.min − E Mott B . Our cal-culation implies NiO is a junction material with large gap andsmall barrier height (comparing with the 1 . . . . Φ NiOB ≈ .
35 eV from the calculatedCBS of the d - d p model of NiO, as indicated in Fig.8. There isnot much reported work on tunnel junction experiments usinglate transition metal monoxides. We found the most relevantwork was done on the Ni-NiO-Ni tunnel junction , wherethe system was found to be a very-low-barrier system with Φ B ≈ . .
22 eVat 4 K to 0 .
19 eV at 295 K. This is in qualitative agreementwith our calculation. Ni was used as metal leads for easy man-ufacture. It is possible to tune the composition of the metaland oxide so as to raise the barrier height slightly.The 4 s band plays an interesting role here. By follow-ing the path: ( Re [ C ] , Im [ C ]) = (0 . , . → (0 . , . → (0 . , − . s band goes from high tolow energy at the Γ point, which is a common feature in Fig.7.NiO presents a di ff erent feature in the complex sector: it firstdecreases to about the Fermi energy and then starts to increaseto higher energy, while the 4 s bands of the other three mate-rials decrease in energy monotonically. In the case of NiO inFig.8, the complex 4 s band crosses the complex extension ofvalence d band at a point very close to E DMFTF , and the real 4 s band has E cond.min ≈ Γ . By using the 4 s band instead ofthe conduction 3 d band, one keeps the Fermi level pinned veryclose to E DMFTF ( i.e. there would be nearly no shift in a tun-nel junction), and Φ B ≈ s band is not responsiblefor correctly determining the Schottky barrier height also be-cause the orbital character should not suddenly change (from3 d to 4 s ) at the branch point. The correct E Mott B in Fig.8 islocated by following the same d orbital character. The sameargument applies to CoO, where the complex 4 s band is justtouching the complex valence d band at a point close to thegap bottom. CoO should also be a very small barrier junc-tion material. The 4 s bands of FeO and MnO fall into an areawhere no other complex band is found. We can expect themto be small-barrier materials based on only complex d bands. FIG. 7: Complex band structures of the d - d p model. The de-cay direction being studied is ˆ k ⊥ = (cid:126) b / | (cid:126) b | , with C = C = Re [ C ] , Im [ C ]) = (0 , − . → (0 , → (0 . , + DMFT Fermi level, E DMFTF , is at zero. -4-2024(0.0,-0.4) (0.0,-0.2) (0.0,0.0) (0.2,0.0) (0.4,0.0) E DMFTF E MottB E d-pB E cond.min Φ B = 0.35 eV NiO
FIG. 8: The CBS of the d - d p model of NiO. x -axis is the pathof ( Re [ C ] , Im [ C ]). E DMFTF is at zero (eV). E Mott B is locatedat the crossing of the imaginary extension of conduction band(small slope arrow) and the imaginary extension of valenceband (large slope arrow). E DMFTF is at zero. E d − pB is foundwithin the imaginary band connecting d - and p -valence bandswhere dE / dk → ∞ . The conduction band minimum E cond.min (short horizontal dash line) is found at Γ at a slightly higherenergy than E Mott B . The Schottky barrier height is measured tobe: Φ B = E cond.min − E Mott B ≈ .
35 eV . It is worthwhile to mention that we do observe a continuoustransition of ”gap” states below E DMFTF , where real p -orbitalbands and d -orbital bands are connected by continuous com-plex bands. This is also shown in Fig.8. The point where dE / dk → ∞ can be located, and we call it the d - p branchpoint. The corresponding energy level is E d − pB . The featureactually exists in the DFT calculation alone, despite the factthat it yields a metallic ground state. The feature largely re-mains after the DMFT calculation. The energy range of thearch-shaped complex band containing E d − pB may decrease dueto the gap opened by DMFT, as seen in the CBS of CoO, FeOand MnO in Fig.7.Numerically some regions of the complex and real bandssmear out, for example near Γ in energy range of [ − , −
3] eVbelow E DMFTF . Such smearing is explained within the quasi-particle picture in which spectral weights have spreads, asseen in the (cid:126) k -resolved spectral functions in Fig.6.In the rest of this section, we turn to the calculation of β at E Mott B for the d - d p model of the four materials and study thedirection dependence of β .The values of the decay rate β should be anisotropicfor crystal structures. As already mentioned in the end ofSec.II.C, we perform the calculation of β in three decay di-rections: (a) ˆ k ⊥ = (1 , , k ⊥ = (1 , , / √ k ⊥ = (1 , , / √
3. The values E Mott B and E d − pB are determinedfor each direction. We found E Mott B , or E d − pB , stays the samefor the three di ff erent decay directions. If the wavefunction propagates perfectly along one direction, then we get the de-cay rate by directly reading o ff Im [ C ] (thus β ≡ | · Im [ C ] | )at the branch point. This has been done for the three direc-tions and for the four materials. The results are summarizedin Table II. We found the values of β at E Mott B are all within therange of [0.29,0.40]. NiO CoO FeO MnO Φ B = E cond.min − E Mott B (eV) 0.35 0.32 0.27 0.30 E Mott B (eV, from E DMFTF ) + + + + β @ E Mott B ˆ k ⊥ : (0,0,1) 0.29 0.31 0.32 0.29ˆ k ⊥ : (0,1,1) 0.33 0.34 0.35 0.33ˆ k ⊥ : (1,1,1) 0.37 0.37 0.40 0.38 E d − pB (eV, from E DMFTF ) -2.5 -2.3 -2.8 -3.0 β @ E d − pB ˆ k ⊥ : (0,0,1) 0.15 0.14 0.16 0.13ˆ k ⊥ : (0,1,1) 0.17 0.16 0.16 0.13ˆ k ⊥ : (1,1,1) 0.18 0.16 0.17 0.14 TABLE II: The calculated values of Φ B and locations of E Mott B and E d − pB , and the decay rates β in di ff erent directions at thecorresponding energy levels. (a) NiO (b) CoO(c) FeO (d) MnO FIG. 9: The distribution of β ( C , C ) at E Mott B for the decaydirection (1,0,0).While Table II gives us an idea of the decay rates at spe-cific single directions, in reality the wavefunctions may notbe perfectly propagating along a given direction. The non-perpendicular incident components can be taken into accountby simply allowing C and C (defined in Eq.(5)) to vary, re-sulting in a distribution of β ( C , C ) for a certain decay direc-tion ˆ k ⊥ . This means we assume the tunnelling barrier system0has translational symmetry in the plane parallel to the inter-face, so that the transmission conserves (cid:126) k (cid:107) . At E Mott B , we haveobtained the β ( C , C ) as shown in Fig.9, Fig.10, and Fig.11. (a) NiO (b) CoO(c) FeO (d) MnO FIG. 10: The distribution of β ( C , C ) at E Mott B for the decaydirection (1,1,0). (a) NiO (b) CoO(c) FeO (d) MnO FIG. 11: The distribution of β ( C , C ) at E Mott B for the decaydirection (1,1,1).At first glance, the symmetry of these distributions is con-sistent with the crystal symmetry, which makes sense becausesingle-site DMFT assumes a completely localized self-energy.The ranges of values of β are quite di ff erent in di ff erent direc-tions. The observation becomes clearer when one introducesthe β -resolved density of states n ( β ). Ideally n ( β ) d β would be the number of ”states” with β values in the infinitesimal in-terval ( β, β + d β ). Numerically, we calculated β on the grid of C × C = ×
80 points within [0,1]x[0,1], and applied linearinterpolation to make the grid denser. The resulting density isshown in Fig.12. It is clear that the decay rates cluster around0.3 for the (1,1,1) direction, and extend to larger and largervalues in (0,1,1) and (0,0,1) directions.FIG. 12: β -resolved density of state n ( β ) for the three direc-tions: (0,0,1) upper, (0,1,1) middle, (1,1,1) down. The β ismeasured at E Mott B . The y -axis scale is same for the three plots.In general, the probability for wavefunction transmissionthrough the tunnel junction can be written as the product of theprobability for transmission across each of the interfaces timesa factor that describes the exponential decay of the electronprobability within the junction material . T B ( (cid:126) k (cid:107) ) = T L ( (cid:126) k (cid:107) ) · T R ( (cid:126) k (cid:107) ) · T evantot ( d ) (7)Here T B is the transmission probability through the junction. T L and T R are the probabilities for an electron to be trans-mitted across the left and right electrode barrier interfaces re-spectively, and T evantot is the evanescent channel contributionto the conductance. We use a simplified model to estimatethe evanescent channel transmission probability T evantot . Thetransmission probability for each (cid:126) k (cid:107) takes the form: T evan (cid:126) k (cid:107) ( d ) = T · exp ( − β (cid:126) k (cid:107) · d ), where d is the thickness (number of layers)of the tunneling barrier and β (cid:126) k (cid:107) is the distribution of β ( C , C )obtained in Fig.9, Fig.10, and Fig.11. The total evanescenttransmission probability for the direction ˆ k ⊥ , at a given en-ergy level, is given by: T evantot ( d ) = N (cid:126) k (cid:107) (cid:88) (cid:126) k (cid:107) T evan (cid:126) k (cid:107) ( d ) (8)We calculate the relative evanescent transmission probabil-ity T evantot ( d ) / T evantot ( d =
0) at E Mott B , as shown in Fig.13. Thecommon feature of exponential decay is clear, and the relativetransmission becomes very small after about 10 layers of theunit cell, for all four materials. The larger values of β in the1(0,0,1) direction results in a slightly faster decay than for theother two directions. The di ff erences in T evantot ( d ) between thethree directions is overall not significant after averaging (cid:126) k (cid:107) . (a) (0,0,1) (b) (0,0,1), log scale(c) (0,1,1) (d) (0,1,1), log scale(e) (1,1,1) (f) (1,1,1), log scale FIG. 13: Relative transmission probability on linear and logscale for the directions: (0,0,1), upper two plots; (0,1,1), mid-dle two plots; (1,1,1), lower two plots;
IV. CONCLUSION
In summary, we have performed DFT plus single-siteDMFT calculation for the four transition metal monoxides,NiO, CoO, FeO and MnO, in the non-spin-polarized phase,and have studied the complex band structures of them by in-cluding the DMFT self-energy in the Green’s function. The d and p orbitals are included in the DMFT calculation, andthe TM 4 s orbital is included in the Green’s function for com-plex band analysis. Both the d p and d - d p models have beenconsidered, and the corresponding screened Coulomb interac-tion parameters are calculated from first-principles using thecRPA method. The resulting spectral functions of the d - d p model present a clear band gap in general agreement with ex-periments, improving upon the gapless DFT ground state cal-culation.By using the full Green’s function that includes s , p , d or-bitals and the DMFT self-energy, we have observed the com-plex band structures of these Mott insulators. Motivated bythe tunnel junction application of these materials, we analyzedthe pinning position of the Fermi level by applying the chargeneutrality condition. In the complex domain, the gap statestransitioning from valence-band character to conduction-bandcharacter displays a jump in dE / dk , rather than the continu-ous transition seen in traditional semiconductors. The branchpoint is found to be at the jump point, and is very close to con- duction band minimum. The TM 4 s band is carefully studied,and we argue that the 4 s band is not responsible for determin-ing the Fermi level pinning position, which is supported byexperiment on NiO tunnel junctions. The calculated resultsare in consistent with experimental observations that NiO haslarge band gap and small Schottky barrier height. The trans-mission decay parameter, β , has been calculated for di ff erentdirections in k -space to give us insight into the evanescentchannel contribution to the conductance. We have investigated β in detail at the Mott branch point within the correlation win-dow. We found that the β parameter has very di ff erent valuesand distributions for di ff erent directions. When the decay di-rection and the incident direction are not the same, β generallybecomes larger and forms non-trivial patterns depending onthe relative direction.The CBS analysis relies on the DMFT self-energy. OurDFT + DMFT calculation based on the modified ELK codeand DCA ++ code is a single shot implementation, not a fullycharge density self-consistent one. The current work assumedthe TM 4 s band is not significantly shifted in a CSC treatment,which is reasonable based on existing studies. We noticed theexisting CSC DFT + DMFT calculations of the same materi-als have shown the necessity of updating charge density whenstructural changes due to pressure and strain are closely tiedto electronic transitions, which raises interesting direction tomotivate our future work.In conclusion, the presented work carries out a fully ab ini-tio study of the charge neutrality level and wavefunction decayrate in the evanescent channel of late transition metal monox-ides in their PM phase under ambient condition, by using thecombination of CBS method and DFT + DMFT calculation.The DFT + DMFT band gap calculations, along with the cRPAcalculations of Coulomb U matrices, are done in a standardway. The main physical features are captured. The newly ob-served feature of the CBS and location of branch point of theseMott insulators are di ff erent from band insulators. In addi-tion, the TM 4 s band of NiO is found to have di ff erent featurein the complex domain than that of the other three materials.The barrier height, values of decay rate and its direction de-pendence can all be obtained from first principle. This numer-ical study could be useful when Mott insulating materials areused for tunnel junction applications. The approach could beapplied to more complicated structures or lower dimensionalcases, as long as the application of singe-site DFT + DMFT isjustified.
V. ACKNOWLEDGEMENTS
This work was supported by the US Department of Energy(DOE), O ffi ce of Basic Energy Sciences (BES), under Con-tract No. DE-FG02-02ER45995. The computation was doneusing the utilities of the National Energy Research ScientificComputing Center (NERSC).2 VI. APPENDIX A
The on-site Coulomb interaction U-matrix is required in-put for the impurity problem. One reliable way to calculatethese parameters from first principles is the constrained Ran-dom Phase Approximation(cRPA) method, which has beenwell described in the literatures . Here we first summa-rize the original idea. Then we explain our implementationbased on the density response function, and list our results ofU-matrices for the two models and the four materials.The cRPA calculation is based on a DFT ground state calcu-lation that includes many empty bands. One aims to get an es-timation of the screened Coulomb interaction for the selectedbands of interest, or an energy window. For this purpose, theparticle-hole polarization between all possible pairs of occu-pied state and unoccupied state are taken into account. Withinthe RPA, the particle-hole polarization is calculated as : P tot ( r , r (cid:48) ; ω ) = occ . (cid:88) i unocc . (cid:88) j ψ ∗ i ( r ) ψ j ( r (cid:48) ) ψ ∗ j ( r ) ψ i ( r (cid:48) ) × ( 1 ω − ε j + ε i + i δ + ω + ε j − ε i − i δ ) (9)where ψ i and ε i are the eigenfunctions and eigenenergies ofthe one-particle Hamiltonian in DFT.The selected bands of interest are often around the Fermilevel, and have a particular orbital character, for example d -like in our case. Following the convention in the literatures,we label the bands of interest or energy window as the d -space. If both the occupied state and the unoccupied stateare within the d -space, then the polarization contributes to P d ( r , r (cid:48) ; ω ). All the other pairs of occupied and unoccupiedstates contribute to P r , where r stands for the rest of thebands. Thus, the total polarization is divided into two parts: P tot = P d + P r . P r is related to the partially screened Coulombinteraction : W r ( ω ) = [1 − ν · P r ( ω )] − · ν (10)In the above equation, ν is the bare Coulomb interaction.It’s obvious that the total polarization, P tot , screens the bareCoulomb interaction, ν , to give the fully screened interac-tion W . With the same logic, P d screens W r to give thefully screened interaction W . Thus, W r is identified as theon-site Coulomb interaction for the bands of interest, i . e . U ( ω ) ≡ W r ( ω ), which has included the screening e ff ect fromthe realistic environment of the material.Our cRPA calculation is based on the partial Kohn-Shamsusceptibility : χ KSr ( r , r (cid:48) ; ω ) = (cid:88) i , j (cid:54)⊂ C.S. ( f i − f j ) ψ ∗ i ( r ) ψ j ( r (cid:48) ) ψ ∗ j ( r ) ψ i ( r (cid:48) ) ω − ε j + ε i + i δ (11)where f i and ε i are the occupancy and energy of the eigenstate ψ i . The summation over band indices runs over all bands ex-cluding the cases where both i and j are inside the correlationsubspace ( C.S. under the summation sign means correlation subspace). In general, the density response function χ ( r , r (cid:48) ; ω )is related to the Kohn-Sham susceptibility by the following in-tegral equation : χ ( r , r (cid:48) ; ω ) = χ KS ( r , r (cid:48) ; ω ) + (cid:34) dr dr χ KS ( r , r ; ω ) × (cid:32) | r − r | + f xc ( r , r ; ω ) (cid:33) × χ ( r , r (cid:48) ; ω ) (12)where the f xc is the functional derivative of the exchange-correlation potential with respect to the charge density, whichis often neglected in the random phase approximation. Sincewe are considering periodic crystal structures, the integralequation is often written in the Fourier-transformed formwhere χ ( r , r (cid:48) ; ω ) becomes χ ( q , q (cid:48) ; ω ) with q and q (cid:48) the recip-rocal lattice vectors. Due to the invariance of the real spaceresponse function with respect to a shift by a lattice vector R : χ ( r + R , r (cid:48) + R ; ω ) = χ ( r , r (cid:48) ; ω ), the χ ( q , q (cid:48) ; ω ) is only nonezerowhen q and q (cid:48) di ff er by a reciprocal lattice vector G . One canreplace q by q + G , replace q (cid:48) by q + G (cid:48) , and restrict q tobe always within the first Brillouin zone. Thus the Fourier-transformed integral equation has the form: χ GG (cid:48) ( q , ω ) = χ KSGG (cid:48) ( q , ω ) + (cid:88) G G χ KSGG ( q , ω ) × (cid:16) v G + q δ G G + f xcG G ( q , ω ) (cid:17) × χ G G (cid:48) ( q , ω ) (13)where v G + q = π/ | G + q | is the expansion coe ffi cient of thebare Coulomb interaction. By using the partial Kohn-Shamsusceptibility χ KSr , GG (cid:48) ( q , ω ) (Fourier transform of Eq.(6)) andneglecting the f xc term in Eq.(7), we reach the equation forthe partial RPA density response function χ RPAr ( q , ω ): χ RPAr , GG (cid:48) ( q , ω ) = χ KSr , GG (cid:48) ( q , ω ) + (cid:88) G v G + q × χ KSr , GG ( q , ω ) × χ RPAr , G G (cid:48) ( q , ω ) (14)Note that the subscript r in Eq.(9) stands for the rest ofthe bands, as same as in Eq.(5). Eq.(9) is first solved for χ RPAr ( q , ω ) in the calculation. The rest calculation is basedon the linear respone theory , where the partially screenedCoulomb interaction W r is related to inverse dielectric func-tion ε − and bare Coulomb interaction ν : W r ( r , r ; ω ) = (cid:82) dr ε − ( r , r ; ω ) ν ( r , r ). The inverse dielectric function ε − isdetermined by ε − ( r , r ; ω ) = + ν · χ RPAr ( r , r ; ω ). Finally thefrequency-dependent screened Coulomb interaction is com-puted from the partial RPA density response function and thebare Coulomb interaction: W r , GG (cid:48) ( q , ω ) = v G + q δ GG (cid:48) + v G + q · χ RPAr , GG (cid:48) ( q , ω ) · v G (cid:48) + q (15)The above calculations have been implemented in theExciting-Plus code (a modified version of ELK code) ,which we used for the U matrix calculations. The Wannierorbitals are constructed by projection to preserve symmetry,and no spatial localization procedure is applied. In addition tothe DFT calculation described in Sec.II.A, 100 empty bandsare included for the cRPA calculation.3NiO, d p model: U σσ mm (cid:48) = .
69 8 .
03 7 .
97 8 .
03 9 . .
03 9 .
69 8 .
82 8 .
03 8 . .
97 8 .
82 10 .
27 8 .
82 8 . .
03 8 .
03 8 .
82 9 .
69 8 . .
10 8 .
25 8 .
19 8 .
25 10 . U σσ mm (cid:48) = .
00 7 .
21 6 .
97 7 .
21 8 . .
21 0 .
00 8 .
24 7 .
21 7 . .
97 8 .
24 0 .
00 8 .
24 7 . .
21 7 .
21 8 .
24 0 .
00 7 . .
66 7 .
39 7 .
15 7 .
39 0 . NiO, d - d p model: U σσ mm (cid:48) = .
29 5 .
67 5 .
42 5 .
67 6 . .
67 7 .
29 6 .
21 5 .
67 5 . .
42 6 .
21 7 .
38 6 .
21 5 . .
67 5 .
67 6 .
21 7 .
29 5 . .
48 5 .
69 5 .
44 5 .
69 7 . U σσ mm (cid:48) = .
00 4 .
85 4 .
44 4 .
85 6 . .
85 0 .
00 5 .
64 4 .
85 4 . .
44 5 .
64 0 .
00 5 .
64 4 . .
85 4 .
85 5 .
64 0 .
00 4 . .
04 4 .
84 4 .
47 4 .
84 0 . CoO, d p model: U σσ mm (cid:48) = .
27 7 .
70 7 .
69 7 .
70 8 . .
70 9 .
27 8 .
47 7 .
70 7 . .
69 8 .
47 9 .
92 8 .
47 7 . .
70 7 .
70 8 .
47 9 .
27 7 . .
73 7 .
95 7 .
93 7 .
95 9 . U σσ mm (cid:48) = .
00 6 .
92 6 .
73 6 .
92 8 . .
92 0 .
00 7 .
91 6 .
92 7 . .
73 7 .
91 0 .
00 7 .
91 6 . .
92 6 .
92 7 .
91 0 .
00 7 . .
30 7 .
12 6 .
93 7 .
12 0 . CoO, d - d p model: U σσ mm (cid:48) = .
47 4 .
94 4 .
71 4 .
94 5 . .
94 6 .
47 5 .
43 4 .
94 4 . .
71 5 .
44 6 .
57 5 .
44 4 . .
94 4 .
94 5 .
44 6 .
47 4 . .
68 4 .
96 4 .
74 4 .
96 6 . U σσ mm (cid:48) = .
00 4 .
16 3 .
78 4 .
16 5 . .
16 0 .
00 4 .
89 4 .
16 4 . .
78 4 .
89 0 .
00 4 .
89 3 . .
16 4 .
16 4 .
89 0 .
00 4 . .
25 4 .
15 3 .
82 4 .
15 0 . FeO, d p model: U σσ mm (cid:48) = .
91 7 .
49 7 .
53 7 .
49 8 . .
49 8 .
91 8 .
18 7 .
49 7 . .
53 8 .
18 9 .
59 8 .
18 7 . .
49 7 .
49 8 .
18 8 .
91 7 . .
40 7 .
75 7 .
79 7 .
75 9 . U σσ mm (cid:48) = .
00 6 .
76 6 .
64 6 .
76 7 . .
76 0 .
00 7 .
65 6 .
76 6 . .
64 7 .
65 0 .
00 7 .
65 6 . .
76 6 .
76 7 .
65 0 .
00 6 . .
98 6 .
98 6 .
89 6 .
98 0 . FeO, d - d p model: U σσ mm (cid:48) = .
89 4 .
50 4 .
33 4 .
50 5 . .
50 5 .
89 4 .
93 4 .
50 4 . .
33 4 .
93 6 .
02 4 .
93 4 . .
50 4 .
50 4 .
93 5 .
89 4 . .
13 4 .
53 4 .
35 4 .
53 6 . U σσ mm (cid:48) = .
00 3 .
78 3 .
47 3 .
78 4 . .
78 0 .
00 4 .
40 3 .
78 3 . .
47 4 .
40 0 .
00 4 .
40 3 . .
78 3 .
78 4 .
40 0 .
00 3 . .
71 3 .
78 3 .
52 3 .
78 0 . MnO, d p model: U σσ mm (cid:48) = .
68 7 .
34 7 .
39 7 .
34 8 . .
34 8 .
68 7 .
99 7 .
34 7 . .
39 7 .
99 9 .
34 7 .
99 7 . .
34 7 .
34 7 .
99 8 .
68 7 . .
18 7 .
59 7 .
65 7 .
59 9 . U σσ mm (cid:48) = .
00 6 .
65 6 .
56 6 .
65 7 . .
65 0 .
00 7 .
47 6 .
65 6 . .
56 7 .
47 0 .
00 7 .
47 6 . .
65 6 .
65 7 .
47 0 .
00 6 . .
78 6 .
87 6 .
80 6 .
87 0 . MnO, d - d p model: U σσ mm (cid:48) = .
59 4 .
28 4 .
13 4 .
28 4 . .
28 5 .
59 4 .
68 4 .
28 4 . .
13 4 .
68 5 .
73 4 .
68 4 . .
28 4 .
28 4 .
68 5 .
59 4 . .
86 4 .
31 4 .
17 4 .
31 5 . U σσ mm (cid:48) = .
00 3 .
60 3 .
33 3 .
60 4 . .
60 0 .
00 4 .
18 3 .
60 3 . .
33 4 .
18 0 .
00 4 .
18 3 . .
60 3 .
60 4 .
18 0 .
00 3 . .
46 3 .
61 3 .
38 3 .
61 0 . The cRPA calculation directly provides the intra-, inter-orbital and exchange interaction parameters. We derived sin-gle value intra-, inter-orbital and exchange interaction param-eters U , U (cid:48) and J from averaging the matrix elements, in order4to compare with the existing calculations. By definition, thediagonal elements of U σσ mm (cid:48) are the intra-orbital interactions,and the average value is: U = T rU σσ mm (cid:48) . The inter-orbitalinteractions are the o ff -diagonal elements of U σσ mm (cid:48) , thus U (cid:48) = (cid:80) m (cid:54) = m (cid:48) U σσ mm (cid:48) . The inter-orbital exchange is the o ff -diagonalelements of U σσ mm (cid:48) − U σσ mm (cid:48) , so J = (cid:80) m (cid:54) = m (cid:48) ( U σσ mm (cid:48) − U σσ mm (cid:48) ). Thevalues of U , U (cid:48) and J for the two models are summarized inTable III. We found general agreement with other cRPA cal-culation of the same materials in the literature .(eV) NiO CoO FeO MnO U , d p U , d - d p U (cid:48) , d p U (cid:48) , d - d p J , d p J , d - d p U , U (cid:48) and J deduced from thecRPA calculation. Values in parenthesis are from Ref. forthe exactly same model construction using a di ff erent codebased on Maximum Localized Wannier Functions.5 ∗ current address: Cognitive Computing and Industry Solutions,IBM Research, Saumerstrasse 4, 8803 Rueschlikon, Switzerland † Contact Email: [email protected]fl.edu; Tel: 352-392-6256 J. K. Tomfohr and O. F. Sankey, Phys. Rev. B , 245105 (2002). A. Ferretti, G. Mallia, L. Martin-Samos, G. Bussi, A. Ruini,B. Montanari, and N. M. Harrison, Phys. Rev. B , 235105(2012). R. O. Jones, Proceedings of the Physical Society , 443 (1966). W. Jauch and M. Reehuis, Phys. Rev. B , 195121 (2004). W. Jauch and M. Reehuis, Phys. Rev. B , 125111 (2002). W. Jauch and M. Reehuis, Phys. Rev. B , 184420 (2003). A. Ayuela, D. Klein, and N. H. March, Croatica Chemica Acta , 463 (2013). S. H. Vosko, L. Wilk, and M. Nusair, Canadian Journal of Physics , 1200 (1980), https: // doi.org / / p80-159. U. von Barth and L. Hedin, Journal of Physics C: Solid StatePhysics , 1629 (1972). Z.-X. Shen, R. S. List, D. S. Dessau, B. O. Wells, O. Jepsen, A. J.Arko, R. Barttlet, C. K. Shih, F. Parmigiani, J. C. Huang, andP. A. P. Lindberg, Phys. Rev. B , 3604 (1991). O. Tjernberg, S. Sderholm, T. Rogelet, U. Karlsson, M. Qvarford,I. Lindau, C.-O. Almbladh, and L. Hellbom, Vacuum , 1215(1995). O. Tjernberg, S. S¨oderholm, G. Chiaia, R. Girard, U. O. Karlsson,H. Nyl´en, and I. Lindau, Phys. Rev. B , 10245 (1996). A. B. Shick, A. I. Liechtenstein, and W. E. Pickett, Phys. Rev. B , 10763 (1999). S. Maekawa, T. Tohyama, S. E. Barnes, S. Ishihara, W. Koshibae,and G. Khaliullin,
Physics of Transition Metal Oxides (Springer,2004). R. W.Yang, Oxford University Press (1994). F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann,and A. I. Lichtenstein, Phys. Rev. B , 195104 (2004). D. Korotin, A. V. Kozhevnikov, S. L. Skornyakov, I. Leonov,N. Binggeli, V. I. Anisimov, and G. Trimarchi, The EuropeanPhysical Journal B , 91 (2008). M. Karolak, G. Ulm, T. Wehling, V. Mazurenko, A. Poteryaev,and A. Lichtenstein, Journal of Electron Spectroscopy and Re-lated Phenomena , 11 (2010), proceedings of InternationalWorkshop on Strong Correlations and Angle-Resolved Photoe-mission Spectroscopy 2009. P. Staar, B. Ydens, A. Kozhevnikov, J.-P. Locquet, andT. Schulthess, Phys. Rev. B , 245114 (2014). A. Kozhevnikov, A. G. Eguiluz, and T. C. Schulthess, 2010ACM / IEEE International Conference for High Performance Com-puting, Networking, Storage and Analysis , 1 (2010). I. Leonov, Phys. Rev. B , 085142 (2015). I. A. Nekrasov, N. S. Pavlov, and M. V. Sadovskii, Journal ofExperimental and Theoretical Physics , 620 (2013). N. F. Mott and R. Peierls, Proceedings of the Physical Society ,72 (1937). N. F. Mott, Proceedings of the Physical Society. Section A , 416(1949). A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev.Mod. Phys. , 13 (1996). W. R. L. Lambrecht and O. K. Andersen, Phys. Rev. B , 2439(1986). N. Marzari and D. Vanderbilt, Phys. Rev. B , 12847 (1997). N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vander- bilt, Rev. Mod. Phys. , 1419 (2012). J. Hubbard, Proceedings of the Royal Society of London A: Math-ematical, Physical and Engineering Sciences , 238 (1963),http: // rspa.royalsocietypublishing.org / content / / / B. A.Altland, Cambridge University Press (2006). E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer,and P. Werner, Rev. Mod. Phys. , 349 (2011). P. Staar, T. Maier, and T. C. Schulthess, Phys. Rev. B , 115101(2013). “ DCA + + : A software framework to solve correlated electronproblems with modern quantum cluster methods (manuscript un-der review),” https://github.com/CompFUSE/DCA (2018). X. Ren, I. Leonov, G. Keller, M. Kollar, I. Nekrasov, and D. Voll-hardt, Phys. Rev. B , 195114 (2006). J. Kuneˇs, V. I. Anisimov, A. V. Lukoyanov, and D. Vollhardt,Phys. Rev. B , 165115 (2007). J. Kuneˇs, V. I. Anisimov, S. L. Skornyakov, A. V. Lukoyanov, andD. Vollhardt, Phys. Rev. Lett. , 156404 (2007). O. Miura and T. Fujiwara, Phys. Rev. B , 195124 (2008). J. c. v. Kolorenˇc, A. I. Poteryaev, and A. I. Lichtenstein, Phys.Rev. B , 235136 (2012). I. A. Nekrasov, V. S. Pavlov, and M. V. Sadovskii, JETP Letters , 581 (2012). P. Thunstr¨om, I. Di Marco, and O. Eriksson, Phys. Rev. Lett. ,186401 (2012). R. Sakuma and F. Aryasetiawan, Phys. Rev. B , 165118 (2013). K. Byczuk, J. Kuneˇs, W. Hofstetter, and D. Vollhardt, Phys. Rev.Lett. , 087004 (2012). J. KuneÅ, A. V. Lukoyanov, V. I. Anisimov, R. T. Scalettar, andW. E. Pickett, Nature Materials , 198 (2008). A. O. Shorikov, Z. V. Pchelkina, V. I. Anisimov, S. L. Skornyakov,and M. A. Korotin, Phys. Rev. B , 195101 (2010). K. Ohta, R. E. Cohen, K. Hirose, K. Haule, K. Shimizu, andY. Ohishi, Phys. Rev. Lett. , 026403 (2012). L. Huang, Y. Wang, and X. Dai, Phys. Rev. B , 245110 (2012). A. A. Dyachenko, A. O. Shorikov, A. V. Lukoyanov, and V. I.Anisimov, JETP Letters , 56 (2012). I. Leonov, L. Pourovskii, A. Georges, and I. A. Abrikosov, Phys.Rev. B , 155135 (2016). O. Grns, I. D. Marco, P. Thunstrm, L. Nordstrm, O. Eriksson,T. Bjrkman, and J. Wills, Computational Materials Science ,295 (2012). H. Park, A. J. Millis, and C. A. Marianetti, Phys. Rev. B ,235103 (2014). S. Bhandary, E. Assmann, M. Aichhorn, and K. Held, Phys. Rev.B , 155131 (2016). I. Leonov, V. I. Anisimov, and D. Vollhardt, Phys. Rev. B ,195115 (2015). E. Greenberg, I. Leonov, S. Layek, Z. Konopkova, M. P. Paster-nak, L. Dubrovinsky, R. Jeanloz, I. A. Abrikosov, and G. K.Rozenberg, Phys. Rev. X , 031059 (2018). R. J. Powell and W. E. Spicer, Phys. Rev. B , 2182 (1970). G. A. Sawatzky and J. W. Allen, Phys. Rev. Lett. , 2339 (1984). G. W. Pratt and R. Coelho, Phys. Rev. , 281 (1959). R. Gillen and J. Robertson, Journal of Physics: Condensed Matter , 165502 (2013). C. Kant, F. Mayr, T. Rudolf, M. Schmidt, F. Schrettle, J. Deisen-hofer, and A. Loidl, The European Physical Journal Special Top-ics , 43 (2009). M. Gvishi and D. Tannhauser, Journal of Physics and Chemistryof Solids , 893 (1972). C. R¨odl, F. Fuchs, J. Furthm¨uller, and F. Bechstedt, Phys. Rev. B , 235114 (2009). H. Bowen, D. Adler, and B. Auker, Journal of Solid State Chem-istry , 355 (1975). H. Shinaoka, M. Troyer, and P. Werner, Phys. Rev. B , 245156(2015). A.-V. Plamad,
A systematic investigation of the double countingcorrections in first principle calculations of strongly correlatedmaterials , Ph.D. thesis, ETH Zurich (2015), dissertation. ETH-Zrich. 2015. Nr. 23141. V. Heine, Phys. Rev. , A1689 (1965). J. Terso ff , Phys. Rev. Lett. , 465 (1984). C. Tejedor, F. Flores, and E. Louis, Journal of Physics C: SolidState Physics , 2163 (1977). F. Flores and J. Ortega, Applied Surface Science , 301(1992). P. C. D. Hobbs, R. B. Laibowitz, and F. R. Libsch, Appl. Opt. ,6813 (2005). P. Mavropoulos, N. Papanikolaou, and P. H. Dederichs, Phys.Rev. Lett. , 1088 (2000). X.-G. Zhang and W. H. Butler, Journal of Physics: CondensedMatter , R1603 (2003). E. Prodan and R. Car, Phys. Rev. B , 035124 (2009). F. Aryasetiawan, K. Karlsson, O. Jepsen, and U. Sch¨onberger,Phys. Rev. B , 125106 (2006). M. Petersilka, U. J. Gossmann, and E. K. U. Gross, Phys. Rev.Lett. , 1212 (1996). L. Vaugier, H. Jiang, and S. Biermann, Phys. Rev. B , 165105(2012). “ELK: an all-electron full-potential linearised augmented-planewave code,” http://elk.sourceforge.net/http://elk.sourceforge.net/