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

TL;DR

This paper develops and applies a combined DFT+DMFT approach to study electronic and structural instabilities in correlated materials, revealing new phases and stability conditions driven by electron correlations.

Contribution

It introduces a formalism for total energy and force calculations within charge self-consistent DFT+DMFT, enabling analysis of correlation-driven phase transitions and instabilities.

Findings

Discovery of spinful exciton condensation in a two-orbital Hubbard model.

Identification of correlation-induced structural instabilities in Fe and FeSe.

Development of a computational scheme for total energies in DFT+DMFT.

Abstract

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.

Figures

Equations

$$\chi_{ij,kl}(\bm{q})=\int_{0}^{\beta}d\tau\sum_{\bm{R}}e^{i\bm{q}\bm{R}}\bigl{(}\langle Tc^{\dagger}_{\bm{R}j}(\tau)c^{\phantom{\dagger}}_{\bm{R}i}(\tau)c^{\dagger}_{\bm{0}k}(0)c^{\phantom{\dagger}}_{\bm{0}l}(0)\rangle-\langle c^{\dagger}_{j}c^{\phantom{\dagger}}_{i}\rangle\langle c^{\dagger}_{k}c^{\phantom{\dagger}}_{l}\rangle\big{)},$$

$${\chi }_{ij,kl}(\mathbf{q})=T\sum _{m,n}{\tilde{\chi }}_{ij,kl}(\mathbf{q};{\omega }_m,{\omega }_n),$$

$$\begin{array}{rl}{\tilde{\chi }}_{ij,kl}(\mathbf{q};{\omega }_1,{\omega }_2)& ={\tilde{\chi }}_{ij,kl}^0(\mathbf{q};{\omega }_1,{\omega }_2)\\ & +T\sum _{{\omega }_3,{\omega }_4}{\tilde{\chi }}_{ij,mn}^0(\mathbf{q};{\omega }_1,{\omega }_3){\Gamma }_{mn,pq}({\omega }_3,{\omega }_4){\tilde{\chi }}_{pq,kl}(\mathbf{q};{\omega }_4,{\omega }_2)\end{array}$$

$$\begin{array}{rl}{\tilde{\chi }}_{ij,kl}({\omega }_1,{\omega }_2)& ={\tilde{\chi }}_{ij,kl}^0({\omega }_1,{\omega }_2)\\ & +T\sum _{{\omega }_3,{\omega }_4}{\tilde{\chi }}_{ij,mn}^0({\omega }_1,{\omega }_3){\Gamma }_{mn,pq}({\omega }_3,{\omega }_4){\tilde{\chi }}_{pq,kl}({\omega }_4,{\omega }_2).\end{array}$$

$$H$$

$$+U\sum _i\left(n_{i\uparrow }^a{n}_{i\downarrow }^a+n_{i\uparrow }^b{n}_{i\downarrow }^b\right)+U^{\prime }\sum _{i,\sigma {\sigma }^{\prime }}{n}_{i\sigma }^a{n}_{i{\sigma }^{\prime }}^b-J\sum _{i\sigma }\left(n_{i\sigma }^a{n}_{i\sigma }^b+a_{i\sigma }^{†}{a}_{i\bar{\sigma }}^{}{b}_{i\bar{\sigma }}^{†}{b}_{i\sigma }^{}\right)$$

$$+J^{\prime }\sum _i\left(a_{i\uparrow }^{†}{a}_{i\downarrow }^{†}{b}_{i\downarrow }^{}{b}_{i\uparrow }^{}+\text{h.c.}\right),$$

$$s_i^{\mu }$$

$$o_i$$

$$d_i^{\mu }$$

$$c_i^{\mu }$$

$$\hat{H}$$

$$[H_{\text{DFT}}^{\sigma }(\mathbf{k}){]}_{\mu \nu }=\sum _i{P}_{i\mu }^{\sigma \ast }(\mathbf{k}){ϵ}_{\mathbf{k}i}^{\sigma }{P}_{i\nu }^{\sigma }(\mathbf{k}),$$

$$\rho (\mathbf{r})=T\sum _{\mathbf{k},i{\omega }_n;ij}{\rho }_{\mathbf{k};ij}{G}_{\mathbf{k};ji}(i{\omega }_n)e^{i{\omega }_{n}0+},$$

$${\rho }_{\mathbf{k};ij}(\mathbf{r})=⟨{\psi }_{\mathbf{k}i}|\mathbf{r}⟩⟨\mathbf{r}|{\psi }_{\mathbf{k}j}⟩+\sum _{I,lm}{Q}_{lm}(\mathbf{r}-{\mathbf{R}}_I)⟨{\psi }_{\mathbf{k}i}|{\beta }_l^I⟩⟨{\beta }_m^I|{\psi }_{\mathbf{k}j}⟩.$$

$$G^{\sigma }(\mathbf{k},i{\omega }_n{)}_{ij}^{-1}=(i{\omega }_n+\mu -{ϵ}_{\mathbf{k}i}^{\sigma }){\delta }_{ij}-{\Sigma }_{ij}^{\sigma }(\mathbf{k},i{\omega }_n),$$

