Multi-orbital two-particle self-consistent approach -- strengths and limitations

a. The two-particle self-consistent approach has proven quite successful in describing the physics of the one-band Hubbard model in two dimensions. It satisfies several exact results. A generalization of this method for general multi-orbital problems is called for. Despite a few previous attempts, the literature has not converged yet towards a unique, benchmarked approach. This paper provides a generalization of the multi-orbital two-particle self-consistent approach, without relying on Dynamical mean filed theory, that it compares with previous works. It constructs 5 different spin-dependent Ansätze and 3 different spin-independent Ansätze that it benchmarks at half-filling.
b. Long derivations that need to be in the literature are presented.
c. They made it possible to close the system of equations for the vertices using sum rules and symmetries.
d. There is an honest declaration of the limitations of the method.
e. The method benchmarked the local vertices with Exact Diagonalisation and Dynamical Mean Field Theory, which are well-known in the field. (See however reservations in the report).

a. It needs further explanations or reference for some statements.
b. Many calculations remain obscure and even probably incorrect as written (typos?).
c. There is a lack details in the benchmark sections on the correspondence with the model parameters and the values of the D,P and C parameters. This makes the results and conclusions hard to validate.
d. There are some statements/conclusions that are unclear.
e. Some of the notation is unclear or undefined and even inconsistent.

This report is for the original submission of this paper. It was discussed with André-Marie Tremblay.
The two-particle self-consistent approach (TPSC) satisfies several exact results, like sum-rules, conservation laws and Pauli principle. It also satisfies the Mermin-Wagner theorem in two dimensions. Other methods, such as RPA violate several of these exact results and others like the renormalization group are much more difficult to implement and require more computing power. This makes TPSC a method of choice for the Hubbard model. Its generalization to multi-orbital models is an important and non-trivial next step that this paper intends to achieve. It is appropriate in principle for Scipost.
Unfortunately, in its present form, this paper cannot be published. It lacks justifications of some steps, presents benchmarks that are not always relevant, contains many typos and lacks some details that would be necessary for reproducibility. I suggest a partial list of changes and verifications of the derivations to make sure possible mistakes are avoided. After the authors respond to the constructive criticism below and after further review by me or some other referee, I hope publication can be recommended.

