Reply to the second review of the Referee 1.

The authors have addressed the criticisms raised by the referees substantively. The revised version of the paper includes additional benchmark data which strengthen confidence in the applicability of the proposed D-TRILEX method. There are now additional explanations that make the paper more readable, and the main points more easily accessible. The additional comments are written concerning

• the meaning of Nimp and the distinction between cluster impurity models and lattices with inequivalent sites

• the possibility and the role of including bosonic hybridization functions

• systematic improvements towards exact theory (in terms of better reference problems and more complete diagrammatic extensions)

• prospects for further generalizations (ordered phases, superconductivity)

The authors also further elaborate on the divergence of bosonic propagators and the strategy to stabilize the calculation. The issue is now presented fairly and in a way that is accessible to a non-specialist. The added Fig. 2 nicely summarizes the workflow of the method, which greatly improves the readability of the paper. Concerning my request from the previous review:

5 - below Eq.2 define →γk,ll′ for the sake of completeness. <<A: This quantity is already defined few lines after Eq.(2).>> Here I meant that one should write the Rashba form of SOC with an equation, rather than to refer to literature. Overall, the revised version of the paper is excellent, and I recommend publication in SciPost Physics.

<<A: We thank the Referee for this nice words and recommendation to accept our work for publication. Concerning the definition of the spin-orbit coupling (SOC): In our work the SOC is described by a momentum-dependent quantity γk,ll′, which is a Fourier transform of effective spin-dependent hopping amplitudes between lattice sites. These effective spin-dependent hoppings can be obtained by diagonalizing a single-particle part of the Hamiltonian that contains the local term ~L·S, which describes the coupling of the orbital angular momentum L to the spin polarization of electrons S=c†σc (for details see, e.g., Ref. [120] – PRB 52, 10239 (1995)). Since the resulting term γk,ll′ is non-local, this form of the SOC is usually referred to as the Rashba form. The precise expression for γk,ll′ depends on the particular system of interest, and for this reason we cannot write down an explicit formula which is valid in general for all the systems that can be studied using our method. To clarify this point, we modified the corresponding sentence in the text of the manuscript. Now it reads: “The non-diagonal contribution in spin space γk,ll′ describes the effect of the external magnetic field and of the spin-orbit coupling (SOC), that is usually expressed in the Rashba form [119]. The latter corresponds to a Fourier transform of effective spin-dependent hopping amplitudes [120].” >>

Reply to the second review of the Referee 2.

The authors added important information that was missing in the previous version showing further limitations of the method, in particular when a non-local interaction is added. However, I am not satisfied with the benchmark of the self-energy results in the case of the dimer. The authors have chosen a set of parameters where the non-local component of the self energy is way smaller than the local one. For this set of parameters, DMFT is already a good approximation for the dimer and the benchmark is useless in this case. The authors should use the same parameters that yields a substantial deviation between ED and DMFT results as shown in Fig.3 when U = t.

<<A: The corresponding results for the non-local self-energy have been added to Fig. 4. These new results show that the non-local part of the D-TRILEX self-energy also follows the trend of the exact result and it is by far the largest contribution to the self-energy for this kind of systems. >>

Also, I am not particularly impressed by the explanation of why this method should work at strong coupling. First, numerical results point to the opposite direction.

<<A: We do not understand what numerical results the Referee is referring to. As we wrote in our reply to a similar question of the Referee in the previous review round: “We want to stress, however, that the [D-TRILEX] 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 [PRB 103, 245123 (2021)], it is shown that the method performs better when the interaction strength is larger than the bandwidth (U=10-12).” As a matter of fact, Fig. 3 of PRB 103, 245123 (2021) clearly demonstrates that the accuracy of D-TRILEX method is worst in the regime between U=6 and U=8 and improves at the interaction strengths starting from U=8. >>

Second, even if this is an expansion around the strong coupling limit, it really does not imply that the method would yield accurate results at the two-particle level at strong coupling. In fact, the papers I have suggested to be referenced address this topic thoroughly: DIFFERENT strong coupling expansions can yield DIFFERENT values of the Neel critical temperature in 3D and sometimes completely destroy antiferromagnetism (AFM). The right approximation schemes are the one that stabilize AFM starting from the atomic limit (where AFM is absent) and yield the right critical temperature. A natural question that arises at this point, would D-Trilex yield the correct value of T-Neel at strong coupling in 3D? Is out there a reference showing this?

