# Unphysical and physical solutions in many-body theories: from weak to strong correlation

Adrian Stan, Pina Romaniello, Santiago Rigamonti, Lucia Reining, Arjan Berger

► To cite this version:

Adrian Stan, Pina Romaniello, Santiago Rigamonti, Lucia Reining, Arjan Berger. Unphysical and physical solutions in many-body theories: from weak to strong correlation. New Journal of Physics, 2015, 17 (9), pp.093045. 10.1088/1367-2630/17/9/093045. hal-01222216

HAL Id: hal-01222216

https://hal.science/hal-01222216

Submitted on 25 Jan 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Unphysical and physical solutions in many-body theories: from weak to strong correlation

This content has been downloaded from IOPscience. Please scroll down to see the full text.

2015 New J. Phys. 17 093045

(http://iopscience.iop.org/1367-2630/17/9/093045)

View the table of contents for this issue, or go to the journal homepage for more

Download details:

IP Address: 134.157.80.136

This content was downloaded on 25/01/2016 at 16:20

Please note that *terms and conditions* apply.
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

**RECEIVED**
14 July 2015

**REVISED**
21 August 2015

**ACCEPTED FOR PUBLICATION**
3 September 2015

**PUBLISHED**
25 September 2015

Content from this work
may be used under the
terms of the *Creative
Commons Attribution 3.0
licence*.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# Unphysical and physical solutions in many-body theories: from weak to strong correlation

Adrian Stan<sup>1,2,3</sup>, Pina Romaniello<sup>4</sup>, Santiago Rigamonti<sup>5</sup>, Lucia Reining<sup>3</sup> and J A Berger<sup>6,7</sup>

¹ Sorbonne Universités, UPMC Université Paris VI, UMR8112, LERMA, F-75005, Paris, France

² LERMA, Observatoire de Paris, PSL Research University, CNRS, UMR8112, F-75014, Paris, France

³ Laboratoire des Solides Irradiés, École Polytechnique, CNRS, CEA-DSM and European Theoretical Spectroscopy Facility (ETSF), F-91128 Palaiseau, France

⁴ Laboratoire de Physique Théorique, CNRS, IRSAMC, Université Toulouse III—Paul Sabatier and European Theoretical Spectroscopy Facility (ETSF), 118 Route de Narbonne, F-31062 Toulouse Cedex, France

⁵ Humboldt-Universität zu Berlin, Institut für Physik and IRIS Adlershof, and European Theoretical Spectroscopy Facility (ETSF), F-12489 Berlin, Germany

⁶ Laboratoire de Chimie et Physique Quantiques, IRSAMC, Université Toulouse III—Paul Sabatier, CNRS and European Theoretical Spectroscopy Facility (ETSF), 118 Route de Narbonne, F-31062 Toulouse Cedex, France

⁷ Author to whom any correspondence should be addressed.

E-mail: arjan.berger@irsamc.ups-tlse.fr

Keywords: multiple solutions, many-body perturbation theory, time-dependent density functional theory, dynamical mean-field theory

## Abstract

Many-body theory is largely based on self-consistent equations that are constructed in terms of the physical quantity of interest itself, for example the density. Therefore, the calculation of important properties such as total energies or photoemission spectra requires the solution of nonlinear equations that have unphysical and physical solutions. In this work we show in which circumstances one runs into an unphysical solution, and we indicate how one can overcome this problem. Moreover, we solve the puzzle of when and why the interacting Green's function does not unambiguously determine the underlying system, given in terms of its potential, or non-interacting Green's function. Our results are general since they originate from the fundamental structure of the equations. The absorption spectrum of lithium fluoride is shown as one illustration, and observations in the literature for some widely used models are explained by our approach. Our findings apply to both the weak and strong-correlation regimes. For the strong-correlation regime we show that one cannot use the expressions that are obtained from standard perturbation theory, and we suggest a different approach that is exact in the limit of strong interaction.

## 1. Introduction

In condensed-matter physics, important formalisms for predicting and understanding material properties, such as density-functional theory (DFT), many-body perturbation theory (MBPT), and dynamical mean-field theory (DMFT), avoid the use of the full many-body wave function and are instead based on simpler quantities, such as densities and Green's functions. The one-body Green's function G, for example, yields expectation values of one-body operators, and the total energy. In particular, G gives access to the spectral function, which is related to photoemission spectra, a widely used experimental technique to determine electronic structure. G is well defined as expectation value of particle addition and removal to the *N*-body ground state or thermal equilibrium. However, to use this definition one should know the many-body wave function, which is out of reach. A framework where the many-body wave function does not appear explicitly is provided by MBPT [1], where the interacting Green's function G is given as a functional of the non-interacting Green's function $G_0$ and the bare Coulomb interaction $v_c$. An important idea of MBPT is to avoid a possibly ill-behaved perturbation expansion of G in terms of $v_c$ and $G_0$ using Dyson equations. These are integral equations that describe the propagation of particles in terms of an effective potential or interaction. For the one-body Green's function, for
example, this effective potential, which is the kernel of the Dyson equation, is the self-energy $\Sigma$. The power of this approach resides in the fact that even a low-order approximation for $\Sigma$ yields contributions to all orders in $v_c$. Following Luttinger and Ward [2], $\Sigma$ is usually expressed as a functional of $G$ instead of $G_0$. However, this makes the Dyson equation nonlinear, which leads to multiple solutions [3, 4]. This is a very fundamental and general problem which, as we will demonstrate, also applies to other quantities such as response functions, from which one can calculate, for example, optical absorption spectra by either solving the Bethe-Salpeter equation (BSE) or the time-dependent DFT (TDDFT) equations. It is different from usual problems of convergence or local minima [5-8]. For example, convergence problems can be readily detected from the oscillatory behavior of the results; see appendix A for an example of such a problem in the case of Hartree-Fock (HF). The appearance of fully converged, but unphysical results, instead, is much more subtle and dangerous, and it has important consequences. It is the main topic of the present work.

We show that the presence of multiple solutions has an impact that reaches far beyond numerical problems, and even points to cases where the currently used strategies to derive approximations break down. We develop our ideas using a model that represents the structure of the exact theory. Calculations for a real system and links to recent observations [9, 10] demonstrate the potential of this approach for analysis and prediction. In particular, we answer several very general and important questions: (i) Is the problem of multiple solutions specific for certain cases, or is it a fundamental problem? Does it only appear for Green's functions or also in the framework of other well-established and widely used methods, like density-functional based approaches?; (ii) Does it impact calculations for real materials?; (iii) How can one detect and avoid unphysical solutions?; (iv) Does the problem depend on the interaction strength, and are there consequences for the many-body theory of strongly interacting materials?

The paper is organized as follows. In section 2 we introduce the one-point model (OPM) and we investigate multiple solutions for the map $G_0 \to G$. We then discuss the consequences of multiple solutions for the calculation of absorption spectra within TDDFT. Finally we investigate the effect of multiple solutions on the inverse map $G_0 \to G$. In section 3 we summarize our findings.

## 2. Theory and discussion

### 2.1. The OPM and the map $G_0 \to G$

To analyze the problem we use the so-called OPM [11-13]. This model is not system specific and can be solved exactly, such that the physical solution is well defined. It represents important structural aspects of the many-body problem, while collapsing all arguments of the Green's functions, self-energy, and the interaction to one point, making the equations scalar. In [3], an approximate version of the OPM was used to discuss multiple solutions within the framework of the GW approximation [14] to the self-energy.

In the present work we use the OPM without approximations, which simulates the full many-body problem. The exact OPM Green's function was derived in [15] from the one-point equivalent of the equation of motion of $G$, expressed as a functional differential equation [16]. The exact solution reads

$$y[y_0, u] = \frac{y_0}{1 + \frac{1}{2}uy_0^2} \quad \text{and} \quad \tilde{s}[y_0, u] = -\frac{1}{2}uy_0, \qquad (1)$$

where $y, y_0$, and $u$ represent $G, G_0$, and $v_c$, respectively. The self-energy $\tilde{s}$ is determined from the Dyson equation $\tilde{s}[y_0, u] = y_0^{-1} - y_0^{-1}[y_0, u]$. In equation (1) $\tilde{s}$ is given as a functional of the bare interaction $u$ and the non-interacting Green's function $y_0$. Usually, however, one works with the self-energy given as a functional of the dressed Green's function, $s[y, u]$. Then the Dyson equation reads

$$y = y_0 + y_0 s[y, u]y. \qquad (2)$$

This is, in general, a nonlinear equation. We first consider the HF self-energy, which in the OPM is<fnref n="8" /> $s^{\text{HF}}[y, u] = -\frac{1}{2}uy$. Let us look at the map $G_0 \to G$, i.e., the usual case, where $y_0$ is set by the system, and one searches $y$. The Dyson equation has two solutions

$$Y_{\text{HF}}^{\pm} = \frac{1}{V} \left[ -1 \pm \sqrt{1 + 2V} \right], \qquad (3)$$

with the rescaled quantities $Y = y/y_0$ and $V = uy_0^2$. Here $Y_{\text{HF}}^{\pm}$ is the physical solution, since it connects smoothly to $Y_0 = 1$ at $V = 0$, and $Y_{\text{HF}}$ is an unphysical solution, that diverges for vanishing interaction. Both are shown in the inset of figure 1. In real problems Dyson equations are solved iteratively. Two possible iteration schemes are:

<footnote n="8">⁸ Note that $s^{\text{HF}}[y, u] = -\frac{1}{2}uy$ is an approximation to the exact self-energy. The latter has the simple functional form $s^{\text{exact}}[y_0, u] = -1/2uy_0$ when expressed in terms of $y_0$, or equation (12) in terms of $y$, as explained below.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1. One-point model (OPM): $Z_0^\pm$ as a function of the interaction V. Squares (red): $Z_0^-$ and solution of scheme (A); circles (blue): $Z_0^+$ and solution of scheme (B); continuous line (orange): the exact solution $Y_0$. Inset: $Y_{HF}^\pm$ as a function of the interaction V.</figcaption></figure>

$$Y^{(n+1)} = \frac{2}{2 + VY^{(n)}} \mathbf{I}; \quad Y^{(n+1)} = \frac{2}{VY^{(n)}} - \frac{2}{V} \mathbf{II}. \qquad (4)$$

While neither of the two schemes has convergence problems when iterating, only scheme (I) converges to the physical solution, whereas scheme (II) converges to the unphysical solution. This happens because the iteration leads to the continued fraction representation of the square root [3] in equation (3)

$$\sqrt{1+x} = 1 + \frac{x/2}{1 + \frac{x/4}{1 + \frac{x/4}{1 + \dots}}} \qquad (5)$$

for $x = 2$ V. The sign of the square root is determined by the continued fraction in the iterative procedure<fnref n="9" />.

## 2.2. TDDFT and the map $\chi_0 \to \chi$

The simple but general structure of the OPM suggests that the same picture should emerge for any Dyson-like equation. For example, optical properties and screening can be calculated by solving the BSE for the two-particle correlation function, using the GW approximation for the self-energy [17]. The screened Coulomb interaction W is calculated from the BSE, and is also part of its kernel. Like in the above HF case, this makes the problem in principle self-consistent (see e.g. [18]). Alternatively, one can use TDDFT, that obeys a similar Dyson-like equation for the reducible polarizability $\chi$ [19]. Here we use TDDFT within the so-called bootstrap approximation proposed in [20]. The corresponding Dyson equation for frequency $\omega$ reads

$$\chi(\omega) = \chi_0(\omega) + \chi_0(\omega) \left[ \nu_c + \frac{1 + \nu_c \chi(\omega = 0)}{\chi_0(\omega = 0)} \right] \chi(\omega), \qquad (6)$$

where $\chi_0$ is the independent-particle polarizability. We evaluate this for a real material with a long-range Coulomb interaction and ab initio band structure. Again, this equation has two solutions that can be obtained by two iteration schemes, analogous to those in equation (4) (see appendix B). Figure 2 shows the absorption spectrum of LiF obtained with the two iteration schemes, as well as the experimental result. The experimental spectrum shows a strongly bound exciton with a binding energy of approximately 1.4 eV [21]. The spectrum obtained from iterating equation (6) with the analogue of scheme (I) is qualitatively correct, since it also shows a strongly bound exciton. The remaining discrepancies with respect to experiment are due to the approximate form of $f_{xc}$ [22]. Instead, iterating equation (6) within scheme (II) makes the exciton disappear completely. This means that, in absence of experimental results, one risks to run into a wrong solution, which would make the theory non-predictive. However, as we showed, the problem can be overcome using the appropriate iteration scheme (I). Note that optical properties can also be calculated from $\epsilon=1-\nu_c P$, where the irreducible polarizability P obeys a Dyson equation with kernel $f_{xc} = 1/[(1 - \nu_c P)\chi_0]$. In this case, the appropriate scheme is scheme (II). This is explained in appendix B.

<footnote n="9">⁹ We note that even using mixing, which is discussed in appendix A, scheme (II) does not yield the physical solution for all V.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 2. Optical absorption spectrum of LiF. Continuous line (red): physical solution; dashed line (blue): unphysical solution; dots (black): experiment [21]. The vertical line indicates the position of the quasiparticle gap.</figcaption></figure>

### 2.3. The map $G_0 \leftarrow G$

So far we have looked at the map $G_0 \to G$ (and $\chi_0 \to \chi$). We now focus on the inverse map $G_0 \leftarrow G$. This map is needed in problems of embedding, where one optimizes an auxiliary quantity $G_0$ in order to produce certain properties of a real system (contained in $G$). The inverse map is also crucial when one wants to express a functional in terms of dressed instead of bare quantities. The most prominent example is the Luttinger-Ward (LW) functional, where the self-energy is given in terms of $G$ instead of $G_0$ [2, 23–25]. For the LW functional to be properly defined, the map $G_0 \leftarrow G$ should be unique.

Within the OPM, consider a system with the bare Green’s function $y_0$, and with the exact, interacting Green’s function $y$ given by equation (1). We now fix $y$ and examine whether the inverse map $z_0 \leftarrow y$ unambiguously leads to $z_0 = y_0$. With the exact self-energy $\tilde{s}[z_0] = -\frac{1}{2}uz_0$ of equation (1), the exact Dyson equation of this problem reads

$$z_0 = y + \frac{1}{2}uyz_0^2, \qquad (7)$$

in which $y$ is known and $z_0$ is to be determined. This equation has again two solutions:

$$z_0^{\pm} = \frac{1}{uy} \left( 1 \pm \sqrt{1 - 2uy^2} \right) \Rightarrow Z_0^{\pm} = \frac{2 + V \pm \sqrt{(2 - V)^2}}{2V}, \qquad (8)$$

where $Z_0 = z_0/y_0$ and we used equation (1). The square root in equation (8) equals the absolute value $|2 - V|$. Because $2 - V$ changes sign at $V = 2$, the physical solution $Y_0 = 1$ is obtained by $Z_0^-$ for $V < 2$ and by $Z_0^+$ for $V > 2$ (see figure 1). In other words neither of the two solutions gives $Z_0 = Y_0$ for all $V$. As a consequence one has to change sign in front of the square root in equation (8) at $V = 2$. This has important consequences for the iterative solution of equation (7): because scheme (I) yields the square root with positive sign, to obtain the map $G_0 \leftarrow G$ we need two different iteration schemes: one for $0 < V < 2$ and the other for $V > 2$. This is different from the map $G_0 \to G$, where one solution gives the physical solution for all $V$, and hence a single iteration scheme suffices.

The need to change iteration scheme is a serious problem. Indeed, Kozik et al [9] pointed out that different iteration schemes, applied to Hubbard and Anderson models, lead to different solutions which cross at a certain interaction. Our OPM results provide the missing explanation: keeping the labels (A) and (B) of [9], the two iteration schemes correspond to

$$\frac{1}{Z_0^{(n+1)}} = 1 + \frac{1}{2}V(1 - Z_0^{(n)}), \qquad (A)$$

$$\frac{1}{Z_0^{(n+1)}} = -1 - \frac{1}{2}V(1 - Z_0^{(n)}) + \frac{2}{Z_0^{(n)}}. \qquad (B)$$

We report the results in figure 1. Scheme (A) converges to the physical solution for $V < 2$ and to the non-physical solution for $V > 2$. Instead, scheme (B) converges to the non-physical solution for $2/3 < V < 2$ and to the physical solution for $10^2 < V < 6$. These results are strictly analogous to those obtained by Kozik et al for

<footnote n="10">¹⁰ For $V < 2/3$ and $V > 6$ scheme (B) does not converge.</footnote>
Hubbard and Anderson models. They can be understood from the fact that scheme (A) creates a continued fraction with positive square root, whereas in scheme (B) the sign of the continued fraction is changed.

This sign problem is a priori a disaster because there is no unique prescription of how to avoid unphysical solutions. The OPM highlights the reducible polarizability [15]

$$ \frac{\chi}{\chi_0} = 2 \frac{2 - V}{(2 + V)^2}, \qquad (11) $$

as critical quantity that changes sign at the crossing¹ V = 2. At the same time, for V > 2 the perturbation expansion of y, in equation (1), diverges. Our result indicates that one should in principle inspect the exact two-particle correlation function as a function of the interaction to detect problems of perturbation theory. This is in line with [10] in which a breakdown of perturbation theory is linked to an eigenvalue of the two-particle correlation function that crosses zero, becoming negative<fnref n="12" />.

Inverting a map between functionals requires a careful definition of their domain [26, 27]. The multiple solutions are the price to pay for the fact that we have not considered this definition in the above discussion. For a system with a non-degenerate ground state this can be understood as follows: if there were two solutions for $G_0 \leftarrow G$, one could obtain the same dressed Green's function from two different $G_0$ and, hence, from two different external potentials<fnref n="13" />. Since the diagonal of G is the density, the Hohenberg-Kohn theorem [28] states that there can only be one external potential, and hence one $G_0$, corresponding to each G. This means that any additional solution $G_0$ is unphysical, in the sense that it cannot be constructed from the solution of a one-body Schrödinger equation. Equivalently, it cannot be written as a sum of simple poles, each with a strength normalized to one. By imposing this condition, one can therefore eliminate unphysical solutions. In the OPM this trivially corresponds to the requirement $Z_0 = 1$, which already implies the solution. In appendix C we work out a more realistic case where the definition of the domain defines the solution when the self-energy is static. A more general discussion on the definition of the domain can be found in [27]. It should be noted that when $G_0$ is an embedding Green's function the discussion is more complicated, because one searches for a fictitious $G_0$ with a frequency dependence that can differ from that of a $G_0$ resulting from a static potential.

With the map $G_0 \leftarrow G$ one can construct the exact self-energy as a functional of G. Using equation (8) in the exact self-energy given in equation (1) we obtain

$$ s^{\pm}[y, u] = -\frac{1}{2y} \left( 1 \pm \sqrt{1 - 2uy^2} \right), \qquad (12) $$

$$ = -\frac{1}{2y} \mp \frac{1}{2y} \pm \frac{1}{2} u \left[ y + \frac{uy^3}{2} + \frac{u^2y^5}{2} + \dots \right]. \qquad (13) $$

We note that the LW functional in equation (12) is unique, but in order to calculate it one has to change the sign. The Dyson equation with the two self-energies of equation (12) leads to two different Green's functions: the physical solution given in equation (1) is obtained from $s^-$ for $V < 2$ and from $s^+$ for $V > 2$, and $y = 0$ is obtained from $s^+$ for any V. Therefore, for weak interaction using the exact self-energy one obtains only one solution, the physical one, contrary to, e.g., the HF approximation. We note that at the point where $s^+$ and $s^-$ meet (at $V = 2$) the derivative $ds^+/dy$ diverges. This could explain the divergence of $\delta\Sigma/\delta G$ observed in [10] in the paramagnetic DMFT solution of a Hubbard model. Note that this divergence occurs at the point where one of the eigenvalues of the polarizability crosses zero.

In equation (13) we Taylor expanded the square root. The convergence radius is infinite, since $0 \le 2uy^2 \le 1$, as can be shown using equation (1). Interestingly, the sum of the first two terms in equation (13) (upper sign) is the first term of an expansion of the self-energy for strong interaction [3]. The remaining terms constitute an expansion in terms of a quantity that is proportional to u and converges for all physical y. This means that one can use perturbation theory over the whole interaction range, but in two different ways for the two different regimes<fnref n="14" />. To lowest order, this corresponds to HF, $s^{HF} = -\frac{1}{2}uy$, for weak interaction, and $s^{SIN-HF} = -\frac{1}{y} - s^{HF}$, for strong interaction. We call this functional strong-interaction HF (SIN-HF)<fnref n="15" />. Both self-energies yield two solutions. We report the physical solution for these two approximations in figure 3. While

<footnote n="11">¹¹ The polarizability is defined as the derivative of G with respect to an external potential. Therefore there is only one branch.</footnote>

<footnote n="12">¹² Negative spectral functions can also be caused by certain approximations [39]. This is a different problem.</footnote>

<footnote n="13">¹³ We note that $G_0$ must satisfy the Kubo-Martin-Schwinger boundary conditions, which means that there is only one solution $G_0$ to the equation $G_0^{-1}G_0 = 1$ [40, 41].</footnote>

<footnote n="14">¹⁴ Note that since $y \to 0$ for $u \to \infty$ the convergence of perturbation theory improves with increasing u.</footnote>

<footnote n="15">¹⁵ Equations (12) and (13) suggest that the exact self-energies in the weak (WIN) and strong interaction (SIN) regimes are related by $\Sigma^{SIN} = -\Sigma^{WIN} - 1/y$, indicating that one could use the knowledge of the weak interaction expansion to obtain good approximations for the strong-interaction case. Note, however, that the generalization of SIN-HF or higher order approximations of the strong-interaction series to real systems is not straightforward.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** One-point model (OPM): Y as a function of the interaction V for 4 < V < 100. Solid line (black): exact solution; dotted line (blue): Hartree–Fock (HF); dashed line (red): strong-interaction HF (SIN–HF). Inset: Y as a function of V for 0 < V < 5.</figcaption></figure>

HF clearly fails for strong interaction, SIN–HF is exact in the strong-interaction limit and performs well for V > 4, while it is worse than HF for V < 4. It is important to note that the physical SIN–HF solution is obtained for V > 1 with the iteration scheme $1/Y^{(n+1)} = 1/Y^{(n)} + \frac{1}{2}VY^{(n)} - 1$. Indeed, as it is also clear from the example of TDDFT, and discussed in appendix B, the appropriate iteration scheme depends on the formulation of the problem. We suggest the OPM as a powerful tool to examine which scheme one should use for a given problem and interaction range.

### 3. Conclusions

We have demonstrated that with a simple but general OPM one can understand and solve structural problems of MBPT. In particular, one can use it sort out the multiple solutions of the nonlinear Dyson equation by choosing the appropriate iteration scheme. We have shown that for the map $G_0 \rightarrow G$ a single iteration scheme suffices to obtain the physical solution for all interaction strenghts. Instead, for the inverse map $G_0 \leftarrow G$ one has to change iteration scheme at the interaction strength at which the reducible polarizability changes sign and perturbation theory of G in terms of $G_0$ starts to diverge. Nonetheless, we have proved that even for strong interaction one can use a perturbative expression for the self-energy in terms of G, which differs from the usual Luttinger-Ward functional. By presenting analogous results for real systems, and by comparing with numerical results in the literature, we have shown that these conclusions go far beyond the OPM. We expect that they will have an impact on many other questions in the domain of many-body physics.

### Acknowledgments

We thank Antoine Georges for stimulating discussions. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007–2013) / ERC Grant Agreement n 320971. SR acknowledges support by the Einstein Foundation Berlin.

