First-principles treatment of Mott insulators: linearized QSGW+DMFT approach

The theoretical understanding of emergent phenomena in quantum materials is one of the greatest challenges in condensed matter physics. In contrast to simple materials such as noble metals and semiconductors, macroscopic properties of quantum materials cannot be predicted by the properties of individual electrons. One of the examples of scientific importance is strongly correlated electron system. Neither localized nor itinerant behaviors of electrons in partially filled 3d, 4f, and 5f orbitals give rise to rich physics such as Mott insulators, high-temperature superconductors, and superior thermoelectricity, but hinder quantitative understanding of low-lying excitation spectrum. Here we present a new first-principles approach to strongly correlated solids. It is Q4 based on a combination of the quasiparticle self-consistent GW approximation and the dynamical mean-field theory. The sole input in this method is the projector to the set of correlated orbitals for which all local Feynman graphs are being evaluated. For that purpose, we choose very localized quasiatomic orbitals spanning large energy window, which contains most strongly hybridized bands, as well as upper and lower Hubbard bands. The self-consistency is carried out on the Matsubara axis. This method enables the first-principles study of Mott insulators in both their paramagnetic and antiferromagnetic phases. We illustrate the method on the archetypical charge transfer correlated insulators La2CuO4 and NiO, and obtain spectral properties and magnetic moments in good agreement with experiments.


