Two-Particle Self-Consistent method for the multi-orbital Hubbard model

One of the most challenging problems in solid state systems is the microscopic analysis of electronic correlations. A paramount minimal model that encodes correlation effects is the Hubbard Hamiltonian, which -- albeit its simplicity -- is exactly solvable only in a few limiting cases and approximate many-body methods are required for its solution. In this review we present an overview on the non-perturbative Two-Particle Self-Consistent method (TPSC) which was originally introduced to describe the electronic properties of the single-band Hubbard model. We introduce here a detailed derivation of the multi-orbital generalization of TPSC and discuss particular features of the method on exemplary interacting models in comparison to dynamical mean-field theory results.


I. INTRODUCTION
In strongly correlated electronic systems the development of many-body techniques is driven by the fact that a description of electronic properties in terms of an independent electron picture fails. Correlation effects result in a plethora of fascinating phenomena such as unconventional superconductivity 1-11, Mott metal-to-insulator transition 12-19, non-Fermi liquid behavior 20-22, or spin liquid phases 23-29 to mention a few. In many materials, correlations originate from a few partially filled orbitals around the Fermi level and, early on, a simplified low-energy description of those orbitals was proposed in terms of the Hubbard model 13, 30? , 31, which maps the electronic part of the full Hamiltonian of the interacting system onto an effective lattice model. This model is expected to capture the correlation effects of the original system but which is still too complex to be solved in the general case except for certain limits. Thus, one requires the development of elaborate approximate manybody methods.
Many promising schemes conceived to describe the electronic properties of correlated materials start from an ab-initio-derived effective non-interacting Hamiltonian where strong correlation effects are then added and treated within an approximate many-body method. Since most recent materials of interest are multi-orbital systems 32-38, an explicit multi-orbital formulation of many-body techniques is required.
Among the large variety of available many-body methods, in this review we will focus on the so-called conserving approximations in the Baym-Kadanoff sense 39 and 40. Those methods are thermodynamically consistent, i.e. thermodynamic expectation values can be obtained as derivatives of the free energy, and preserve many important physical constraints like the Ward identities for the collective modes and conservation laws for singleparticle properties. Still, they can differ from one another * zantout@itp.uni-frankfurt.de in how far they fulfill other physical constraints like local spin and charge sum rules or the Mermin-Wagner theorem 41. In what follows we shortly review some of these methods.
A very powerful and successful combination of ab-initio and many-body techniques is density functional theory plus dynamical mean field theory (DFT+DMFT) [42][43][44] where DFT provides a reasonable starting point for the electronic structure of the system and DMFT [45][46][47] introduces all correlation effects that appear in terms of a dynamical (frequency-dependent) but local selfenergy Σ(ω). Approximating the full momentum-and frequency-dependent self-energy by a local dynamical quantity amounts to restricting all correlation effects beyond DFT to a single site. While this local approximation has been very successful in explaining many experimentally observable properties of strongly correlated systems 47-55, there are materials where non-local correlations are non-negligible. This is the case, for instance, close to a phase transition, or in phenomena like the pseudo-gap physics in high-Tc cuprates 56-60. A straightforward way of including non-local correlation effects in the DMFT framework are the cluster DMFT method (CDMFT)47, 61-65 or the dynamical cluster approximation (DCA)62, 63, 66? -69, which explicitly treat short-range correlations between neighbouring sites by enlarging the unit cell to comprise multiple atoms of the same type, but are restricted in spatial resolution due to the large computational cost.
In general, many full momentum-dependent approximations that directly operate in the thermodynamic limit are available, both perturbative and non-perturbative in different quantities 62, 70, and 71. Most straightforward weak-coupling perturbative expansions in the electronelectron interaction approximate the one-particle irreducible vertex, i.e. the self-energy Σ. This one-particle vertex describes the renormalization of an electron due to the electron-electron interaction in the background of all other electrons arising from scattering processes. Such an expansion can also be done in other quantities like the screened interaction W . This is the case in the GW approximation 72-75 where the Dyson equation relates the unrenormalized single-particle Green's function G 0 and the single-particle vertex Σ with the renormalized Green's function G.
Another approach is to expand Σ non-perturbatively in the interaction, but perturbatively in the locality of the diagrams 70. The DMFT approximation is then the lowest order term in the sense of locality, since it approximates the one-particle vertex Σ to be local, but generates it from a summation of all diagrams that can be obtained from local propagators.
An alternative procedure is to approximate two-particle quantities like the irreducible vertex Γ, from which one-particle quantities can be derived that usually contain a richer structure of correlation effects. On the two-particle level, the irreducible vertex Γ contains information about two-particle scattering processes. Here, the Bethe-Salpeter equation represents a two-particle analogue of the single-particle Dyson equation, relating two-particle Green's functions like the bare and renormalized generalized susceptibilities with the twoparticle irreducible vertex 76 and 77. Methods like the Random Phase Approximation (RPA) or the Fluctuation Exchange Approximation (FLEX) 78-83 sum certain subsets of diagrams to approximate the twoparticle vertex, while DΓA 70, 84, and 85 approximates the vertex as a dynamical but local quantity, including all local diagrams. Further two-particle extensions for the vertex are for example TRILEX 86 and 87, QUADRILEX 88, dual boson 89 and dual fermion techniques 90 or GW+DMFT 91 and 92 which use the local DMFT solution and vertex as a starting point for a perturbation series to generate non-local diagrams.
In this review we focus on the Two-Particle Self-Consistent approach (TPSC). This is a method developed within the Baym-Kadanoff scheme 39, 40, 77? ? that retains the dynamical and non-local nature of electronic correlations, while using physical sum rules to obtain consistent values for all the quantities that are approximated. As described in Section II, one approximates the two-particle irreducible vertex Γ, usually depending on three frequency-and momentum-indices 71? ? ? ? ? ? , to be frequency and momentum independent, i.e. the two-particle irreducible vertex Γ is a mere constant. The vertex Γ is then determined by requiring the spin and charge susceptibilities to obey physical summation rules. This is in contrast to many approaches where the (approximate) irreducible vertices are obtained by solving complicated Parquet equations 70 and 71. From the equations of motion derived for the Green's function (also called Schwinger-Dyson equation) the self-energy can then be obtained from the bare interaction, bare Green's function, two-particle irreducible vertices and generalized susceptibilities. This local and static approximation of the vertex Γ in the two-particle sector results in a nonperturbative, fully frequency and momentum dependent single-particle self-energy Σ, which has been shown to be able to describe many electronic correlation effects including pseudo-gap physics and superconductivity63 and is by definition the sum of all closed two-particle irreducible skeleton diagrams. In the diagrammatic representation the bold lines correspond to full Green's functions while the single-wiggled lines are interaction vertices. In the TPSC approximation one assumes that all diagrams can be approximated by a diagram of first order where the interaction vertex is replaced by an effective irreducible interaction vertex (double-wiggled line).
93 as will be discussed in section III. In this review we present a detailed derivation of an extension of the standard single-orbital Two-Particle Self-Consistent method to the multi-orbital case as presented in 94 and discuss applications to model systems, as well as benchmarks with other methods.

