Unphysical and Physical Solutions in Many-Body Theories: from Weak to Strong Correlation

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 non-linear 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.


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 Σ. The power of this approach resides in the fact that even a low-order approximation for Σ yields contributions to all orders in v c . Following Luttinger and Ward [2], Σ 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][6][7][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 G. 0 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 In section 3 we summarize our findings.

Theory and discussion
2.1. The OPM and the map To analyze the problem we use the so-called OPM [11][12][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 manybody 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]. i.e., the usual case, where y 0 is set by the system, and one searches y. The Dyson equation has two solutions with the rescaled quantities Y = y/y 0 and = V uy .
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) The sign of the square root is determined by the continued fraction in the iterative procedure 9 .

TDDFT and the map c c  0
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 χ [19]. Here we use TDDFT within the so-called bootstrap approximation proposed in [20]. The corresponding Dyson equation for frequency ω reads where χ 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 ò=1−v c P, where the irreducible polarizability P obeys a Dyson equation with kernel c xc 0 In this case, the appropriate scheme is scheme (II). This is explained in appendix B. 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][24][25]. For the LW functional to be properly defined, the map ¬ G G 0 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 y 0 unambiguously leads to = z y. . Because -V 2 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 Y 0 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 G 0 we need two different iteration schemes: one for 0 < V < 2 and the other for V > 2. This is different from the map 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  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] as critical quantity that changes sign at the crossing 11 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 twoparticle 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 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 G, 0 one could obtain the same dressed Green's function from two different G 0 and, hence, from two different external potentials 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 G 0 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 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 d d S 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   uy 0 2 1, 2 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 14 . To lowest order, this corresponds to HF, =s uy, HF for strong interaction. We call this functional strong-interaction HF (SIN-HF) 15 . Both self-energies yield two solutions. We report the physical solution for these two approximations in figure 3.While 11 The polarizability is defined as the derivative of G with respect to an external potential. Therefore there is only one branch. 12 Negative spectral functions can also be caused by certain approximations [39]. This is a different problem. 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 G 1 0 1 0 [40,41]. 14 Note that since  y 0 for  ¥ u the convergence of perturbation theory improves with increasing u. 15 Equations (12) and (13) suggest that the exact self-energies in the weak (WIN) and strong interaction (SIN) regimes are related by S = -S y 1 ,

SIN WIN
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. 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 = + - 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.

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 G 0 a single iteration scheme suffices to obtain the physical solution for all interaction strenghts. Instead, for the inverse map ¬ G G 0 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. 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][31][32][33][34][35] where   a 0 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 g = ---| | Ue i j i j 0.5 (with i, j running over the sites). Without mixing (α = 0.0) the HF total energy is linear in the interaction, but only for  U 5. For U > 5 the scheme oscillates between