Reply to Anonymous Report 1

Report

The authors give a thorough account of the D-TRILEX formalism, developed for a rather general, multi-band lattice model. Details of the implementation are discussed, and the method is tested on several examples, and compared directly to other methods (DMFT, DiagMC@DB). In cases where no benchmark data is available, the physical picture emerging in the tests is analyzed in some detail and is found to be consistent with the findings of previous works.

The subject of the paper is important. The development of methods that capture non-local correlations in 2D lattice models falls under central tasks in condensed matter theory. Added to that, the ability to treat multiband models may prove essential for the description of numerous materials of interest. In my view, the work presented in this paper is very well motivated.

The numerical results appear very convincing: Figs.2-4 demonstrate that D-TRILEX performs consistently better than DMFT and gives results which are comparable to (I imagine) much more expensive DiagMC@DB. The amount of data presented also indicates that the computations are, indeed, feasible and relatively cheap.

<<A: We thank the Referee for the in-depth analysis of our work, as well as for the useful suggestions and comments.>>

However, I find that some important points are not sufficiently emphasized. As far as I understand, one of the big advantages of the D-TRILEX method is that it does not require retarded interactions in the impurity problem (unlike TRILEX). There are no mentions of this in neither the introduction or the abstract. In the conclusions, it is stated that even frequency-dependent interactions in the original model would not change this desirable property of the D-TRILEX method. This should be more elaborated as it may not be immediately obvious to a casual reader.

<<A: We think that a more correct statement would be that one of the advantages of the D-TRILEX method is that it can be formulated on the basis of an arbitrary interacting reference problem. The reason is that in our approach the splitting between the impurity part and the rest of the lattice action can be introduced arbitrarily by adding and subtracting the fermionic and bosonic hybridization functions, as highlighted in the derivation of the partially bosonized theory (Appendix A.1). For this reason, even if the original problem contains the non-local interaction, the reference problem does not necessarily have to (but, in principle, can) contain the corresponding frequency-dependent bosonic hybridization. To our knowledge, our ideology if different from the one of GW+EDMFT approach, where the frequency-dependent bosonic hybridization function is always considered when taking into account the non-local interaction. The TRILEX approach is derived is a similar way as the GW+EDMFT method. For instance, in Ref.[PRB 97, 155145 (2018)] the authors say ``This simplified TRILEX version can be seen as a GW+EDMFT-like scheme which, however, can treat simultaneously both charge and spin fluctuations''. Nevertheless, we don't know whether taking into account the bosonic hybridization is necessary in TRILEX when the non-local interaction is considered. This question should be addressed to the authors of the TRILEX method. We comment on the use of the bosonic hybridization function in the paragraph after Eq. (4) when discussing different reference systems.>>

I also find that the benchmark in the case of the extended Hubbard model (Section 5.2) is out of place. There have been previous implementations of the D-TRILEX method which could have treated this same model (e.g. Ref. 91). I do not understand why the multiband generalization is needed for this benchmark. It would have made more sense to focus on actual multiband models in the main text. At least, the authors should have emphasized that Fig.4 does not actually test any of the new generality of D-TRILEX developed in this work.

<<A: Previous works on D-TRILEX [PRB 100, 205115 (2019), PRB 103, 245123 (2021), npj Comput. Mater. 8, 118 (2022)] indeed suggest that the method is able to account for the effect of the non-local interactions. However, no thorough benchmarking of the results for the extended Hubbard model has been performed so far. For this reason, in this work we decided also to investigate the performance of D-TRILEX in the case of a single-orbital extended Hubbard model on a square lattice with the local U and the nearest-neighbor V interactions between electronic densities. The aim of the current work is to provide a comprehensive discussion of the D-TRILEX method and show how the method performs against exact methods. For this reason, we believe that this Section fits well within the scope of the manuscript, especially since the numerical calculations in the entire work including the ones for the single-band extended Hubbard model have been performed using the code implementation discussed in this manuscript. We have clarified this point at the beginning of the Section 5.2 of the revised manuscript.>>

I also find a general lack of discussion regarding the systematic errors in the theory. One of the main touted strengths of D-TRILEX is that the expression for the self-energy and polarization in D-TRILEX satisfy symmetries that are not satisfied in standard TRILEX. In standard TRILEX, the expressions Σ=λGWΛ and P=λGGΛ (featuring only one renormalized vertex Λ) allow for a straightforward formulation of systematic corrections to the theory: In the cluster TRILEX, the expressions for Σ and P are of the same form as in the single-site theory, and as Nimp→∞, one recoveres the exact result . It is not immediately clear that the single-site D-TRILEX can be systematically improved towards the exact theory. What is the Nimp→∞ limit of D-TRILEX? I think that in the introductions and conclusion, the properties of D-TRILEX should be discussed in terms of control parameters and sources of systematic error. Is there a series of systematic improvements to D-TRILEX which leads to the exact solution of a model?

