First principles vs second principles: Role of charge self-consistency in strongly correlated systems

First principles approaches have been successful in solving many-body Hamiltonians for real materials to an extent when correlations are weak or moderate. As the electronic correlations become stronger often embedding methods based on first principles approaches are used to better treat the correlations by solving a suitably chosen many-body Hamiltonian with a higher level theory. Such combined methods are often referred to as second principles approaches. At such level of the theory the self energy, i.e. the functional that embodies the stronger electronic correlations, is either a function of energy or momentum or both. The success of such theories is commonly measured by the quality of the self energy functional. However, self-consistency in the self-energy should, in principle, also change the real space charge distribution in a correlated material and be able to modify the electronic eigenfunctions, which is often undermined in second principles approaches. Here we study the impact of charge self-consistency within two example cases: TiSe$_{2}$, a three-dimensional charge-density-wave candidate material, and CrBr$_{3}$, a two-dimensional ferromagnet, and show how real space charge re-distribution due to correlation effects taken into account within a first principles Green's function based many-body perturbative approach is key in driving qualitative changes to the final electronic structure of these materials.

Density functional theory [1][2][3] has been the workhorse for material specific electronic structure calculations for the last half of the century.Despite enormous success in many respects, it has however some intrinsic limitations.First of all, although the Hohenberg-Kohn theorem [1] guarantees the existence of some density functional providing an exact ground state energy at a given charge density distribution ρ, its exact form is unknown.In practice, this functional is considered as being local or almost local (generalized gradient corrections), which is generally speaking an uncontrollable approximation (for detailed discussions see the review [3]).Next, and even more importantly, the Kohn-Sham quasiparticles [2] are, generally speaking, just auxiliary quantities to calculate the total energy and their direct comparison with experimental spectroscopic information is hardly justifiable.Although this is regularly done with partial excellent agreement, there are numerous counterexamples starting from the famous "gap problem" in semiconductors [4].
An alternative approach is based on the concept of Green's function functionals.Luttinger-Ward [5] and Baym-Kadanoff [6] theorems respectively prove the existence of such functionals in-and out-of-equilibrium.Conceptually, this way is more attractive since the knowledge of an exact single-and two-particle Green's functions guarantees an accurate description of spectroscopic properties of solids [7].On the other hand, again, an exact form of this functional is practically unknown and we have just its formal definition in terms of infinite sums of skeleton free-leg diagrams [5,6].If we are interested in a description of subtle phenomena such as, e.g., the Kondo effect [8] or nonquasiparticle states in half-metallic ferromagnets [9], the necessary sequence of diagrams seems to be too complicated to be practically taken into account for a complete first-principle realization.
Therefore, alternative embedding approaches were introduced which combine first-principle calculations with model treatments to describe the strong correlations within some low-energy subspace.This way, weakly correlated states at high energies are described within a low-level theory, while the strongly correlated sub-space is treated in higher level approaches.This is popularly done by mapping the low-energy space to multi-band generalized Hubbard models, which are afterwards often solved, e.g., using dynamical mean-field theory (DMFT) [10], a program suggested and called "LDA ++ " in Ref. 11 and which we refer to in the following as "second principles".In many cases, this leads to a dramatic improvement of description of strong correlation effects in real materials with itinerant-electron magnets [12] and heavy-fermion compounds [13] being two major successful examples (for detailed reviews see [9,14,15]).Of course, DMFT [16] is a local approximation which takes only the energy dependence of electronic self energy into account and completely neglects its momentum dependence.However, the latter can be taken into account via various beyond-DMFT diagrammatic approaches [17], rendering it a technical problem rather than a fundamental one.Also the way how one can map the first-principle electronic structure onto efficient Hamiltonians can be, in principle, improved.The contemporary way is based on the so-called constrained RPA (cRPA) approach [18,19] but there are no principle obstacles to improve it further if necessary.
A key impediment to second principles approaches is, however, that multiple energy scales are operative: the highenergy scales controlling low-energy fluctuations cannot be integrated out without model assumptions.High-energy scales contain information about chemistry and disorder specific to real materials.Yet while first principles theories contain this information, their application to strongly correlated systems has been limited.In weakly correlated materials, "first principles" approaches tend to predominate because they rely on a minimum of model assumptions, and are often predictive.This is not the case when correlations are strong because standard methods, usually based on extensions of density functional theory (DFT), lack the sophistication to encapsulate the strong spin and charge fluctuations, or the fidelity to characterize one-particle properties near the Fermi level (which are essential to capture low-energy excitations characteristic of correlated systems); nor are they adequately equipped to generate the (twoparticle) suceptibilities.Even in weakly correlated cases, dynamical screening (perhaps the most important many-body effect [20]) is not well treated by such standard one-body descriptions [21,22].The difficulties are even more severe for spin fluctuations, where the characteristic energy scale for excitations can be very small.Models such as the Hubbard Hamiltonian do indeed contain, within some region of parameter space, key manybody effects such as the metal-insulator transition, pseudogap phases, quantum criticality, and both conventional and unconventional superconductivity.Thus, applications of this model has become the canonical approach to characterizing such phenomena.However, the limits to such an approach become apparent when the high-energy scales that control parameters for the low-energy ones are nontrivial.Furthermore, within any model Hamiltonian, there are only two possibilities for the correlation effects to modify the electronic structure, namely via the energy-and/or the momentum-dependencies of the self-energy Σ.In first-principle approaches there is, however, an additional possibility in form of the charge self-consistency.That is, if the correlation effects can essentially modify the charge distribution ρ, than one needs to recalculate the model Hamiltonian at every iteration which makes the mapping procedure very cumbersome or even practically useless.The physical question of fundamental importance is the following: can a real-space charge redistribution due to correlation effects be qualitatively important leading not just to a moderate renormalization of the model parameters but also to a reconstruction of the electronic structure beyond any purely We show in the following how different levels of theory significantly modify the effective one-body potential through changes in the electron density.To this end, we employ three different levels of theory: the local-density approximation (LDA), quasiparticle self-consistent GW theory [23][24][25] (QSGW ), which, in contrast to conventional GW, modifies the charge density and is determined by a variational principle [26], and finally an extension of QSGW, where the polarizability needed to construct W is computed including vertex corrections (ladder diagrams) by solving a Bethe-Salpeter equation (BSE) for the two-particle Hamiltonian [27].We denote the latter QSG W , with the substitution W → W signifying that a BSE was solved to compute W.These first-principles approaches allow us to carefully analyze the impact of the full charge self consistency taking correlation effects with increasing diagrammatic precision into account.
In terms of diagram classes taken into account QSGW and QSG W represent the forefront of currently available first-principle approaches.As we show, it is essential that the first-principles starting point is of sufficiently high fidelity to capture physics the second-principles scheme cannot reach.First-principles schemes are too cumbersome to handle more than a limited class of diagrams, and it may still be true in general that second-principles schemes may still be needed to capture physics outside the reach of the first-principles scheme.Low-energy spin fluctuations, Kondo effect and non-quasi-particle states in weakly doped Mott insulators and half-metallic ferromagnets seem to be the archetypal example of this.For TiSe 2 and CrBr 3 , QSGW / QSG W appears, however, to adequately describe most physical observables, obviating the need for second-principles schemes.TiSe 2 : TiSe 2 is a layered diselenide compound with space group P 3m1 (Fig. 1).Below 200K, it undergoes a phase transition to a charge-density wave (CDW), forming a commensurate 2×2×2 superlattice (P 3c1) of the original structure.At the transition there is a softening of the zone boundary phonon accompanied by changes in the transport properties [28,29].
A number of works have tried to determine the energy dispersion around the Fermi level in both phases with a special focus on the overlap/gap between the Se-4p valence band at Γ and the Ti-3d conduction band at L, sometimes with discordant results.In the high-temperature P 3m1 phase, reports range from predicting a semimetal (overlap <120 meV) between these bands, to an insulator with <60 meV gap, depending on the study [30][31][32][33][34][35][36][37].In the CDW P 3c1 phase there is a greater consensus, namely that the gap is small and positive.In brief, upon cooling the system, the CDW transition induces a distortion that either (slightly) increases the existing gap or leads to a gap opening between these bands with the gap for the CDW phase being ∼100-150 meV [30,[33][34][35]37].
At high temperature, the positive or negative (i.e., overlap) indirect gap between Se-4p and Ti-3d bands is larger than the negative gap obtained with standard DFT calculations.In fact, DFT is not helpful because it predicts a negative gap in both the undistorted and CDW phases.Bianco et al. [38] finds with LDA+U (U =3.9 eV) a gap of 14 meV in the P 3m1 phase and 200 meV in the CDW phase.LDA+U is a kind of "second principles" method because the answer depends U , which is not known.To check whether the negative gap is merely an artifact of the model, Cazzaniga et al. [39] considered a G 0 W 0 calculation based on the LDA, and found a gap of ∼200 meV in the high-temperature P 3m1 structure.
We show here that while GW does indeed modify the quasiparticle spectrum, the true situation is more complex.This is because not only the eigenvalues but the density is significantly renormalized relative to the LDA.This induces a corresponding change in the effective potential through the inverse of the susceptibility, χ −1 (x 1 , x 2 ) = δV (x 1 )/δn(x 2 ).What appears to be special about TiSe 2 is that χ −1 (1, 2) is large, and the correction to the LDA density modifies the effective potential V in such a way as to reduce the splitting between occupied and unoccupied levels.This is very unusual: it has long been established that in the vast majority of cases, GW based on the LDA continues to underestimate the gap in semiconductors, albeit less so than the LDA [40].
These effects can only be found through self-consistency.QSGW is ideally suited for this case, as its excitation spectra are generally superior to fully self-consistent GW [41][42][43][44][45].We find that the P 3m1 is indeed semimetallic, as is the case with DFT, but for different reasons.We first revisit the GW calculation of the undistorted P 3c1 structure, but with some modifications: • we did not include a Z factor.There are various justifications for this, most notably as an approximate way to incorporate self-consistency in G with fixed W ; see Appendix in Ref. 40.Omission of Z tends to widen bandgaps.
• the full matrix G LDA W LDA is used, in the QSGW sense [24]: Panel (a) of Fig. 2 shows LDA and G LDA W LDA bands similar to the G LDA W LDA calculation of Ref. 39. Focusing on the LDA bands, the highest occupied state at Γ turns red very close to Γ, indicating the penetration of the Ti-derived conduction band into the valence band (resembling a "negative" bandgap).This is an artifact of the LDA's well known tendency to underestimate splittings between occupied and unoccupied states, and as standard G LDA W LDA increases this separation (blue dashed lines).The (indirect) G LDA W LDA gap of 300 meV is slightly larger than Ref. 39, in line with the unit Z factor used in the present calculation.
Fig. 2(b) shows that self-consistency is crucially important in TiSe 2 .The off-diagonal elements of Σ 0 ij modify the density n(r) and thus V (r).A simple way to estimate ∆V is to make an ansatz that the LDA adequately yields χ −1 = δV /δn.For a modified n the potential becomes V (n) = Σ 0 − V xc (n LDA ) + V xc (n).n can be determined selfconsistently in the usual manner by adding a fixed external potential Σ 0 − V xc (n LDA ) to the LDA Hamiltonian and allowing it to go self-consistent.Remarkably, the gap becomes negative again, as shown by the blue dashed lines in Fig. 2(b), but the dispersion is very different from the LDA.In particular the inverted gap character at Γ disappears, which is topologically essential for a gap to form at Γ.The quality of the ansatz can be checked by carrying out a complete QSGW calculation.This is shown as solid lines in Fig. 2(b), and it demonstrates the ansatz is reasonable.As we will show elsewhere, the observed low-temperature gap forms as a consequence of the charge density-wave instability.CrBr 3 : Monolayer (1L) of CrBr 3 is a two-dimensional ferromagnetic (FM) insulator where the magnetic moments of monolayer CrBr 3 align normal to the plane (see Fig. 1 for the crysal structure).The spontaneous magnetization persists in monolayer CrBr 3 with a Curie temperature of 34 K [46].Within a purely atomic picture, fully determined by the crystal field environment and the Hund's multiplet structure [47], the low energy properties of the materials and the magnetism should be entirely governed by Cr-d electrons.However, the ligands, their masses and the number of core states they have, play an integral role in determining the low energy properties of CrBr 3 .In a separate work we discuss the role of the ligands like Cl, Br and I in determining the crucial low energy properties of the entire class of 1L Chromium trihalides [48].The Br-p states strongly hybridize with Cr-d states in CrBr 3 .In the present work we show how charge self-consistency at different levels of the theory controls the nature of the eigenfunctions and the Br component in the valence band manifold of 1L CrBr 3 .
We simulate the free standing 1L of ferromagnetic CrBr 3 within LDA, QSGW and QSG W .We also perform a rigorous check for vacuum correction to all relevant observables by increasing the vacuum size from 20 Å to 80 Å [48].We check for convergence and scaling of band gap and the dielectric constant ∞ with vacuum size as discussed in a separate work [48].We observe that FM-1L CrBr 3 is an insulator with 1.3 eV of electronic band gap in LDA, which is significantly lower than in QSGW yielding a gap of 5.7 eV.The large enhancement in QSGW band gaps relative to the LDA is standard in polar compounds [23].Nevertheless, within the random phase approximation (RPA), it has long been known that W is universally too large [49,50], which is reflected in an underestimate of the static dielectric constant ∞ .Empirically, ∞ seems to be underestimated in QSGW by a nearly universal factor of 0.8 [51], for a wide range of insulators [52,53] resulting in slightly overestimated [23] band gaps.This can be corrected by extending the RPA screening to introduce an electron-hole attraction in virtual excitations.These extra (ladder) diagrams are solved by a BSE, and they significantly improves on the optics, largely eliminating the discrepancy in ∞ [54].When ladder diagrams are also added to improve W in the GW cycle (W → W ), it significantly improves the one-particle gap as well [27,55].In detail, our QSG W implementation is self-consistent in the sense that the updated W also updates Σ and hence the cycle continues until W , Σ and G converge iteratively.This scenario is played out in CrBr 3 : the QSGW bandgap is slightly larger than QSG W bandgap, as seen in Table I.Also we converge the observables like band-gap and ∞ by increasing the size of the two-particle Hamiltonian within our self-consistent BSE implementation.We find that for CrBr 3 to converge both observables we find it necessary to include 24 valence bands and 24 conduction bands in two-particle Hamiltonian that we solve within BSE.Larger sizes of the two-particle Hamiltonian did not lead to  any changes in these observables.The convergence with respect to the number of states entering into the two-particle is significantly slower than in simple sp semiconductors.Once the two-particle Hamiltonian size is converged, we converge the observables in terms of vacuum size.Next we examine independent variations of the Hartree potential, via the density ρ, and the self-energy Σ 0 .In the GW case, Σ 0 denotes the quasiparticlized version of the dynamical self-energy Σ(ω); for the LDA it denotes the LDA exchange-correlation potential.Unless stated otherwise, results are presented with ρ and Σ 0 internally self-consistent.Considering this case at first, there is a remarkable difference between the LDA and QSGW electronic band structures.Within LDA (see Fig. 3(a)), the valence band maximum falls at the M point, while within QSGW it shifts to the Γ point (see Fig. 3(b)).The eigenfunctions are also quite different: the Br contribution to the low energy valence band manifold is significantly larger within QSGW.At a still higher level of theory replacing W → W , a portion of the strong perturbation of the LDA band structure is partially undone (Fig. 3(c)); shifts Br contribution to the valence eigenfunctions in the direction of the LDA (see table I).This is readily understood as a softening of W by the ladder diagrams, as noted above.With the QSG W , the bandgap in CrBr 3 is reduced slightly to 4.65 eV.The top most valence band in QSG W has a shape similar to LDA but the band gap is approx.three times as large and the valence bandwidth gets renormalized by a factor of ∼ 2. The observation that QSG W more closely resembles the LDA than QSGW is remarkable and calls for further analysis.
To this end, we consider independent variations of ρ and Σ 0 in the following senses: • ρ from LDA and Σ 0 from QSGW, which we denote as Σ QSGW [ρ(LDA)].This scheme produces a valence band structure similar to LDA (see Fig. 3(d)), but with 6.0 eV electronic gap, close to the QSGW gap of 5.7 eV.This clearly establishes the important role of the density in determining the effective one-body hamiltonian.As in the case of simple sp tetrahedral semiconductors where the LDA density is already rather good, the gap change is mainly controlled by the nonlocality in the self-energy which the LDA misses [56,57].
• ρ from QSG W and Σ from QSGW, which we denote as Σ QSGW [ρ(QSG W )] (see Fig. 3(e)).Now the valence band structure is much closer to the QSGW band structure, although the top most valence band is significantly narrowed.Also, the gap is similar to QSGW gap (5.7 eV).This tells us that the Hartree and many-body contributions cannot be decoupled.The addition of electron-hole ladder diagrams should considerably improve on the RPA's known inadequacy in describing short-ranged correlations [58], and here we see that it affects both Hartree and exchangecorrelation parts.
• ρ from QSGW and Σ from QSG W , which we denote as Σ QSG W [ρ(QSGW )] (bottom right panel of Fig. 3(f)).This shows in a different way how the Hartree and correlation contributions to the potential are interwined.
To further probe the role of the charge density, we plot ρ in the planes passing through the Cr and Br atoms at different levels of theory (see Fig. 4).The density is plotted in real space, and the abscissa and ordinate are defined by the the inverse transpose of the 2×2 matrix composed of b 1 and b 2 (see Fig. 1) with x and y defined by aligning b 2 parallel to y.In this notation the M point is on the b 2 line, or the y axis.On formation of the 2D crystal charge is augmented on the Cr-Cr and Br-Br bonds, taking it away from the atoms.QSGW accentuates this tendency (see Fig. 4(a,b)), as does QSG W , but to a relatively lesser extent (see Fig. 4(c,d)).However, although the structure of the top most valence band seems similar within LDA and QSG W and different within QSGW and QSG W , we show in Fig. 4(e,f) that the real-space ρ(QSG W ) is much closer to ρ(QSGW ) compared to LDA.In short, QSG W weakly modifies and slightly localizes charges in comparison to QSGW.
In conclusion, using a self-consistent first principles Green's function approach we show how correlations induce large changes in both the one-body (Hartree) and many-body contribution to the potentials, and that the two are inherently intertwined.To demonstrate the effect we considered two currently popular materials systems: a three dimensional charge-density-wave candidate TiSe 2 and a two-dimensional ferromagnet CrBr 3 .Such changes to the electronic wavefunction go way beyond any weak renormalization of the parameters that determine the electronic structure within a second principles approach and thus calls for development of better first principles approaches that solve many-body Hamiltonians for real materials with better approximations.MIK, ANR and SA are supported by the ERC Synergy Grant, project 854843 FASTCORR (Ultrafast dynamics of correlated electrons in solids).MvS and DP are supported by the National Renewable Energy Laboratories.We acknowledge PRACE for awarding us access to Irene-Rome hosted by TGCC, France and Juwels Booster and Clusters, Germany, STFC Scientific Computing Department's SCARF cluster.