II. SINGLE-BAND TPSC METHOD
Before we outline the history of the TPSC method and present a detailed derivation of the multi-orbital TPSC we would like to sketch in this section the main ideas of the single-band TPSC 41 and 93 that is formulated for a Hubbard model with on-site interaction U . Since TPSC is derived within the Baym-Kadanoff scheme 39, 40, 77? ? , we start with the description of the interacting system in terms of a Luttinger-Ward functional Φ[G] 95-97, which is a scalar functional of the dressed many-body Green's function G. Specifically, Φ[G] incorporates all closed two-particle irreducible skeleton diagrams constructed from G and the on-site interaction U . In general, however, Φ[G] cannot be evaluated explicitly but one can approximate it which is the idea of conserving approximations. In TPSC one assumes that at intermediate interaction strengths one can absorb the effect of diagrams of all orders into an effective irreducible interaction four-point vertex Γ that is local and static and appears only in the first order diagram (see fig. 1).
The effective irreducible interaction vertex Γ can be decomposed into a spin vertex Γ sp and a charge vertex Γ ch . These two vertices are then determined from the spin and charge susceptibilities χ sp/ch ; more specifically from the so-called local spin and charge sum rules, where n is the filling and n ↑ n ↓ is the double occupation. The single-band TPSC approach needs an ansatz for the calculation of the spin vertex Γ sp or, equivalently, of the double occupancy n ↑ n ↓ given by Having determined the spin and charge vertices Γ sp , Γ ch one uses the Bethe-Salpeter equation 76 and 77 for the two-particle Green's function to compute where χ 0 is the particle-hole bubble diagram −G 0 * G 0 calculated from the non-interacting Green's function G 0 . By construction, the self-energy Σ is computed from δG which equals some constant in the case of the TPSC Luttinger-Ward functional. This constant can be absorbed into the chemical potential and therefore no single-particle renormalizations take place. In the framework of conserving approximations one can further improve on this result by using the Bethe-Salpeter equation for the self-energy 41 and 77, where the non-interaction Green's function G 0 appears instead of the dressed Green's function G to preserve consistency with the TPSC Luttinger-Ward functional and where Σ HF is the Hartree-Fock self-energy and * denotes the convolution. This improvement yields an approximation that is not conserving in the strict Kadanoff-Baym sense 39, 40, and 77 but still retains conservation laws to a high degree 41. Additionally, one can further improve the self-energy by taking crossing symmetry of the two-particle irreducible vertex Γ into account 98 which gives Finally, one uses the Dyson equation to obtain the full Green's function G.