### Appendix A. Convergence problems: the example of HF

As an example of convergence problems in many-body theory we consider the HF approximation. The task is to solve the Dyson equation $G^{-1} = G_0^{-1} - \Sigma^{\text{HF}}[G]$. Usually this is done iteratively as, for example, $G^{-1(n+1)} = G_0^{-1} - \Sigma^{\text{HF}}[G^{(n)}]$. This iteration scheme can converge to a HF solution that is not necessarily the ground-state, or it can oscillate between two states, neither of which is a HF solution [29]. There are several approaches to overcome convergence problems (see [30–35]). For example, one can mix the results of previous iteration steps, as $G^{-1(n+1)} = G_0^{-1} - \Sigma[(1 - \alpha)G^{(n)} + \alpha G^{(n-1)}]$, where $0 \le \alpha \le 1$ is the mixing parameter. In figure 4 we plot the HF total energy as a function of the interaction strength U for a chain of ten sites with a Coulomb-like interaction $\gamma = Ue^{-0.5|i-j|}/|i-j|$ (with i, j running over the sites). Without mixing ($\alpha = 0.0$) the HF total energy is linear in the interaction, but only for $U \le 5$. For $U > 5$ the scheme oscillates between
<figure><img src="image_7.png" /><figcaption>Figure 4. The Hartree–Fock (HF) total energy of a ten sites chain with a Coulomb-like interaction, as a function of the interaction U. Squares (red): results without mixing ($\alpha = 0.0$); circles (blue): results with mixing ($\alpha = 0.5$).</figcaption></figure>

