Multiloop functional renormalization group for the two-dimensional Hubbard model: Loop convergence of the response functions

(1) Manuscript develops and discusses a very advanced many-body method
(2) This method is applied to a complicated model (the 2D Hubbard model)
(3) Many convergence and benchmark checks are performed
(4) The developed method has much potential for future applications

(1) Some technical aspects could be better illustrated (see below)
(2) The violation of the Mermin-Wagner theorem could be better discussed
(3) Possible qualitatively new results for the 2D Hubbard model could be better highlighted (see below)

This manuscript is about the development and application of a fRG algorithm for two-dimensional Hubbard models, improving previous implementations in various different aspects. Firstly, the full frequency and momentum dependence of the vertex functions are explicitly accounted for and, secondly, the feedback of the self energy into the two-particle vertex flow is included. Furthermore, the multi-loop fRG scheme is implemented which successively adds higher loop orders to the diagrammatic approximation and allows to systematically improve the scheme. The implementation of the multi-loop fRG also represents the main advancement compared to previous fRG studies of two-dimensional Hubbard models. Another focus lies on response functions which are computed in two different ways, via a post-processing of the vertex functions and via explicitly solving the flow equations for the susceptibilities. In the latter case, the authors show how these flow equations can be generalized to a multi-loop structure. As a test system, the authors apply their technique to the two dimensional Hubbard model on the square lattice at half filling. Extensive convergence and benchmark checks concerning the number of fermionic frequencies and momentum patching points, the impact of the self-energy insertion, and the number of loops are performed. As a main result, the authors find a rapid convergence of the susceptibility when increasing the number of loops. This convergence also applies to the aforementioned two methods of computing the susceptibility.

Overall, I find it quite impressive how this work improves the fRG on so many different fronts. Even though the paper has a clear methodological focus I am also impressed that the authors manage to apply their technique to a quite non-trivial system such as the two dimensional Hubbard model. Indeed, what they call their "test model" is already a very complex system. Finally, the comprehensive and careful convergence checks add a significant amount of trust in the validity of their analysis. I see significant potential of this method for future applications to other (possibly more complicated) Hubbard models. For these reasons I recommend the publication of this work after the authors have considered the points/questions listed below.

(1) The Section 2.2 about the derivation of the multi-loop fRG equations for the susceptibility is rather technical and I had problems understanding how this multi-loop extension works. The section would profit much from a diagrammatic illustration of the scheme. In Fig. 1 an illustration is already given, but I didn't find it very enlightening.

(2) As shown in Fig. 6 the real part of the self energy vanishes at two points at the Fermi surface; hence at these points the Fermi surface remains unchanged by the self energy feedback. Does this apply to the whole Fermi surface and for higher loop orders as well? In other words, does the perfect nesting effect always remain intact?

(3) The authors argue that in the limit of large loop numbers, the Mermin-Wagner theorem should be fulfilled. On the other hand, their numerical results in Fig. 10 still show a significant violation of Mermin-Wagner even at l=8. The authors provide arguments why it is very challenging to suppress the pseudo-critical temperature in the absence of extremely long-range fluctuations. However, given the fact that the fRG is often criticized because of the violation of Mermin-Wagner, it would be desirable to see a more detailed discussion here. What I mean is, if a method that is known to fulfill Mermin-Wagner still finds a finite ordering temperature, then this temperature should have a physical meaning. For example, it should be related to the extent of correlations taken into account and it should maybe even be possible to calculate/estimate it without solving the fRG equations. Can the authors comment on this?

(4) Besides quantitative changes in the results when applying their extended fRG scheme to the two-dimensional Hubbard model, do the authors also find qualitatively new properties of this system? If yes, I think it should be better highlighted.

(5) It would be very interesting to know about the numerical efforts of these calculations. How long did they take and on how many CPUs did they run? Do these calculations require massive parallelization and supercomputing facilities?

(6) The numbering of the subsections in Sec. 4 is very strange as it starts with a zeroth subsection.