III. PREVIOUS FORMULATIONS OF THE TPSC METHOD
In this section we provide a brief overview on past developments of TPSC.
Early on, prior to the formulation of the two-particle self-consistent method, Vilk et al. introduced in 99 a simple self-consistent way of obtaining approximate spin and charge susceptibilities of the two-dimensional oneband Hubbard model without adjustable parameters. The ansatz was motivated by the local field approximation 100 and 101. This approach provided results comparable to quantum Monte Carlo simulations at weak to intermediate coupling strength and respected the Mermin-Wagner theorem that prohibits in two dimensions a spontaneous breaking of the SO(3) spin symmetry at finite temperature. A few advantages of this ansatz, which are also present in the TPSC equations, are the inclusion of short-range quantum fluctuations and finite temperature effects and the ability to reproduce Kanamori-Brueckner (KB) screening, which describes the reduction of the bare interaction in the spin channel 31. First applications to the single-band Hubbard model far from van Hove singularities revealed valuable insights into the spin and charge fluctuations 99, 102-104.
The first complete formulation of single-band TPSC was introduced in 105 and a very extensive and thorough presentation of it can be found in 41, 93, and 106.
The TPSC method was developed with the aim of fulfilling essential physical properties such as the local spin and charge sum rules, which is, for instance, violated by the FLEX approximation, and conservation laws of the spin and charge susceptibilities, like χ sp/ch (q = 0, iq n ) = 0 for q n = 0, which are also not fulfilled in the FLEX approximation. Actually, TPSC not only fulfills the Mermin-Wagner theorem and local spin and charge sum rules exactly, but it also satisfies Luttinger's theorem and the f-sum rule to a high degree (the deviation is of the order of a few percent) 41. Extensions of TPSC including the transversal particle-hole channel contributions to the electronic self-energy 98 and 106 preserve also the crossing symmetry of the irreducible four-point vertex Γ, while obeying spin rotational invariance.
The single-band TPSC method has been successfully applied to a multitude of physical phenomena described by the single-band Hubbard model 41, 105, 107-112. For instance, mostly in the framework of high-Tc cuprate superconductors, studies on the precursor antiferromagnetic bands, the pseudogap phase, superconducting transition temperatures, spectral and dynamical properties 113-121 and extent of quantum criticality 122 and 123 provided further insight on these materials and model systems. Specifically, this approach was used to study universal critical behavior, where TPSC was shown to be in the same universality class as the n → ∞ limit of the O(n) classical vector model 108 while n = 3 would correspond to the correct limit for a magnetic transition. The proof in the referenced publication is cutoff independent and has a lower critical dimension d = 2, where the correlation length grows exponentially. TPSC was used to estimate the extent of quantum criticality in the Hubbard model phase diagram 122 with applications to high-Tc cuprate superconductors 111. Interestingly, the same critical behavior is also observed ? in a special formulation of the DΓA which might be due to a possible similarity of both methods, where TPSC (static vertex function) is some limiting case of DΓA (dynamical vertex function).
However, TPSC fails to be a good approximation in the strong-coupling regime where the frequency and momentum dependence of the irreducible four-point vertex Γ becomes important while the TPSC vertex is, by construction, constant. The method is therefore not able to describe Hubbard satellites, but only precursors of antiferromagnetic bands and pseudogaps 41. Only very recently, the so-called TPSC+ approach was introduced in 71 where effective frequency-dependent vertex corrections were included in the TPSC equations to improve the results in the intermediate coupling regime.
The efficiency and reliability of TPSC has been discussed in comparison to a few other approximations like paramagnon and perturbation theories 124, the dynamical cluster approximation 114 and the FLEX approximation 125. For more detailed comparisons of TPSC to other many-body method we point to 70 and 71. In this context, in Ref. 126 a measure of the self-energy dispersion was defined emphasizing the role of local versus non-local correlations taking as reference TPSC versus DMFT calculations. Based on this concept, nonlocal correlation effects on topological properties in the Haldane-Hubbard model 127 and the bilayer "twistronic" 1T-TaSe 2 128 have recently been investigated.
The single-band TPSC has been also applied to the attractive Hubbard model, where spin and charge fluctuations are replaced by pairing fluctuations. The general scheme of the derivation remains identical, and the conserving and respective single-and two-particle properties are retained 129-133. The testing ground for this TPSC implementation was found in underdoped cuprates, with a focus on the interplay between superconductivity and pseudogap physics 130 and 134 and good agreement to QMC data in the weak to intermediate coupling regime was reported 134.
Further extension of TPSC included the effect of nearest-neighbour interactions V in the extended singleband Hubbard model 135 and the effect of pair-correlation functions derivatives which had been neglected so far 136 and 137.This formulation yields not only good agreement with QMC when V is small, but also in the limit where charge fluctuations are the main contribution, in contrast to FLEX and mean-field calculations.
Extensions to multi-site TPSC were introduced in Ref. 138 and applied to study the semimetal to antiferromagnetic phase transition on a honeycomb lattice 138 and, the metal to superconducting transition in organic superconductors 139-141.
The first extension of TPSC to a multi-orbital formulation was introduced by Miyahara et al. 142 where superconducting critical temperatures and gap symmetries with additional comparisons to RPA and FLEX were performed for high-Tc superconductors.
The authors of this review introduced in 94 an alternative multi-orbital formulation of TPSC that is presented in the next section. We will discuss the differences between both multi-orbital formulations in the next section. An extended study 143 on a larger class of iron-based superconductors focussed on a comparison of multi-orbital TPSC with FLEX and RPA and underlined the importance of non-local effects in those materials.

IV. DERIVATION OF MULTI-ORBITAL TPSC
In this section we present a detailed derivation of multiorbital TPSC 94. The aim is to formulate an approximation to solve the multi-orbital Hubbard Hamiltonian that does not violate conservation laws. To do so one starts from the Luttinger-Ward functional 95 and 96 and applies approximations on the four-point vertex as in fig. 1. In a conserving approximation one restricts all possible closed skeleton two-particle diagrams, i.e. diagrams that contain fully dressed Green's functions without explicit self-energy lines, to a subset of diagrams. From functional differentiation it is then possible to compute a selfenergy that is consistent with the chosen set of skeleton diagrams. If the chosen subset of diagrams leads to convergent series with a physical solution, one can be sure to have a conserving approximation in the Baym-Kadanoff sense. All expressions are given in Planck units where = k B = 1.

A. Definitions
The full Hamiltonian that is considered in this work is given by where t Ri− Rj αβ are all hoppings concerning orbitals α and β that are connected by lattice vectors R i − R j and we dropped the spin index since we assume a paramagnetic state without breaking of time-reversal symmetry. U αβ denote the onsite orbital-dependent Hubbard interactions and J αβ are the onsite inter-orbital Hund's couplings. The operator c α,σ ( R i , τ ) destroys an electron with spin σ in the α-orbital at unit cell position R i at the imaginary time τ and c † β,σ ( R j , τ ) creates an electron with spin σ in the β-orbital at unit cell position R j at τ and n α, . Note that we dropped the time dependence of the Hamilton operator since it is not explicitly time dependent. The multi-orbital Green's function for a lattice system is defined as Due to space-time translational invariance one can rewrite Taking advantage of the spatial translational invariance the Fourier transform of the expression above yields the Green's function in reciprocal space G µν,σ ( k, τ − τ ), where k is a reciprocal lattice vector. Furthermore, one uses the anti-periodicity of the Green's function (11) to define the Matsubara Green's function where the fermionic Matsubara frequencies are defined by ω n := (2n+1)π β . The kinetic part of the Hamilton operator can be diagonalized via a Fourier transformation which leads to a set of eigenvectors a b ( k) and corresponding where N orb is the number of orbitals and Thus, the non-interacting Matsubara Green's function reads From the Matsubara Green's function we can calculate the filling n via To abbreviate the notation we introduce the following convention for space-time vectors and functions: The multi-orbital non-interacting susceptibility is defined as where the bosonic Matsubara frequencies are given by q n := 2πn/β. The non-interacting susceptibilities are 4index-tensors. We have dropped the spin index since the non-interacting Hamiltonian is spin-rotational invariant.
B. Functional Derivative Approach for the longitudinal channel in the multi-orbital Hubbard model In order to be able to include the orbital-dependent interactions U µν , J µν in the preceding formalism one introduces an artificial scalar field φ µν,σ (1) that depends on the orbitals µ, ν and the position and imaginary time coordinate 1 (see eqs. (18)) as developed by Kadanoff and Baym in 39 and 40. This field is supposed to couple to the field operators of equal spin (longitudinal channel) like in the generalized partition function in the grand canonical ensemble where we use a short-hand notation for integrations and summation, i.e. symbols with a bar are summed or integrated over. Such a coupling that leads to spin and charge fluctuations in the particle-hole channel only are the most important modes in the repulsive Hubbard model 98 and 106. Although an inclusion of the transversal channel (φ couples to opposite spins) leads also to quantum fluctuations, it was not considered so far 94 and 142 (see section VI). Obviously, one can recover the usual physics by setting the field φ to zero. The advantages of such a field become visible in the following steps. First, a generalized Green's function G σ (1, 2; [φ]) =: In the framework of conserving approximation also higher-order derivative with respect to the field φ, are considered. This expression will be an important ingredient in the formulation of a self-consistency.