two values. Instead, using $\alpha = 0.5$ leads to a total energy that is linear in $U$ over a large range. In other words, this problem can be easily detected and solved. It does not appear in calculations that we have performed on Hubbard clusters, irrespective of the interaction strength.

## Appendix B. Unphysical solutions in TDDFT

Optical absorption spectra can be obtained from the imaginary part of the macroscopic dielectric function $\epsilon_M$ that is related to the reducible polarizability $\chi$ according to

$$ \epsilon_M(\omega) = 1/\left[1 + \nu_c\chi(\omega)\right]. \qquad (B.1) $$

Here we neglected crystal local-field effects. The reducible polarizability $\chi(\omega)$ can be obtained from a Dyson-like equation:

$$ \chi(\omega) = \chi_0(\omega) + \chi_0(\omega)\left[\nu_c + f_{xc}\right]\chi(\omega). \qquad (B.2) $$

The bootstrap approximation to the exchange-correlation kernel $f_{xc}$ reads [20]:

$$ f_{xc} = \frac{1 + \nu_c \chi(\omega = 0)}{\chi_0(\omega = 0)}. \qquad (B.3) $$

This kernel is static by definition.

Since the bootstrap kernel depends on $\chi$, equations (B.2) and (B.3) were solved iteratively in [20]<fnref n="16" />. However, in [22] it was shown that these equations can be solved analytically and have two solutions, of which only one is physical. It turns out that the iteration scheme used in [20] gives the physical solution. Here we show that this iteration scheme is equivalent to scheme (I). Moreover, we show an iteration scheme, equivalent to scheme (II) that leads to well converged, but unphysical optical spectra.

