LDA+DMFT approach to ordering phenomena and the structural stability of correlated materials

Materials with correlated electrons often respond very strongly to external or internal influences, leading to instabilities and states of matter with broken symmetry. This behavior can be studied theoretically either by evaluating the linear response characteristics, or by simulating the ordered phases of the materials under investigation. We developed the necessary tools within the dynamical mean-field theory (DMFT) to search for electronic instabilities in materials close to spin-state crossovers and to analyze the properties of the corresponding ordered states. This investigation, motivated by the physics of LaCoO$_3$, led to a discovery of condensation of spinful excitons in the two-orbital Hubbard model with a surprisingly rich phase diagram. The results are reviewed in the first part of the article. Electronic correlations can also be the driving force behind structural transformations of materials. To be able to investigate correlation-induced phase instabilities we developed and implemented a formalism for the computation of total energies and forces within a fully charge self-consistent combination of density functional theory and DMFT. Applications of this scheme to the study of structural instabilities of selected correlated electron materials such as Fe and FeSe are reviewed in the second part of the paper.


Introduction
The theoretical understanding of complex materials with strongly interacting electrons is one of the most challenging areas of current research in condensed matter physics. Experimental studies of such materials have often revealed rich phase diagrams originating from the interplay between electronic and lattice degrees of freedom [1][2][3]. These compounds are therefore particularly interesting in view of possible a e-mail: kunes@ifp.tuwien.ac.at arXiv:1705.00875v2 [cond-mat.str-el] 14 Aug 2017 technological applications. Namely, the great sensitivity of many correlated electron materials to changes of external parameters such as temperature, pressure, magnetic and/or electric fields, doping, etc., can be employed to construct materials with useful functionalities.
The electronic properties of materials can be computed from first principles by density functional theory (DFT), e.g., in the local density approximation (LDA) [4], the generalized gradient approximation (GGA) [5], or using the so-called LDA+U method [6]. Applications of these approaches describe the phase diagrams of many simple elements and semiconductors, and of some insulators quite accurately. Moreover, they often allow to make correct qualitative predictions of the magnetic, orbital, and crystal structures of solids where the equilibrium (thermodynamic) structures are determined by simultaneous optimization of the electron and lattice systems [7][8][9]. However, these methods usually fail to describe the correct electronic and structural properties of electronically correlated paramagnetic materials. Hence the computation of electronic, magnetic, and structural properties of strongly correlated paramagnetic materials remains a great challenge. Here the computational scheme obtained by the combination of DFT and dynamical mean-field theory (DMFT) [10][11][12][13], usually referred to as DFT+DMFT (or more explicitly as LDA+DMFT, GGA+DMFT, etc.) [14][15][16][17][18][19][20][21], provides a powerful new method for the calculation of the electronic, magnetic, and structural properties of correlated materials from first principles. The DFT+DMFT approach is able to give a good qualitative and quantitative description of the properties of correlated solids in both their paramagnetic and magnetically ordered states and demonstrates the crucial importance of electronic correlations in determining the properties of these multiband materials. In fact, by supplementing the DFT+DMFT approach with methods of quantum information theory it even becomes possible to quantify correlations in materials and to compare the correlation strength of different materials, e.g., transition-metal monoxides [22,23].
The DFT+DMFT scheme makes it possible to describe and explain the effect of finite temperatures, including thermally driven phase transitions, in real materials. By overcoming the limitations of conventional band-structure methods it opens the way for fully microscopic investigations of the structural properties of strongly correlated systems. So far the DFT+DMFT approach has been mainly used to calculate oneparticle and local two-particle observables needed to explain experimental results obtained, for example, by photoemission and x-ray absorption spectroscopies. The method has been widely applied to d and f transition metals and their oxides or Fe-based superconducting materials.
Employing a novel implementation of DFT+DMFT in combination with planewave pseudopotentials it was recently demonstrated that, by performing total-energy calculations for correlated materials, it is possible to compute atomic displacements, perform structural optimizations, determine the lattice dynamics, and explore structural transitions caused by electronic correlations [24][25][26][27][28][29]. Thereby the DFT+DMFT computational scheme is able to treat electronic and structural properties of strongly correlated materials on the same footing.
By combining DFT+DMFT with linear response techniques the scheme can be employed to investigate both electronic and lattice instabilities of correlated materials. In particular, it is possible to search for electronic instabilities in an unbiased manner, i.e., without prior assumptions about the broken symmetry of the ordered phase. Lattice instabilities can be explored by evaluating a complete set of the forces acting on the atoms or the phonon dispersion relations.
In this article we review some of our results obtained during the funding period of the DFG Research Unit FOR 1346 for electronic and lattice instabilities. The paper is organized as follows: In Sec. 2 we discuss the DMFT treatment of electronic (particle-hole) instabilities and states with spontaneously broken symmetry. As a pilot problem we study materials close to a spin-state crossover. This already leads to a rich phase diagram with several competing instabilities. After a general introduction to the formalism in Sec. 2.2.1, we present results for the spin-state crossover material LaCoO 3 in Sec. 2.3. In Sec. 2.4 we introduce a minimal two-orbital model of spinstate crossover and investigate its instabilities as well as broken-symmetry phases in Sec. 2.5.
In 3.1 we review the DFT+DMFT computational scheme for the computation of the electronic structure and phase stability of correlated materials. In particular, we present a detailed formulation of the fully charge self-consistent DFT+DMFT scheme implemented with plane-wave pseudopotentials. This method allows one to explore structural transformations (e.g., structural phase stability) caused by electronic correlations. We discuss applications of the DFT+DMFT scheme to the study of the electronic structure and phases stability of correlated materials, such as elemental iron in Sec. 3.2 and the parent compound of the Fe-based superconductors, FeSe, in Sec. 3.3.
In Sec. 3.4, we provide a detailed derivation of the DFT+DMFT method implemented within linear-response theory with respect to atomic displacements. It has been shown that this approach allows one to evaluate a full set of forces acting on the nuclei and thereby perform a structural optimization of the lattice. We apply this new technique to study the structural phase stability of a particularly simple test material, namely solid hydrogen. An outlook is provided in Sec. 4.

Electronic instabilities
There are two general theoretical methods to explore the existence of long-range ordered phases with spontaneously broken symmetries: Method 1 : investigate a particular kind of long-range order which is assumed from the beginning, and Method 2 : investigate the instabilities of the normal phase.
The first method has the obvious drawback that the broken symmetry of the ordered phase needs to be anticipated. This concerns in particular the translational symmetry, which is assumed in all calculations for extended systems. A unit cell which can host the ordered state is required, which may be difficult to achieve, or even impossible in the case of incommensurate order. Apart from the translational symmetry, care may also be needed to remove symmetries that might be present implicitly in the computational scheme or in the initial conditions for iterative methods such as DMFT.
In particular, the two-band Hubbard model proved to be simple enough to be computationally tractable, while exhibiting surprisingly rich physics. Kuneš and collaborators found various ordered phases of this model which can be described either as spin-state ordered states or excitonic condensates [59][60][61][62][63]. An exciton is a bound state between an electron and a hole with total charge zero. In this article we focus on excitons with spin S = 1. Due to their bosonic nature excitons can Bose condense (exciton condensation), which gives rise to anomalous matrix elements of the one-particle density matrix. In the case of excitons with spin S = 1 the condensation breaks spin isotropy, while time-reversal symmetry is not necessarily broken. Related results were reported by Kaneko et al. using the variational cluster approximation [64,65] and Vanhala et al. using cellular DMFT [52]. Other examples of unconventional ordered states obtained within DMFT are the breaking of the channel symmetry in the twochannel Kondo lattice reported by Hoshino et al. [66] and the spin/orbital freezing scenario for cuprate superconductors discussed by Werner et al. [67].
While the DMFT formalism needed for the study of broken-symmetry phases is not fundamentally different from that for normal (e.g., paramagnetic) phases, specific implementations, impurity solvers and post-processing tools need modifications. For example, the calculations in the exciton condensate reported below require the possibility of off-diagonal hybridization functions. For this purpose we modified the segment CT-HYB QMC code [68,69] to allow for a real off-diagonal hybridization.
As a minor methodological development we briefly mention a trick adapting the maximum entropy formalism [70] for indefinite (off-diagonal) spectral densities. It consists in writing the off-diagonal spectral function in the form A(ω) = A + (ω) − A − (ω), where A ± (ω) ≥ 0, which requires only minor code modifications. While this decomposition is highly non-unique we found the stability of A(ω) to be satisfactory. An alternative approach working with an indefinite A(ω) was proposed in Ref. [71].