The self-energy and Dyson equation
In this framework the self-energy Σ can be defined implicitly via and leads to the Dyson equation 2. A self-consistent formulation for δG δφ So far, we only have an implicit equation for the selfenergy (eq. (23)) but the right-hand side of the equation is already known from eq. (22) and, therefore, one can rewrite where the notation 1 + and 1 ++ is used to make the expression well-defined: The ordering of the operators within the functions is defined via τ ++ 1 > τ + 1 > τ 1 with infinitesimal small differences between them.
The difficult part now is to evaluate the variational differentiation on the right-hand side of equation (25) which can be further simplified with the Bethe-Salpeter equation Before this step we present relationships between δG/δφ and spin and charge susceptibilities.

C. Spin and charge susceptibilities
One defines the charge susceptibility as the linear response to charge perturbations where generalized density operators are defined by and it was shown that certain sum rules that are called local charge sum rules can be defined in the paramagnetic phase, and for µ = ν, where the expectation values can be calculated from the Green's function: Those equalities can be directly proven from the definition (eq. (27)) and the Pauli principle. The importance of those strictly valid equations for TPSC will be explained later. Taking a closer look at eq. (22) one identifies which motivates the definition of a generalized three-point charge susceptibility that reproduces the previously defined charge susceptibility (eq. (27)) in the limit 3 → 1 + . This can be further evaluated with the self-consistent equation for δG δφ (eq. (26)) and leads to where one defined the (irreducible) charge vertex as Γ ch γβλρ (4, 5; 6, 7) := σ δΣ βγ,σ (4, 5) φ δG ρλ,↑ (6, 7) φ φ=0 .
Now, all equations can be put together to derive the self-energy as a function of the generalized three-point spin and charge susceptibilities. Starting from the implicit equation for the self-energy (eq. (25)) and multiplying with G −1 γδ,σ (2, 2) results in Making use of the equations of motion for the three-point susceptibilities (eqs. (45), (36)) and dropping the φ subscript to reduce notation yields to an expression of the where one defines and 1. The vertices Γ sp/ch,0 : Differences between RPA and TPSC Comparing the interaction vertices with RPA results 144 one observes a difference in the non-interacting vertices and This is because in RPA the spin and charge vertices δΣ δG are calculated while discarding every higher order contribution in eq. (25), i.e. the RPA vertices are constructed from functionally differentiating the Hartree-Fock selfenergy. which is a functional of the fully dressed single-particle Green's function G. It is the sum of all closed two-particle irreducible skeleton diagrams that can be constructed from G and the on-site interactions U, J. The first and second order functional derivatives of Φ[G] are related to the single-particle self-energy Σ and the irreducible vertex function Γ, respectively Any approximation to the Luttinger-Ward functional that consists of a diagrammatic truncation by taking only a certain convergent subset of diagrams into account can be shown to be conserving in the Baym-Kadanoff sense. The most simple example is the Hartree-Fock approximation, which only includes the first-order diagrams contributing to the Luttinger-Ward functional (see fig. 1 top line). The well-known FLEX approximation consists of taking a subset of diagrams that can be summed via the geometrical series 78 and 80. Another widely used approximation is the dynamical mean-field approximation, which approximates the Luttinger-Ward functional by taking into account only diagrams generated from local propagators Φ[G] ≈ Φ[G loc ] summed up to all orders by solving a local impurity model. In the multi-orbital TPSC formulation of 94 and 142 one proceeds as in the single-band TPSC 41 approach where the Luttinger-Ward functional is approximated by Φ = GΓG where Γ is static and local. In the same fashion as in the single-orbital case this leads to a constant selfenergy Σ(1, 2) = Γn(1)δ(1 − 2) and local and static (but orbital-dependent) spin and charge vertices Γ sp/ch αβγδ (1, 2; 3, 4) = Γ sp/ch The TPSC self-energy is closely related to the Hartree-Fock self-energy but with a renormalized effective interaction Γ. This is very different to the Hartree-Fock approximation, which equates Γ with the bare interactions U, J while in TPSC one a priori does not impose a limitation on the value of Γ. The constant self-energy contribution can be assumed to be already included in the input from density functional theory. To further improve the self-energy within TPSC one reinserts the bare Green's function into the self-energy equation eq. (48). This gives an improved self-energy where the collective modes enter while keeping the level of appoximation, i.e. susceptibilities are computed from G 0 and the input Green's function for the self-energy is also G 0 . This can be motivated by the fact that collective modes influence single-particle properties but the opposite effect is much smaller. Obviously, as argued in the single-band TPSC 41 one has to keep all equations at this level of iteration, i.e. singleshot calculations, since reiterating the susceptibilities or self-energy with the full Green's function would not only violate the local spin and charge sum rules but also be in contradiction to the assumption of static and local irreducible vertices Γ sp/ch . Since in 94 and 142 one uses input from DFT it is assumed that the static Hartree-Fock terms are already accounted for and therefore one drops them in the self-energy expression of eq. (48). Thus, equation (48)) and equation (55) lead to the final expression for the self-energy, A simple Fourier transformation leads to The resulting self-energy has the shape from paramagnon theories and allows for interpretations where the coupling to bosonic modes in the spin and charge channel is dressed with renormalized vertices due to higher order correlation functions. Since the derivation of the selfenergy was done only in the longitudinal channel the four-point vertex Γ does not fulfill crossing symmetry (see also Sec. VI).
One continues with the susceptibilities (eqs. (45), (36)) by inserting the TPSC spin and charge vertices from eq. (55) which simplifies the expressions to and thus By taking advantage of the index combination (αβγδ) → (αβ), (γ, δ) to reduce tensor equations to matrix equations it was shown that a Fourier transformation yields Analogously, the equation in the charge channel is The high-frequency behavior of the self-energy can be calculated via the local spin and charge sum rules but does not coincide with the exact high frequency result145 and 146. The reason lies within the constant vertices Γ sp/ch since a proper tail in the longitudinal particlehole channel can only occur if the interaction vertices renormalize to the bare interaction in the high-frequency limit41. The same statement is not true for the transversal particle-hole channel although averaging the TPSC self-energy expressions gave improved high-frequency behavior 98 which leads us to the assumption that this improvement is an error cancellation effect.