<<A: We acknowledge that this is an important point to address. First of all, as we write at the beginning of Section 3, ``the introduced effective fermion-boson action (5) allows for the calculation of observable quantities by using diagrammatic techniques. This goal can be achieved either by performing approximated diagrammatic expansions or by applying exact numerical methods such as the diagrammatic Monte Carlo (DiagMC) scheme. In this section we discuss the simplest diagrammatic approximation, namely the D-TRILEX method, that represents a feasible approach for actual calculations in the multi-band case.'' For this reason, the developed approach allows one to systematically improve the theory by including more diagrammatic contributions. For instance, one can consider much more complex diagrammatic structures that account for vertical (transverse) insertions of momentum- and frequency-dependent bosonic fluctuations, which are present in the DiagMC@DB approach but are not considered in ladder-like dual approximations including the D-TRILEX approach (this sentence is added to the revised main text). In this case, the polarization operator of D-TRILEX will also take the form P=λGGΛ, where Λ will be a vertex renormalized by the transverse fluctuations. However, it should be pointed out that these diagrammatic contributions can be considered only by sacrificing the simplicity of the computational scheme proposed here. In fact, the D-TRILEX approach has one of the most complex diagrammatic structures that still allow for feasible numerical calculations in the multi-orbital case. Another systematic improvement can be achieved by ``tuning the parameters of the reference system, especially in situations where the DMFT impurity problem does not provide a good starting point for the D-TRILEX diagrammatic expansion. For instance, the reference system can be improved by considering a cluster of suitable size, or by introducing more appropriate fermionic and bosonic hybridization functions.'' This discussion has been added to Conclusions. Regarding the N→∞ limit, if we consider a cluster of size Nl (with some requirements on symmetry) as the reference system, taking Nl →∞ also leads to the exact solution in our method. The reason behind this statement is that D-TRILEX is derived as a diagrammatic expansion in terms of dual Green's functions Gd=G-g on top of an exactly solved reference system. Here, G is a cluster DMFT Green's function and g is the Green's function of the reference system. In the limit of Nl →∞, the reference system coincides with the lattice problem. This implies that G=g and also means that in the considered limit all diagrammatic contributions in D-TRILEX vanish.>>

Can the extended Hubbard model (discussed in Section 5.2) be rewritten in a 2-site supercell notation and then solved straightforwardly by D-TRILEX as a multiband model? If possible, such benchmark would certainly be of great interest for the readers, and would greatly strenghthen the paper. In fact, in the introduction it is stated that the authors formulate D-TRILEX "in a multi-orbital and multi-site framework", yet no actual benchmarks with Nimp>1 were given, to my understanding.

<<A: Yes, the extended Hubbard model on a square lattice can be rewritten in a supercell notation and solved in the same way as a multi-band system. However, we think that these calculations would not add any new results. In fact, a similar calculation has already been performed for a bilayer square lattice (Section 5.4). In that case we considered two equivalent representations for the model, namely with either one (Nimp = 1) or two (Nimp = 2) lattice sites in the unit cell, and showed that the results for these two identically coincide. As to the second part of the question, we realized that we did not clearly state what we mean by the expression ``multi-site'', so we decided to add the following clarification after Eq.(28): ``In this context, Nimp is the number of independent impurities in the unit cell of the reference problem. Note that the case of Nimp > 1 corresponds to a collection of impurities, as explained in Ref.[PhysRevB.97.115150], and not to a cluster of Nimp sites. If the impurities are all identical, then the reference system reduces to a single site impurity problem. If some of them are different, it is sufficient to solve an impurity problem only for the non-equivalent ones. In the multi-impurity case, fluctuations between the impurities are taken into account diagrammatically in the framework of D-TRILEX approach. On the other hand, a cluster reference system corresponds to a multi-orbital problem with Nimp =1. In this case, Nl is the total number of orbitals and sites of the considered cluster. The separation between orbitals and sites that we introduce is useful to reduce the computational complexity when addressing problems with several atoms in the unit cells.''>>