<<A: It is not surprising that different approximations lead to different results for the Néel critical temperature (or for any other quantity of interest) in some intermediate regime between the weak-coupling and the atomic limit. Therefore, we do not fully understand what the Referee means with “the right Néel temperature” and “right approximation schemes”. If we use DMFT as a reference point, we are guaranteed that our expansion will work in some vicinity of the atomic limit. The same holds true for the other choices of the interacting reference problem that are exact in the atomic limit. As any approximation can in principle give a different value for the Néel temperature, we think that comparing D-TRILEX results for the susceptibility with non-exact results presented in the two works [PRB 99, 235106 (2019) and PRB 104, 235128 (2021)] mentioned by the Referee does not add any value to our work. Indeed, in these two papers the Néel temperatures are obtained by calculating the DMFT-like susceptibilities. However, in the 3D case considered there, the DMFT susceptibility does not correspond by any means to the exact susceptibility of the problem. Additionally, we would like to stress again that no approximation for the three-point vertex functions is utilized in D-TRILEX, contrarily to what is done in the works [PRB 99, 235106 (2019) and PRB 104, 235128 (2021)].

That said, we would like to provide some arguments on why D-TRILEX provides reasonably accurate results also for two-particle quantities. First of all, in our approach, single- and two- particle quantities are obtained self-consistently from the same functional introduced in Eq. (75). As a consequence, one can expect the accuracy of the D-TRILEX method to be comparable between the single- and two-particle levels. Secondly, as explicitly demonstrated in Ref. [PRB 103, 245123 (2021)], the D‑TRILEX diagrammatic expansion effectively takes into account only longitudinal fluctuations in the two-particle Bethe-Salpeter equation (BSE) considered in dual theories and also when calculating the DMFT susceptibility (see Fig. 1 in Ref. [PRB 103, 245123 (2021)]). These longitudinal fluctuations represent the leading contribution to BSE as we explicitly demonstrated in Ref. [PRB 103, 245123 (2021)] by means of exact diagrammatic Monte Carlo calculations. Therefore, D-TRILEX should be in a good agreement with the other dual theories and the DMFT result for the susceptibility when non-longitudinal fluctuations can be neglected.

Based on the same argument, we expect our method to give accurate results for the Néel temperature for the single-band Hubbard model on a cubic lattice. In order to support this claim, we performed calculations for the Néel temperature for the half-filled single-band Hubbard model on a cubic lattice at t=1 and U=8. We found a Néel temperature TNéel = 0.33, which is in very good agreement with exact results obtained using Determinant Diagrammatic Monte Carlo (DDMC) and Quantum Monte carlo (QMC) approaches, as well as using the ladder dual fermion method. All these results are nicely collected in a recent work on DiagMC [arXiv:2112.15209 (2021)]. A more thorough analysis of the Néel transition in a cubic lattice will be addressed in a separate work and does not fit within our discussion of the method beyond the single-band Hubbard model. >>

Also, even if D-Trilex does not use the four-point vertex function, DB technique does, therefore having reliable approximate schemes to calculate the four-point vertex function would reduce the computational complexity for numerical calculation of strongly correlated systems in regime where D-Trilex does not work. I think these are subtle and important aspects that the authors missed and could leave for future work, but I truly believe that mentioning them in the discussion section and adding the suggested references would enrich their work.

<<A: The current work does not address the DB theory, except as an intermediate step to derive the D-TRILEX method. As a consequence, we believe that the discussion of approximations for the calculation of the four-point vertex is not strictly related with the content of our work. Going beyond the simple D-TRILEX diagrammatics would indeed necessarily involve dealing with the four-point vertex functions, which will unavoidably ruin the simplicity of the computational scheme proposed in our work. However, a detailed discussion of the multi-band DB theory and/or much more complex extensions of D-TRILEX theory is beyond the scope of this work. >>

1. We clarified the meaning of the Rashba form of the SOC in the introduction.
2. We introduced two new panels in Fig.4, showing the local and non-local self-energy as a function of interaction and frequency.
3. We added a description of the new content of Fig.4 in Sec. 5.1.