The analytical solution of the system of equations (B.2) and (B.3) for $\omega = 0$ gives:

$$ \begin{aligned} \chi^{\pm} &= -\frac{\chi_0}{2} \pm \sqrt{\left(\frac{\chi_0}{2}\right)^2 - \frac{\chi_0}{\nu_c}}, \\ &= -\frac{\chi_0}{2} \pm \left|\frac{\chi_0}{2}\right| \sqrt{1 - \frac{4}{\nu_c \chi_0}}. \end{aligned} $$

The sign of the square root corresponding to the physical solution can be found by analyzing the limit of large screening $\nu_c \chi_0 \to -\infty$. For $\chi^-$, the correct result $\chi^- \to -1/\nu_c$, i.e., the large screening random-phase approximation (RPA) result, is obtained. Note that in this limit, equation (B.3) also gives the expected RPA result, that is $f_{xc} = 0$. For $\chi^+$, the unphysical limit $\chi^+ \to -\chi_0$ is obtained. In practice, one would iterate equations (B.2) and (B.3) up to self-consistency. Two iteration schemes emerge, which are given by setting $Y = -\nu_c \chi$ and $V = -2/\nu_c \chi_0$ in schemes (I) and (II). The physical solution is obtained using scheme (I), which is equivalent to the iterative solution proposed in [20].