In general, I would suggest that the authors expand the introduction and conclusions so as to rely less on the reader's ability to sift through very long formalism. Main points should be accessible in a quick reading of the paper.

Requested changes

1 - in the introduction and conclusions, emphasize more and elaborate on the properties of the impurity problem. If no retarded interactions are needed in the impurity problem, why is that? Under what conditions is this no longer true?

<<A: As we explain above, our method relies on the arbitrary splitting between the reference problem and the rest of the action. In this case, considering the frequency dependent hybridization function is not mandatory, but still can be done if needed for a particular calculation. As we demonstrate in Section 5.2, ``the D-TRILEX approach can accurately treat the non-local interactions on the basis of the DMFT impurity problem. The latter does not contain the bosonic hybridization function and thus is easier to solve numerically.'' At the same time, performing additional calculations for the susceptibility for the dimer problem presented in Section 5.1 we find that the largest mismatch between the D-TRILEX and the exact results appears only for a very strong non-local interaction V>0.3U. We argue there that ``the inclusion of the bosonic hybridization function in the spirit of EDMFT could improve the result, because in this case some contributions of the non-local interaction would be taken into account in the impurity problem via the bosonic hybridization.''>>

2 - explain better the motivation for the benchmark in the case of the extended Hubbard model. What does this have to do with the multiband generalization of D-TRILEX? What about the multi-site (Nimp > 1)solution for the extended Hubbard model?

<<A: This question has been addressed above and the corresponding changes to the text have been added.>>

3 - discuss more the systematic errors, and possibility for systematic improvements. How does the formalism change if we improve the dual boson starting point? Explain what is the Nimp→∞ limit of the method.

<<A: This question is also addressed above.>>

4 - Discuss more the relationship of D-TRILEX to TRILEX and other methods. Citation of Phys. Rev. B 96, 104504 (2017) in the context of TRILEX and GW+DMFT is due, as it showcases the versatility of those methods. Can D-TRILEX be formulated in Nambu space to treat superconductivity (it is, in a way, a 2-orbital calculation)?

<<A: We thank the Referee for pointing to the missing citation. The corresponding discussion has been added to Conclusions and reads: ``Even though not included in our current computational scheme, one can also address the problem of superconductivity in the framework of the developed approach. In order to find the transition temperature to a superconducting state one can look at the divergence of the lattice susceptibility (72) in the corresponding particle-particle channel. To this aim one can either consider a suitable reference problem related to the desired superconducting order parameter (see, e.g., Refs.PRB 94, 125133 (2016), PRB 100, 024510 (2019), PRB 101, 045119 (2020), npj Quantum Mater. 7, 50 (2022)]) and use the D-TRILEX form for the polarization operator in the particle-particle channel, or to additionally account for the scattering on the transverse momentum- and frequency-dependent bosonic fluctuations in the polarization operator in the particle-particle channel in the case of a single-site reference system [PRB 90, 235132 (2014), PRB 96, 104504 (2017), PRB 99, 041115 (2019)]. Going inside the superconducting phase would require to introduce an anomalous component of the Green's function working in the Nambu space formalism similarly to what has been proposed in the framework of the TRILEX approach [PRB 96, 104504 (2017)].''>>

5 - Discuss more the stability of the Dyson equation in the second paragraph of page 11. Have you encountered problems with diverging bosons? Is there a general problem in reaching low temperatures? Can one suppress ordering instabilities? Can one apply D-TRILEX in ordered (say, AFM) phases?

<<A: Yes, we have issues with the convergence of this equation when we are close to a phase boundary. In the current implementation of the method we cannot go to dynamically broken symmetry phases, so we experience a divergence at the transition point. We added a more detailed discussion of strategies to avoid the divergence in low-dimensional systems, where Mermin-Wagner theorem forbids certain forms of symmetry breaking. We expect to being able to partially cure this issue by considering cluster or a multi-impurity reference problems. We would like to point out that this issue is not restricted to D-TRILEX. Most of the diagrammatic methods exhibit this issue in one way or another.>>

5 - below Eq.2 define →γk,ll′ for the sake of completeness.

<<A: This quantity is already defined few lines after Eq.(2).>>

6 - last paragraph before Section 6, fix ``the complex physics in the considered this model''.

<<A: The misprint has been corrected.>>

Reply to Anonymous Report 2