$$E=E_{\text{DFT}}[\rho (\mathbf{r})]+⟨{\hat{H}}_{\text{DFT}}⟩-\sum _{m,\mathbf{k}}{ϵ}_{m\mathbf{k}}^{\text{DFT}}+⟨{\hat{H}}_U⟩-E_{\text{DC}},$$

$$\sum _{m,\mathbf{k}}{ϵ}_{m,\mathbf{k}}^{\text{DFT}}=T\sum _{i{\omega }_n,\mathbf{k}}\text{Tr}[H_{\text{DFT}}(\mathbf{k})G_{\mathbf{k}}^{\text{DFT}}(i{\omega }_n)]e^{i{\omega }_n{0}^{+}}.$$

$$⟨{\hat{H}}_{\text{DFT}}⟩$$

$$F_s$$

$${\delta }_s⟨{\hat{H}}_{\mathrm{DFT}}⟩$$

$${\delta }_s{\hat{G}}_{\mathbf{k}}(\omega )={\hat{G}}_{\mathbf{k}}(\omega )[{\delta }_s{\hat{H}}_{\mathrm{DFT}}^{\mathbf{k}}+{\delta }_s\hat{\Sigma }(\omega )-{\delta }_s\mu ]{\hat{G}}_{\mathbf{k}}(\omega ).$$

$$({\delta }_s{\hat{H}}_{\mathrm{DFT}}^{\mathbf{k}}{)}_{nm}$$

$${\delta }_s{V}_{i\mathbf{k}}^{\mathrm{KS}}$$

$${\delta }_s\hat{G}({\tau }_1-{\tau }_2)=-\hat{\chi }({\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4){\delta }_s{\hat{\mathcal{G}}}^{-1}({\tau }_3,{\tau }_4)$$

$$\hat{\chi }({\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4)=⟨{\mathcal{T}}_{\tau }\hat{c}({\tau }_1){\hat{c}}^{†}({\tau }_2){\hat{c}}^{†}({\tau }_3)\hat{c}({\tau }_4)⟩-⟨{\mathcal{T}}_{\tau }\hat{c}({\tau }_1){\hat{c}}^{†}({\tau }_2)⟩⟨{\mathcal{T}}_{\tau }{\hat{c}}^{†}({\tau }_3)\hat{c}({\tau }_4)⟩,$$

$${\delta }_s\hat{\Sigma }\simeq U{\delta }_s\hat{N}$$

Full text

1 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 Rev1 ; Rev2 ; Rev3 . These compounds are therefore particularly interesting in view of possible 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) LDA , the generalized gradient approximation (GGA) PB96 , or using the so-called LDA+U method Anisimov91 . 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 Baroni01 ; RevWIEN2k ; RevVASP . 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) metzner89 ; georges96 ; Kotliar04 ; P2:springerbook , usually referred to as DFT+DMFT (or more explicitly as LDA+DMFT, GGA+DMFT, etc.) Anisimov97 ; Lichtenstein98 ; P2:LDADMFT2002b ; Held06 ; Kotliar06 ; P2:Held07 ; P2:Katsnelson_RMP08 ; P2:FOR1346-Proceedings-2011 , 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 Byczuk12 ; Byczuk12_ .

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 one-particle 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 plane-wave 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 Leonov08 ; Leonov08_ ; Leonov08__ ; Leonov11 ; Leonov12 ; Leonov15a . 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 LaCoO3 in Sec. 2.3. In Sec. 2.4 we introduce a minimal two-orbital model of spin-state 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.

2 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.

2.1 Method 1

Numerous forms of magnetic order of strongly correlated electron systems have been studied intensively by DMFT-based techniques in basic models chitra99 ; zitzler02 ; sangiovanni06 ; hoshino10 and materials ranging from oxides of $3d$ kunes09 ; kunes12 , $4d$ jakobi11 ; kim15 ; dang15 or $5d$ arita12 elements, pnictides yin11 ; aichhorn11 all the way to elemental metals katsnelson99 ; lichtenstein01 . Apart from superconductivity in the attractive keller01 ; capone02b ; capone02 ; toschi05 ; bauer09 ; bauer09b ; koga10 and multiband koga15 ; vanhala15 ; hoshino15 ; hoshino16 Hubbard models DMFT investigations of long-range ordered states other than magnetically ordered phases are rare. We note that superconductivity in repulsive Hubbard model was extensively studied maier00 ; maier05 ; haule07b with cluster extensions of DMFT dca .

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 kunes11a ; kunes14b ; kunes14c ; kunes16 ; sotnikov16 . 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 kaneko14 ; kaneko15 and Vanhala et al. using cellular DMFT vanhala15 . Other examples of unconventional ordered states obtained within DMFT are the breaking of the channel symmetry in the two-channel Kondo lattice reported by Hoshino et al. hoshino11 and the spin/orbital freezing scenario for cuprate superconductors discussed by Werner et al. werner17a .

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 werner06 ; alps1 to allow for a real off-diagonal hybridization.

As a minor methodological development we briefly mention a trick adapting the maximum entropy formalism maxent for indefinite (off-diagonal) spectral densities. It consists in writing the off-diagonal spectral function in the form $A(\omega)=A_{+}(\omega)-A_{-}(\omega)$, where $A_{\pm}(\omega)\geq 0$, which requires only minor code modifications. While this decomposition is highly non-unique we found the stability of $A(\omega)$ to be satisfactory. An alternative approach working with an indefinite $A(\omega)$ was proposed in Ref. otsuki17 .