N.B. The changes in notation as one progresses through the paper make it difficult to find out what is a typo and what is a mistake. I have rederived most equations and indicated where I have questions, but I had to stop at some point.
1. In eq 1 : The interaction term should be 1/4 if you are assuming that the matrix element is antisymmetric, as is shown in eq 57 of Salmhofer and Honerkamp, that you reference. I personally think one should immediately start from the unsymmetrized form of the matrix element, and keep ½.
2. Also, in Eq.(1), why are you keeping such a general form of the interaction while you are assuming, without saying it explicitly, that the two-body interaction called U here is spin-independent (as follows from Eq.(2)).
3. After eq. 1 : You define the creation (annihilation) operator with position and imaginary time indices, but in eq. 1, all operators are not explicitly imaginary time-dependent since the Hamiltonian, as usual, is imaginary-time independent. Also, you did not define $\tau_1$ as imaginary time. I know you need to define creation (annihilation) operators with all the position, spin, imaginary time and orbital dependences since you need it later.
4. In eq.3, the term C contains $S \cdot S$ which includes $S_z S_z$, so it does not contain only spin flips that are given by $S_x S_x+S_yS_y = (S_+S_- + S_- S_+)/2$. Shouldn’t it be impossible to have same spin for all the creation (annihilation) operators, because then it isn’t really a spin-flip? The name seems inappropriate.
5. Throughout the derivations (eq 4 and after) : Why did the order of indices change? Please keep the same order of numbered indices as in previous equations, otherwise it makes the derivation hard to follow.
6. You also use two different types of notation : an explicit one (e.g. $\Sigma_{o_1 s_1}(r_1, \tau_1)$) and an abbreviated one (e.g $
\Sigma_1$). First, I would define what $\Sigma_1$ stands for. Second: I would suggest that you keep the same notation throughout, at least within the same eq. I do know how long those equations are, so you try to keep them more concise, but the most important thing is to explain when and how you use a short-hand. (For example, in the appendix, you use both at the same time, and when there are deltas, it is not clear whether these are delta functions for all variables, or just for position and time or otherwise.)
7. In eq. 6 you seem to change notation from position and imaginary time dependent creation (annihilation) operators to non-explicitly dependent ones: Is this to make the eqs. simpler and easier to follow? Please comment on that.
8. Going from eqs. 5 to 6 : I am not sure how you handled the switch in operators and the indices of U. Could you elaborate on that? When I try to do it, I switch the two annihilation operators of $o'_{3} s_1$ and $o'_{2} s'_{4}$, and because they anticommute, only the sign changes. Then I switch $o'_2$ and $o'_3$ for both U and the operators. That leaves me with two identical four-point function but a sum of $U_{o_2' o_3' o_4' o_1}$ and $U_{o_3' o_2' o_1 o_4'}$. Although, because of the symmetry used, they are equivalent like you said and it cancels the 1/2, but I don’t have any $\delta$ function or $(1-\delta)$. I think they should normally come from the model. And, by the way, it would have been easier to start from an interaction matrix element that was not antisymmetrized.
9. Definition of $\bar{\delta} $: “$\bar{\delta}$ is a short hand for $1-\delta$”: In eq. 6 I read this as :
$\bar{\delta}_{s_1,s_4'} \bar{\delta}_{o_2',o_3'} = (1 - \bar{\delta}_{s_1,s_4'} ) (1- \bar{\delta}_{o_2',o_3'} )$ : Which cancels all the aligned spin four-point functions, even when the orbitals are different. Thus, no Hund’s coupling of aligned spins on different orbitals. In the benchmark, this is the term on the right in eq 50. Hence, if this above definition is the one you mean, I think it makes the benchmark a lot less accurate, and that this would explain why the case of the Kanamori Hamiltonian works better, because the other terms in the Kanamori model are not neglected by this. Although, when I follow your derivation, I think you mean it in this way $\bar{\delta}_{s_1,s_4'} \bar{\delta}_{o_2',o_3'} \rightarrow (1 - \delta_{s_1,s_4'} \delta_{o_2',o_3'} )$, so that it respects the Pauli principle, but this should come from the model directly, as in eq 50. The only zero term is the case where two electrons are destroyed on the same orbital with same spin. Was this your intention? Could you clarify, or add a section in the appendix with the whole derivation in details?
10. Another way to see that the product of delta functions that appear in eq.(7) are incorrect (unless redefined as above) is to look at the numerator of eq.(9). If there is $\bar{\delta}_{o_1,o_4}$ as a prefactor of D in eq.(7), it must be a prefactor for the whole numerator of eq.(9). The numerator of eq.(9) is correct I think, but not the delta functions in eq.(7).
11. In eq. 7 : If we follow your delta, why are there $1- \delta$ for all terms except for the P term? Could you clarify? When I redid this calculation, nothing mathematical would have removed a $1-\delta$ over orbitals.
12. In eq. 7: I think there is a typo : There is a sum over spin s_4', but primed variables already means that they are summed over in your notation.
13. Eqs 9 through 14 : I do realize the ambiguity in the spin or spin-independent Ansätze, I have three interrogations upon this realization: the first one is that the SU(2) Symmetry on the interaction U allowed to remove the spin-dependency of the interaction term, why is it ok to reinsert it? The second one is how can you differentiate the two different interaction terms of eq 50 with a spin independent interaction? The third one is that while trying to recover this Ansätze, I came upon the realization that we have to assume that the Ansätze are all decoupled, for D, C and P what is your insight on that?
14. In eq.8 the global sign for P is incorrect if eq. 7 is correct.
15. Why are the indices of D interchanged in going from Eq.(8) to Eq.(9)?
16. Eq. 15: Please comment on why you did not put any $(r,\tau)$ dependence.
17. Eqs. 15 and 16 : Only label one number for the whole eq.
18. Eqs. 21 and 22, you analyze the functional derivative contribution of the Ansätze, but only for one term of D and one of C. What about the P coefficient? Also, is it the same result in the case where the Ansätze is spin-independent?
19. What assumptions are involved in claiming that the functional derivatives on the D Ansätze coefficients of eq. 15 vanish in eq. 21.
20. Eq. 21 : I think the left-hand side of the eq. needs to have a sum over s_4? Also, I think that normally, there should be a sum of the orbital indice $o_4$ of $\tilde{D}_{o_4, o_1}^{s_4, \sigma}$ , which is not the same index as the one in the Green’s function of the functional derivative.
21. Eqs. 23 and 24: Following your expression for the self-energy in eqs. 15 and 16 and the definition for Vertex in eq. 20, the sign in both eqs. 23 and 24 should be the opposite.
22. In eq. 34 through 38, it took me a while to understand that 1/2 meant either 1 or 2 and not spin half. Could you maybe just mention a clarification? Once I understood, the notation made sense, but it would be easier for a first-time reader if you put a small explanation or example.
23. In Section II.A : You write “We include both the particle-hole symmetrized and the usual version of the Ansätze in our implementation.“ Does usual mean the one you developed above? What is the particle-hole symmetrized Ansätze? Did you take it from another reference? If so, please add the reference and the eq. number, if you symmetrized your Ansätze, please show the final eq. of those.
24. In Section II.A : You write “Apart from these differences, the sum rules utilized in Ref. [87] and Ref. [88] are a subset of the sum rules employed here, and we checked that by constraining the eqs. we can reproduce the previous results.”. Could you add an appendix on those verifications, what are the constraints that you added? I obtained the same sum rules as you did, but not the ones in the reference you mention.
25. After eq 50 : You talk about how U=U', but you did not mention U' before, please define it or refer to its definition.
26. Between eqs. 50 and 51, about half-filling : What is the definition of your half-filling? One electron per orbital, or per site? Also, could you explain a little more on how you obtain eq. 52 by explaining the values of $n_{o_i,o_j}^{s_i,s_j}$ in the half-filling definition? It could be added in an appendix.
27. Eq 49 : It looks like the numbered orbitals are mixed up, because the right-hand side of the eq. does not close a bubble if we follow your definitions in the appendix.
28. Eq. 57 has two $c^\dagger$ in the hopping term. Is this intended? Is it Hermitian?
29. Section III : You did not mention what was the size of the real-space and reciprocal space used to do TPSC calculations. In figure 3, you seem to be doing TPSC for two site, but even in the single-orbital case, TPSC is not valid below a 64 site x 64 site, and normally it is more used in 128x128 or 256x256, and even 512x512 in the low temperatures regimes. Could you please mention what is the system size that you are using? If it is different for different benchmarks, please mention each time. Also, if it is a 2-site only, then it is expected that TPSC fails, yes.
30. Section III: Benchmarks and models : It could be a good thing to recall the different interaction terms of eq. 7 and give the values you used in the calculations. For both eqs. 56 and 50. Also, in the case of spin-independent Ansätze, how did you differentiate the (U-2J) and (U-3J) in eq. 56?
31. Section III.A, page 8, second paragraph : I think you should add a reference on the results of TPSC for the single orbital model.
32. Section III.A, page 9, first paragraph : The sentence : “As an illustration, let us considering a dimer at half-filling with a Hubbard- Kanamori interaction.” I think you mean “[…] let us consider a dimer […]”.
33. Section III.A, page 9, first paragraph : About the illustration of a dimer at half-filling with a Hubbard Kanamori interaction, you mean the $H_{int}$ of equation 56? If so, I think it would be best if you moved that to section III.B. where eq. 56 is already introduced. If not, could you add the details (eq. number): which interactions exactly you are referencing to? Also, could you show what are the different TPEVs? Which one of the TPEVs is represented by the Hartree-Fock decoupling? TPSC is not expected to be valid in the atomic limit, t=0. So your comment about the fact that the basics ideas of TPSC breaks down is correct, but this was also expected, you could add a reference, from many TPSC papers.
34. In section III.A : after eq 50 : you write “This interaction in combination with the kinetic term leads to a vanishing Σ0” what is your $\Sigma^0$? The non-interacting self-energy? The non-interacting self-energy should be zero, yes, but because the non-interacting green’s function is the one that solves the non-interacting Hamiltonian, which solves exactly, thus there is no self-energy in its expression. Please explain what is your $\Sigma^0$.
35. In section III.B., in the second paragraph, you explain why the higher values of J works better for TPSC while looking at the different groups of spins. You talk about the eigenvalues, but there is not further explanation on how to obtain those eigenvalues. I guess it would come from the diagonalization of the Interaction part only of the Kanamori Hamiltonian, but it would be nice to have more description about that, in an appendix perhaps, or maybe you could add a reference where it is well explained?
36. In section III.B. third paragraph, you write “ This also indicates that this specific part of TPSC performs better at lower temperatures (since it is basically a ground state targeting approach) - however, the local and static approximation of the vertex becomes more and more inappropriate at lower temperatures which is why TPSC typically fails at low temperatures.”
TPSC is not “basically a ground state targeting approach”. It is actually far from it. It is not valid in 2D deep in the renormalized classical regime, which is a region of the phase diagram in temperature and density that starts at finite temperature and goes to zero temperature.
Also, each part of the sentence separated by “-” are two contradicting statements. In general, I would say that TPSC can reach lower values of temperature than DMFT, but one needs to be careful in the grid definition, because the correlation length should be smaller than the lattice size to capture long wavelenght correlations. A lot of those explanations can be found in [Gauvin, Phys Rev B 108, 7 (2023). DOI : 10.1103/PhysRevB108.075144].