Method 2
The investigation of instabilities of the normal phase removes the bias towards assuming a broken symmetry of the ordered phase. The price to be paid is that only instabilities, i.e. tendencies to ordering, can be identified, while the detailed properties of the ordered phases and their possible instabilities are not accessible. Furthermore, when several susceptibility modes diverge simultaneously, the symmetry of the ordered state is not uniquely determined by the susceptibility alone. It may be determined by higher-order terms in the expansion of the free energy, see Ref. [61] for an example. Calculations in linear-response are more demanding since correlation functions need to be evaluated which are not needed for the DMFT self-consistent cycle.

Computation of susceptibilities within DMFT
The response to a weak external field is one of the most important characteristics of any physical system. For example, it yields information about the system's stability and possible phase transitions. For fields which couple to single-particle operators the response of the many-body system is described by susceptibilities, i.e., two-particle correlation functions. Their evaluation is a difficult and in general unsolved problem. Computations in the framework of DMFT are simplified by two facts: (i) susceptibilities do not enter the self-consistent solution of the DMFT equations, i.e., can be computed after a converged DMFT-solution is obtained, (ii) the essential building block, the vertex function, is local. Nevertheless, even with these simplifications the calculation of susceptibilities remains a formidable numerical task, which has been performed only in special cases so far. One example is the evaluation of equal-times correlations such as the double occupancy, which is a standard part of DMFT implementations. However, these quantities do not contain information about the dynamics and the excitations of the system. Studies of the two-particle dynamics are usually limited to the local (impurity) spin and charge susceptibilities, which can be determined directly from the solution of the impurity problem [72]. Since they describe the response to a local applied field, they do not provide information about instabilities towards long-range order. Such information is contained in q-dependent susceptibilities.
When searching for (particle-hole) instabilities within the DMFT formalism one is interested in static susceptibilities of the type where q is a vector from the first Brillouin zone. In general, only the correlation functions with orbitals i and j belonging to the same atom (and similarly for k and l) are relevant. Due to the locality of the two-particle irreducible vertex [11] such functions appear on both sides of the Bethe-Salpeter equation and may cause a divergence. The DMFT order parameter is therefore local, but can exhibit a modulation in space (divergence at q = 0). The susceptibility (1) is obtained from the two-particle correlation functionχ ij,kl (q; ω m , ω n ) which is the solution of a pair of Bethe-Salpeter equations [11,73] Hereχ pq,kl (ω, ω ) is the impurity two-particle correlation function, andχ 0 ij,kl (ω, ω ) andχ 0 ij,kl (q; ω, ω ) are the local and lattice "bubbles", respectively, which are constructed from the lattice Green functions. Truncated summations over the Matsubara frequencies and the treatment of the high-frequency tails are discussed in Ref. [74].
DMFT calculations within linear response were performed for one-orbital models in the early days of DMFT [73,75,76]. It is only quite recently that linear response has been used for an unbiased search for instabilities in two-and three-orbital Hubbard models [53,54,67,[77][78][79]. The generalization of the above formalisms to dynamical susceptibilities is straightforward, but involves an analytic continuation of the bosonic frequency. The computation of dynamical susceptibilities within DMFT was reported for simple models in Ref. [80] and, with a simplified RPA-like vertex, even for real materials [81]. The degeneracy of atomic states is a common cause for electronic instabilities in strongly correlated systems. The spin degeneracy of a ground-state atomic multiplet in magnetic insulators is a typical example. The vicinity of a spin-state crossover, where several atomic multiplets become quasi-degenerate, may give rise to complex ordering phenomena. Their physics is the subject of the following sections.
A crossover between the low-spin state (LS) and the high-spin state (HS) ("spinstate crossover") is essentially an atomic effect -a consequence of the competition between the Hund's rule coupling J and the crystal field (CF) splitting ∆. Varying these parameters may lead to a level crossing [82] and thus change the single-ion ground state. The CF splitting is typically controlled by external or chemical pressure. Spin-state crossovers were observed in many oxides of transition metals from the middle of the periodic table such as MnO, Fe 2 O 3 , FeO, CoO, which were theoretically studied with the LDA+DMFT approach [34,[83][84][85][86][87]. In bulk materials it usually gives rise to a smooth crossover or a first order transition, and often involves a sizeable change of the specific volume. For most materials pressures in the range of tens of GPa are required to induce a spin-state crossover, which corresponds to variations of the CF on the scale of several 100 meV.
LaCoO 3 is an interesting exception. In this material the parameters are finetuned in such a way that a (partial) spin-state crossover can be studied by varying the temperature. In fact, the strongly temperature dependent magnetic and transport properties of LaCoO 3 and related compounds have attracted much attention and have been studied already for half a century [88][89][90]. Below 50 K LaCoO 3 appears to be a band insulator. However, above 100 K it exhibits a magnetic response typical for local moments while the charge gap continuously disappears between 450 and 600 K. This suggests that the material is much more than an ordinary band insulator.
Traditionally several different approaches have been employed to explain the physics of LaCoO 3 : (i) the single-ion picture of a LS, S = 0, ground state of the Co 3+ ion, with intermediate spin (IS), S = 1, or HS, S = 2, excitations augmented by spin-exchange between these states on the lattice, (ii) band structure approaches with electrons interacting via a static mean-field, and (iii) a combination of both in terms of DMFT. The technically simplest approach, (i), has been used extensively to interpret the experimental T -dependence of the magnetic susceptibility and specific heat. It describes the physics of a thermally induced statistical mixture of different atomic multiplets, but does not capture the extended nature of a bulk material. Models including IS, HS or both can be found in the literature [91][92][93][94]. Generally, it is not possible to describe LaCoO 3 by a single-ion model with constant parameters. Band structure methods, (ii), have been used to describe the one-particle spectra in the low-T regime [95][96][97]. A Hartree-Fock-like LDA+U method has been employed to study stable spin states and their possible ordered patterns [98,99]. However, these methods cannot capture the temporal fluctuations between atomic multiplets and are by construction limited to electrons at T = 0. The physical temperature enters the calculations only via the thermal expansion of the lattice. The LDA+DMFT approach, (iii), includes the best features of the approaches (i) and (ii), and combines the extended nature of the systems with local electronic correlations [100][101][102][103].
Křápek et al. [101] used LDA+DMFT to study the temperature dependence of the atomic spin state and the one-particle spectra. In Fig. 1 their results are compared with the experimental photoemission spectra. As expected the results are quite sensitive to the choice of the parameters used in the computation, in particular the Hund's exchange J. The thermal expansion of the lattice has a non-negligible effect on the electronic properties, but by itself it cannot explain the experimental susceptibility. On the other hand, even with a rigid lattice one could observe the effect of thermal population of the spinful atomic states. Perhaps the most important result is the negligible contribution of the IS state, irrespective of the values of the interaction and double-counting parameters. This result is consistent with single atom calculations [82,105], which yield the LS-HS or HS-LS sequence of the multiplet states in the vicinity of the spin-state crossover, but never a LS-IS sequence with a small gap.  [104]. The LS solution corresponds to T = 580 K and the HS solution was obtained for T = 1160 K. The measurements were taken at 65 K (denoted as LS) and 300 K (denoted as HS+LS, since the temperature is not high enough for the full spin-state crossover). Adopted from Ref. [101] with the permission of the authors. Recently, Sotnikov and Kuneš [63] proposed a mechanism in which the IS states still play an important role in LaCoO 3 . Starting from a global LS ground state, the IS states can be viewed as tightly bound excitons which carry spin S = 1 and appear in three orbital flavors. A perturbation expansion in the hopping [63] shows that the excitons have a rather high anisotropic mobility, i.e., each orbital flavor moves predominantly in one of the three cubic planes. The HS states are then viewed as tightly bound bi-excitons with total spin S = 2 formed by two IS states with different orbital flavors. The estimated mobility of HS bi-excitons is much smaller than that of IS excitons. The proposed excitation spectrum of LaCoO 3 at low temperature is shown in Fig. 2. The energies of atomic multiplets follow the LS-HS-IS sequence. However, the dispersion of IS excitons causes the lowest excitation on the lattice to be an IS wave with a specific crystal momentum q.
This picture is clearly beyond the capabilities of the single-ion model, (i), as well as band structure approaches, (ii). The question is then whether DMFT is able to capture the proposed behavior. Although the answer turns out to be negative, it can provide useful insights. Roughly speaking DMFT sees the IS excitation at its atomic q-averaged energy, i.e., substantially above the HS state. It does not take into account that in parts of the Brillouin zone propagating excitons have substantially lower energy. This explains the results of Křápek et al. [101]. A linear response calculation within DMFT performed on top of the low-temperature LS solution should be able to identify the dispersive excitonic branch -at least its low-energy part that can be distinguished from the particle-hole continuum. Such a calculation, while possible in principle, will be quite demanding and has yet to be done. Altogether one may say that DMFT can describe the dispersion of a single IS excitation on top a LS state, but cannot describe the state in which such (strongly interacting) excitations are thermally populated.
Experimentally, an observation of the low-energy excitonic branch in the lowtemperature (LS) state of LaCoO 3 based on two-particle spectroscopy, e.g., resonant x-ray scattering, offers a way to test the proposed scenario. It is interesting to note that on the one-particle level stoichiometric LaCoO 3 at low temperature appears to be an uncorrelated band insulator with sharp quasi-particle bands [101]. The "hidden" correlations in the material become manifest in the one-particle spectra only upon heating or doping.

Two-orbital Hubbard model
In the investigation of electronic instabilities in the vicinity of a spin-state crossover, such as in cobaltites, there are two reasons not to start with realistic models that include the full d shell. The first one is purely technical. Linear-response calculations are very demanding and their feasibility for the full d shell with the present (still developing) codes is questionable. The second reason has to do with the interpretation of the results. As shown by recent LDA+U calculations [106], the excitonic order parameter in d 6 systems has 18 components, and numerous almost degenerate phases are possible. We find it therefore necessary to start with simpler systems and gradually increase the complexity.
The two-orbital (i.e., two-band) Hubbard model (2BHM) provides such a minimal microscopic model to study the competition between CF splitting and Hund's rule exchange without the orbital degeneracy of realistic d shell. The Hamiltonian reads with the notationσ = −σ. It describes electrons with two orbital (a and b) and two spin (σ =↑, ↓) flavors moving on a lattice and interacting via an on-site interaction H int . Here, a † iσ (a iσ ) is an operator creating (annihilating) a fermion with orbital flavor a and spin σ on the lattice site i, n a iσ = a † iσ a iσ is the corresponding local density operator, and analogously for the b fermions. The sum ij runs over the nearest neighbour (nn) bonds. Here we consider not only the hopping between orbitals with the same flavour (hopping amplitude t), but explicitly include also the cross-hopping between orbitals with different flavors (hopping amplitude V ). The parameter ∆ describes the CF splitting. The standard Slater-Kanamori interactions are parametrized by U and J (with U = U −2J, J = J). Numerical quantum Monte-Carlo calculations are greatly simplified when the interaction is approximated by a density-density interaction, whereby spin-flips (the second contribution in the J term) and pair-hopping J are neglected. This approximation will be used here.
In the following we will discuss the 2BHM at, and close to, half-filling. Let us first summarize the basic properties of its normal phase. The physics at strong and intermediate coupling is controlled by the competition between the Hund's coupling J and the crystal-field splitting ∆ [107,108]. Large ∆ favors the singlet LS state, while large J favors the triplet HS state. When spontaneous symmetry breaking is excluded the ∆ − U phase diagram at fixed J/U can be divided into three regions: HS Mott insulator, LS band insulator, and a metallic state ("metal") [107]. The firstorder metal-insulator transition turns into a crossover at higher temperatures. The low-energy physics deep in the Mott phase is described by the S = 1 Heisenberg model with antiferromagnetic interaction. The LS band insulator far from the phase boundaries is a global singlet with a gapped excitation spectrum. In the vicinity of the HS-LS crossover both LS and HS states have to be taken into account. The near degeneracy of the atomic multiplets in this region gives rise to several instabilities.
Kuneš and Křápek [59] reported spin-state order (a checker-board pattern of HS and LS states) on a square lattice in the 2BHM with very asymmetric (narrow/wide) bands. An interesting feature of this transition, which can be well explained by the classical Blume-Emery-Griffiths model [109], is its reentrant character in parts of the phase diagram [110,111]. The above (classical) order is not the only instability of the model. Magnetic ordering of the HS states and condensation of spin-triplet excitons [54,64,77,[112][113][114][115][116] are competing instabilities. In the following, we focus on the latter one.

Spin-triplet exciton condensation
Long-range order in materials with singlet atomic ground states is rare. If the energy of an atomic excitation is comparable to its amplitude to propagate to neighbouring atoms, spontaneous symmetry breaking may take place. This mechanism, called exciton condensation or excitonic magnetism, was recently proposed to be realized in 4d 4 materials such as Ca 2 RuO 4 [117,118] or 3d 6 cubic materials such as the Pr 0.5 Ca 0.5 CoO 3 family [60]. Excitonic magnetism has been studied in simple lattice models such a 2BHM or effective pseudospin models [60,64,65,77,117,119]. In the following we review the DMFT results for the spin-triplet exciton condensate.

Linear response
The linear-response formalism described in Sec. 2.2.1 (Method 2) allows for an unbiased search for instabilities. Recently, several authors [54,77,120] applied this approach to the 2BHM in the parameter range close to the triple point between the metal and the HS and LS insulators. In Fig. 3 we show a typical result of such a calculation. The generalized particle-hole susceptibility exhibits three distinct modes with sizeable amplitudes, which correspond to the spin susceptibility, orbital susceptibility, and spin-triplet excitonic susceptibilities, respectively [121]. These describe  the response of the system to the fields with site indices i, spin indices α, β =↑, ↓ and Cartesian indices µ = x, y, z of the Pauli matrices σ µ . Note that in the actual calculations with the density-density interaction, it is only the spin s z i and the excitonic c x i , c y i , d x i , d y i components that lead to a finite response. The U (1) spin rotation symmetry implies a degeneracy of the x and y excitonic modes. In the absence of cross-hopping (and pair-hopping) the spin-density d µ i and spin-current-density c µ i modes are degenerate. In Fig. 3 both the orbital and the excitonic susceptibility diverge at q = (π, π) when the temperature is decreased, with the orbital instability leading. This divergence is a signature of the checker board instability discussed in the previous section. Varying the ratio of the bandwidths one can switch between the orbital and excitonic instability. The orbital susceptibility is favored by strongly asymmetric bands, while for a moderate asymmetry the exciton condensation is the leading instability.

Excitonic condensate
In the remainder of this section we study the properties of the exciton condensate (Method 1). We choose the model parameters in the region dominated by the excitonic instability (arrow in Fig. 3), but unlike in Fig. 3 we choose, for computational Energy (eV) Next, we address the Pr 3+ →Pr 4+ valence transition and the fact that experime is observed in a doped system. An isostructural valence transition points to a ne of the two charge states of the Pr 4f shell, which therefore acts as a charge reservo CoO 3 subsystem at a fixed chemical potential rather than fixed particle density.
have repeated the calculations at fixed chemical potential. In Fig. 2 we show th dependences of the order parameter |φ| and the particle density n, which reflects t valence in the real material. Doping the system away from half filling leads to re While n(T ) is almost constant above T c , below T c the system draws particles from t approach the half filling in a process controlled by gain in the condensation energy a cost of adding electrons. This agrees well with the experimental behavior of the Pr v in experiment and in the model the average Pr valence changes rapidly, but continuou While theoretical n(T ) starts growing precisely at T c its experimental counterpart convenience, t a t b < 0 to get a uniform condensate (q = 0) [111]. The condensate is characterized by a finite expectation value φ i = αβ σ αβ a † iα b iβ = d i +ic i . It is instructive to consider the condensate wave function in the strong-coupling limit, which can be approximated as a product i (s|LS i + ξ|HS i ) of coherent superpositions of LS and HS states. Depending on the HS component we can obtain condensates with distinct symmetries and physical properties. In particular, HS states with non-zero expectation value of the spin operator S i = 0, e.g., |S z = 1 give rise to a ferromagnetic condensate (FMEC) with finite, ferromagnetically ordered spin moments, while HS states with S i = 0, e.g., |S z = 0 or 1 √ 2 (|S z = 1 + |S z = −1 ), lead to a condensate with vanishing ordered moments which we call polar excitonic condensate (PEC). Which of these states is realized depends on the details of the model in Eq. (5). While FMEC and PEC states have many different properties they can be distinguished solely by the structure of the corresponding order parameter φ i [61].
In addition to the bosonic physics of the condensate, the 2BHM describes also fermionic excitations and their interaction with the condensate. In Fig. 4 we show the evolution of the one-particle spectral function for the half-filled 2BHM when, starting from a semi-metallic normal state, the system is cooled through the EC instability. For the stoichiometric filling the system adopts the PEC phase, a preference that originates in the anti-ferromagnetic super-exchange between the HS states [111]. Below the transition temperature a sizeable gap ( k B T ) opens. The similarity to an s-wave superconductor is not accidental. There is a mapping between the excitonic insulator and a superconductor, which dates back to the 1960's [54,112] and involves a similar algebra for both phenomena. Another notable property of the PEC phase is a temperature independent spin susceptibility [60].

Doping and cross-hopping
Next, we consider the effect of doping on the exciton condensate. The DMFT phase diagram of a doped excitonic insulator is shown in Fig. 5. A destructive effect of doping on the exciton condensation was observed experimentally [122] and can be traced back to the competition between the condensation and the kinetic energy. In the present model the excitonic instability vanishes for more that 0.17 holes per atom. However, before that happens a transition to the FMEC phase takes place. At low temperatures the transition is of first order and is accompanied by charge separation, while at higher temperatures it proceeds via two continuous phase transitions. The microscopic origin of the FMEC phase has been explained in terms of a generalized double-exchange mechanism [62]. The spin of the doped carrier couples to the net spin polarization of the condensate. At low doping, the antiferromagnetic super-exchange between HS states prevails and the system settles in the PEC state. At higher doping, the spin polarized condensate is more transparent for the doped carriers and the gain in the kinetic energy stabilizes the FMEC state. Finally, we discuss the effect of moderate cross-hopping |V ab V ba | |t a t b |. We consider two possibilities for the cross-hopping: V ab = V ba (even) and V ab = −V ba (odd) [62]. The corresponding phase diagrams, shown in Fig. 6, exhibit the same basic features as their "parent" system in Fig. 5 with several notable differences. First, the PEC order parameter, which has an arbitrary complex phase (at least on the meanfield level) for V ab = 0, becomes purely real (V ab V ba < 0) or imaginary (V ab V ba > 0). The corresponding states are characterized by a local spin polarization with vanishing dipole but finite multipole moment (referred to as spin-density-wave (SDW) phase in the terminology of [112]), or by a local pattern of spin currents (referred to as spincurrent-density-wave (SCDW) phase). Second, a continuous phase transition from the normal phase directly to the FMEC phase is possible only via an intermediate PEC phase (primed phases in Fig. 6).
All the SDW, SCDW, SDW , and SCDW phases of Fig. 6 are instances of the PEC state. Nevertheless, the primed and unprimed phases have distinct properties. The primed phases were shown to host spin textures in the reciprocal space [62]. In particular, textures with odd k-parity are interesting since they do not break time reversal

Structural stability of correlated materials
In this section we present a detailed formulation of a variant of the DFT+DMFT computational scheme [10-12, 14, 15, 17, 18], which allows one to compute atomic displacements and thereby detect correlation-induced structural transformations, determine the phase stability, calculate the lattice dynamics, and perform a structural optimization of correlated electron materials [24][25][26][27][28][29].
We also present a detailed formulation of a fully charge self-consistent DFT+DMFT scheme, which includes the effect of electron correlations on the charge density [123][124][125][126][127]. The correlations lead to a charge redistribution between the correlated 3d or 4f orbitals and non-correlated s, p states, which has been shown to play an important role, e.g., near the Mott metal-insulator transition [29].
The DFT+DMFT approach overcomes the limitations of standard band-structure techniques and allows one to perform microscopic investigations of the electronic structure and lattice properties of correlated electron materials. It has been widely employed in recent investigations of the electronic and structural properties of strongly correlated electron materials [22-29, 34, 43, 83-86, 128-159].
In particular, this method was used to study the electronic and structural properties of elemental Fe [27,28,135] and the iron chalcogenide FeSe [145], and to explain the cooperative Jahn-Teller effect in paramagnetic KCuF 3 and LaMnO 3 [24][25][26]. Recently it was employed to explore the electronic properties and the phase stability of paramagnetic V 2 O 3 [29] and the transition metal monoxides MnO, FeO, CoO, and NiO near the pressure-induced Mott-Hubbard metal-insulator transition [87].
We will demonstrate that the DFT+DMFT approach provides a qualitative and quantitative description of the electronic properties and phase stability of all these materials, in spite of their chemical, structural, and electronic differences. The scheme is robust and makes it possible to address, on the same footing, electronic, magnetic, and structural properties of strongly correlated materials.

Methodological developments: Total-energy functional and a fully charge self-consistent DFT+DMFT scheme
In order to compute the electronic properties and determine the phase stability of correlated electron materials, we implemented [29] a fully charge self-consistent DFT+DMFT approach [123][124][125][126][127], which is formulated in terms of plane-wave pseudopotentials [7,160]. Following Anisimov et al. [161] and Trimarchi et al. [162] we construct a basis set of atomic-centered symmetry-constrained Wannier functions [161][162][163][164], which allows us to build up an effective low-energy HamiltonianĤ DFT for the partially occupied correlated orbitals, e.g., d of f orbitals of transition-metal ions. The HamiltonianĤ DFT , which provides a realistic description of the singleelectron band structure of a material, is further supplemented by on-site Coulomb interactions for the correlated orbitals. In the density-density approximation for the electronic interaction (see Sec. 2.4) the many-body Hubbard Hamiltonian takes the formĤ wheren mσ =ĉ † mσĉmσ is the local density operator for the orbital m and spin σ. Here U σ,σ m,m and U σ,σ m,m are the reduced interaction matrices for equal and opposite spins, respectively. The interaction matrices are expressed in terms of the Slater integrals F 0 , F 2 , and F 4 , which for the case of 3d electrons are related to the local Coulomb and Hund's rule coupling as U = F 0 , J = (F 2 + F 4 )/14, and F 2 /F 4 = 0.625. The termĤ DC is a double-counting correction which accounts for the electronic interactions already described by DFT. The Coulomb repulsion U and Hund's coupling J can be evaluated using constrained DFT and/or constrained random phase approximation (RPA) methods within a Wannier-functions formalism [164]. The U and J values are related to the matrix notation as U = {m} U σσ mm /N 2 and where N is the number of correlated orbitals. We implement the fully charge self-consistent DFT+DMFT approach in terms of plane-wave ultrasoft pseudopotentials [7,160]. In this scheme, the matrix elements of H DFT are determined as where P σ iν (k) = ψ σ ki |Ŝ|φ σ kν are the matrix elements of the orthonormal projection operators expressed in the basis of local orbitals φ σ kν . The charge density ρ(r) is calculated as where T is the temperature. The charge density matrix elements ρ k;ij in the basis of the Kohn-Sham wave functions ψ ki are computed as Here, the index I refers to an atomic site, the functions β I m (r) are computed in an atomic calculation and differ for different atomic species. Furthermore, Q I lm (r) is the augmentation function that is strictly localized in the core region, and G k;ij (iω n ) is the lattice Green's function in the KS wave functions basis at a given k-point. The (inverse) lattice Green's function is evaluated as where σ ki are the Kohn-Sham eigenvalues. Σ σ ij (k, iω n ) is the self-energy obtained from the solution of the DMFT impurity problem by "upfolding" of the impurity self-energy from the localized Wannier to the Kohn-Sham basis. We note that for the non-correlated states which do not hybridize with the correlated orbitals the density matrix collapses to the Fermi function for the state k and band i, i.e., T iωn G k;ij (iω n )e iωn0+ = f ki δ ij . In practice, it is convenient to compute the charge density as ρ(r) = ρ DFT (r) + ∆ρ(r), i.e., to split the DFT contribution and the correlation-induced difference ∆ρ(r). Full charge self-consistency is achieved when both the charge density ρ(r) and the local Green's function are converged.
Within the DFT+DMFT calculation the total energy is then evaluated using the expression [24][25][26]131] where E DFT [ρ(r)] is the DFT total energy obtained for the self-consistent charge density ρ(r). The third term on the right-hand side of Eq. (15) is the sum of the DFT valence-state eigenvalues which is evaluated as the thermal average of the DFT Hamiltonian with the non-interacting DFT Green's function G DFT k (iω n ): Ĥ DFT is evaluated in a similar way but with the full Green's function including the self-energy. To calculate these two contributions, the summation is performed over the Matsubara frequencies iω n , taking into account an analytically evaluated asymptotic correction. Thus, for Ĥ DFT one has where the first moment m k 1 is computed as m k 1 = H DFT (k) + Σ(i∞) − µ. The asymptotic part of the self-energy Σ(i∞) is calculated as the average of Σ(iω n ) over the last several iω n points. The interaction energy Ĥ U is computed from the double occupancy matrix. The double-counting correction E DC is evaluated as the average Coulomb repulsion between the N d correlated electrons in the Wannier orbitals.
This DFT+DMFT approach allows us to determine correlation-induced structural transformations of correlated materials together with the corresponding change of the atomic coordinates and the unit cell shape. It can also be used to explain experimentally observed structural data and to predict structural properties of real correlated materials.
In particular, we performed fully charge self-consistent computations of the electronic structure and phase stability of paramagnetic V 2 O 3 across the Mott-Hubbard metal-insulator transition (MIT) [29]; see also Refs. [143,144]. Namely, in charge nonself-consistent calculations the a 1g -e π g crystal-field splitting of paramagnetic V 2 O 3 is found to be strongly enhanced [138][139][140][141][142]. This leads to a substantial redistribution of the charge density and thereby influences the lattice structure due to electron-lattice coupling. For that reason full charge self-consistency turns out to be crucial to obtain a more realistic description of the physical properties of V 2 O 3 [29] near the Mott metal-insulator transition. For a detailed discussion we refer to Ref. [29].
In the following we report our results obtained by the DFT+DMFT scheme for the electronic and structural properties of elemental Fe [27,28,135] and the iron chalcogenide FeSe [145].

Lattice dynamical properties of paramagnetic Fe
Elemental iron is an exceptionally important material even for present-day technology. Iron exhibits a rich phase diagram with at least four allotropic forms [165,166]. At ambient conditions it is ferromagnetic and has a bcc crystal structure (α iron). Upon heating above the Curie temperature T C ∼ 1043 K, α iron becomes paramagnetic, but remains in its bcc crystal structure. Only upon further increase of the temperature above T struct ∼ 1185 K α iron exhibits a structural phase transition to a fcc structure (γ phase). Under pressure α iron makes a transition to a paramagnetic hcp structure ( phase) at ∼ 11 GPa.

Lattice stability and phonon spectra near the α-γ structural transition
State-of-the-art band structure methods provide a qualitatively correct description of various electronic and structural properties of iron [167][168][169][170][171][172][173][174][175]. For example, these methods provide a good quantitative understanding of the equilibrium crystal structure and the lattice dynamical properties of the ferromagnetic α phase. However, applications of these techniques to describe, e.g., the α-γ phase transition in iron, do not lead to satisfactory results. They predict a simultaneous transition of the structure and the magnetic state at the bcc-fcc phase transition while, in fact, the bcc-to-fcc phase transition occurs only about 150 K above T C . Moreover, the elastic and dynamical stability of the bcc phase is found to depend sensitively on the value of the magnetization. For example, in the absence of the magnetization, standard band-structure methods predict bcc iron to be unstable [176]. We now understand that this is due to the presence of local moments above T C which cannot be treated realistically by conventional band structure techniques. This problem has been overcome by employing the DFT+DMFT approach which allows one to study correlated materials both in to the long-range ordered and paramagnetic state [27,28,43,[133][134][135][136][137]. The DFT+DMFT method naturally accounts for  [28]. The DFT+DMFT result is further interpolated using a Born-von Kármán model with interactions expanded up to the 5-th nearest-neighbor shell. The results are compared with neutron inelastic scattering measurements at 1173 K [177] and 1428 K [178] for the bcc and the fcc Fe phase, respectively. From Ref. [28] with permission of the authors. the existence of local moments above T C and has shown to provide a good quantitative description of the properties of α iron. Moreover, applications of DFT+DMFT to study the equilibrium crystal structure and phase stability of iron at the α-γ phase transition reveal that the bcc-to-fcc phase transition takes place at a temperature of about 1.3 T C , i.e., well above the magnetic transition, in agreement with experiment [27].
We determine the structural phase stability and lattice dynamics of paramagnetic iron at finite temperatures by employing the DFT+DMFT approach implemented with the frozen-phonon method [28,135]. The approach is implemented with planewave pseudopotentials [7,160] which allows us to compute lattice transformation effects caused by electronic correlations [24][25][26]. We employ this technique to study the temperature dependent phonon dispersion relations and phonon spectra of paramagnetic iron at the bcc-fcc phase transition. Our presentation follows the discussion of Ref. [28].
It was previously shown that band structure calculations within the nonmagnetic generalized gradient approximation (GGA) cannot explain the experimentally observed phase stability of paramagnetic iron at the bcc-fcc phase transition, since they do not describe electronic correlations adequately. We now include the effect of electronic correlations by constructing an effective low-energy Hamiltonian for the partially filled Fe s, d orbitals based on the results of the nonmagnetic GGA. We construct a basis of atomic-centered symmetry-constrained Wannier functions for the Fe s, d orbitals [161,162,164], with U = 1.8 eV and J = 0.9 eV as obtained by previous theoretical and experimental estimations. To solve the realistic many-body problem we employ the Hirsch-Fye algorithm without charge self-consistency.
The phase stability and lattice dynamical properties of iron near the bcc-fcc phase transition are computed within the DFT+DMFT approach implemented with the frozen-phonon method [28]. The phonon frequencies are calculated by introducing a small set of displacements in the corresponding supercells of the equilibrium lattice which results in an energy difference with respect to the undistorted structure. We first focus on the lattice dynamical properties of iron near the bcc-to-fcc phase transition and perform calculations at temperatures T = 1.2 T C and 1.4 T C , which are below and above the temperature T struct ∼ 1.3 T C where the structural phase transition occurs.
We present our results for the phonon dispersion relations and phonon spectra in Fig. 7. The computations are performed for the equilibrium volume calculated at this particular temperature (lattice constant a = 2.883Å, which almost coincides with experiment). To evaluate the phonon frequencies for arbitrary wave vectors in the Brillouin zone, we performed lattice dynamical calculations on the basis of a Bornvon Kármán model with interactions expanded up to the 5-th nearest-neighbor shell. The calculated phonon dispersions of the bcc phase of iron show the typical behavior of a bcc metal with an effective Debye temperature ∼ 458 K. The phonon frequencies are overall positive, implying mechanical stability of the bcc lattice structure at ∼ 1.2 T C , i.e., well above the Curie temperature, in agreement with experiment. This corrects the results obtained with the non-magnetic GGA which finds the bcc lattice to be dynamically unstable even for the equilibrium lattice constant a = 2.883Å. Most importantly, our calculations clearly demonstrate the crucial importance of electronic correlations to explain both the thermodynamic and the lattice dynamical stability of the paramagnetic bcc phase of iron. For details we refer to Ref. [28]. Overall, the structural phase stability, equilibrium lattice constant, and phonon frequencies of bcc iron obtained by DFT+DMFT are in remarkably good agreement with the experimental data which were taken at nearly the same reduced temperature T /T C [177].
We also study the lattice dynamical properties of the paramagnetic fcc phase which is found to become energetically favorable above ∼ 1.3 T C . To prove the mechanical stability of the fcc phase we evaluate the lattice dynamics of the fcc phase at T ∼ 1.4 T C . Our results for the phonon dispersion relations and phonon spectra, which were calculated for the equilibrium lattice constant a = 3.605Å are also shown in Fig. 7. The effective Debye temperature is obtained as ∼ 349 K. The phonon frequencies are overall positive, implying mechanical stability of the fcc lattice structure at T ∼ 1.4 T C . Although nonmagnetic GGA calculations also find the fcc lattice structure to be mechanically stable (for the GGA equilibrium volume), the GGA energy for fcc iron is higher than that for the close-packed hcp structure. By contrast, DFT+DMFT calculation find the simultaneous thermodynamic and lattice dynamical stability of the paramagnetic fcc phase of iron at ∼ 1.4 T C , in agreement with experiment. Altogether our results for the structural phase stability, equilibrium lattice constant, and phonon frequencies agree remarkably well with the available experimental data taken at nearly the same reduced temperature T ∼ 1.4 T C [178]. Again it should be noted that the application of the nonmagnetic GGA to fcc iron finds phonon frequencies which differ considerably from experiment. These findings clearly demonstrate the importance of electronic correlations for the lattice dynamical properties of fcc iron.

Origin of the phase stability of δ iron
At even higher temperatures, above 1670 K, iron is known to make a structural transition to the δ phase, which is stable up to the melting curve [165,166]. To explain the phase stability of δ iron, we calculate the temperature evolution of the phonon dispersions near the α-to-γ and γ-to-δ phase transformations within DFT+DMFT [135]. The results for the phonon dispersion curves at temperatures near the magnetic phase transition temperature T C are compared in Fig. 8 (left). Phonon dispersion curves obtained in the temperature range of 1.4 -1.8 T C are also compared in Fig. 8 (right).
The phonon frequencies are seen to be overall positive, both in the ferromagnetic and paramagnetic phase at around T C . This implies the structural stability of the bcc phase in accordance with experiment. In the given temperature range the calculated phonon dispersions show only a rather weak temperature dependence. We notice a good quantitative agreement between the calculated phonon dispersions of bcc iron The results are compared with neutron inelastic scattering measurements at 1043 K [177]. Right panel: Calculated phonon dispersions of paramagnetic bcc iron near the α-to-γ and γ-to-δ phase transitions for different temperatures. From Ref. [135] with permission of the authors. and the experimental data taken at T ∼ T C [177]. In particular, we find a remarkable anomaly in the transverse (110) acoustic mode (T 1 mode) along the Γ -N branch near the bcc-fcc phase transition temperature at ∼ 1.2 T C . As will be discussed below this behavior is a dynamical precursor of the bcc-to-fcc phase transition [179] which occurs at ∼ 1.3 T C .
Upon further temperature increase our results, which were obtained in the harmonic approximation, clearly show that the bcc phase becomes dynamically unstable at T ∼ 1.4 T C , i.e., above the α-to-γ phase transition temperature [135]. The origin of the instability lies in the T 1 mode in the Γ -N direction discussed above, whose frequency becomes imaginary near the N -point. At the same time other phonon dispersion branches remain stable and are only weakly temperature dependent near the transition point. At even higher temperatures (T ∼ 1.8 T C ) the bcc structure of the δ phase is dynamically unstable in the whole Γ -N branch, again due to the T 1 mode, with an additional anomaly near the ( 2 3 , 2 3 , 2 3 ) point. Therefore we expect that the T 1 transverse (110) phonon mode continuously softens upon temperature increase near the α-γ and γ-δ phase transition temperature -a feature which is typical for simple nonmagnetic metals. By contrast, at temperatures T < T struct the phonon dispersion modes of bcc iron are relatively rigid, and resemble those of non-polymorphic bcc metals such as Cr, Mo, W. This should be contrasted with the phonon spectra of the denser fcc (γ) phase of iron which shows only a rather weak temperature dependence for all temperatures.
An estimate of the lattice free energy, which is based on the quasi-harmonic method described in Ref. [180], shows that the strong anharmonicity due to the T 1 mode at high temperatures is sufficient to overcome the total energy difference between the bcc and the fcc phases, resulting in the phase stability of the bcc (δ) phase [135]. On the basis of these results we therefore conclude that the hightemperature bcc (δ) phase of iron is stabilized by the lattice entropy, which gradually increases upon heating due to the increasingly anharmonic behavior of the T 1 phonon mode.

Correlation-induced topological Fermi surface transition in FeSe
As a second application of the DFT+DMFT scheme discussed in this section we address the electronic and structural properties of FeSe. This compound can be regarded as the parent compound of the Fe-based superconductors [181][182][183][184][185][186][187]. It has been recently shown that the critical temperature T c of FeSe, which is about T c ∼ 8 K at ambient pressure [188], depends very sensitively on changes of the lattice structure of FeSe due to pressure or chemical doping. In particular, T c is found to increase up to ∼ 37 K under hydrostatic pressure of about 7 GPa and to ∼ 14 K upon chemical (isovalent) substitution with Te [189][190][191][192][193][194][195]. The angle-resolved photoemission and band structure calculations reveal that FeSe has the same Fermi surface topology as the pnictides [196][197][198][199][200][201][202][203][204][205]. It is characterized by an in-plane magnetic nesting wave vector Q m = (π, π), consistent with s ± pairing symmetry [206,207]. Moreover, both pnictides and chalcogenides display a strong enhancement of short-range spin fluctuations near T c , with a resonance at (π, π) in the spin excitation spectra [208][209][210][211], suggesting a common origin of superconductivity in pnictides and chalcogenides, e.g., due to spin fluctuations [206,207]. On the other hand, the related isoelectronic compound FeTe exhibits no superconductivity [193,[212][213][214], but shows long-range antiferromagnetic (π, 0) order below T N ∼ 70 K, suggesting that the solid solution Fe(Se,Te) should have a remarkable crossover from the (π, π) to (π, 0) magnetic behavior upon substitution Se with Te. In addition, FeTe exhibits a remarkable phase transition under pressure, from a tetragonal to a collapsed-tetragonal phase, with a simultaneous collapse of local moments, indicating that the solid solution Fe(Se,Te) is close to an electronic and/or lattice transition [215,216].
We will now discuss the origin of this surprising behavior as well as the properties of the Fe(Se,Te) solid solution, following the presentation by Leonov et al. [145]. Most recently these calculations were generalized to take into account full charge self-consistency [217]. It should be noted, however, that charge self-consistency is not essential here, i.e., qualitatively similar results are obtained for FeSe already within the non-charge self-consistent calculations reported earlier [145].
The electronic structure and phase stability of FeSe is computed as a function of lattice volume, employing the DFT+DMFT approach implemented with plane-wave pseudopotentials [24][25][26]. For the partially filled Fe 3d and Se 4p orbitals we construct a basis set of atomic-centered symmetry-constrained Wannier functions [161,162,164]. To solve the realistic many-body problem, we employ the continuous-time hybridization-expansion quantum Monte-Carlo algorithm [68,218]. The calculations are performed at three different temperatures: T = 290 K, 390 K, and 1160 K. In these calculations we use the average Coulomb interactionŪ = 3.5 eV and Hund's exchange J = 0.85 eV, in accord with previous estimates for pnictides and chalcogenides [40,41,126,[219][220][221][222][223][224][225][226][227][228]. We employ the fully-localized double-counting correction, evaluated from the self-consistently determined local occupancies, to account for the electronic interactions already described by GGA. To investigate the phase stability, we take a tetragonal crystal structure (space group P 4/mmm) with the lattice parameter ratio c/a = 1.458 and Se position z = 0.266, and calculate the total energy as a function of volume.

Total energy and fluctuating moments
Our results are presented in Fig. 9. In particular, we find the equilibrium lattice constant a = 7.07 a.u., which is within 1% of the experimental value [189,190]. The calculated bulk modulus is B ∼ 70 GPa, which is comparable with that for iron pnictides. Most importantly, our DFT+DMFT calculations predict a structural transition of FeSe upon ∼ 10 % expansion of the lattice volume. This result is unexpected and, in fact, is very different from that obtained with the non-spin-polarized GGA, where the bulk modulus comes out much higher. Once again the results obtained within DFT+DMFT demonstrate the crucial importance of electronic correlations for the 7 Fig. 9. Total energy (left) and fluctuating local moment (right) of paramagnetic FeSe as a function of lattice constant calculated within DFT+DMFT at a temperature of 490 K with (sc) and without (nsc) charge self-consistency. The total energy is compared with that obtained within the non-magnetic GGA ("nm GGA"). After Ref. [217]. electronic structure and phase stability of FeSe. Namely, the repulsive interaction leads to an increase of the unit cell volume and hence results in a reduction of the bulk modulus.
We interpret the structural transition as a collapsed-tetragonal (low-volume) to tetragonal (high-volume) phase transformation upon expansion of the lattice volume. The phase transition is accompanied by a strong increase of the fluctuating local moment m 2 z , which grows monotonically upon expansion of the lattice. The highvolume phase has a much larger local moment and a softer lattice with a much lower bulk modulus of 35 GPa.

Spectral function and Fermi surface structure
The results for the integrated spectral function are shown in Fig. 10 (left). In agreement with previous studies [41], we notice a remarkable reduction of the Fe 3d bandwidth near the Fermi energy caused by electronic correlations. The lower Hubbard band is located at about −1.5 eV for both phases. Upon expansion of the lattice, we find a substantial spectral weight transfer. In particular, the spectral function for the low-volume phase exhibits a well-defined quasiparticle peak located below the Fermi level at −0.19 eV, which is absent in the high-volume phase. We note that the peak originates from the van Hove singularity of the Fe xz/yz and xy bands at the M -point. Moreover, we find a substantial qualitative change in the self-energy upon expansion of the lattice, resulting in a significant orbital-selective renormalization of the Fe 3d bands (not shown here). The xz/yz and xy bands exhibit significantly stronger correlations than the z 2 and x 2 − y 2 bands. While in the low-volume phase the quasiparticle mass enhancement is moderate, ∼ 2.0 -2.5, our calculations for the high-volume phase yield a strong renormalization ∼ 4 for the xz/yz orbitals and ∼ 6 for the xy orbitals. This shows in particular that the effect of orbital-selective correlations increases upon expansion of the lattice.
Our results for the Fermi surface presented in Fig. 10 (right) reveal a complete reconstruction of the electronic structure upon expansion of the lattice, resulting in a dramatic change of the Fermi surface topology ("Lifshitz transition") [145]. In particular, the Fermi surface at the M -point collapses, leading to a large square-like hole pocket around the M -point in the high-volume phase, in surprising analogy with the cuprates. In addition, the hole pockets around the Γ -point transform into incoherent spectral weight at the Fermi level along the Γ -X direction. The reconstruction of the Fermi surface topology leads to a corresponding change of the magnetic correlations in FeSe. We find in-plane nesting with Q m = (π, π), connecting hole and electron parts of the Fermi surface, to be dominant in the low-volume phase. Upon expansion of the lattice, the Lifshitz transition sets in, resulting in the (π, 0)-type magnetic correlations in the high-volume phase.
Our findings suggest that the proximity of a van Hove singularity to the Fermi level strongly influences, or even induces, (unconventional) superconductivity in the chalcogenide FeSe 1−x Te x series [145].

First-principles calculation of atomic forces and structural distortions
So far the phase stability of strongly correlated materials was investigated by performing total-energy calculations within DFT+DMFT. These calculations are very demanding even for simple materials, since they require the minimization of the total energy as a function of all atomic displacements. The computational effort therefore increases exponentially, which strongly limits the applicability of total-energy based techniques. This obstacle can be overcome by computing the complete set of interatomic forces using the Hellmann-Feynman theorem, whereby it becomes possible to compute the lattice structure even of complex materials. We conclude this section by formulating the DFT+DMFT approach for the calculation of interatomic forces and structural distortions in correlated materials based on the implementation of DFT+DMFT within the linear-response formalism [149]; see also a recent formulation based on a stationary implementation of the DFT+DMFT functional [229,230]. The calculation of forces makes it possible to compute atomic displacements and equilibrium atomic positions and, hence, to explain the origin of lattice transformations caused by electronic correlations. Moreover, it allows one to determine the equilibrium lattice structure of correlated systems even in the vicinity of a Mott metal-insulator transition -a computation which was not feasible up to now.
We discuss the DFT+DMFT scheme for the calculation of interatomic forces and structural distortions in correlated materials following Ref. [149]. For illustrative purposes we restrict our presentation to a discussion of structural transitions in a correlated model system, elemental (solid) hydrogen. By choosing different values of the local Coulomb interaction parameter U , we are able to explore the properties of solid hydrogen near a Mott metal-insulator phase transition.
The interatomic force acting on the atom s is calculated as the first-order derivative of the total energy functional (15): Here δ s ≡ d/dR s denotes the derivative with respect to the atomic position R s , and F s DFT is the force on the atom s calculated within DFT. Furthermore, δ s Ĥ DFT is the thermal average of the force operator δ sĤDFT , which leads to the Hellmann-Feynman contribution given by the first-order changes of the DFT Wannier HamiltonianĤ DFT , plus the term arising from the explicit dependence of the local Green function on the atomic positions: The derivative of the local Green function is given by To compute the interatomic forces caused by the Coulomb interaction (the 4th term on the right-hand side of Eq. (18)), we make use of the first-order derivative of the Galitskii-Migdal formula [231,232] for the interaction energy F s Here we assume that the average Coulomb interactionŪ and Hund's rule coupling J retain their values when the atomic positions change. The force operator δ sĤLDA and the first-order change of the self-energy δ sΣ (ω) are the two independent variables in the force functional [Eq. (18)] which need to be evaluated to compute the interatomic forces [233].
To determine δ sĤDFT , we generalize the projection scheme used to evaluate the DFT Wannier Hamiltonian [161,162]. The former is based on the projection of the set of site-centered atomic-like trial-orbitals |φ n on the Bloch functions |ψ ik of the selected bands with indices ranging from N a to N b . In this way the force operator may be written as where δ s V KS ik and δ s V Hxc ik denote the first-order changes in the LDA Kohn-Sham and the Hartree and exchange-correlation potentials, respectively [234]. Within the planewave pseudopotential approach [7,160] the Kohn-Sham contribution δ s V KS ik can be calculated as where V KS s (G, G ) is the Kohn-Sham potential for atom s (for details see Refs. [235,236]). The contribution δ s V Hxc ik is obtained from linear-response DFT calculations [235,236].
To evaluate the change of the self-energy δ sΣ (ω) we perform a functional derivative of the impurity Green function (we suppress the spin/orbital indices and assume a summation over repeated indices) and make use of the first-order derivative of the local Green function [Eq. (20)]. Eqs. (20) and (23) are solved self-consistently by employing δ sĜ −1 = δ sĜ −1 +δ sΣ and the two-particle correlation function, i.e., the generalized susceptibility, χ(τ 1 , τ 2 , τ 3 , τ 4 ) calculated within DMFT. As a starting point, we evaluate an initial guess for δ sΣ employing the static Hartree approximation This gives us the initial guess for δ sĜ −1 which we then insert into Eq. (23) to evaluate δ sĜ . From δ sĜ and δ sĜ we obtain a new estimate for δ sΣ which allows us to compute the new δ sĜ using the k-integrated Dyson equation [Eq. (20)]. This scheme is iterated until self-consistency over the δ sΣ is reached.

Applications to elemental solid hydrogen
To test the formalism proposed here, we perform a series of calculations for the simplest correlated electron problem, namely elemental (solid) hydrogen assuming a cubic structure with lattice constant a = 8 atomic units (a.u.). We then compare our results for the total energy computed as a function of atomic displacement with those obtained by the numerical integration of the corresponding forces [149]. Furthermore, by varying the local Coulomb interaction U , we explore the structural properties of solid hydrogen near a Mott metal-insulator phase transition. Thereby we can test whether our approach is able to determine structural transformations in the vicinity of the Mott-Hubbard transition, which is still a challenging problem of present-day solid state physics. The nonmagnetic LDA calculations for cubic hydrogen yield a metallic solution with a half-filled hydrogen s band of 3 eV width located at the Fermi level. To evaluate the force, we consider a supercell with two hydrogen atoms, in which one of the atoms is displaced by a distance δ with respect to its crystallographic position. In Fig. 11 we show the results for the total energy obtained by LDA as a function of δ. The nonmagnetic LDA calculations find the cubic lattice of hydrogen to be unstable since the total energy decreases with δ. Now we take into account the electronic correlations by calculating the properties of paramagnetic hydrogen using the DFT+DMFT method. For the partially filled hydrogen s orbitals a basis of atomic-centered symmetry constrained Wannier functions is constructed. The calculations are performed for the U values in the range of 1 -4 eV at a temperature T = 0.1 eV.
In Fig. 11 we also present our results for the total energy calculated by DFT+DMFT for paramagnetic hydrogen as a function of the displacement δ. By changing the U values, we were able to check the accuracy of our method by calculating the kinetic and interaction contributions, respectively, to the total force. By integrating the corresponding force with respect to δ, we find excellent agreement (within 1 -2 meV) between the force-based and the total energy calculations -even for large displacements δ (up to ∼ 10 % of the lattice constant a). Most interestingly, by increasing U , the cubic lattice (more precisely, the investigated displacive mode) becomes (meta-) stable for U ≥ 4 eV. These results clearly demonstrate the crucial importance of electronic correlations for the lattice stability of correlated materials.
Further applications of the DFT+DMFT scheme implemented within the linearresponse formalism to SrVO 3 and KCuF 3 are discussed in Ref. [149]. Our results for solid hydrogen, the correlated metal SrVO 3 , and the correlated Mott-Hubbard insu-lator KCuF 3 , demonstrate that the DFT+DMFT linear-response method presented here provides a robust computational tool for the study atomic displacements caused by electronic correlations. In particular, it allows one to determine the structural phase stability of both metallic and insulating correlated materials. The approach opens the way to calculate forces and thereby explore lattice instabilities induced by electronic correlations. Lattice dynamical properties of correlated electron materials can also be calculated by implementing the approach with, for example, the so-called small displacements method (see, e.g., Ref. [237]).

Conclusions and Outlook
In this article we reviewed the main results of our research into multi-band correlation phenomena during the funding period of the DFG Research Unit FOR 1346. In the first part of the paper we discussed the DMFT treatment of spontaneous symmetry breaking. We outlined the two general approaches to this problem: first, monitoring the divergent susceptibilities in the normal phase that indicate an instability towards long range ordering and, second, calculations in the phases with broken symmetry. As a pilot problem we studied materials close to the spin-state crossover. The investigation of the two-band Hubbard model, which provides a minimal description of this phenomenon, revealed a rich phase diagram with numerous phases with diverse properties. The linear-response approach to the calculation of susceptibilities turned out to be very useful for the identification of continuous phase transitions without any prior knowledge of the type of symmetry breaking. Nevertheless calculations in the ordered phases are indispensable in order to identify first-order transitions or phase separation as well as to investigate the properties of the ordered states.
Besides serving as a non-trivial test ground for the implemented formalism, these model calculations had a very specific materials based motivation, namely the longstanding problem of spin-state crossover in compounds from the LaCoO 3 family. We performed a series of DFT+DMFT studies of various aspects of the LaCoO 3 physics, the results of which led us to propose a new picture of this material based on the concept of spinful excitons and their condensation. Experimental efforts to test the proposal are under way.
The competition of spin states in the two-band Hubbard model and the resulting long-range order, including excitonic magnetism and superconductivity, have been studied so far only on simple lattices, such as the Bethe and hypercubic lattice. Other lattice geometries and hopping patterns will undoubtedly lead to new types of order, e.g., the formation of a supersolid or magnetically ordered supersolid, which is of great theoretical interest. Experimental realizations are ultimately necessary to test these results. A first step is to establish experimental techniques capable of unambiguously identifying these exotic states of matter. Scattering experiments, e.g., employing resonant inelastic x-ray scattering (RIXS), and their dependence on the polarization may have this capability. Theoretical simulations of RIXS spectra in the different states are therefore necessary.
In the second part of the paper we reviewed the DFT+DMFT scheme for the computation of the electronic structure and phase stability of correlated materials. This method allows one to explore structural transformations of the ionic lattice caused by correlations among the electrons, and to investigate lattice dynamical properties of correlated materials. We employed the DFT+DMFT scheme to explain the phase stability and lattice dynamical properties of elemental paramagnetic iron near the bcc-fcc phase transition. Furthermore, we discussed the origin of the phase stability of the high-temperature δ-phase. Our calculations clearly demonstrated the crucial importance of electronic correlations for the thermodynamic and lattice dynamical stability of the paramagnetic bcc phase of iron.
We also investigated the parent compound of the Fe-based superconductors, FeSe. In particular, we computed the electronic structure and phase stability of FeSe as a function of lattice volume and predicted a topological change (Lifshitz transition) of the Fermi surface upon expansion of the lattice as can be achieved experimentally by substituting Se by Te. This reconstruction is accompanied by a sharp increase of the local moments.
The present implementation of DFT+DMFT employs an adiabatic (Born-Oppenheimer) approximation, where it is assumed that the electrons move in an effective potential produced by the heavy nuclei, i.e., the electrons and nuclei are treated independently. This simplification excludes a dynamic coupling between the electronic and lattice degrees of freedom, which is of importance for the investigation of, for example, dynamical polaron effects. In principle, appropriate generalizations of the DFT+DMFT scheme can be achieved by combining it with linear-response techniques to compute the electron-phonon coupling and by the implementation of non-equilibrium DFT+DMFT methods employing the Keldysh formalism for the electrons with (quantum) molecular dynamics simulation techniques [238]. We leave these projects for the future.
The field of DMFT has developed into many directions in the period covered here. Among them is the growing interest in two-particle quantities which, after all, describe particularly important physical properties and also the technologically most exploited features of materials. New impurity solvers were implemented [239][240][241][242][243], and new algorithms developed [244][245][246][247] that greatly improved the computation of two-particle correlation functions. Along with the development of diagrammatic approaches beyond DMFT [248][249][250], fundamental questions about two-particle consistency [251], the consequences of its violation in DMFT, and their possible remedies were raised and are expected to give fresh impetus. The availability of two-particle correlators for multi-orbital quantum impurities, which allows for the application of linear response to realistic systems, e.g., with full d-shells, calls for further developments in the implementation of this data-intensive formalism.