F. Ansatz equations for the irreducible vertices
So far, we only described how a local and static fourpoint vertex can simplify the self-energy and susceptibility expressions. In this subsection it is explained how to determine its value.

The spin vertex
In order to determine the renormalized vertices Γ sp within TPSC one makes use of the local spin sum rules (eqs. (42), (43) and 44). Unfortunately, these sum rules also include unknowns, namely the double occupations n µ,σ n ν,σ and the system of equations is therefore under-determined. Therefore, more information is needed to fix both the double occupations and the vertex Γ sp . The simplest ansatz is to do a Hartree-Fock decoupling for the right-hand side of ΣG (eq. (23)) for each expectation value and write a prefactor A, B in front to recover the result for equal time/position/orbital evaluation, i.e. 2 → 1 and α = γ: To recover now the original result for equal time/position we set Substituting A, B back into the ansatz equation (eq. (62)) and multiplying with G −1 γν (2, 2) one obtains To get now the renormalized vertices one performs the functional derivatives from eq. (46). This leads to the ansatz Thus, equation (65) together with the local spin sum rules gives us a set of equations to uniquely determine the spin vertex Γ sp . Moreover, in 94 one further sets Γ sp µννµ = Γ sp µνµν = Γ sp µµνν due to the symmetry of Γ sp,0 (eq. (50)). This is in contrast to 142 which is discussed in sec. IV I. Since the ansatz equations in (63) are not particle-hole symmetric it was shown in Ref. 94 that this can be enforced by symmetrizing those expressions, n µ,σ n µ,−σ n µ,σ n µ,−σ + n α,σ n µ,σ n α,σ n µ,σ + This enforcement of particle-hole symmetry is only one way to deal with the breaking of particle-hole symmetry in eq. (65). Alternatively, one can do a particle-hole transformation in the case of electron doping and keep eq. (65) in the case of hole doping as explained in 41. Note, that the ansatz fails for J = 0, U = 0 since in that case the ansatz for Γ sp µµνν renormalizes to zero and one ends up with a non-interacting double-occupancy n µ n ν = n µ n ν for all values of U .
It was shown in the original single-band TPSC that such an ansatz can be also motivated from the local-field approach for the electron gas and reproduces Kanamori-Bruecker screening 99. In principle, one would have to use the occupations n α,σ from the interacting system in eq. (66) but in Ref. 94 one uses the occupations of the non-interacting system and assumes that those are close to the occupations of the interacting system. For the specific material study in Ref. 94 one can show that this is a rather good approximation and consistent with the idea of single-band TPSC that spin and charge fluctuations are effectively calculated from the bare Green's function. This additional approximation is employed in order to avoid an additional self-consistency loop between the spin and charge susceptibilities and the interacting Green's function. Alternatively, one could use n µ,σ n ν,σ as an input from some other method like DMFT and solve the equations of the spin and charge sum rules directly (eq. (32), (43), (44)) without the need of the ansatz in eq. (65) and (66).

The charge vertex
In addition to the remarks on the self-energy tail at the end of section IV E the value of Γ ch µµνν , µ = ν, might be negative in TPSC if one would enforce the local charge sum rules without further constraints, since it allows for positive Σ(iω n ) which is unphysical and gives e.g. negative spectral weight.
To avoid this one restricts the elements Γ ch µµνν , µ = ν to be non-negative. In Ref. 94 one optimizes the charge sum rules (eq. (32)) by minimizing the difference between right-hand and left-hand side while keeping the charge vertex Γ ch positive while in Ref. 142 the negative contributions of the charge vertex are set to zero. Both schemes yield small errors in the local charge sum rules as is discussed in section V.

G. Flow diagram of multi-orbital TPSC
For a better overview we present in fig. 10 a flow diagram of the multi-orbital TPSC scheme as introduced in Ref. 94.