<footnote n="16">¹⁶ The iterations are done only for $\omega = 0$.</footnote>
If, instead of equation (B.2), the Dyson equation for the irreducible polarizability $P = \chi_0 + \chi_0 f_{xc} P$ is used, and the optical spectrum is obtained from $\epsilon_M = 1 - \nu_c P$, one arrives again at schemes (I) and (II), with $Y = \nu_c P$ and $V = -2/\nu_c \chi_0$. In this case, however, the physical solution is given by scheme (II), which yields the correct large screening limit $P = \chi_0 \to -\infty$. This can be understood from the structure of the equations. It shows that the iteration scheme leading to the physical solution depends on the way in which the Dyson equation is formulated and that the OPM is a powerful tool to determine which iteration scheme has to be used.

Ab initio calculations of the optical spectrum of LiF were performed with the exciting code [36] using a grid of 4096 non-equivalent k-points. The LFE have been neglected. The Kohn–Sham eigenvalues from a ground-state calculation within the local-density approximation (Perdew–Wang) [37] were corrected through a scissors shift of 5.22 eV in order to mimic the experimental quasiparticle gap $E_g^{LiF} \approx 14.2$ eV [38].

## Appendix C. Definition of the domain

Here we give an illustration for the restriction of the number of solutions for $G_0 \leftarrow G$ to only one in the case of a static self-energy by restricting the domain of $G_0$ [26, 27]. Any physical Green’s function should obey several physical constraints. For example, a physical $G_0$ can be expressed in its single-particle form