37. Appendix A, eq. A23 : Normally, the two derivatives w.r.t. $\tau_a$ and $\tau_b$ for $G(\tau_a, \tau_b) $give a different order for $\Sigma G$ : one is $\Sigma G$ and the other is $G \Sigma$.
38. In eq. A14, what does $S^0$ superscript mean?
39. The definition of susceptibilities in eq. (A11) does not lead to the matrix structure in eq. 25.

Ask for major revision

In this manuscript, the authors present a revised formulation of the multi-orbital Two-Particle-Self-Consistent (TPSC) method, addressing inconsistencies in a previous attempt. They compare the numerical results of the new formulation with those from the original approach and other many-body methods, including DMFT, ED, and D-TRILEX.

While the manuscript may warrant publication in some form, there are several open issues that, in my opinion, need to be addressed before I can recommend it for publication.

Following the Editor advice, I have revised the v2 version that can be found on this link https://arxiv.org/abs/2410.00962.

1.) The authors propose two different formulations of the multi-orbital TPSC method, referred to as TPSC3 and TPSC5. I am curious why this distinction arises. If one starts with the assumption that the two-particle Green's function is constructed as its disconnected part multiplied by a time-independent constant—chosen to reproduce the exact same time and position average—shouldn’t this lead exclusively to TPSC5? If this is not the case, could the authors clarify the source of this ambiguity? Additionally, a diagrammatic representation of the two formulations in terms of vertices could be helpful for less-expert readers to better visualise and distinguish the different approaches.