H. Internal accuracy check
The TPSC approach provides an internal accuracy check by combining the equation that relates the product of the self-energy and Green's function (eq. (23) ) to correlation functions that can be obtained from the local spin and charge sum rules. Evaluating the left-hand side of eq. (23) at equal orbital and time/position one gets In the case of TPSC this equation is fulfilled for G = G 0 since the self-energy is from a single shot approach 41 but in our multi-orbital case where the charge vertex cannot be calculated to fulfill the charge sum rules exactly we will also see a small deviation here (see sec. V).
I. Comparison to the multi-orbital TPSC formulation from Ref. 142 While the overall form of the TPSC equations presented in 94 (and also in this review) and the ones in 142 are exactly the same we discuss in this section the specific differences between both implementations and comment the implications on the final results. First, one observes that in the formulation of Ref. 94 the charge vertex Γ ch is determined such that its components are non-negative and minimize the error in the local charge sum rules. This is in contrast to the scheme in 142, where the negative contribution to the charge vertex component Γ ch ααββ , that eventually leads to unphysical negative spectral weight, is set to zero by hand. Since the optimization procedure is less invasive we expect that the charge sum rules in Ref. 94 are fulfilled to a higher degree than in Ref. 142, although in general the optimization procedure converges to values of Γ ch ααββ that are very close to zero compared to the other nontrivial elements of Γ ch (see e.g. fig. 4). For this reason we assume that the treatment of the charge channel is still very similar in both schemes.
Next, we find that Miyahara et al. employ the spin vertex ansatz in eq. (63) while we use the particle-hole symmetrized ansatz in eq. (66). It was already pointed out in 142 that the particle-hole transformation only leads to marginal changes for the high-Tc cuprate superconductor studied but in general this will strongly depend on the material. In any case, one needs to consider some kind of particle-hole symmetrization to meet the original ideas of this ansatz 41 and 100.
Another important difference between both formulations is the choice of the local spin and charge sum rules that are used to determine the spin and charge vertex Γ sp/ch . The sum rules presented here were chosen such that the right-hand side is only dependent on the double occupations and orbital occupations which allows for easier numerical root search but the outcome of both schemes should be the same.
Finally, an important difference is the form of the bare vertex functions Γ sp/ch,0 . While in the derivation we reproduced here the expressions explicitly follow from the Bethe-Salpeter equation and the equation of motion of the self-energy, Miyahara et al. employed the RPA expressions of those vertices. This leads to a difference in the vertex components αβαβ for α = β. Since those vertices enter the self-energy expression in eq. (57) which fulfills the self-consistency relation tr(ΣG 0 ) exactly (see eq. (67)) we conclude that the formulation presented here will be more accurate in terms of this sum rule.