2.2 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. kunes14c 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.

2.2.1 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 mar08 . 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 $\bm{q}$-dependent susceptibilities.

When searching for (particle-hole) instabilities within the DMFT formalism one is interested in static susceptibilities of the type

[TABLE]

where $\bm{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 georges96 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 $\bm{q}\neq 0$). The susceptibility (1) is obtained from the two-particle correlation function $\tilde{\chi}_{ij,kl}(\bm{q};\omega_{m},\omega_{n})$

[TABLE]

which is the solution of a pair of Bethe-Salpeter equations jarrell92 ; georges96

[TABLE]

Here $\tilde{\chi}_{pq,kl}(\omega,\omega^{\prime})$ is the impurity two-particle correlation function, and $\tilde{\chi}^{0}_{ij,kl}(\omega,\omega^{\prime})$ and $\tilde{\chi}^{0}_{ij,kl}(\bm{q};\omega,\omega^{\prime})$ 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. kunes11b .

DMFT calculations within linear response were performed for one-orbital models in the early days of DMFT jarrell92 ; jarrell95 ; ulmke95 . It is only quite recently that linear response has been used for an unbiased search for instabilities in two- and three-orbital Hubbard models kunes14a ; boehnke14 ; hoshino15 ; hoshino16 ; steiner16 ; werner17a . 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. boehnke12 and, with a simplified RPA-like vertex, even for real materials park11 .

2.3 Spin-state crossover: The case of LaCoO3

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) (“spin-state crossover”) is essentially an atomic effect — a consequence of the competition between the Hund’s rule coupling $J$ and the crystal field (CF) splitting $\Delta$. Varying these parameters may lead to a level crossing tanabe 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, Fe2O3, FeO, CoO, which were theoretically studied with the LDA+DMFT approach kunes08 ; kunes09 ; shorikov10 ; ohta12 ; leonov15 ; leonov16 . 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.

LaCoO3 is an interesting exception. In this material the parameters are fine-tuned 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 LaCoO3 and related compounds have attracted much attention and have been studied already for half a century heikes64 ; naiman65 ; raccah67 . Below 50 K LaCoO3 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 LaCoO3: (i) the single-ion picture of a LS, $S$ $=$ [math], ground state of the Co3+ 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 yamaguchi96 ; zobel02 ; kyomen03 ; kyomen05 . Generally, it is not possible to describe LaCoO3 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 abbate93 ; korotin96 ; knizek05 . A Hartree-Fock-like LDA+U method has been employed to study stable spin states and their possible ordered patterns hsu09 ; knizek09 . 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 eder10 ; krapek12 ; zhang12 ; karolak15 .

Křápek et al. krapek12 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 tanabe ; haverkort06 , 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.