2.) Another point of concern is that the authors compare the TPSC results on the square lattice with methods that are not exact in two dimensions. Could the authors elaborate on the specific limitations of DMFT and D-TRILEX in capturing physics in two-dimensional systems? Additionally, could quantum corrections in two dimensions further reduce the double occupancy values? Providing clarification on these aspects would strengthen the context of the comparisons and the interpretation of the results.

3.) An interesting point raised by the authors concerns the limitations of the theory in the half-filled SU(4)-symmetric case at J =0 and with zero inter-orbital hopping. This is a peculiar case where the TPSC formulation clearly breaks down, and the authors acknowledge and discuss this, which is an important critical observation. However, this breakdown occurs in a very specific limit. I believe that for the SU(4) case at different integer fillings (i.e., n=1,3), variations in the chemical potential relative to the half-filled case could lead to differences from the single-band model. How might these differences affect the comparison between single-orbital and multi-orbital TPSC?

4) One of the strengths of TPSC is its ability to compute a non-local and frequency-dependent self-energy with reasonable computational cost. However, the authors do not discuss how the new formulations might alter this quantity compared to the previous approach. I believe the authors should address this by calculating the self-energy numerically and comparing it with other many-body methods capable of evaluating this quantity such as DMFT and D-TRILEX.

5) I am uncertain whether comparing TPSC against ED provides meaningful quantitative insights for its potential applications. TPSC is intended for systems in the thermodynamic limit within the weak-to-intermediate coupling regime. A dimerized state may emerge in one-dimensional systems or when interaction strengths exceed the range where TPSC is effective. Additionally, a two-site model lacks a bulk structure, meaning there is no distinction between open-boundary (OBC) and periodic-boundary conditions (PBC). Since TPSC is typically applied to models with PBC, using a cluster (e.g. four site with PBC in 2D) that respects these conditions might offer a more appropriate basis for quantitative comparison.

Here there are some minor aspects that the authors could address for improving the presentation of their paper.

(i) Why does Figure 2 appear before Figure 3, despite the first reference to Figure 2 occurring after the first reference to Figure 3 in the text?

(ii) The first three panels on the left in Figure 3 depict the same quantities. I suggest displaying only one panel and noting in the caption that the three averages should be identical.

(iii) In the legend of figure 4 Ref. 88 is indicated. Did the authors mean Ref.90?

Ask for major revision