The authors present an extension of the D-TRILEX numerical method for many body systems introduced in Ref. [91] for a single-band model to the more generic case of multi-orbital systems. After introducing the formalism in Sections 2 and 3 (with the support of the Appendices) and discussing the details of the numerical implementation in Sec. 4, the authors showed applications of the method to several cases, benchmarking it when possible to exact solutions or less approximated methods. Although there are no conceptual advances with respect to the single orbital case, the article reports technical advances due to numerical difficulties for the theoretical treatment of more complex multi-orbital systems.

<<A: We thank the Referee for reviewing our manuscript as well as for the useful comments and suggestions.>>

The benchmark with the Kanamori dimer is particularly interesting. In fact, given the extreme simplicity of it, the dimer has an exact solution that can be used to validate the method. In section 5.1., it is shown how DMFT fails to predict the correct value of the imaginary part of the local self-energy at low frequency while the D-TRILEX gives the correct result. However, the benchmark with the exact solution is not complete. In fact, the authors should show how the non-local components of the self-energy, that are missing in DMFT, compare with the exact solution.

<<A: We thank the Referee for suggesting this addition to our manuscript. Indeed, the self-energy directly carries the information about the local part included in DMFT and the neglected non-local part. We added two panels in Figure 3 (old Figure 2) showing the corresponding results for the self-energy.>>

At the same time, it would be important to check that also the susceptibilities calculated with D-TRILEX compare well with the exact solution. Furthermore, since the method can treat non-local interactions, it would be instructive to show a comparison with the dimer in presence of non-local interactions as well.

<<A: Based on this suggestion, we added the results for the static charge and spin susceptibilities obtained for different values of the local U and the non-local V interactions to the revised manuscript. They are presented in Figure 5. We find that the D-TRILEX susceptibility matches very well with the exact result for small-to-moderate non-local interaction strength and breaks down at large non-local interactions.>>

In section 5.2 the authors apply the D-TRILEX to the extended Hubbard model on a square lattice and compare the results to the ones obtained using DIAGMC@DB in Ref.[82]. While for weak-to-intermediate values of the onsite repulsion U the agreement between the two methods is good, their predictions start to deviate for larger values of the Hubbard U. While the authors report this problem they do not really discuss in too much detail why it is happening and how this could be overcome.

<<A: The reason for the deviation from the exact result at larger values of U is that in this regime the system lies close to a magnetic instability, where magnetic fluctuations become strongly non-linear (see, e.g., Ref. [PRB 102, 224423 (2020)]). This non-linear behavior originates from the mutual interplay between different bosonic modes as well as from an anharmonic fluctuation of the single mode itself. The description of these effects requires to consider much more complex diagrammatic structures that account for vertical (transverse) insertions of momentum- and frequency-dependent bosonic fluctuations, which are present in the DiagMC@DB approach but are not considered in ladder-like dual approximations including the D-TRILEX approach. This point is clarified in the revised manuscript. We want to stress, however, that the method is exact in the large-U regime, as it is an expansion around an impurity that represents the exact solution in that limit. Indeed, in work mentioned by the Referee, it is shown that the method performs better when the interaction strength is larger than the bandwidth (U=10-12).>>

In these works Phys. Rev. B 99, 235106 (2019) and Phys. Rev. B 104, 235128 (2021) efficient ways for calculating vertex functions in the strong coupling limit have been put forward reproducing the asymptotic power law behavior of the Neel Temperature in three dimensions for large values of U. The authors should cite these papers and comment if it is possible to use approximate forms for the vertices to obtain more accurate results at stronger couplings.

<<A: The weaknesses of the D-TRILEX method and the mismatches with the exact results are not related to the calculation of impurity vertices. These vertices are computed exactly using continuous-time quantum Monte Carlo method within the w2dynamics software package. Furthermore, these papers seem not related in any way to our work as they rely on the calculation of the four-point vertex function, which is completely absent in our approach. In D-TRILEX, the four-point vertex functions are eliminated from the theory by using a special path-integral (Hubbard-Stratonovich) transformation that generates a partially bosonized approximation for the four-point vertex function. This approximation was first introduced by some of us in a set of works in a single-orbital case [PRL 121, 037204 (2018), PRB 99, 115124 (2019), PRB 100, 205115 (2019), PRB 103, 245123 (2021)], and in the current work was extended to a multi-orbital case.>>

As a suggestion to improve the readability of the paper, I recommend the authors to add a figure where they show schematically the workflow of the D-TRILEX discussed in section 4.

<<A: The corresponding figure has been added to the text as a new Figure 2.>>

I find the manuscript interesting enough and I recommend its publication after my comments have been properly addressed.

<<A: We thank again the Referee for the useful suggestions. >>