Recently, Sotnikov and Kuneš sotnikov16 proposed a mechanism in which the IS states still play an important role in LaCoO3. 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 sotnikov16 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 LaCoO3 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 $\bm{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 $\bm{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. krapek12 . 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 low-temperature (LS) state of LaCoO3 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 LaCoO3 at low temperature appears to be an uncorrelated band insulator with sharp quasi-particle bands krapek12 . The “hidden” correlations in the material become manifest in the one-particle spectra only upon heating or doping.

2.4 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 afonso16 , 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

[TABLE]

with the notation $\bar{\sigma}=-\sigma$. It describes electrons with two orbital ($a$ and $b$) and two spin ($\sigma=\uparrow,\downarrow$) flavors moving on a lattice and interacting via an on-site interaction $H_{\text{int}}$. Here, $a_{i\sigma}^{\dagger}$ ($a^{\phantom{\dagger}}_{i\sigma}$) is an operator creating (annihilating) a fermion with orbital flavor $a$ and spin $\sigma$ on the lattice site $i$, $n_{i\sigma}^{a}=a_{i\sigma}^{\dagger}a^{\phantom{\dagger}}_{i\sigma}$ is the corresponding local density operator, and analogously for the $b$ fermions. The sum $\sum_{\langle ij\rangle}$ 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 $\Delta$ describes the CF splitting. The standard Slater-Kanamori interactions are parametrized by $U$ and $J$ (with $U^{\prime}=U-2J$, $J^{\prime}=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^{\prime}$ 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 $\Delta$ werner07 ; suzuki09 . Large $\Delta$ favors the singlet LS state, while large $J$ favors the triplet HS state. When spontaneous symmetry breaking is excluded the $\Delta-U$ phase diagram at fixed $J/U$ can be divided into three regions: HS Mott insulator, LS band insulator, and a metallic state (“metal”) werner07 . The first-order 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 kunes11a 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 BEG , is its reentrant character in parts of the phase diagram hoston91 ; kunes15 . The above (classical) order is not the only instability of the model. Magnetic ordering of the HS states and condensation of spin-triplet excitons halperin68b ; zocher11 ; brydon09a ; brydon09b ; kaneko13 ; kaneko14 ; kunes14a ; hoshino16 are competing instabilities. In the following, we focus on the latter one.

2.5 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 Ca2RuO4 khaliullin13 ; jain15 or $3d^{6}$ cubic materials such as the Pr0.5Ca0.5CoO3 family kunes14b . Excitonic magnetism has been studied in simple lattice models such a 2BHM or effective pseudospin models khaliullin13 ; kunes14a ; kunes14b ; kaneko14 ; kaneko15 ; chaloupka16 . In the following we review the DMFT results for the spin-triplet exciton condensate.

2.5.1 Linear response

The linear-response formalism described in Sec. 2.2.1 (Method 2) allows for an unbiased search for instabilities. Recently, several authors kunes14a ; boehnke15 ; hoshino16 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 correspondence . These describe the response of the system to the fields

[TABLE]

with site indices $i$, spin indices ${\alpha,\beta=\uparrow,\downarrow}$ and Cartesian indices $\mu=x,y,z$ of the Pauli matrices $\sigma^{\mu}$.

Note that in the actual calculations with the density-density interaction, it is only the spin $s_{i}^{z}$ and the excitonic $c_{i}^{x},c_{i}^{y},d_{i}^{x},d_{i}^{y}$ 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}^{\mu}$ and spin-current-density $c_{i}^{\mu}$ modes are degenerate. In Fig. 3 both the orbital and the excitonic susceptibility diverge at $\bm{q}=(\pi,\pi)$ 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.

2.5.2 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 convenience, $t_{a}t_{b}<0$ to get a uniform condensate ($\bm{q}=0$) kunes15 . The condensate is characterized by a finite expectation value $\bm{\phi}_{i}=\sum_{\alpha\beta}\bm{\sigma}_{\alpha\beta}\langle a_{i\alpha}^{\dagger}b_{i\beta}^{\phantom{\dagger}}\rangle=\langle\bm{d}_{i}+i\bm{c}_{i}\rangle$. It is instructive to consider the condensate wave function in the strong-coupling limit, which can be approximated as a product $\prod_{i}\left(s|\text{LS}\rangle_{i}+\xi|\text{HS}\rangle_{i}\right)$ 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 $\langle\bm{S}_{i}\rangle\neq 0$, e.g., $|S_{z}=1\rangle$ give rise to a ferromagnetic condensate (FMEC) with finite, ferromagnetically ordered spin moments, while HS states with $\langle\bm{S}_{i}\rangle=0$, e.g., $|S_{z}=0\rangle$ or $\tfrac{1}{\sqrt{2}}(|S_{z}=1\rangle+|S_{z}=-1\rangle)$, 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 $\bm{\phi}_{i}$ kunes14c .

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 kunes15 . Below the transition temperature a sizeable gap ($\gg 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 halperin68b ; hoshino16 and involves a similar algebra for both phenomena. Another notable property of the PEC phase is a temperature independent spin susceptibility kunes14b .

2.5.3 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 eisenstein04 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 kunes16 . 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}|\ll|t_{a}t_{b}|$. We consider two possibilities for the cross-hopping: $V_{ab}=V_{ba}$ (even) and $V_{ab}=-V_{ba}$ (odd) kunes16 . 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 mean-field 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 halperin68b ), or by a local pattern of spin currents (referred to as spin-current-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 kunes16 . In particular, textures with odd $\bm{k}$-parity are interesting since they do not break time reversal symmetry and give rise to an effective Hamiltonian similar to Rashba-Dresselhaus spin-orbit coupling. Microscopically the spin textures can be traced to an effective spin dependent hopping of the doped carriers on the condensate background kunes16 .

3 Structural stability of correlated materials

In this section we present a detailed formulation of a variant of the DFT+DMFT computational scheme metzner89 ; georges96 ; Kotliar04 ; Kotliar06 ; Anisimov97 ; Lichtenstein98 ; Held06 , 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 Leonov08 ; Leonov08_ ; Leonov08__ ; Leonov11 ; Leonov12 ; Leonov15a .

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 Pourovskii07 ; Haule07 ; Amadon08 ; Aichhorn09 ; Park14 . 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 Leonov15a .

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 Savrasov01 ; Savrasov01_ ; McMahan03 ; Amadon06 ; Bieder14 ; lichtenstein01 ; Katanin10 ; Leonov11 ; Leonov12 ; Belozerov13 ; Leonov14a ; Glazyrin13 ; Glazyrin13_ ; Held01 ; Held01_ ; Held01__ ; Held01___ ; Poteryaev07 ; Grieger12 ; Grieger12_ ; leonov15 ; Leonov15a ; Leonov15c ; Pavarini04 ; Biermann05 ; kunes09 ; Poteryaev04 ; Leonov08 ; Leonov08_ ; Leonov08__ ; Leonov14b ; Pavarini08 ; Pavarini10 ; kunes08 ; Kunes08_ ; shorikov10 ; ohta12 ; Huang12 ; Huang12_ ; Ren06 ; Ren06_ ; Ren06__ ; Ren06___ ; Ren06____ ; Byczuk12 ; Byczuk12_ .

In particular, this method was used to study the electronic and structural properties of elemental Fe Leonov11 ; Leonov12 ; Leonov14a and the iron chalcogenide FeSe Leonov15c , and to explain the cooperative Jahn-Teller effect in paramagnetic KCuF3 and LaMnO3 Leonov08 ; Leonov08_ ; Leonov08__ . Recently it was employed to explore the electronic properties and the phase stability of paramagnetic V2O3 Leonov15a and the transition metal monoxides MnO, FeO, CoO, and NiO near the pressure-induced Mott-Hubbard metal-insulator transition leonov16 .

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.

3.1 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 Leonov15a a fully charge self-consistent DFT+DMFT approach Pourovskii07 ; Haule07 ; Amadon08 ; Aichhorn09 ; Park14 , which is formulated in terms of plane-wave pseudopotentials Baroni01 ; Giannozzi09 . Following Anisimov et al. Anisimov05 and Trimarchi et al. Trimarchi08 we construct a basis set of atomic-centered symmetry-constrained Wannier functions Marzari97 ; Anisimov05 ; Trimarchi08 ; Korotin08 , which allows us to build up an effective low-energy Hamiltonian $\hat{H}_{\text{DFT}}$ for the partially occupied correlated orbitals, e.g., $d$ of $f$ orbitals of transition-metal ions. The Hamiltonian $\hat{H}_{\text{DFT}}$, which provides a realistic description of the single-electron 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

[TABLE]

where $\hat{n}_{m\sigma}=\hat{c}^{\dagger}_{m\sigma}\hat{c}_{m\sigma}$ is the local density operator for the orbital $m$ and spin $\sigma$. Here $U^{\sigma,\sigma}_{m,m^{\prime}}$ and $U^{\sigma,\bar{\sigma}}_{m,m^{\prime}}$ 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 $\hat{H}_{\text{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 Korotin08 . The $U$ and $J$ values are related to the matrix notation as $U=\sum_{\{m\}}U^{\sigma\bar{\sigma}}_{mm^{\prime}}/N^{2}$ and $J=U-\sum_{\{m\neq m^{\prime}\}}U^{\sigma\sigma}_{mm^{\prime}}/N(N-1)$, 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 Baroni01 ; Giannozzi09 . In this scheme, the matrix elements of $\hat{H}_{\text{DFT}}$ are determined as

[TABLE]

where $P^{\sigma}_{i\nu}(\bm{k})=\langle\psi^{\sigma}_{\bm{k}i}|\hat{S}|\phi^{\sigma}_{\bm{k}\nu}\rangle$ are the matrix elements of the orthonormal projection operators expressed in the basis of local orbitals $\phi^{\sigma}_{\bm{k}\nu}$. The charge density $\rho(\bm{r})$ is calculated as

[TABLE]

where $T$ is the temperature. The charge density matrix elements $\rho_{\bm{k};ij}$ in the basis of the Kohn-Sham wave functions $\psi_{\bm{k}i}$ are computed as

[TABLE]

Here, the index $I$ refers to an atomic site, the functions $\beta^{I}_{m}(\bm{r})$ are computed in an atomic calculation and differ for different atomic species. Furthermore, $Q^{I}_{lm}(\bm{r})$ is the augmentation function that is strictly localized in the core region, and $G_{\bm{k};ij}(i\omega_{n})$ is the lattice Green’s function in the KS wave functions basis at a given $\bm{k}$-point. The (inverse) lattice Green’s function is evaluated as

[TABLE]

where $\epsilon^{\sigma}_{\bm{k}i}$ are the Kohn-Sham eigenvalues. $\Sigma^{\sigma}_{ij}(\bm{k},i\omega_{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 $\bm{k}$ and band $i$, i.e., $T\sum_{i\omega_{n}}G_{\bm{k};ij}(i\omega_{n})e^{i\omega_{n}0+}=f_{\bm{k}i}\delta_{ij}$. In practice, it is convenient to compute the charge density as $\rho(\bm{r})=\rho_{\text{DFT}}(\bm{r})+\Delta\rho(\bm{r})$, i.e., to split the DFT contribution and the correlation-induced difference $\Delta\rho(\bm{r})$. Full charge self-consistency is achieved when both the charge density $\rho(\bm{r})$ and the local Green’s function are converged.

Within the DFT+DMFT calculation the total energy is then evaluated using the expression Amadon06 ; Leonov08 ; Leonov08_ ; Leonov08__

[TABLE]

where $E_{\text{DFT}}[\rho(\bm{r})]$ is the DFT total energy obtained for the self-consistent charge density $\rho(\bm{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^{\text{DFT}}_{\bm{k}}(i\omega_{n})$:

[TABLE]

$\langle\hat{H}_{\text{DFT}}\rangle$ 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\omega_{n}$, taking into account an analytically evaluated asymptotic correction. Thus, for $\langle\hat{H}_{\text{DFT}}\rangle$ one has

[TABLE]

where the first moment $m^{\bm{k}}_{1}$ is computed as $m^{\bm{k}}_{1}=H_{\text{DFT}}(\bm{k})+\Sigma(i\infty)-\mu$. The asymptotic part of the self-energy $\Sigma(i\infty)$ is calculated as the average of $\Sigma(i\omega_{n})$ over the last several $i\omega_{n}$ points. The interaction energy $\langle\hat{H}_{U}\rangle$ is computed from the double occupancy matrix. The double-counting correction $E_{\text{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 V2O3 across the Mott-Hubbard metal-insulator transition (MIT) Leonov15a ; see also Refs. Grieger12 ; Grieger12_ . Namely, in charge non-self-consistent calculations the $a_{1g}$–$e_{g}^{\pi}$ crystal-field splitting of paramagnetic V2O3 is found to be strongly enhanced Held01 ; Held01_ ; Held01__ ; Held01___ ; Poteryaev07 . 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 V2O3 Leonov15a near the Mott metal-insulator transition. For a detailed discussion we refer to Ref. Leonov15a .

In the following we report our results obtained by the DFT+DMFT scheme for the electronic and structural properties of elemental Fe Leonov11 ; Leonov12 ; Leonov14a and the iron chalcogenide FeSe Leonov15c .

3.2 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 Basinski55 ; Wyckoff63 . At ambient conditions it is ferromagnetic and has a bcc crystal structure ($\alpha$ iron). Upon heating above the Curie temperature $T_{C}\sim$ 1043 K, $\alpha$ iron becomes paramagnetic, but remains in its bcc crystal structure. Only upon further increase of the temperature above $T_{\text{struct}}\sim$ 1185 K $\alpha$ iron exhibits a structural phase transition to a fcc structure ($\gamma$ phase). Under pressure $\alpha$ iron makes a transition to a paramagnetic hcp structure ($\epsilon$ phase) at $\sim$ 11 GPa.

3.2.1 Lattice stability and phonon spectra near the $\alpha$-$\gamma$ structural transition

State-of-the-art band structure methods provide a qualitatively correct description of various electronic and structural properties of iron Hasegawa83 ; Singh91 ; Singh91_ ; Singh91__ ; DalCorso00 ; Singh91____ ; Sha06 ; Kormann12 ; Ruban12 . For example, these methods provide a good quantitative understanding of the equilibrium crystal structure and the lattice dynamical properties of the ferromagnetic $\alpha$ phase. However, applications of these techniques to describe, e.g., the $\alpha$-$\gamma$ 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 Hsueh02 . 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 lichtenstein01 ; Katanin10 ; Leonov11 ; Leonov12 ; Belozerov13 ; Leonov14a ; Glazyrin13 ; Glazyrin13_ . The DFT+DMFT method naturally accounts for the existence of local moments above $T_{C}$ and has shown to provide a good quantitative description of the properties of $\alpha$ iron. Moreover, applications of DFT+DMFT to study the equilibrium crystal structure and phase stability of iron at the $\alpha$-$\gamma$ 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 Leonov11 .

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 Leonov12 ; Leonov14a . The approach is implemented with plane-wave pseudopotentials Baroni01 ; Giannozzi09 which allows us to compute lattice transformation effects caused by electronic correlations Leonov08 ; Leonov08_ ; Leonov08__ . 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. Leonov12 .

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 Anisimov05 ; Trimarchi08 ; Korotin08 , 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 Leonov12 . 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_{\text{struct}}\sim 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 Born-von 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 $\sim$ 458 K. The phonon frequencies are overall positive, implying mechanical stability of the bcc lattice structure at $\sim$ 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. Leonov12 . 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}$ Neuhaus97 .

We also study the lattice dynamical properties of the paramagnetic fcc phase which is found to become energetically favorable above $\sim 1.3~{}T_{C}$. To prove the mechanical stability of the fcc phase we evaluate the lattice dynamics of the fcc phase at $T\sim 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 $\sim$ 349 K. The phonon frequencies are overall positive, implying mechanical stability of the fcc lattice structure at $T\sim 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 $\sim 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\sim 1.4~{}T_{C}$ Zarestky87 . 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.

3.2.2 Origin of the phase stability of $\delta$ iron

At even higher temperatures, above 1670 K, iron is known to make a structural transition to the $\delta$ phase, which is stable up to the melting curve Basinski55 ; Wyckoff63 . To explain the phase stability of $\delta$ iron, we calculate the temperature evolution of the phonon dispersions near the $\alpha$-to-$\gamma$ and $\gamma$-to-$\delta$ phase transformations within DFT+DMFT Leonov14a . 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 and the experimental data taken at $T\sim T_{C}$ Neuhaus97 . In particular, we find a remarkable anomaly in the transverse (110) acoustic mode ($T_{1}$ mode) along the $\Gamma$-$N$ branch near the bcc-fcc phase transition temperature at $\sim$ 1.2 $T_{C}$. As will be discussed below this behavior is a dynamical precursor of the bcc-to-fcc phase transition Petry95 which occurs at $\sim$ 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\sim 1.4~{}T_{C}$, i.e., above the $\alpha$-to-$\gamma$ phase transition temperature Leonov14a . The origin of the instability lies in the $T_{1}$ mode in the $\Gamma$-$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\sim$ 1.8 $T_{C}$) the bcc structure of the $\delta$ phase is dynamically unstable in the whole $\Gamma$-$N$ branch, again due to the $T_{1}$ mode, with an additional anomaly near the $(\frac{2}{3},\frac{2}{3},\frac{2}{3})$ point. Therefore we expect that the $T_{1}$ transverse (110) phonon mode continuously softens upon temperature increase near the $\alpha$-$\gamma$ and $\gamma$-$\delta$ phase transition temperature – a feature which is typical for simple nonmagnetic metals. By contrast, at temperatures $T<T_{\text{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 ($\gamma$) 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. Drummond02 , 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 $(\delta)$ phase Leonov14a . On the basis of these results we therefore conclude that the high-temperature bcc $(\delta)$ 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.

3.3 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 kamihara08 ; Kamihara08_ ; Kamihara08__ ; Paglione10 ; Paglione10_ ; stewart11 ; Paglione10___ . It has been recently shown that the critical temperature $T_{c}$ of FeSe, which is about $T_{c}\sim$ 8 K at ambient pressure Hsu08 , 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 $\sim 37$ K under hydrostatic pressure of about 7 GPa and to $\sim 14$ K upon chemical (isovalent) substitution with Te Margandonna08 ; Margandonna08_ ; Mizuguchi08 ; Mizuguchi08_ ; Medvedev09 ; Sales09 ; Sales09_ . The angle-resolved photoemission and band structure calculations reveal that FeSe has the same Fermi surface topology as the pnictides Xia09 ; Tamai10 ; Tamai10_ ; Tamai10__ ; Nakayama10 ; Nakayama10_ ; Nakayama10__ ; Nakayama10___ ; Nakayama10____ ; Nakayama10_____ . It is characterized by an in-plane magnetic nesting wave vector $Q_{m}=(\pi,\pi)$, consistent with $s_{\pm}$ pairing symmetry Mazin08 ; chubukov08 . Moreover, both pnictides and chalcogenides display a strong enhancement of short-range spin fluctuations near $T_{c}$, with a resonance at $(\pi,\pi)$ in the spin excitation spectra Christianson08 ; Christianson08_ ; Christianson08__ ; Christianson08___ , suggesting a common origin of superconductivity in pnictides and chalcogenides, e.g., due to spin fluctuations Mazin08 ; chubukov08 . On the other hand, the related isoelectronic compound FeTe exhibits no superconductivity Medvedev09 ; Bao09 ; Bao09_ ; Bao09__ , but shows long-range antiferromagnetic $(\pi,0)$ order below $T_{N}\sim$ 70 K, suggesting that the solid solution Fe(Se,Te) should have a remarkable crossover from the $(\pi,\pi)$ to $(\pi,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 Li09 ; Li09_ .

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. Leonov15c . Most recently these calculations were generalized to take into account full charge self-consistency Skornyakov17 . 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 Leonov15c .

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 Leonov08 ; Leonov08_ ; Leonov08__ . For the partially filled Fe $3d$ and Se $4p$ orbitals we construct a basis set of atomic-centered symmetry-constrained Wannier functions Anisimov05 ; Trimarchi08 ; Korotin08 . To solve the realistic many-body problem, we employ the continuous-time hybridization-expansion quantum Monte-Carlo algorithm werner06 ; Gull11 . 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 $\bar{U}=$ 3.5 eV and Hund’s exchange $J=$ 0.85 eV, in accord with previous estimates for pnictides and chalcogenides Haule08 ; Haule08_ ; Aichhorn09 ; Haule08___ ; Haule08____ ; Haule08_____ ; yin11 ; aichhorn11 ; Haule08________ ; Haule08_________ ; Haule08__________ ; Haule08___________ ; Haule08____________ . 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 $P4/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.

3.3.1 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 Margandonna08 ; Margandonna08_ .

The calculated bulk modulus is $B\sim$ 70 GPa, which is comparable with that for iron pnictides. Most importantly, our DFT+DMFT calculations predict a structural transition of FeSe upon $\sim$ 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 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 $\sqrt{\langle m^{2}_{z}\rangle}$, which grows monotonically upon expansion of the lattice. The high-volume phase has a much larger local moment and a softer lattice with a much lower bulk modulus of 35 GPa.

3.3.2 Spectral function and Fermi surface structure

The results for the integrated spectral function are shown in Fig. 10 (left). In agreement with previous studies aichhorn11 , 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, $\sim$ 2.0 – 2.5, our calculations for the high-volume phase yield a strong renormalization $\sim$ 4 for the $xz/yz$ orbitals and $\sim$ 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”) Leonov15c . 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 $\Gamma$-point transform into incoherent spectral weight at the Fermi level along the $\Gamma$-$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}=(\pi,\pi)$, 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 $(\pi,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 FeSe1-xTex series Leonov15c .

3.4 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 Leonov14b ; see also a recent formulation based on a stationary implementation of the DFT+DMFT functional haule2016f ; haule2016g . 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. Leonov14b . 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):

[TABLE]

Here $\delta_{s}\equiv d/d\bf{R}_{s}$ denotes the derivative with respect to the atomic position $\bf{R}_{s}$, and $F^{s}_{\mathrm{DFT}}$ is the force on the atom $s$ calculated within DFT. Furthermore, $\delta_{s}\langle{\hat{H}_{\mathrm{DFT}}}\rangle$ is the thermal average of the force operator $\delta_{s}\hat{H}_{\text{DFT}}$, which leads to the Hellmann-Feynman contribution given by the first-order changes of the DFT Wannier Hamiltonian $\hat{H}_{\text{DFT}}$, plus the term arising from the explicit dependence of the local Green function on the atomic positions:

[TABLE]

The derivative of the local Green function is given by

[TABLE]

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 Migdal ; Galitskii58 for the interaction energy $F^{s}_{\mathrm{U}}=-\frac{1}{2}\text{Tr}\sum_{i\omega_{n}}[\delta_{s}\hat{\Sigma}(i\omega_{n})\hat{G}(i\omega_{n})+\hat{\Sigma}(i\omega_{n})\delta_{s}\hat{G}(i\omega_{n})]e^{i\omega_{n}0+}$. Here we assume that the average Coulomb interaction $\bar{U}$ and Hund’s rule coupling $J$ retain their values when the atomic positions change. The force operator $\delta_{s}\hat{H}_{\text{LDA}}$ and the first-order change of the self-energy $\delta_{s}\hat{\Sigma}(\omega)$ are the two independent variables in the force functional [Eq. (18)] which need to be evaluated to compute the interatomic forces chem_pot .

To determine $\delta_{s}\hat{H}_{\text{DFT}}$, we generalize the projection scheme used to evaluate the DFT Wannier Hamiltonian Anisimov05 ; Trimarchi08 . The former is based on the projection of the set of site-centered atomic-like trial-orbitals $|\phi_{n}\rangle$ on the Bloch functions $|\psi_{ik}\rangle$ of the selected bands with indices ranging from $N_{a}$ to $N_{b}$. In this way the force operator may be written as

[TABLE]

where $\delta_{s}V_{i\bf{k}}^{\mathrm{KS}}$ and $\delta_{s}V_{i\bf{k}}^{\mathrm{Hxc}}$ denote the first-order changes in the LDA Kohn-Sham and the Hartree and exchange-correlation potentials, respectively basis . Within the plane-wave pseudopotential approach Baroni01 ; Giannozzi09 the Kohn-Sham contribution $\delta_{s}V_{i\bf{k}}^{\mathrm{KS}}$ can be calculated as

[TABLE]

where $V_{s}^{\mathrm{KS}}(\bm{G},\bm{G^{\prime}})$ is the Kohn-Sham potential for atom $s$ (for details see Refs. DalCorso93 ; DalCorso93_ ). The contribution $\delta_{s}V_{i\bf{k}}^{\mathrm{Hxc}}$ is obtained from linear-response DFT calculations DalCorso93 ; DalCorso93_ .

To evaluate the change of the self-energy $\delta_{s}\hat{\Sigma}(\omega)$ we perform a functional derivative of the impurity Green function (we suppress the spin/orbital indices and assume a summation over repeated indices)

[TABLE]

with

[TABLE]

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 $\delta_{s}\hat{\mathcal{G}}^{-1}=\delta_{s}\hat{G}^{-1}+\delta_{s}\hat{\Sigma}$ and the two-particle correlation function, i.e., the generalized susceptibility, $\chi(\tau_{1},\tau_{2},\tau_{3},\tau_{4})$ calculated within DMFT. As a starting point, we evaluate an initial guess for $\delta_{s}\hat{\Sigma}$ employing the static Hartree approximation

[TABLE]

This gives us the initial guess for $\delta_{s}\hat{\mathcal{G}}^{-1}$ which we then insert into Eq. (23) to evaluate $\delta_{s}\hat{G}$. From $\delta_{s}\hat{G}$ and $\delta_{s}\hat{\mathcal{G}}$ we obtain a new estimate for $\delta_{s}\hat{\Sigma}$ which allows us to compute the new $\delta_{s}\hat{G}$ using the $\bm{k}$-integrated Dyson equation [Eq. (20)]. This scheme is iterated until self-consistency over the $\delta_{s}\hat{\Sigma}$ is reached.

3.4.1 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 Leonov14b . 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 $\delta$ with respect to its crystallographic position. In Fig. 11 we show the results for the total energy obtained by LDA as a function of $\delta$. The nonmagnetic LDA calculations find the cubic lattice of hydrogen to be unstable since the total energy decreases with $\delta$.

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 $\delta$. 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 $\delta$, we find excellent agreement (within 1 – 2 meV) between the force-based and the total energy calculations — even for large displacements $\delta$ (up to $\sim$ 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$ $\geq$ $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 linear-response formalism to SrVO3 and KCuF3 are discussed in Ref. Leonov14b . Our results for solid hydrogen, the correlated metal SrVO3, and the correlated Mott-Hubbard insulator KCuF3, 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. Alfe09 ).

4 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 long-standing problem of spin-state crossover in compounds from the LaCoO3 family. We performed a series of DFT+DMFT studies of various aspects of the LaCoO3 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 $\delta$-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 aoki14a . 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 alps2 ; triqs ; parragh12 ; lu14 ; wolf15 , and new algorithms developed hafermann12 ; boehnke11 ; gunacker15 ; gunacker16 that greatly improved the computation of two-particle correlation functions. Along with the development of diagrammatic approaches beyond DMFT toschi07 ; rubtsov08 ; rubtsov12 , fundamental questions about two-particle consistency vanloon16 , 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.

Acknowledgments

We gratefully acknowledge support by the Deutsche Forschungsgemeinschaft through FOR 1346.