METHOD
Single particle calculations (DFT, and energy band calculations with the static quasiparticlized QSGW self-energy Σ 0 (k)) were performed on a 16×16×16 (TiSe 2 ) and 16×16×1 (CrBr 3 ) k -mesh while the (relatively smooth) dynamical self-energy Σ(k) was constructed using a 8×8×8 (TiSe 2 ) and 6×6×1 (CrBr 3 ) k -mesh and Σ 0 (k) extracted from it.For each iteration in the QSGW self-consistency cycle, the charge density was made self-consistent.The QSGW cycle was iterated until the RMS change in Σ 0 reached 10 −5 Ry.Thus the calculation was self-consistent in both Σ 0 (k) and the density.Numerous checks were made to verify that the self-consistent Σ 0 (k) was independent of starting point, for both QSGW and QSG W calculations; e.g. using LDA or Hartee-Fock self-energy as the initial self energy for QSGW and using LDA or QSGW as the initial self-energy for QSG W .
For the present work, the electron-hole two-particle correlations are incorporated within a self-consistent ladder-BSE implementation [27,54] with Tamm-Dancoff approximation [59,60].The effective interaction W is calculated with ladder-BSE corrections and the self energy, using a static vertex in the BSE.G and W are updated iteratively until all of them converge and this is what we call QSG W . Ladders increase the screening of W, reducing the gap besides softening the LDA→QSGW corrections noted for the valence bands.
For CrBr 3 , we checked the convergence in the QSG W band gap by increasing the size of the two-particle Hamiltonian.We increase the number of valence and conduction states that are included in the two-particle Hamiltonian.We observe that for all materials the QSG W band gap stops changing once 24 valence and 24 conduction states are included in the two-particle Hamiltonian.While the gap is most sensitive to the number of valence states, 14 conducting states produces results within 2% error of the converged results from 24 conduction states.