Introduction.
The first principles description of strongly-correlated materials is currently regarded as one of the greatest challenges in condensed matter physics. The interplay between localized electrons in open d-or fshell and itinerant band states gives rise to rich physics that makes these materials attractive for a wide range of applications such as oxide electronics, high temperature superconductors and spintronic devices. Various theoretical approaches are currently being pursued [1]. One of the most successful approaches is the dynamical mean field theory (DMFT) [2]. In combination with density functional theory [3,4], it has described many features of strongly-correlated materials successfully and highlighted the surprising accuracy of treating correlations local to a small subset of orbitals exactly, while treating the reminder of the problem in a static mean field manner. [5,6].
Correlations are built into the full electron Green's function (G), through electron self-energy (Σ) G(r, r ′ , iω n ) = r| iω n 1 −Ĥ H −Σ(iω n ) where r is position vector, ω n is Matsubara frequency andĤ H is the Hartree Hamiltonian. The successes of the DMFT method suggest that in certain energy range a good expression for the self-energy Σ is provided by Σ(r, r ′ , iω n ) = Rαβ P * Rα (r)Σ α,β loc (iω n )P Rβ (r ′ ) + V (r, r ′ ) (2) with Σ loc described by an impurity model, where R is lattice vector and P Rα (r) is the projector to a set of local orbitals α centered at R.
The numerous successes of DMFT in different classes of correlated materials revived the interest in the long sought goal of achieving a diagrammatically controlled approach to the quantum many body problem of solids. The free energy functional can be expressed in terms of the Green's function G and the screened Coulomb interactions W [7,8]. The lowest order perturbation theory in this functional gives rise to the GW approximation [9] while the local approximation applied to the most correlated orbitals gives rise to an extended DMFT approach to the electronic structure problem [8]. The addition of the GW and DMFT graphs was proposed and implemented in model Hamiltonian studies [10] and in realistic electronic structure [11,12]. There is now intense activity in this area with many recent publications [13][14][15][16] triggered by advances in the quality of the impurity solvers [17][18][19], insights into the analytic form of the high frequency behavior of the self-energy [20] and improved electronic structure codes.
Several conceptual issues remain to be clarified before the long sought goal of a robust electronic structure method for solids is attained. The first issue is the choice of local orbitals (P Rα (r) in Eq. (2)) on which to perform the DMFT method (summation of all local Feynman graphs). The second issue is the level of self-consistency that should be used in the calculation of various parts of the diagrams included in the treatment ( free or bare Green's function G 0 vs self-consistent interacting Green's functions G).
The self-consistency issue appears already at the lowest order, namely the GW level, and it has been debated over time. The corresponding issue in GW+DMFT is expected to be at least as important, but has not been explored, except for model Hamiltonians [21,22]. At the GW level, it is now well established that Hedin's fully self-consistent formulation [9], while producing good total energies in solids [23], atoms and molecules [24,25], does not produce a good approximation to the spectra of even simple semiconductors or weakly correlated metals. Instead, using a free (quasiparticle) Green's function in the evaluation of the polarization graph of the GW method gives much better results for spectral functions. This is the basis of the one-shot quasiparticle (QP) GW, starting from LDA [26] or from others [27,28]. Unfortunately, the answer depends on the starting point. A solution for this problem is to impose a self-consistency equation to determine G 0 . This method, called the quasiparticle self-consistent GW (QSGW) [29], is very success- ful reproducing the spectra of many systems [30][31][32].
Previous GW+DMFT studies typically used one shot QPGW and projectors spanning a relatively small energy window [13][14][15][16]. In this work, we propose a different approach to the two issues: the level of self-consistency and the choice of the DMFT orbital. First we use the quasiparticle self-consistent GW in combination with DMFT to address the self-consistency issue. Next, we choose a very localized orbital for the summation of the higher order Feynman diagrams in DMFT, therefore the hybridization spans much larger energy window.
In the LDA+DMFT context, the choice of very localized orbitals has provided a great deal of universality since the interactions do not vary much among compounds of the same family. This has been demonstrated in the studies of iron pnictides [33] and transition metal oxides [34]. This choice results in a second advantage as we will show below, namely the frequency dependence of the interaction matrix can be safely ignored. After the orbitals are chosen, all the parameters are selfconsistently determined. This is the first ab initio quasiparticle self-consistent GW+DMFT implementation and the first study on a paramagnetic Mott insulator within the GW+DMFT method.
Methods. Our approach is carried it out entirely on the Matsubara axis, which requires a different approach to the quasiparticle self-consistency in GW [35], called Matsubara Quasiparticle Self-consistent GW (MQSGW), where the quasiparticle Hamiltonian is constructed by linearizing the self-energy and renormalization factor  [36]. Working on the Matsubara axis, is numerically very stable, provide a natural interface with advanced DMFT solvers such as continuous-time quantum Monte-Carlo (CTQMC) [17][18][19] and has very good scaling in system size as in the space-time method [37]. (see Supplemental Material [38] for details).
For DMFT, it is essential to obtain bandstructure in a fine enough crystal momentum (k) mesh to attain desired frequency resolution of physical quantities. To achieve such momentum resolution, we use a Wannier-interpolated MQSGW bandstructure in a large energy window using Maximally localized Wannier function (MLWF) [39], and than constructed local projector in a fine momentum mesh. In contrast to SrVO 3 [13][14][15][16] where a set of t 2g states is reasonably well separated from the other bands, correlated 3d orbitals in La 2 CuO 4 are strongly hybridized with other itinerant bands. In this case, it is necessary to construct local projectors from states in a wide enough energy windows to make projectors localized near the correlated atoms. We constructed local projectors in the energy window E F ±10eV  Fig. 1(b). The Dashed lines in (b) and (c) represent electronic bandstructures within non spin-polarized MQSGW and spin-polarized MQSGW, respectively in which there are ∼82 bands at each k point, where E F is the Fermi level. Then we confirmed that absolute value of its overlap to the muffin-tin orbital (of which radial function is determined to maximize electron occupation in it) is more than 95%. Our choice of energy window is justified by the Cu-3d spectra being entirely contained in this window. Using constructed MLWFs in large energy window, we defined our local-projector is quasiparticle wavefunction with an index n, and N k is the number of k points in the first Brillouin zone.
Static U d and J H are evaluated by a modification of the constrained RPA method [40], which avoids screening by the strongly hybridized bands. This screening by hybridization is included in our large energy window DMFT. For details, see Supplemental Material [38]. We divide dynamic polarizability within MQSGW approximation χ QP into two parts, χ QP = χ low QP + χ high QP . Here, χ low QP is defined by all transitions between the states in the energy window accounted for by the DMFT method (E F ±10eV ). Using χ high QP , we evaluate partially screened Coulomb interaction U −1 (r, r ′ , k, iω n ) = V −1 (r, r ′ , k) − χ high QP (r, r ′ , k, iω n ) and parametrize static U d and J H by Slater's integrals [41,42], where V is bare Coulomb interaction.
The Feynman graphs included in both MQSGW and DMFT (double-counting) are the local Hartree and the local GW diagram. They are computed using the local projection of the MQSGW Green's function (Ĝ QP ) Finally, for the stable numerics, we approximated Σ DC (iω n ) ≃Σ DC (iω n = 0) since these low order diagrams are dominated by the Hartree-Fock contribution.
Results. Fig. 2(a) shows the frequency dependence of real and imaginary parts of U d . It is calculated on an imaginary frequency axis and analytically continued using a Pade approximant [44]. We also plot the fully screened Coulomb interaction W d for comparison. Static U d is 12.5 eV and U d remains almost constant up to 10 eV. In contrast, in W d , there are several peaks due to low-energy collective excitations below 10 eV. At very high energy, U d approaches the bare coulomb interaction of 28 eV. Calculated J H is 1.4 eV and has negligible frequency dependence. By contrast, conventional constrained-RPA, in which 10 bands of mostly Cu-3d character are excluded from screening, results in static U d = 7.6 eV, which is too small to open the Mott gap, and which is also inconsistent with photoemission experiments on CuO charge transfer insulators [45].
We also computed the static U d and J H by requiring that the calculated excitation spectra of MQSGW+DMFT with (local) GW as the impurity solver matches the spin-polarized MQSGW spectra. Here we used non spin-polarized MQSGW band structure and allowed spontaneous magnetic long range order by embedding impurity self energy, which is function U d and J H , within spin-polarized GW approximation. In Fig.  2(b), we allowed U d to vary between 8-13 eV (at fixed J H = 1.4 eV) and we plot the size of the indirect gap. The gap size of this method matches the gap of spin-polarized MQSGW when U d ≈ 12 eV. If the choice of U d and J H is correct, the resulting spectra must be similar to the prediction of spin-polarized MQSGW method. We show this comparison in Fig. 2(c) to confirm a good match. In addition, the relative position of Cu-d band (the lowest en- ergy conduction band at S) to the La-d band (the lowest energy conduction band at Y) is also well matched justifying the approximation ofΣ DC (iω n ) ≃Σ DC (iω n = 0). Σ DC (iω n = 0) for Cu-d x 2 −y 2 orbital differs from nominal double counting energy [46] by only 1%, highligting again the advantages of using a broad window and narrow orbitals.
We now discuss the magnetic moment associated with Cu and the electronic excitation spectra by using MQSGW+DMFT (with U d = 12.5eV , J H = 1.4eV ) in which the impurity is solved by the numerically exact CTQMC [17,18] and compare them with other methods. LSDA does not have a magnetic solution. In contrast, spin-polarized MQSGW, QSGW [29], and MQSGW+DMFT predict 0.7 µ B , 0.7 µ B , and 0.8 µ B , respectively. This is consistent with experimental measurements, although the later span quite large range 0.4µ B ∼ 0.8µ B [47][48][49]. We find that MQSGW opens too large gap in the antiferromagnetic (AFM) phase, while it remains metallic in the paramagnetic (PM) phase. LDA+DMFT predict too small excitation energy of Laf levels. We show that these deficiencies can be remedied by adding all local Feynman diagrams for the Cu-d orbitals using the DMFT and treating the other states within GW approximation.
In the low-energy spectrum, LSDA does not have a insulating solution; there is a single non-magnetic solution with zero energy gap as shown in the bandstructure (Fig.  3(a)) and total density of states ( Fig. 4(a)). The non spin-polarized MQSGW also predicts metal as shown in Fig. 4(a), but the two bands of primarily Cu-d x 2 -y 2 character near the Fermi level are here well separated from the rest of the bands (dashed lines in Fig. 3(b)). Spinpolarized MQSGW calculation (dashed lines in Fig. 3(c)) yields qualitative different results from LSDA and non spin-polarized MQSGW calculation. The two Cu-d x 2 -y 2 bands are now well separated from each other with a bandgap of 3.4 eV. Spin-polarized QSGW [29] also yields insulating phase with a gap of 4.0 eV. In the experiment, the larger direct gap, as measured by optics, is ∼ 2eV [50,51].
We show that these deficiencies of LDA, QSGW and MQSGW in the low-energy spectra can be remedied by adding all local Feynman diagrams for the Cu-d orbitals using the DMFT. The LDA+DMFT calculation in Fig. 4(a), carried out by the all-electron LDA+DMFT method [34,46], predicts reasonable gap of 1.5 eV and 1.8 eV in PM and AFM phases, in good agreement with experiment and previous LDA+DMFT studies [34,[52][53][54][55]. Within MQSGW+DMFT, we find gaps of 1.5 eV and 1.6 eV in PM and AFM phases, respectively. The excitation spectra of MQSGW+DMFT in PM and AFM phase as shown in Fig. 3(b) and 3(c) are very similar as both are insulating with well separated Cu-d x 2 -y 2 bands, which is now also substantially broadened due to large scattering rate in Hubbard-like bands. The projected density of states of Cu-d x 2 -y 2 orbital in Fig. 4(b) shows the size of the gap more apparently. This gap is of correlated type as it results from the non-perturbative pole in the impurity self-energy near zero frequency. This pole is connected to spin fluctuations, which at low temperature results in an ordered AFM phase. The Zhang-Rice peak appears around ∼−2eV and strongly overlapps with O-p states .
In the high energy region, the most distinctive difference is the position of La-f peak. It appears at ∼ 3eV within LDA and LDA+DMFT, but at around ∼ 9eV , in the inverse-photoemission spectra (cyan dotted line in Fig. 4(a)) [43]. By treating La-f within GW approximation, it appears at ∼ 10eV within MQSGW and MQSGW+DMFT.
In summary, we introduced a new methodology within MQSGW+DMFT and tested it in the classic charge transfer insulator La 2 CuO 4 . Our methodology predicts a Mott-insulating gap in the PM phase, thus overcoming the limitation of LDA and QSGW. It yields more precise peak positions of the La-f states, thus improving the results of LDA+DMFT. The method should be useful in understanding electronic excitation spectrum of other strongly-correlated materials, in particular those where precise position of both the itinerant and correlated states is important.

COMPUTATIONAL DETAILS: BASIS SET
All calculations are performed using our relativistic spin-polarized, full-potential, linearizedaugmented-plane-wave package (RSPFLAPW) [1,2], which is based on full-potential linearized augmented plane wave plus local orbital (FPLAPW+lo) method. For this particular calculation, experimental lattice constants and atomic positions at the low-temperature orthorhombic phase [3]

MATSUBARA QSGW CALCULATIONS
The electron self-energy can be systematically expanded in terms of the dressed Green's function G and the screened Coulomb interaction W . Within GW approximation, we keep only the first term of the series expansion of self-energy in W : where r is the position vector in a unit cell, k is the crystal momentum, R is the lattice vector, ω n is Matsubara frequency, τ is Matsubara time, and s is spin index. Within Matsubara quasiparticle self-consistent GW (MQSGW) approximation, we calculate the dynamic polarizability and the electron self-energy by using quasiparticle Green's function (G QP ) instead of Full GW Green's function. First, we construct QP Green's function using the quasiparticle Hamiltonian (H QP ) by For the first iteration, we regard Hamiltonian within local density approximation as H QP .
Next, the screened Coulomb interaction W QP is evaluated using the dynamic polarizability χ QP (r, r ′ , k, ′ iω n ) by W −1 QP (r, r ′ , k, iω n ) = V −1 (r, r ′ , k) − χ QP (r, r ′ , k, iω n ) within random phase approximation (RPA): where V is bare Coulomb interaction. Then, we calculate MQSGW electron self-energy Then, we constructed quasiparticle Hamiltonian with linearized self-energy and renormal- whereĤ H (k) is Hartree Hamiltonian. This process is repeated until self-consistency is attained.
The static U d and J H are evaluated by a method similar to constrained RPA [4], but here we avoid screening by the strongly hybridized bands, for which the screening by hybridization is included in our large energy window DMFT. We divide dynamic polarizability χ QP into two parts, χ QP = χ low QP + χ high QP . Here, χ low QP is defined by all transitions between the states in the energy window accounted for by the DMFT method (E F ± 10eV ): Using χ high QP , we evaluate partially screeend Coulomb interaction by U −1 (r, r ′ , k, iω n ) = V −1 (r, r ′ , k) − χ high QP (r, r ′ , k, iω n ). Then, We then parametrize static U = F 0 and J H = (F 2 + F 4 )/14 by Slater's integrals [5], drdr ′ U(r, r ′ , R = 0, iω n = 0)W * R=0,m 1 (r)W * R=0,m 2 (r ′ )W R=0,m 3 (r ′ )W R=0,m 4 (r), where W R,i (r) is Wannier function centered at R with index m and Y lm is spherical har- Here, P i,n (k) = R < W Ri |ψ nk > e ik·R / √ N k is projector to the correlated subspace at each k and ψ nk (r) is quasiparticle wavefunction with an index n. N k is the number of k points in the first Brillouin zone.Ê imp and∆ imp are impurity level energy and hybridization function, which are inputs to impurity solver.Σ embeded = P † (k) Σ imp (s, iω n ) −Σ DC (iω n ) P (k) is embeded self-energy with impurity self-energy (Σ imp ) and double-counting correction (Σ DC ).