$$G_0(x, x', \omega) = \sum_s \frac{\phi_s^*(x) \phi_s^*(x')}{\omega - \epsilon_s}, \qquad (C.1)$$

where $\phi_s(x)$ and $\epsilon_s$ are the single-particle wave functions and energies, respectively. Let us now suppose that the Dyson equation given by

$$G_0^{-1} = G^{-1} + \Sigma[G]$$

with $G^{-1}$ fixed, has two solutions of the form (C.1), $G_0$ and $G'_0$. The difference of the two corresponding Dyson equations can be written as

$$[G'_0]^{-1} = G_0^{-1} - \Sigma[G_0] + \Sigma[G'_0]. \qquad (C.2)$$

Now we write this equation in the basis labelled $\ell$ in which $G_0$ is diagonal. This basis exists and is the same for all frequencies because $G_0$ has the form given in equation (C.1). The right-hand side of equation (C.2) in the basis reads $\delta_{\ell\ell'}(\omega - \epsilon_\ell) - \Sigma_{\ell\ell'}[G_0] + \Sigma_{\ell\ell'}[G'_0]$. In the simple case of a static HF self-energy this means that the off-diagonal elements of $[G'_0]^{-1}$ must be static and therefore equal to their $\omega \to \infty$ limit. In this limit, for a spectrum $\epsilon_s$ that is bounded, these elements are zero and as a consequence $G'_0$ must be diagonal. This means that its wave functions are the same as those of $G_0$ and therefore $\Sigma_{HF}[G_0] = \Sigma_{HF}[G'_0]$. We finally obtain that $G'_0 = G_0$. We conclude that in this example the constraint (C.1) is sufficient to guarantee that there is only one solution. Hence, the problem of multiple solutions is related to the definition of the domain of representability of $G_0$.

## References

[1] Fetter A and Walecka J D 2003 *Quantum Theory of Many-Particle Systems* (New York: Dover)

[2] Luttinger J M and Ward J C 1960 *Phys. Rev.* **118** 1417

[3] Lani G, Romaniello P and Reining L 2012 *New J. Phys.* **14** 013056

[4] Tandetzky F, Dewhurst J K, Sharma S and Gross E K U 2015 *Phys. Rev.* B **92** 115125

[5] Fukutome H 1971 *Prog. Theor. Phys.* **45** 1382

[6] Seeger R and Pople J A 1977 *J. Chem. Phys.* **66** 3045

[7] Bauernschmitt R and Ahlrichs R 1996 *J. Chem. Phys.* **104** 9047

[8] Thom A J W and Head-Gordon M 2008 *Phys. Rev. Lett.* **101** 193001

[9] Kozik E, Ferrero M and Georges A 2015 *Phys. Rev. Lett.* **114** 156402

[10] Schäfer T, Rohringer G, Gunnarsson O, Ciuchi S, Sangiovanni G and Toschi A 2013 *Phys. Rev. Lett.* **110** 246405

[11] Molinari L G 2005 *Phys. Rev.* B **71** 113102

[12] Molinari L G and Manini N 2006 *Eur. Phys. J.* B **51** 331

[13] Pavlyukh Y and Hübner W 2007 *J. Math. Phys.* **48** 052109

[14] Hedin L 1965 *Phys. Rev.* **139** A796

[15] Berger J A, Romaniello P, Tandetzky F, Mendoza B S, Brouder C and Reining L 2014 *New J. Phys.* **16** 113025

[16] Martin P C and Schwinger J 1959 *Phys. Rev.* **115** 1342

[17] Onida G, Reining L and Rubio A 2002 *Rev. Mod. Phys.* **74** 601

[18] Chang E K, Rohlfing M and Louie S G 2000 *Phys. Rev. Lett.* **85** 2613

[19] Gross E K U, Dobson J F and Petersilka M 1996 Topics in current chemistry *Density Functional Theory* II vol 181 ed R F Nalewajski (Berlin: Springer) pp 81–172

[20] Sharma S, Dewhurst J K, Sanna A and Gross E K U 2011 *Phys. Rev. Lett.* **107** 186401

[21] Roessler D M and Walker W C 1967 *J. Opt. Soc. Am.* **57** 835

[22] Rigamonti S, Botti S, Veniard V, Draxl C, Reining L and Sottile F 2015 *Phys. Rev. Lett.* **114** 146402

[23] Chitra R and Kotliar G 2001 *Phys. Rev.* B **63** 115110

[24] Potthoff M 2006 *Condens. Matter Phys.* **9** 557
[25] Blöchl P E, Pruschke T and Potthoff M 2013 *Phys. Rev. B* **88** 205139

[26] Potthoff M 2003 *Eur. Phys. J. B* **32** 429

[27] Eder R 2014 arXiv:1407.6599

[28] Hohenberg P and Kohn W 1964 *Phys. Rev.* **136** B864

[29] Eric C and Claude L B 2000 *ESAIM: M2AN* **34** 749

[30] Hartree D R 1957 *The calculation of Atomic Structures* (New York: Wiley)

[31] Saunders V R and Hillier I H 1973 *Int. J. Quantum Chem.* **7** 699

[32] Pulay P 1980 *Chem. Phys. Lett.* **73** 393

[33] Pulay P 1982 *J. Comput. Chem.* **3** 556

[34] Cancès E and le Bris C 2000 *Int. J. Quantum Chem.* **79** 82

[35] Kudin K N, Scuseria G E and Cancès E 2002 *J. Chem. Phys.* **116** 8255

[36] Gulans A, Kontur S, Meisenbichler C, Nabok D, Pavone P, Rigamonti S, Sagmeister S, Werner U and Draxl C 2014 *J. Phys.: Condens. Matter.* **26** 363202

[37] Perdew J P and Wang Y 1986 *Phys. Rev. B* **33** 8800

[38] Piacentini M, Lynch D W and Olson C G 1976 *Phys. Rev. B* **13** 5530

[39] Uimonen A-M, Stefanucci G, Pavlyukh Y and van Leeuwen R 2015 *Phys. Rev. B* **91** 115104

[40] Kubo R 1957 *J. Phys. Soc. Japan* **12** 570

[41] Martin P C and Schwinger J 1959 *Phys. Rev.* **115** 1342