V. APPLICATIONS TO A MODEL SYSTEM
In this section we present calculations on the twoorbital Hubbard model on the square lattice with nearestneighbor hopping t in order to show some general features of the method and point out where the limitations are. The interaction matrices U µν , J µν are in Hubbard-Kanamori form 31 and are thus derived from the interaction values U, J, i.e. U µν = U if µ = ν and U µν = U − 2J else while J µν = J. All presented results are calculated at k B T /t = 0.5 and half filling if not mentioned differently.
First, we present the double occupancy n µ,σ n ν,σ in fig. 2 where TPSC results are compared to DMFT. fig. 2(a) shows that the double occupations n µ,σ n µ,−σ decrease when increasing the interaction values U, J as expected from the local on-site interactions. This effect is further enhanced by enlarging the Hund's coupling J which favors equal spin states in different orbitals. Moreover, in fig. 2(b) and (c) the double occupations n µ,σ n ν,σ drop by increasing U/t when σ = σ . This effect is again due to the on-site repulsion that penalizes double occupation. On the other hand, the equal spin double occupation n µ,σ n ν,σ increases with U/t and is thus enhanced by increasing the Hund's coupling J. This is the counterpart to the decrease in n µ,σ n ν,−σ . The results between both methods are comparable, except at large U/J (see fig. 2(b) and (c)). Specifically, while DMFT shows a drop in the equal spin double occupation n µ,σ n ν,σ at low U and J it is always enhanced compared to the noninteracting value 0.25 in TPSC. Since this drop in DMFT occurs at small values of the Hund's coupling we conclude that the tendency to form a highspin configuration is not yet dominating over the direct intraorbital interaction U ≈ U in the itinerant phase. As U , and correspondingly J, is increased, the system becomes more localized and a high spin state developes, as evident by the subsequent increase in n µ,σ n ν,σ . This effect is not captured in TPSC, which we attribute to the Hartree-Fock-like decoupling in the ansatz equations for the spin vertex (see eq. (62)).
Next, we investigate the renormalized interaction matrices in the spin and charge channel Γ sp and Γ ch . We start with Γ sp since it is determined first in the TPSC procedure. Note that under the constraint of intermediate interaction values U/t and J/t to ensure numerical stability, it is always possible to find Γ sp such that the local spin sum rules (eqs. (42), (43) and (44)) can be fulfilled exactly. The spin vertex Γ sp is shown in fig. 3. Similar to single-orbital TPSC, fig. 3 shows the Kanamori-Brueckner screening, namely the saturation of the spin vertex Γ sp with increasing U/t. The screening of Γ sp µµµµ is stronger (compare U/J = 3 and U/J = 6) if the Hund's coupling is larger while the opposite is true for Γ sp µµνν , µ = ν. This can be understood from the fact that larger Hund's coupling favors double occupation of equal spins in different orbitals (see fig. 2) while it suppresses double occupation of opposite spins in the same orbital. This effect reduces the screening (see eq. (65)) in Γ sp µµνν if µ = ν and suppresses the screening if µ = ν.
Having obtained the spin vertex Γ sp and the double occupations n µ,σ n ν,σ the next step in the TPSC procedure is to determine the charge vertex from the local charge sum rules (eqs. (32)). The calculation of Γ ch from the local charge sum rules is doable in a systematic way but, unfortunately, is not straightforward and this might have several reasons as we will elucidate in the following discussion.  because they are equal to Γ sp µµνν , µ = ν. As is in singleorbital TPSC, the (Kanamori-Brueckner) screening, i.e. the saturation of the spin vertex with increasing U/t, is observed. Further, it is evident that larger Hund's couplings suppress the matrix elements Γ sp µµµµ while the opposite is the case for Γ sp µµνν , µ = ν.
have two important issues in the cases of large Hund's coupling J and in the case of large Hubbard interactions U . In the first case we observe that Γ ch µνµν , µ = ν, has the tendency of strong divergence which could possibly make the TPSC procedure numerically unstable due to overflow or precision errors when handling diverging matrix elements. In this study the diverging charge vertex components went up to values ∼ 10 11 and we had no numerical problems but we cannot exclude the previously mentioned errors in the general case. On the other hand fig. 4 (b)-(d) shows that for lower J/U values the matrix element Γ ch µµνν , ν = µ, has the tendency to go to negative values. This gives a negative contribution to the spectral function A( k, ω) and must be avoided. To do this, one imposes numerically the constraint of positive solutions and computes Γ ch such that the difference between left-hand and right-hand side of equations (32) is minimal. The result of this restricted charge vertex calculation is shown as filled symbols in fig. 4. As one can see this has in general only an effect if U/t is large and J/U is small. In the other cases the unconstrained determination of Γ ch gives purely non-negative results -except for the point at U/t = 3.8, 4.0 in fig. 4(a) which is due to the divergence of Γ ch 0101 . Moreover, fig. 4 shows that the effect of the restriction has a small impact on the matrix elements Γ ch 0101 and Γ ch 0000 . As in 142 we also find that in the cases of small Hund's coupling J compared to U the matrix element Γ ch 0011 is small compared to the other matrix elements and this further corroborates the procedure to set those matrix elements to zero.  show that for large Hund's coupling J/U the matrix element Γ ch 0101 diverges. On the other hand for small Hund's coupling the matrix element Γ ch 0011 has the tendency to converge to negative solutions which contribute to a negative spectral function A( k, ω) that is unphysical. Therefore, we restrict Γ ch to positive matrix elements and minimize the difference between right-hand side and left-hand side of the local charge sum rule equations (eqs. 32). The results for this constrained calculations are the filled symbols. The general tendency of the charge vertex is similar to the single-orbital TPSC with a diverging behavior which leads to the suppression of charge fluctuations in the system. Finally, we remark the general tendency of the charge vertex to diverge as a function of increasing U/t. This leads -as in the single-orbital case -to the freezing of charge fluctuations in the system. This non-perturbative divergence in the charge channel was originally interpreted as a precursor effect of the Mott transition 41, that is not only observed in TPSC but also in DMFT ? and DCA ? . More recent studies ? were able to show that such divergences can also occur independent of the Mott transition and that the formation of local magnetic moments can also lead to such divergences ? . Indeed, the divergence of the charge vertex can lead to interesting physical consequences such as phase-separation instabilities as observed in DMFT ? ? . Since those charge vertex divergences are an interesting and active field of research we give a sketch of how they appear within TPSC. For simplicity we go back to the single-band TPSC where one can see from the Bethe-Salpeter equation (eq. (4)), the local spin sum rule (eq. (1)) and the ansatz eq. (3) that the double occupation n ↑ n ↓ is determined as the root of the equation Obviously, an increasing value of Hubbard interaction U will lead to an increase on the left-hand side of the local spin sum rule hence the double occupation n ↑ n ↓ has to decrease in order to readjust equality because the non-interacting susceptibility is positive and real. If we consider the decreasing double occupation in the local charge sum rule, 1 βN q q 2χ 0 (q) 1 + χ 0 (q)Γ ch = n + 2 n ↑ n ↓ − n 2 , (69) we easily conclude that the charge vertex has to increase to compensate for the loss on the right-hand side. In the case of half-filling and the limit where the double occupation goes to zero, the right-hand side of the local charge sum rules also vanishes and thus necessitates a divergence of the charge vertex. Within TPSC the feature of the diverging charge vertex is therefore intimately connected with the Pauli principle from which the local sum rules were derived.
Since the restriction to positive matrix elements of Γ ch will necessarily lead to an error in the local charge sum rule it is worth to take a look at this error. The relative error of the local charge sum rule is defined as the sum of differences between right-hand side and left-hand side of the equations (32) divided by the respective right-hand sides. Those errors are shown in fig. 5. Clearly, the plot shows that for small interactions U/t the constrained and unconstrained determination yield the same (non-negative) Γ ch that fulfills the local charge sum rules exactly. On the one hand, if the rel. error in charge sum rules U/t Figure 5. The relative error in the charge sum rules for a constrained determination of Γ ch is shown in dependence of U/t. In the case of small interaction values U/t, Γ ch can be determined such that the local charge sum rules are fulfilled exactly. Only in the cases of large Hubbard interaction U/t in combination with small Hund's coupling to Hubbard interaction ratio J/U we observe deviations up to 50% or in the case where Γ ch µµνν undergoes a sign change (see fig. 4) a clear jump in the error happens at U/t = 2.6.
Hubbard interaction is large compared to the hopping amplitude t (U/t ≥ 3.5) and to the Hund's coupling J (U/J ≥ 3), large deviations appear that rise close to 50%. The largest deviations, as can be already guessed from fig. 4, happen where the restriction induces largest differences, namely Γ ch µµνν , µ = ν and therefore the largest deviations appear in the µµνν-charge sum rule. This also explains the jump in the error at U/t = 2.6 where Γ ch µµνν changes sign and the local charge sum rules cannot be fulfilled anymore. The fact that even the unconstrained solutions for Γ ch can lead to an error of up to 15% as for U/J = 3 at U/t = 3.8, 4 is due to the numerical instability of the diverging charge vertex. Although these numbers seem rather discomforting one has to keep in mind that they only occur at large interaction values U/t (U/t ≥ 3.5) and in those regions it is the spin susceptibility that dominates the contribution in the self-energy equation (57). One can therefore expect to obtain a self-energy that is still qualitatively and even quantitatively accurate as long as the original assumption of this method, namely that Γ sp/ch are constant, is still a good approximation.
To conclude the discussion of the spin and charge vertices Γ sp/ch , it is worthwhile to investigate the degree of local spin and charge sum rules' violation if one does not  1 βN q q,iqn χ sp + χ ch µµµµ ( q, iq n ) = 2 n µ − n µ 2 = 1, where we have used that our calculations are at half filling, n µ = 1. The result of the comparison is shown in fig. 6. TPSC fulfills the equal orbital sum rule by construction up to the largest values of U/t considered, namely U/t = 4. At U/t > 2.6 a small deviation becomes visible which is due to the fact that in this regime one has to restrict Γ ch which leads to small deviations in all local charge sum rules and consequently also in the equal orbital sum rule (eq. 70). On the other hand, this same calculation once more shows that the largest deviations due to the constrained determination of Γ ch happen in the µµνν-sum rule where µ = ν which is not part of the equal orbital sum rule and therefore the error is still very small compared to the ones presented in fig. 5. Nonetheless, fig. 6 demonstrates that no renormalization of the spin and charge vertex (RPA) can lead to strong . Quasi-particle spectral weight Z(kx, ky) for U/t = 2.5 and U/J = 4, n = 0.8 and kBT /t = 0.03. The quasiparticle weight shows a strong momentum dependence where the minima are located along in the vicinity of the four hot spots along the lines kx = 0 and ky = 0 lines. Z was obtained by linear extrapolation of the imaginary part of Σ(k, iωn) to iωn = 0. In black we show the Fermi surface of the system. violation of the equal orbital sum rule and therefore to a specific aspect of the Pauli principle. It was shown in Ref. 41 that this deviation is quadratic when the interaction parameters are small and we reproduce the same result in the multi-orbital case.
Next, we consider the results of the self-energy where we have picked k B T /t = 0.03 and n = 0.8 for all calculations. Note that we dropped the orbital index of all objects in the following since the self-energy and the Green's function are diagonal and orbital degenerate in orbital-index space. This means that all presented functions are always computed from orbital averages of the diagonal matrix elements. For all the following results the deviation from left-hand to right-hand side of the internal accuracy check equation (67) was always lower than 2.5% for both G and G 0 . For U/t = 2.5 and U/J = 4 the spin fluctuations in the system are already large and the momentum-dependence of the self-energy leads to strongest suppression of quasiparticle spectral weight Z( k) at the hot spots where the number drops down to 0.75 (see fig. 7). To compare the momentum dependent quasi-particle spectral weight Z( k) to the momentum independent result of DMFT we go into the local limit of TPSC by replacing Σ( k, iω n ) by its momentum average Σ(iω n ) := 1 . Local quasi-particle spectral weight Z is obtained from the local limit of the TPSC self-energy Σ(iωn) (see main text for definition). TPSC and DMFT agree well in the parameter range considered. of TPSC and DMFT is shown in fig. 8. Moreover, we compared the self-energy for U/t = 2.5 and U/J = 4 of DMFT and the local limit of TPSC in fig. 9. We observe that the momentum averaged TPSC and DMFT are in good agreement and show the same qualitative trend. The largest discrepancy is in the imaginary part of the self-energy at large Matsubara frequencies ω n which is due to the fact that the renormalized spin and charge vertices Γ sp , Γ ch do not rescale to the bare values -since the vertices are constant -to give the right constant in the 1/iω n expansion of the self-energy.