FIG. 2 .
FIG. 2. Energy bands of the undistorted P 3c1 structure.(a): solid lines are LDA results, with red and green depicting a projection onto Ti and Se orbital character, respectively.Blue dashed line shows shifts calculated in the GW approximation based on the LDA, as described in the text.(b): blue dashed line shows results from GW based on the LDA (same Σ as in panel (a)), with an extra potential ∆V LDA deriving from a ρ shift computed from the rotation of the LDA eigenvectors.Solid lines are QSGW results, with the same color scheme as in panel (a).

FIG. 3 .
FIG. 3. Electronic Band Structure of CrBr3 : (a) LDA, (b) QSGW and (c) QSG W . W , we mean the polarizability that enters into W is calculated beyond the time-dependent Hartree approximation (RPA) but also includes ladder diagrams computed from a four-point Bethe-Salpeter equation (BSE).This significantly improves W as seen by comparing the macroscopic dielectric function to experiment.Colors correspond to Br-px+py (red), Br-pz (green), Cr-d (blue).A fourth color (white) selects the Cr-eg manifold, and washes out the color to the extent it is present.Topmost valence bands of LDA and BSE have similar shape, apart from a strong narrowing the bandwidth.The shape changes in QSGW.with the VBM shifting to Γ.(d) Σ QSGW [ρ(LDA)], (e) Σ QSGW [ρ(QSG W)], and (f) Σ QSG W [ρ(QSGW )].

FIG. 4 .
FIG. 4. Charge densities (ρ) of CrBr3: of the highest valence band state with the abscissa and ordinate x and y.All the left panels (a,c,e) pass through a Cr plane and right panels (b,d,f) pass through a Br plane.Panels (a,b) display constant-amplitude contours for QSGW ρ after subtracting out the LDA ρ.Contours are taken in half-decade increments in ρ, with a factor of 300 between highest contour (red) and lowest (blue).Panels (c,d) show the change in ρ passing from LDA to QSG W eigenfunctions; panels (e,f) show the corresponding change passing from QSGW to QSG W .In the bottom four panels (c,d,e,f), blue→red has a similar meaning as in the top panels (increasing positive δρ), while contours of negative δρ are depicted by increasing strength in the change blue→green.

TABLE I .
Fraction of spectral weight that the Halogen (Br) contributes to the total within an energy window of 0 (Fermi energy) to 0.6 eV below the Fermi energy (bound states).The Cr-d occupancies, and the electronic band gaps with different choices of self consistent ρ and Σ are shown.