VI. SUMMARY AND OUTLOOK
We reviewed the multi-orbital extensions of the Two-Particle Self-Consistent approach that was originally formulated by Vilk and Tremblay 41 for the single-orbital case. This approach is able to enforce physical properties of the system such as the Pauli principle, the Kanamori-Brueckner screening, the Mermin-Wagner theorem and many conservation laws that are inherent to conserving approximations. However, deviations can occur by increasing the interaction strength U/t since by construction one considers the spin and charge vertices Γ sp , Γ ch to be orbital-dependent constants and this is not compatible with all local charge and spin sum rules (see eqs. (32)), (42), (43), and (44)). . Real and imaginary part of the local limit of the self-energy Σ(iωn) from TPSC and DMFT. Except for the large ωn-limit of the imaginary part of the self-energy TPSC and DMFT agree well. This discrepancy is due to the missing frequency dependencies of the spin and charge vertices Γ sp and Γ ch .
Although this constraint leads to deviations of up to 50% in the local charge sum rules, the latter has usually a small effect on the self-energy since the spin fluctuations are the major contributions in repulsive Hubbard models. All self-energy calculations in this review showed a deviation in the internal accuracy equation (67) of at most 1.5%. Most importantly for the self-energy calculations we saw that indeed the quasi-particle weight Z( k) shows momentum dependence where the strongest suppression of spectral weight occurs close to the X and Y point. TPSC is therefore a close relative of conserving approximations that not only respects the local spin and charge sum rules and the tr(ΣG) sum rule to high degree but does also not violate the Mermin-Wagner theorem. Moreover, due to the convolution expressions (eqs. (19) and (57)) it is very fast to perform TPSC calculations via Fast Fourier Transform and the convolution theorem. In addition, new numerical techniques were developed to further improve numerical efficiency 147? , 148. The method provides momentum-and frequency-dependent self-energies in the regime of weak to intermediate coupling strength and was successfully applied to multi-band high-Tc superconductors 94, 142, and 143. As for future work it might be interesting to extend multi-orbital TPSC to study further neighbor interactions as was already done for the single-orbital TPSC 135-137 or superconductivity as was already started in 142. From the method point of view it would be also worthwhile to incorporate at least some kind of frequency dependence to the spin and charge vertices Γ sp , Γ ch to be able to fulfill the local spin and charge sum rules to a higher degree and get the right high frequency behavior of the imaginary part of the self-energy. This would also improve the method in the sense that Hubbard physics would be visible. This is especially desirable because the feature of diverging charge vertices has attracted a lot of attention in recent years and contains itself interesting physics ? ? ? ? ? ? . Moreover, one could include self-energy contributions from the transversal channel like in 98 and 106 which was found numerically to improve on the tail of the selfenergy. Unfortunately, the form of Γ sp,0 does not have the simple form of equation (50) and one cannot straightforwardly set Γ sp µµνν = Γ sp µνµν = Γ sp µννµ as we did for the longitudinal channel. This also means that there are not enough local spin sum rules and ansatz equations to determine Γ sp .
Alternatively, it may be worthwhile to use double occupations from precise methods and calculate Γ sp from the local spin sum rules without need of an ansatz equation. This procedure was already suggested and used in the single-orbital case in 113.