Tracking the Footprints of Spin Fluctuations: A Multi-Method, Multi-Messenger Study of the Two-Dimensional Hubbard Model

The Hubbard model represents the fundamental model for interacting quantum systems and electronic correlations. Using the two-dimensional half-filled Hubbard model at weak coupling as testing grounds, we perform a comparative study of a comprehensive set of state of the art quantum many-body methods. Upon cooling into its insulating antiferromagnetic ground-state, the model hosts a rich sequence of distinct physical regimes with crossovers between a high-temperature incoherent regime, an intermediate temperature metallic regime and a low-temperature insulating regime with a pseudogap created by antiferromagnetic fluctuations. We assess the ability of each method to properly address these physical regimes and crossovers through the computation of several observables probing both quasiparticle properties and magnetic correlations, with two numerically exact methods (diagrammatic and determinantal quantum Monte Carlo) serving as a benchmark. By combining computational results and analytical insights, we elucidate the nature and role of spin fluctuations in each of these regimes and explain, in particular, how quasiparticles can coexist with increasingly long-range antiferromagnetic correlations in the metallic regime. We also critically discuss whether imaginary time methods are able to capture the non-Fermi liquid singularities of this fully nested system.

The Hubbard model represents the fundamental model for interacting quantum systems and electronic correlations. Using the two-dimensional half-filled Hubbard model at weak coupling as testing grounds, we perform a comparative study of a comprehensive set of state of the art quantum manybody methods. Upon cooling into its insulating antiferromagnetic ground-state, the model hosts a rich sequence of distinct physical regimes with crossovers between a high-temperature incoherent regime, an intermediate temperature metallic regime and a low-temperature insulating regime with a pseudogap created by antiferromagnetic fluctuations. We assess the ability of each method to properly address these physical regimes and crossovers through the computation of several observables probing both quasiparticle properties and magnetic correlations, with two numerically exact methods (diagrammatic and determinantal quantum Monte Carlo) serving as a benchmark. By combining computational results and analytical insights, we elucidate the nature and role of spin fluctuations in each of these regimes and explain, in particular, how quasiparticles can coexist with increasingly long-range antiferromagnetic correlations in the metallic regime. We also critically discuss whether imaginary time methods are able to capture the non-Fermi liquid singularities of this fully nested system. The Hubbard model [1][2][3][4] has, for interacting quantum systems, a similar status as the Ising model for classical phase transitions and magnetism. It is the simplest possible model that can be considered, which nonetheless captures essential aspects of the physical phenomena of interest. In relation to materials with strong electronic correlations, the Hubbard model in its simplest form (especially, with a single band) is at best an approximation to reality. However, experimental progress in the field of cold atomic gases in optical lattices now yield rather accurate physical realizations of this simple model in the laboratory [5][6][7].
In contrast to the Ising model however, our current understanding of the Hubbard model is still lacunary. A thorough understanding can be reached in the limiting cases of one dimension [8] and infinite dimensions (infinite lattice connectivity) [9][10][11], thanks to efficient analytical and computational methods available in these limits. In contrast, the two-dimensional case relevant to both cuprate superconductors [12][13][14][15] and cold atomic gases [5][6][7] still holds many open questions, both in relation to its phase diagram as a function of interaction, particle density and temperature, and regarding the nature of excited states as well as response functions and transport properties. There is a broad consensus in the community that progress on these outstanding issues is essential, even if addressed through the deceptively simplelooking Hubbard model.
In recent years, a number of computational methods have been developed to this end. In this context, it is of crucial importance to interrogate these methods regarding their respective ability to address regimes of physical interest. Furthermore, increasing emphasis is being put on establishing definite results with controlled computational methods, which can then serve as benchmarks [16][17][18] for approximate, often more flexible and computationally efficient methods.
Here we focus on a regime of the Hubbard model which is simple at first sight but, as we shall see, deceptively so: small interaction values (often referred to as "weak coupling") and half-filling on the square lattice with nearest-neighbour hopping. The main purpose of this article is to assess the ability of state-of-the-art computational methods to address the finite-temperature physics of the model in this regime. We provide an extensive comparison between basically all methods currently available for this purpose, with two distinct Monte Carlo methods serving as reference benchmarks. Despite the apparent simplicity of this regime, we shall see that all methods face rather severe limitations, especially regarding the lowest temperature that can currently be reached. Our study also interrogates the model through a set of different physical observables, spanning thermodynamic properties, single-particle correlation functions (the Green function and associated self-energy) as well as two-particle correlations (the spin correlation function and correlation length). Because of this wide spectrum of both methods and observables, we have borrowed terminology from the astrophysics community in designating our work as a 'multi-method, multi messenger study' [19].
Despite the deceptive simplicity of the two-dimensional Hubbard model in this parameter regime, the physics is quite rich and non-trivial. As is well established, the ground-state is an antiferromagnetic insulator which can be qualitatively understood using Slater's classic description [20]. However, finite-temperature properties display a rich sequence of interesting crossovers between physically distinct regimes as the system is cooled down towards its antiferromagnetic insulating ground-state. Two key features make these finite-temperature properties and crossovers non-trivial: (i) the fact that, despite perfect nesting, fluctuations destroy antiferromagnetic longrange order at any non-zero temperature -while the correlation length being exponentially large (Mermin-Wagner theorem [21,22]) and (ii) the van Hove singularity present at the 'antinodal' points of the Fermi surface, which further suppresses coherence of single-particle excitations near these points. Furthermore, the perfect nesting of the Fermi surface leads to departure from Fermi liquid behaviour in the metallic regime.
Another important goal of this article is therefore to discuss and characterize these different physical regimes and crossovers in details, with a particular focus on assessing the ability of the different methods to capture their physics properly. By combining computational results and analytical insights, we elucidate in particular the role and nature of spin fluctuations in the different regimes. As we shall see, this reveals some unexpected features of the metallic regime which are likely to have broader implications for materials with electronic correlations. Furthermore, we discuss whether the imaginary time computational methods considered in this article are able to probe the subtle non-Fermi liquid singularities of the metallic state caused by the perfect nesting of the Fermi surface. This rich physics thus makes this regime of the Hubbard model a perfect opportunity for systematic benchmarks, as well as useful testing grounds for future work on more complex regimes of parameters as well as real materials.

B. Overview of the methods assessed
We categorize the algorithms considered in this article into the following groups: • Benchmark methods. We consider two very different Monte Carlo methods. The first one is determinantal quantum Monte Carlo (DQMC, [23]), the other one is the diagrammatic Monte Carlo method (further referred to as DiagMC) [24] (in its recent connected determinant implementation [CDet] for connected one-particle reducible quantities [25] and ΣDDMC for one-particle irreducible quantities [26,27], respectively). The reason for their application is twofold: (i) The methods are controlled and, hence, in the regime where they can be applied and converged, these methods are numerically exact. Therefore they are serving as a benchmark for the other approximate methods considered in this article. (ii) The regimes where the methods actually break down will be assessed -a crucial information in relation to their application to more challenging regimes. For both benchmark methods, we show error bars in the figures (which may, however, be smaller than the respective symbol sizes). All data points in this manuscript obtained by these methods are numerically exact.
• Mean field methods. In Sec. III, we discuss dynamical mean-field theory (DMFT, [28]) as a reference point beyond which spatial fluctuations must be included to properly address the twodimensional model. In that section, we also briefly discuss simple static mean-field theory (MFT).
• Cluster extensions of DMFT. The dynamical cluster approximation (DCA), and cellular DMFT (CDMFT) provide one possible route to systematically include spatial correlations within the DMFT framework, beyond the single-site approximation [29][30][31][32]. Note that cluster-based methods (like CDMFT, DCA and also cluster-TRILEX, see next paragraph) are controlled methods with the control parameter being the size of the cluster. However, in some regimes shown in this paper (low temperatures) these algorithms could not be converged as a function of cluster size (for reasons explained in App. D 5 and D 6).
• Vertex based extensions of DMFT. Another route for including spatial correlations beyond single-site DMFT is relying on higher order Green functions (vertex functions). In the main text we present results from the dynamical vertex approximation (DΓA, ladder version), the triply irreducible local expansion (TRILEX) in various flavors, the dual fermion (DF, ladder version) and the dual boson (DB, single-shot) approach [33].
This list covers the vast majority of currently available computational methods able to address finite temperature properties. One notable exception is the minimally entangled typical thermal state method (METTS) and related approaches, which combine together tensor network representations and stochastic sampling [40][41][42]. A systematic exploration of this method as applied to the Hubbard model is currently being actively pursued by several groups and comparisons with the present methods will have to be performed in future work. The basic principles of each of these methods, useful references for further reading and results from slightly differing implementations of the respective methods and algorithms are summarized in App. D. Throughout the paper we consider the interaction value of U = 2t. Let us stress from the outset that, despite this rather moderate interaction value, each of these methods encounters limitations in their regime of applicability. These limitations stem either from (i) the approximation performed or (ii) algorithmic obstacles.
We find that the lowest reachable temperature for the DiagMC algorithm is T DiagMC min ≈ 0.06t. In this case, reaching lower temperatures is hindered by the difficulty in summing the perturbative series. Interestingly, we find that the limitation of the DQMC algorithm is similar, T DQMC min ≈ 0.06t. In that case, the limitations originate from the exponentially growing correlation length which would require the simulation of prohibitively large systems at lower T . DMFT, in contrast, can be converged to very low temperatures and also at T = 0. Self-consistent methods (e.g. TRILEX) suffer from convergence problems at low-T , whereas calculations involving a 'singleshot' correction beyond DMFT without self-consistency such as DΓA, DB or DF can be performed as long as the correlation length can be accurately resolved (DΓA, DB) or, for DF, as long as the starting point -paramagnetically restricted DMFT -remains reasonably accurate. The finite momentum grid also limits the application of fRG and PA, and, to a lesser degree, TPSC and TPSC+. An intrinsic limitation of TPSC occurs in the renormalized classical regime (see App. D 11) leading to a rather severe overestimate of the onset temperature of the pseudogap. TPSC+ has been proposed to remedy this: in the present paper, the first application of TPSC+ is actually presented, but its applicability has yet to be explored more widely. An obvious limitation of quantum cluster theories are the cluster sizes which they can reach, which have to be compared to the correlation length -a very demanding criterion in the present case, as will be shown later. The Fermi surface of the half-filled system (µ = 0) is diamondshaped (bold black), the black arrows indicate the nesting vectors, mutually connecting Fermi surface points. b) The corresponding (particle-hole symmetric) density of states (DOS) as a function of energy ρ0(ε). ε = 0 corresponds to half-filling. c) The value of the dispersion relation along a high-symmetry path exhibits a plateau around (π, 0), leading to a vanishing Fermi velocity vF.
C. Definition of the model, the role of the van Hove singularity and nesting We consider the single-band Hubbard model defined by the following Hamiltonian: where t is the (nearest-neighbor) hopping amplitude, ij denotes summation over nearest-neighbor lattice sites, σ ∈ {↑, ↓} the electron's spin, U the strength of the (purely local) Coulomb interaction and n iσ = c † iσ c iσ the spin resolved number operator. Throughout the paper all energies are given in units of t = 1. Furthermore we seth = 1 and k B = 1. We consider the case of U = 2 (usually regarded as "weak coupling") at half-filling n = n ↑ + n ↓ = 1, corresponding to a chemical potential of µ = U/2 = 1 and the simple square lattice, resulting in the following dispersion relation for the electrons (lattice constant a = 1): The particular form of the dispersion and the case of half-filling leads to a very peculiar diamond-shaped Fermi surface, already resulting in an interesting behavior with- -dimensional Hubbard model on a simple square lattice in the weak-coupling regime around U = 2t. T N QP and T AN QP denote the onset of coherence at the nodal k = (π/2, π/2) (triangle) and antinodal k = (π, 0) (dot) Brillouin zone point respectively. The onset of the (pseudo-)gap at antinode and node is denoted by T AN * and T N * respectively. Right: Qualitative sketches of imaginary parts of the self-energy on the Matsubara axis (extracted from DΓA calculations) for temperatures corresponding to the colors given in the phase diagram. out interactions present: (i) it exhibits a ("perfect") nesting by the momentum vector Q = (π, π), which connects every Fermi surface point to another respective one (see Fig. 1), leading to an enhanced susceptibility at q = Q and (ii) the momenta around the (stationary) antinodal Fermi surface point k AN = (π, 0) imply a logarithmic divergence in the density of states ρ 0 (ε) at the Fermi level (van Hove singularity, Fig. 1), leading to a larger scattering phase space than at the nodal point k N = (π/2, π/2). Furthermore, because we consider only nearest-neighbor hopping, the diamond-shaped Fermi surface displays perfect nesting by the whole family of wave vectors of the form (q x , ±q x ), with consequences for the nature of the metallic regime.

D. Organization of this article
This article is organized as follows: In Sec. II, we describe the different physical regimes encountered in this model as a function of temperature, using results from our two benchmarks methods (DiagMC and DQMC). In Sec. III, we discuss dynamical mean-field theory, which serves as a starting point of several approximate methods considered in this article. In Sec. IV, we discuss the calculation of single-particle properties using all the different methods introduced above. In Sec. V, we discuss the T -dependence of the double occupancy and its physical significance. In Sec. VI, we discuss two-particle response functions and the T -dependence of the magnetic correlation length. In Sec. VII, we discuss the implications of our computational results for the physics of spin fluctuations in this model. Finally, a discussion and conclusions are provided in Sec. VIII. A number of appendices present more technical points as well as details of the different methods.

II. QUALITATIVE DESCRIPTION: PHYSICAL REGIMES AND CROSSOVERS
Before presenting detailed results from a variety of many-body approaches in Sec. III-VI, in this Section we sketch the overall physical picture which emerges from this study in Fig. 2 (see also [54][55][56]). The left panel indicates in a schematic manner the key crossover scales that delimit different physical regions as a function of temperature T , for a given value of U . Our quantitative study focuses on U = 2 but the qualitative statements made here are expected to apply throughout the weak to intermediate coupling regime (see e.g. [57] for a study of the evolution to higher couplings). As the system is cooled down from high temperature, we observe several regimes with qualitatively different physical properties (a quantitative criterion for the onset of these scales will be given at the end of the section).
At high temperature, thermal fluctuations prevent the formation of long-lived quasiparticles: this regime can be thought as an 'incoherent soup' of fermions above their degeneracy temperature and is depicted as the red shaded area 1 in Fig. 2. Cooling the system progressively extinguishes these thermal fluctuations, leading to increased coherence in the single-particle spectrum and the appearance of long-lived quasiparticles. Here and below, we use the term 'quasiparticle' in a general and somewhat loose sense of a dispersing single-particle excitation with a 'long enough' lifetime. For the specific model at hand, because of perfect nesting, the 'quasiparticles' do not obey Landau Fermi liquid theory however: this is discussed in more details in Sec. VII C 3. At the node,  this quasiparticle coherence scale T N QP corresponds to the temperature at which the thermal de Broglie wavelength v * F /(πT ) along the nodal direction becomes larger than the lattice spacing, with v * F being the effective Fermi velocity renormalized by interactions. The metallic regime is depicted as region 3 (light blue) in Fig. 2.
The crossover scale associated with the passage from region 1 to region 3 is not the same all along the Fermi surface however. Because of the van Hove singularity stemming from the antinodal points of the Fermi surface such as (π, 0) (see Sec. I C), the coherence temper-ature T N QP at the nodal point k N = (π/2, π/2) is higher than the coherence temperature at the antinodal point T AN QP < T N QP . This defines an extended crossover region 2 in which the system is coherent near the nodes but still incoherent near the antinodes (orange shaded area in Fig. 2). Although further lowering the temperature in the metallic regime 3 initially results in freezing out thermal fluctuations and hence in an increase of the quasiparticle lifetime, this does not persist down to the lowest temperature. Indeed, antiferromagnetic correlations develop as T is lowered, with an exponentially growing correlation length, eventually diverging at T = 0 when the groundstate with antiferromagnetic long-range order is reached.
In this low-T regime, long-wavelength antiferromagnetic fluctuations (Slater paramagnons) lead to an enhancement of the quasiparticle scattering rate upon cooling and to the formation of a pseudogap in the singleparticle spectrum, which evolves into a sharp gap in the Slater-like insulator at T = 0 [20]. Once again, the crossover temperature T * corresponding to the suppression of coherence and the opening of the pseudogap is not uniform along the Fermi surface: it is larger at the antinodes where the destruction of coherence occurs first upon cooling, and smaller at the nodes: T N * < T AN * . Hence, in the dark blue shaded area 4 where T N * < T < T AN * , one has a partially (pseudo-)gapped Fermi surface. Eventually, all states of the Fermi surface are suppressed by antiferromagnetic fluctuations for T < T N * resulting in a fully open pseudogap everywhere on the Fermi surface (purple shaded area 5 ). Let us stress again that longrange antiferromagnetic order and a true gap only sets in at T = T Néel = 0 as a consequence of the Mermin-Wagner theorem [21,22].
Since all the temperature scales described above correspond to crossovers, an appropriate criterion must be defined to identify and quantify them. These scales mostly refer to the presence or absence of characteristic spectral features in the single-particle properties, and hence an obvious observable would be the momentum-and energyresolved spectral function: A(k, ω) = − 1 π Im G(k, ω + i0 + ) and corresponding self-energy as a function of real frequency ω (see also App. B and [44]). However, as all the methods considered in the following are formulated on the Matsubara (imaginary) frequency axis, a much more practical criterion can be obtained via the imaginary frequency dependence of the (imaginary part of the) momentum-resolved self-energy Im Σ(k, iω n ), consistently with previous work [54,55,58]. Representative results for this quantity in the five different regimes discussed above are displayed on the right panel of Fig. 2. At high temperature, the thermal fluctuations lead to a divergent behavior of Im Σ(k, iω n ) at low frequencies. T QP can thus be defined as the temperature where this divergent behavior is eased, i.e. when the slope between the first and second Matsubara frequency changes sign and becomes negative [54,55,58]. In the metallic regime (region 3 ), the behavior of Im Σ(k, iω n ) over the lowest Matsubara frequencies can be approximated by a Taylor series Im Σ(k, In a standard Fermi liquid, Z k and γ k are the quasiparticle spectral weight (Z k < 1) and inverse quasiparticle lifetime, respectively (see App. B for details and a critical discussion). Likewise, the onset of insulating/pseudogapped behaviour at T N,AN * is signalled by an additional change of slope, which becomes positive again at lower temperatures. The actual quantitative criterion for extracting the data of Tab. I using the quasiparticle weight Z k is also discussed in App. B.
To illustrate the onset of these different regimes at the one particle level, we display in Fig. 3 results for the imaginary part of the self-energy on the Matsubara axis calculated with the two numerically exact methods used in this paper, namely diagrammatic Monte Carlo (DiagMC, left panels, see App. D 1) and determinantal quantum Monte Carlo (DQMC, right panels, for technical details, including finite size scaling and Trotter extrapolation see App. D 2, for preceding works see also [54,59,60]). Both exact methods shown in Fig. 3 exhibit T QP and T * with a momentum differentiation between antinode (upper panels) and node (lower panels) with T N QP ≈ 0.42 and T AN QP ≈ 0.35 for the onset of the quasiparticle coherence and T AN * ≈ 0.065 and T N * ≈ 0.0625 for the onset of the insulating behavior (more intermediate temperatures have been calculated for the extraction, see Fig. 17 and App. B). DiagMC and DQMC are in agreement within error bars and (for the one-particle quantities) are able to be converged until T ≈ 0.063. With the aim of giving an overview and for further reference, we summarize in Table I the results of each of the different methods investigated below for the crossover temperatures delimiting the five distinct regimes discussed above.

III. THE DYNAMICAL MEAN-FIELD REFERENCE POINT AND THE ROLE OF FLUCTUATIONS
Several methods considered in the following use dynamical mean-field theory (DMFT) as a starting point. Within DMFT, local fluctuations are taken into account (i.e. quantum and thermal fluctuations between the four possible local states on each site), but spatial fluctuations are not. In the two-dimensional model considered here, these fluctuations are strong and DMFT should be viewed as a zeroth-order approximation, which needs to be extended in order to take these fluctuations into account. DMFT is exact in the formal limit of infinite dimensions or infinite lattice connectivity, in which spatial fluctuations do become negligible [9,10].

A. Self-energies and quasiparticle coherence
Because non-local fluctuations are neglected, the selfenergy is approximated within DMFT by a function which is local in real space, i.e. independent of momentum. In Fig. 4, we display the DMFT self-energy for several temperatures, and compare it to the benchmark (DiagMC) result for both the antinodal and nodal points (left and central panel, respectively). It is seen that the DMFT approximation is accurate at very high temperatures, where indeed correlations are very local (as shown below, the magnetic correlation length is only a couple of lattice spacing down to T 0.2). Deviations between the DMFT self-energy and the DiagMC benchmark are already apparent at T 0.33: at this temperature, these deviations are small at the nodal point but already significant at the antinodal point due to the proximity of the van Hove singularity and the resulting momentumdependence of the self-energy. It is interesting to note, however, that the local component of the self-energy (i.e on-site or, equivalently, momentum integrated) is in good agreement with DMFT down to a much lower temperature. As shown in Fig. 5, the local self-energy obtained from DiagMC is still on top of the DMFT result for temperatures as low as T = 0.1, close to the DMFT ordering temperature. This may come as a surprise in view of the fact that the antiferromagnetic correlation length (Sec. VI) is sizeable at this temperature, of order five lattice spacings. In Sec. VII, we provide an explanation to this observation and discuss further its physical significance.
From the results in Fig. 4 (and more temperature points, not shown) one observes that, within DMFT, the onset of quasiparticle coherence associated with the crossover from the high-T incoherent regime into the lower-T metallic regime (regime 1 to regime 2 in Fig. 2) occurs at T DMFT QP 0.45. This is in reasonable agreement with the QP coherence scale at the node from our benchmark calculations T N QP 0.42, but somewhat higher than the antinodal value T AN QP 0.35, again due to the lack of momentum dependence of the self-energy, which is es-sential to account for nodal-antinodal differentiation. We have calculated the quasiparticle coherence temperature within DMFT as a function of coupling U , and display the result in Fig. 6 (red line in upper panel). As expected, this scale decreases as the coupling is increased.
B. Crossover scales: the DMFT viewpoint When used as an approximation for the twodimensional half-filled Hubbard model of interest here, DMFT yields a symmetry breaking phase transition into an insulating phase with antiferromagnetic (AF) longrange order at finite temperatures for any value of U . This is expected because DMFT does not take into account spatial fluctuations which destroy finite-T AF ordering in two dimensions (and, hence, does not satisfy the Mermin-Wagner theorem [21,22]). The Néel temperature T Néel obtained within DMFT is displayed in Fig. 6 as a function of U (blue curve in upper panel), as well as the corresponding staggered magnetization m as a function of temperature (lower panel). As expected, the Néel temperature displays a maximum, which signals the crossover between the weak-coupling regime in which the Néel temperature and magnetization are exponentially small, and the strong coupling regime. For very large U , the charge gap is of order U and the magnetization saturates. In that regime, for T U , spin degrees of freedom are described by an effective Heisenberg model with superexchange J = 4t 2 /U , and the DMFT Néel temperature is proportional to J. For reference, we also display in Fig. 6 the result of the standard static mean-field theory (i.e. the Hartree mean-field for the transition into the spin-density wave phase, denoted below by 'MFT'). MFT considerably overestimates the Néel temperature as well as m for most values of U . This is particularly pronounced at large U where static MFT does not correctly separate spin and charge degrees of freedom so that the ordering temperature incorrectly coincides with the charge gap U m (∼ U at large U ) (for a comparison of MFT and DMFT see [61]).
The key observation is that, taken together, the DMFT quasiparticle coherence scale (red line in Fig. 6) and the DMFT Néel temperature (blue line) can be taken as semiquantitative mean-field estimates for the crossover lines separating the different regimes discussed in Fig. 2 above. The DMFT quasiparticle coherence scale is an estimate for the crossover into the metallic regime 3 in Fig. 2 (with T DMFT QP 0.45 comparable within error bars to the nodal value from our benchmark, while the antinodal one is about 12% smaller). In turn, since in this weak-coupling regime the insulating gap is associated with magnetic quasi-long range order, the DMFT Néel temperature is an estimate for the onset of the insulating behavior -regime 5 in within both MFT and DMFT. In both cases this is achieved by the appropriate Bethe-Salpeter equations for the correlation function: with a static vertex equal to the bare U within MFT (the calculation then reduces to the random phase approximation, RPA), and with a fully frequency-dependent but spatially local vertex within DMFT, see [28]. It is often overlooked that mean-field methods (MFT or DMFT) provide us with a determination of correlation functions and indeed it naively appears as somewhat paradoxical that spatially dependent correlation functions can be obtained from a mean-field ansatz which is local at the one-particle level. This is actually a classic question in statistical mechanics, whose resolution lies in a careful interpretation of linear-response theory [62]. The mean-field correlation lengths are displayed in Fig. 7 in comparison to the DiagMC benchmark. It is apparent that the DMFT correlation length is in excellent agreement with the benchmark down to T 0.1. The figure also clearly illustrates the connection between long-range magnetic correlations and insulating behavior.

C. The insulating regime
We finally comment on the comparison between the DMFT self-energies and the true solution in the low-T insulating regime. Within DMFT, the insulating gap in this weak-coupling regime corresponds to a Slater mechanism associated with AF long-range order [20]. It is associated with spin polarization and a non-zero value of the real part of the self-energy ReΣ ↑ (0) = −ReΣ ↓ (0) = 0 such that the quasiparticle equation ω + µ − ReΣ(ω) = 0 has solutions only for ω outside the gap. The low-frequency behavior of the DMFT selfenergy is non-singular, and in particular it has (using also particle-hole symmetry) a spin-independent linear term ImΣ(iω n ) ∝ (1 − 1/Z)ω n + · · · , corresponding to ReΣ σ (ω) = ReΣ σ (0) + (1 − 1/Z)ω + · · · on the real fre- quency axis, which is similar in structure in the metallic phase and in the AF-ordered insulating phase. This is clearly apparent in Fig. 4. In contrast, because a spin polarization is absent at any non-zero T in the true solution, the self-energy must be much more singular at low frequencies to open the insulating gap, as becomes clear from Fig. 3 and Fig. 4. The precise nature of this singularity is discussed in Sec. VII. At T = 0, however, the model does order and the general structure of the self-energy is expected to be similar to the Slater (and DMFT) one discussed above.
We finally display in Fig. 8 the value of the spectral function extrapolated to zero frequency, as obtained from DMFT as well as from the DiagMC benchmark at the node and antinode. Interestingly, there is rather good agreement between these different methods for this quantity, despite the different low-frequency behavior of the self-energy. In particular, we note that the crossover into the insulating regime is very sharp when seen from this physical observable, which is rather well approximated by the DMFT solution that has long range order in the insulating regime and hence displays a singularity at the Néel temperature. Detailed analysis close to the crossover temperature would of course reveal differences.
Summarizing, we see that DMFT provides a reasonable approximate description of the key crossovers encountered as a function of temperature. As a mean-field theory, it of course mimicks the crossover into an insulating phase with a large AF correlation length as a phase transition into a phase with long-range AF order. Including fluctuations beyond DMFT is therefore especially crucial in two dimensions, in which long-range order exists only at T = 0. This is the purpose of the cluster and diagrammatic (vertex-based) extensions of DMFT discussed in the following sections.

IV. INCLUDING FLUCTUATIONS BEYOND MEAN-FIELD: SINGLE-PARTICLE PROPERTIES
In this section, we show how the inclusion of (non-Gaussian bosonic) spatial fluctuations influences the oneparticle properties (self-energies) beyond the mean field picture. This can be achieved by several methods: cluster extensions (dynamical cluster approximation DCA, cellular DMFT CDMFT, Sec. IV A), and diagrammatic extensions of DMFT (Sec. IV B) as well as with other approaches (Sec. IV C). In this section, the results for the self-energy are presented and compared to the exact (benchmark) methods.
A. Cluster extensions of DMFT: DCA and CDMFT Quantum cluster methods [29] are quite natural extensions of DMFT in which the non-local components of the self-energy are computed by considering a cluster of N c sites self-consistently embedded in the lattice. This provides a controlled sequence of approximations which converge to the exact result when N c → ∞. Depending on the regime of parameters, convergence in this asymptotic limit may or may not be attained in practice however.
There are several flavors available within the broad family of quantum cluster theories, depending on how exactly the interacting cluster is described and how the embedding is performed. First we consider the dynamical cluster approximation (DCA, [63][64][65]), where the Brillouin zone is paved with patches, and the self-energy is approximated as a piece-wise constant function over these patches whose components are calculated from a cluster of N c sites with periodic boundary conditions, hence defining the discrete momenta associated with each patch. For methodological details see [29] and App. D 5. Fig. 9 (left) shows the self-energy on the Matsubara axis of a DCA calculation with N c = 128 momentum patches for several temperatures and the benchmark Di-agMC data as reference. For these calculations symmetry breaking to an AF long-range ordered phase was not allowed. We see that the improvement brought over singlesite DMFT by these calculations is that the momentumdependence of the quasiparticle coherence scale along the Fermi surface and nodal-antinodal differentiation is quantitatively reproduced (regimes 1 -3 ).
However, the (non-magnetic) DCA with N c = 128 is not able to open the pseudogap at low temperatures (shown in Fig. 9 until T = 0.067, i.e. β = 1/T = 16). The reason for this failure is easily understood by noting that the correlation length, previously displayed in Fig. 7, is ξ ≈ 15 lattice sites at the lowest T available (β = 14) and even larger at lower T. A periodic cluster of finite size N c cannot resolve a correlation length larger than ξ max DCA ∼ √ N c /2 which is 5.7a for N c = 128. From this argument, one would expect that a cluster of order 1000 sites would be required to capture the opening of the insulating pseudogap with DCA in the case at hand. As discussed in Appendix D 5, this is beyond practical reach of the algorithm that we use here in the context of DCA.
A different strategy in the family of cluster extensions of DMFT is to consider a real-space embedding and hence a cluster with open boundary conditions, as in the cellular DMFT (CDMFT) approach [32,66]. Fig. 9 also displays data from CDMFT with N c = 64 = 8 × 8. As apparent from these data, the real-space embedding method again captures the momentum differentiation of the coherence scale, but fails in opening the pseudogap for the same reasons as DCA (ξ max CDMFT ≈ 8/2 = 4). From a quantitative perspective, a recently introduced extrapolation scheme (center focused extrapolation, CDMFT+CFE [67], App. D 6) improves the comparison with the benchmark, however for the U = 2 case also fails to open a gap. Hence, the size of the clusters has to be increased much beyond the values of N c considered here in order to reproduce the pseudogap at this small interactions strength (see also [55]).
Summarizing, we conclude that extensions of DMFT based on cluster embedding methods succeed in reproducing the momentum dependence of the quasiparticle coherence scale (nodal-antinodal differentiation) in the metallic regime, but that much larger cluster sizes would be necessary in order to capture the opening of the insulating pseudogap, due to the very large correlation length in this weak coupling regime. Cluster embedding approaches perform much better in regimes with a smaller correlation length, such as larger values of U (i.e. in the strong coupling regime) [68] or (disregarding the signproblem) when doping away from half-filling, as documented by the success of these methods in capturing the strong coupling pseudogap of the doped Hubbard model [31,[69][70][71][72][73][74][75][76][77][78][79][80][81][82][83][84][85][86][87]. Let us also note here that, analogously to the spirit of Fig. 8, allowing for AF order within these cluster methods -and hence mimicking the large correla-tion length regime as an ordered state -would most likely improve the agreement, see also [88].

B. DΓA, TRILEX, and dual fermions/bosons
In view of the limitations of cluster embedding theories for this U = 2 weak coupling regime, we now turn to an alternative way of treating long-ranged correlations more efficiently here, namely diagrammatic extensions of DMFT.
We display in Fig. 10 the results of several such diagrammatic extensions for the frequency-dependence of the self-energy at the antinode at different temperatures (apart from the diagrammatic extensions of DMFT also other approximation methods are considered there, see next subsection). Fig. 11 displays the results for the nodal point. Please note that the lowest temperature displayed is not always the same for the different methods, and it is useful to refer to Table I as a reminder of the important crossover scales. For the sake of comparison, the first panel shows again the data from the DiagMC benchmark.
We observe that all of the diagrammatic extensions of DMFT presented here [ladder-DΓA with a Moriya λcorrection in the spin channel (App. D 7), TRILEX Λ 2 (App. D 8), ladder DF (App. D 9) and single-shot DB (please note that, in the absence of the nonlocal interaction, the fully self-consistent DB theory would coincide with the DF approach when the bosonic hybridization function is discarded, see App. D 10)] are able to correctly reproduce the crossover from the incoherent to the metallic regime. Indeed, all methods display incoherent behavior (region 1 ) at high temperatures, before the onset of quasiparticles becomes visible first for the nodal point (region 2 ) and then, at lower temperatures, for the antinode (region 3 ). The temperatures of this onset T QP , if at all, only slightly deviate from each other and the benchmarks within the numerical accuracy. Larger deviations, both on a qualitative and quantitative level become visible, however, when lowering the temperature into the insulating pseudogap regime associated with growing magnetic correlations. Let us remind the reader that this crossover is signalled by a second change of slope in the self-energies -first at the antinode (region 4 ) and then at the node (region 5 ) -corresponding to a scattering rate that grows upon cooling.
Whereas DΓA and DF correctly reproduce these crossovers into the pseudogap regime, TRILEX does not exhibit these changes of slope, down to the lowest temperatures where we could converge the method. The DF method also succeeds rather quantitatively, both at the node and antinode, while the DB method appears to perform better at the antinode than at the node (but does not open the gap at the accessible temperatures). From a more quantitative point of view, DF and DB slightly underestimate the scattering rate at the node with respect to DiagMC whereas DΓA seems to slightly overestimate the scattering rate at the antinode and simultaneously exhibits a slightly lower T N * than the benchmark. Summarizing, we conclude that among the diagrammatic extensions of DMFT presented here, the DΓA and the DF method appear to be best at capturing the dif-ferent crossover regimes for the self-energy. In terms of the practical ability of performing calculations in this parameter regime, we must point out that all methods suffer from convergence problems when going down to lower and lower temperatures. The reason for these problems varies from method to method. For the benchmark methods: in DiagMC the series cannot be summed at low-T and the DQMC suffers from the exponentially growing correlation length for T < T min ≈ 0.063. In the case of the DΓA (T min ≈ 0.05), lower temperatures can be reached if one is able to converge in the internal momentum grids. The same is true for TRILEX (T min ≈ 0.05), DF (T min ≈ 0.05) and DB (T min ≈ 0.063). Please also note that within DiagMC the lowest reachable temperature is different for node and antinode (1/T AN min = 18 vs. 1/T N min = 16). and antinode. However, there are significant deviations from the benchmark regarding the onset of the insulating pseudogap behavior. TPSC is one of the first methods in which a detailed understanding of the mechanism responsible for the weak-coupling pseudogap was achieved early on (see Refs. [34,89,90] and Sec. VII). As seen from Figs. 10 and 11, the change of slope in the self energies associated with the pseudogap opening is indeed qualitatively captured by TPSC, but the onset temperatures T AN,N * are severely overestimated. As discussed in Sec. VI, this is due to an overestimation of spin fluctuations in this method. A recent improvement of the method, TPSC+ [36], leads to a definite improvement in this respect as shown on the figures. TPSC+ partially feeds back the self-energy into the fluctuation propagators, mimicking frequency-dependent vertex corrections.
The PA appears to eventually capture insulating behavior at the antinode, although at lower temperatures T < 0.05 in comparison to DiagMC, but doesn't open a pseudogap at the node at this temperature.
The fRG calculations are possible only down to a "pseudocritical" temperature scale T 0.07 at which the running coupling constants diverge (see also the discussion in Ref. [50]). Down to this temperature, however, fRG is in qualitative agreement with the benchmark and shows a non-metallic behavior at the antinode (regime 4 ).

V. DOUBLE OCCUPANCY AND POMERANCHUK EFFECT
In view of its physical significance discussed below, we present in this section the temperature-dependence of the double occupancy: It is displayed as a function of temperature in Fig. 12, as obtained from different methods. We see (left panel, DQMC and DiagMC benchmarks) that three regimes are found: at high T (down to about T 1) D(T ) decreases upon cooling and then reaches a minimum, at intermediate temperatures D(T ) actually increases upon cooling and, finally, D(T ) sharply drops when entering the gapped regime. The high T regime is expected and easy to understand: as T is raised, an increasing number of high-energy doubly occupied configurations are thermally populated. This is apparent from the simple ex- pression of D(T ) in the atomic (zero hopping) limit: which approaches D 0 = 0.25 for U/T → 0, i.e. at very high temperatures or, alternatively, in the noninteracting limit. The intermediate regime in which D(T ) increases upon cooling is more interesting -note that this regime includes in particular large parts of the metallic region 3 . As observed early on in Ref. [91], this apparently counter-intuitive non-monotonic behaviour of D(T ) can be understood qualitatively from entropy considerations. Observing that the entropy S is obtained from the free energy F as S = − ∂F ∂T and that D = ∂F ∂U , one obtains the thermodynamic Maxwell relation: Increasing U at fixed temperature in the metallic regime leads to an increase in entropy. Indeed, in this regime, the entropy is linear in T with a slope related to the effective mass which grows as U is increased. Hence, the right-hand side of Eq. (5) is negative and thus D(T ) must increase upon cooling/decrease upon heating in this regime. This is the same phenomenon as the famous 'Pomeranchuk effect' in liquid 3 He (Clausius-Clapeyron equation): upon heating the system, a tendency to increased localization (smaller D) is found, because localization leads to a higher (spin) entropy and is thus favorable thermodynamically (until thermal population of doubly occupied sites kicks in at higher T ). Note that the sign of ∂D/∂T also directly determines the shape of the isentropy lines in the (U, T ) plane, which are defined by S(T i (U ), U ) = const. and hence obey [92,93]: with c = T ∂S ∂T the specific heat per lattice site. Hence cooling can in principle be achieved by increasing the interaction strength adiabatically in the regime where ∂D/∂T < 0, as initially suggested in [92], further discussed in [93], and experimentally realized in cold atomic systems with extended SU (6)-symmetry [94].
Finally, the low T behaviour in which D sharply decreases again upon cooling corresponds to temperatures around regime 4 in which the system behaves as an antiferromagnetc insulator. We note that the total energy of the Hubbard model is given by: The drop in D can be understood from the fact that at small U the crossover into the (Slater) antiferromagnetic correlation regime corresponds to a gain in potential energy [44,54,88]. It was shown in Ref. [88] that at strong coupling, in contrast, D increases when entering the (Heisenberg) antiferromagnetic correlations regime, corresponding to a gain in kinetic energy. We see (left panel of Fig. 12) that DMFT reproduces all three regimes qualitatively, but overestimates the amplitude of the 'Pomeranchuk effect' by approximately a factor of two, as discussed in Ref. [93]. The reason is that the spin entropy in the localized state at high-T is overestimated in DMFT, due to its neglect of spatial correlations. Obviously, in DMFT, the sharp drop at low T corresponds to the phase transition into an ordered phase -a mean-field description of the actual crossover (see above). This can also be seen as an increase of the magnetic local moment, since via S 2 z = 1 − 2D the double occupancy is related to the increase of AF correlations.
All the methods displayed in Fig. 12 qualitatively reproduce the non-monotonous behaviour of D(T ) but the sharp drop at low T is, as expected, only present in the methods which can describe the low T antiferromagnetic insulator regime. It is, for this reason, absent in the TRILEX based methods for the range of temperatures studied (middle panel). Single-site TRILEX also overestimates D significantly, but cluster extensions of TRILEX are closer to the benchmark at intermediate T .
DF similarly overestimates the double occupancy. For details, including differences between different variants of these 'dual' methods, see the discussion about selfconsistency in App. D 10. In contrast, DΓA and DB (single-shot) follow the benchmark quite accurately. Also, the parquet approximation (PA) appears to capture all three described regimes, while TPSC and TPSC+ qualitatively predict the correct physical picture albeit less accurate quantitatively than PA in comparison to the benchmarks. We do not display fRG data here, since too few frequencies in the self-energy could be calculated reliably (see Sec. IV).
We end this section with a technical remark on the actual calculation of D in the investigated methods. In DQMC and DiagMC D can be directly calculated as an equal-time correlation function. For the other methods the calculation of the double occupancy is possible either via the (lattice) Green's function and self-energy, utilizing the equation of motion technique based on the Galitskii-Migdal formula [95]: (used in the DΓA, DF, PA and TPSC/TPSC+, for details on the latter see App. D 11) or, when allowed by the algorithm used for the impurity solver, from a direct computation on the impurity (for DMFT and TRILEX, see also [96]). The DB result was obtained via the local part of the lattice susceptibility as discussed in Ref. [96]. For a comparison to self-consistent DB we refer to App. D 10.
We note in passing that, when applying Eq. (8), treating the high frequency tails is very important to obtain accurate results.

VI. INCLUDING FLUCTUATIONS BEYOND MEAN-FIELD: MAGNETIC CORRELATIONS
In this section, we turn to two-particle correlation functions and compare the different methods for two important observables probing the magnetic correlations: the static antiferromagnetic spin susceptibility χ sp (q = (π, π), iΩ n = 0) and the magnetic correlation length ξ.

A. Antiferromagnetic static spin susceptibility
The spin susceptibility (or spin correlation function) is given by: where we define here and β = 1/T . Note that we omit for simplicity an additional prefactor of 1/2 for the spin operator and that n ↑ = n ↓ = 0.5 at half-filling in the paramagnetic phase (T = 0). The temperature dependence of the static spin susceptibility χ sp (q, iΩ n = 0) is of particular interest for the case of the half-filled Hubbard model on the square lattice, as the perfect nesting (see Fig. 1) leads to a strong enhancement at Q = (π, π). The temperature dependence of χ sp (q = Q, iΩ n = 0) reflects the increasing dominance of antiferromagnetic spin fluctuations upon cooling. One starts at high T with almost independent fluctuating moments and a Curie law χ sp (q = Q, iΩ n = 0) ∝ T −1 (bosonic mean-field behavior). Approaching the T = 0 ground-state with antiferromagnetic long-range order, the range of spin correlations grows and nonlocal spin fluctuations in the paramagnetic phase (antiferromagnetic paramagnons) develop. At low-T , a regime with an exponentially growing correlation length is found (see below). Fig. 13 displays χ sp (q = Q, iΩ n = 0) for various methods as a funtion of (inverse) temperatures on a logarithmic scale. We start our analysis from the left panel. As already described in the analysis of the meanfield picture in Sec. III both MFT (T MFT Néel ≈ 0.21, orange triangles) and DMFT (T DMFT Néel ≈ 0.08, gray squares), due to their mean-field nature, incorrectly predict finite Néel temperatures: the crossover is mimicked as a true (continuous) phase transition (left panel). The thermodynamic transition manifests itself as a divergence of the susceptibility at the corresponding wave vector.
Of course, the Mermin-Wagner theorem prohibits finite temperature ordering in 2D, which is reflected in the data of both benchmark methods DiagMC (black triangles) and DQMC (black circles), whose susceptibilities do not diverge at finite T . Rather, after following the mean-field curves at high temperatures, the benchmark data enter an intermediate regime in which χ(Q) appears to increase approximately exponentially. This regime largely coincides with the metallic regime 3 , eventually crossing over into a second exponential regime at low-T (see also Sec. VI C), which sets in at a temperature close to the DMFT Néel ordering. This low T exponential regime is to be expected, since there, the charge degrees of freedom are frozen out by the gap and the system enters the insulating regime 5 , as indicated by the black dashed line T AN,DiagMC * . In this regime, the effective spin dynamics is expected to be described by a nonlinear sigma model and this exponential growth is typical of the lower critical dimension d = 2 [34,[97][98][99][100]. The first exponential observed in the metallic regime is more surprising and will be discussed in more details below. Let us stress that, due to the reasons already mentioned for single-particle quantities, both benchmark methods are limited in terms of the temperatures they can reach (T min ≈ 0.07). Please note that the minimal reachable temperature differs from the one for the self-energy.
Turning to the diagrammatic extensions of DMFT, one can see that DΓA (left panel), which respects the Mermin-Wagner theorem [43][44][45]101], captures well the different regimes of the benchmark (Curie law at high temperatures and the two exponential regimes). The small quantitative underestimation observed here may potentially be cured by an improved version of the Moriya λ-correction [44] or a more thorough treatment of the asymptotics of the vertex function as a function of frequency [102][103][104]. DF and single-shot DB agree well with the benchmark except for the lowest temperature point where the algorithm can be reliably converged.
In the case of TRILEX (central panel) we present results for different variants of the method, all of which seem to capture a low temperature exponential scaling, however, with different degrees of accuracy: whereas single-site TRILEX (red circles, solid line) and cluster TRILEX with two cluster sites (TRILEX N c = 2, orange squares) largely underestimate the values of the susceptibility, the result is significantly improved by increasing the cluster size to four (TRILEX N c = 4, yellow circles). Remarkably, the TRILEX variant with the electron-boson vertex inserted on both sides (TRILEX Λ 2 , red triangles, see App. D 8) seems to be on top of the cluster N c = 4 results. Investigating how these values eventually converge to the exact results with N c → ∞ is left to future studies.
We finally turn to the other methods (right panel): Both PA and fRG capture quantitatively the high-T Curie regime. The PA appears to systematically underestimate χ sp from 1/T = 10 on, whereas fRG is over-estimating it. Let us just comment here that although the presented fRG scheme does not fulfill the Mermin-Wagner theorem, its recent multiloop extension [47][48][49] and the PA does [51,53] (see also [46]). Both TPSC and TPSC+ are in agreement with the benchmarks at high temperature, but, as already mentioned in the previous section about the self-energy, the spin fluctuation intensity is significantly overestimated in TPSC with respect to the benchmark (right panel). This overestimation is improved again in the TPSC+.
Our results also demonstrate that, remarkably, neither the fact that the Mermin-Wagner theorem is respected by a theory nor that a theory uses self-consistent interacting Green functions, guarantees a quantitatively better result, as can be inferred by a comparison of the TPSC or PA result with the benchmark, respectively.

B. Magnetic correlation length
The spatial range of antiferromagnetic spin correlations can be quantified by the magnetic correlation length ξ, which in practice can be extracted by a fitting procedure using the Ornstein-Zernike (OZ) form of the bosonic propagator [105]: The results presented in this section were obtained using a slightly modified form of the OZ expression, appropriate for a model on a lattice, as described in App. A. Note that we neglect (small) deviations from the OZ form beyond mean-field. Within the OZ form, the AF static susceptibility discussed previously obeys: χ(Q, iΩ n = 0) ∼ Aξ 2 .
Fig. 14 shows the correlation length ξ obtained via such a fitting procedure applied to the susceptibility data from the different methods. One can immediately see that the curves are, to a large extent, qualitatively similar to the susceptibility curves of Fig. 13. Interestingly, the intermediate-T exponential behaviour in the metallic regime is clearly visible in the DiagMC benchmark data. The lowest-T exponential regime in the insulator is hard to reach with DiagMC, but is obtained in DΓA (and, less clearly, in PA and TRILEX) which can be used down to lower temperature than the benchmark methods. We note, however, that sizeable quantitative differences do exist between the different methods at low-T , and hence we conclude that the precise determination of the correlation length in the low-T regime, where it becomes exponentially large, is a challenge for all state of the art computational methods currently available. Summarizing, we observe three successive regimes for magnetic correlations as the temperature is lowered. At high-T ( > ∼ 1/5), a Curie mean-field behaviour is found with a small correlation length ( < ∼ 2), and all methods (except static MFT) are basically in quantitative agreement in this regime. At intermediate temperatures , we find a correlation length that appears to increase exponentially to a good approximation, with values ξ(T = 1/5) 2 reaching ξ(T = 1/13) 10. This is somewhat surprising since this regime is metallic with quasiparticles that become more coherent as temperature is lowered, at least in the nodal region. We discuss further the physical significance of this finding in Sec. VII by comparing to spin-fluctuation theories. All methods reproduce this intermediate exponential regime qualitatively. On the quantitative level, excellent agreement is found between the benchmark Di-agMC and especially DF and DB throughout this regime. We also note very good agreement for the determination of the correlation length with DMFT and TPSC+ down to T 1/10. Most other methods (DΓA, the various variants of TRILEX and the PA) also provide a satisfactory determination of the susceptibility and correlation length in this intermediate regime.
The low-temperature insulating regime T < ∼ 1/12.5 T DMFT Néel can hardly be probed with current state of the art benchmark methods. The DΓA results down to T = 1/20 are consistent with the exponential growth of the correlation length [44,54], expected from the low-energy description of the spin degrees of freedom by a non-linear sigma model once a charge (pseudo)gap opens up. We note significant discrepancies between the different approximate methods available in this regime however. The growth rate of this exponential regime appears to be different (and faster) than the one in the intermediate-T metallic regime, a finding which will have to be confirmed in future work by exact computational methods when they become capable of reaching lower temperatures.

VII. INSIGHTS INTO THE NATURE AND ROLE OF SPIN FLUCTUATIONS
In this section, we provide understanding into the physical nature of the different regimes highlighted above, especially the metallic regime and the (pseudo)gapped insulating regime. We focus in particular on the nature and role of spin fluctuations. We explore whether spin fluctuation theory of the weakcoupling type is able to describe qualitatively the various regimes: this analysis also provides analytical insights into the physics of these different regimes. We focus in particular on the following key questions: (i) What is the physical mechanism for the opening of the pseudogap at low temperature?
(ii) What are the implications of the growing antiferromagnetic correlation length upon cooling for the coherence of quasiparticles in the metallic regime?
(iii) To what extent can this metallic regime be characterized as a Fermi liquid?
A. Weak-coupling spin fluctuation theory In this subsection, we interrogate weak-coupling spinfluctuation theory and ask whether it provides a satisfactory description of the behaviour of the self-energy, at least on a qualitative level, in the different regimes of temperature. This will also provide guidance throughout the rest of this section in identifying which fluctuations contribute most to the self-energy.
In the simplest version of spin fluctuation theory the self-energy can be expressed quite generally as (omitting the Hartree term): In this expression, χ(q, iΩ n ) is the momentum-and frequency-dependent spin susceptibility and G(k, iω n ) is the (lattice) Green function. In the following we always choose G = G 0 to be the non-interacting Green function, see the end of Sec. VII B for comments about the drawbacks of self-consistent schemes.
The coupling constant g in Eq. (12) characterizes the coupling between electrons and the spin collective modes. It is renormalized as compared to its bare value as U increases. In the following, we are more interested in asking whether weak coupling spin fluctuation theory qualitatively captures the different regimes than in quantitative statements. However, when a comparison is attempted, a value of g has to be chosen. For example, in the version of TPSC considered in the present article, g is given by (focusing on the spin channel only): g 2 TPSC = 3 8 U U sp , with U sp = U n ↑ n ↓ / n ↑ n ↓ . The factor 3/8 ensures rotational invariance by properly accounting for the relative contributions of longitudinal and transverse spin fluctuations [106], while the expression of U sp is key to the TPSC scheme, see also App. D 11. In the following we shall, for simplicity, use g 2 = 3U 2 /8 when performing quantitative comparisons. Fig. 15 presents a comparison between the DΓA results for the self-energy at the antinode k = (π, 0) and the spin-fluctuation expression Eq. (12) in which different choices are made for the spin susceptibility χ(q, iΩ n ). We perform the comparison using the susceptibility data from DΓA due to its good agreement with the benchmark, the availability of the q-resolved susceptibility and its ability to enter the highly insulating regime 5 . This comparison is performed for three different temperatures, corresponding to the incoherent regime 1 , the metallic regime 3 and the pseudogapped insulating regime 5 .
We begin the discussion with the result obtained by using the full (solid blue curve, square markers) DΓA susceptibility. We see that, at a qualitative level, the spin-fluctuation approximation succeeds in reproducing the characteristic low-frequency dependence of the selfenergy associated with all three regimes. There is furthermore good quantitative agreement in the incoherent regime 1 . With the chosen value of g, the selfenergy is overestimated in both the metallic and pseudogap regimes. Using g = g TPSC together with the value of the double occupancy calculated in Fig. 12 largely remedies this overestimation in a large part of the metallic regime.
In order to better understand which spin fluctuations dominate the different regimes, we have also applied Eq. (12) by restricting the integration over the momentum transfer q to the vicinity of the antiferromagnetic wave vector Q = (π, π) according to |q i − Q i | < 2ξ −1 , i ∈ {x, y} (dashed blue curves in Fig. 15). We see that in the pseudogap regime 5 the qualitative frequency dependence of the antinodal self-energy is not modified by this restriction (although obviously, as a result of the restriction in the integration, the self-energy is underestimated by an amount which is approximately frequency independent). This indicates that, in this regime, the self-energy is indeed dominated by the antiferromagnetic collective modes associated with the vicinity of Q = (π, π). In contrast, a completely different situation is found in the metallic regime 3 : restricting the momentum sum to the vicinity of Q yields a spin-fluctuation self-energy which has a qualitatively different frequency dependence than the actual one, inconsistent with the properties of a metal (and actually closer to the shape associated with a pseudogap). It also underestimates the overall order of magnitude of the self-energy by more than a factor of four. We thus conclude that the single-particle properties in the metallic regime are not dominated by the spin fluctuations associated with the antiferromagnetic wave vector: taking into account spin fluctuations at all wave vectors is crucial. We will come back in detail to this observation in Sec. VII C. We finally note that the incoherent regime is not affected by the momentum restriction, simply because the correlation length is so small there that our criterion on the transfer momentum actually does not restrict the integration domain significantly.
Finally, we consider whether the spin susceptibility in Eq. (12) can be approximated, as often done in spin fluctuation theories, by an Ornstein-Zernike form emphasizing long-wavelength antiferromagnetic fluctuations: As detailed in App. A, we use a lattice generalization of this expression and perform a fit of the momentum and frequency dependence of the DΓA susceptibility in order to determine the amplitude A, correlation length ξ (Sec. VI) and Landau damping coefficient γ as a function of temperature. The antinodal self-energy obtained by inserting expression (13) into the spin-fluctuation expression (12) is displayed in Fig. 15 (purple data points). This approximation is quantitatively accurate in the high-temperature incoherent regime. It captures qualitatively the characteristic frequency dependence of a pseudogap in the low-temperature regime, as also detailed analytically in the following section, but considerably overestimates its magnitude. In contrast, in the metallic regime, it fails quite severely both qualitatively and quantitatively, indicating that the Ornstein-Zernike approximation fails to capture important spectral weight in χ(q, iΩ n ) at q far from (π, π) which is crucial in the metallic regime as emphasized above and further detailed in Sec. VII C.

B. Analytical insights into the pseudogap insulating regime
We provide here some analytical insights into the low-T pseudogap regime, based on the weak coupling spin fluctuation theory and the Ornstein-Zernike approximation for the spin susceptibility in Eq. (13). As seen above, these approximations are qualitatively reasonable (although not quantitatively accurate) in this regime.
In this 'renormalized classical' regime, only the lowest Matsubara frequency needs to be retained in Eq. (12) leading to [34,85,90]: In this expression, the prefactor A of the spin susceptibility has been combined with g 2 and numerical prefactors to yield a couplingg ∝ Ag 2 with the dimension of an energy. In Eq. (14), only the dominant term has been retained, k F is a point on the Fermi surface away from the van Hove singularity (i.e. from the antinode) so that the Fermi velocity v F is non-vanishing, and ω c designates the important low-energy scale: which explicitly depends on the Fermi velocity (i.e. k F ).
In the last expression, we have postulated an exponential growth of the correlation length, characteristic of the low-T gapped regime, with ρ s as the spin stiffness. Expression (14) has the following behaviour at low and high (imaginary) frequency: These expressions shed light on the dependence of the self-energy on Matsubara frequency reported above in the low-T regime. The lowest Matsubara frequency in this regime is always larger than the tiny low-energy scale ω c : in fact the condition πT ω c defines the range of low temperatures in which the insulating/pseudogap behaviour is observed [34,90]. In this regime, the Matsubara frequency self-energy seemingly displays a downwards divergent behaviour as in Eq. (16). However, if a calculation for low-ω < ∼ ω c was possible, it would actually display a saturation to a finite value as in Eq. (17).
For completeness, we discuss the implications of this expression for the low-frequency spectral function and self-energy on the real frequency axis. Analytically continuing the above expression yields (for k = k AN ): with the following asymptotic forms: |ω| with a slope obviously having a sign opposite to that of a Fermi liquid ('Z > 1'). The corresponding self-energy and spectral function on the Fermi surface A(k, ω) = −Im[ω + µ − ε k − Σ] −1 /π are plotted in Fig. 16. The two prominent peaks at the edge of the pseudogap in the spectral function can be understood by using the asymptotic form of ReΣ for ω > ω c in the quasiparti-cle equation ω − ReΣ(ω) = 0. We see that the poles of the Green function are located at a frequency ω such that ω 2 = 2gT ln ω/ω c . The dominant term on the r.h.s. of this equation is the −2gT ln ω c one, and we finally obtain the locations of the main peaks at: where the corrections are linear in temperature. Hence, the location of these peaks become temperatureindependent in the low-T limit. As T → 0 these two peaks evolve continuously into the gap-edge peaks defining the single-particle excitations across the insulating gap due to long-range antiferromagnetic order. Note that, within mean-field theory the T = 0 gap is given by U m s with m s as the staggered magnetization, and that close to a quantum critical point the scaling m s ∝ √ ρ s indeed holds, establishing consistency with Eq. (22) (in the present case of d = 2 additional logarithmic corrections intervene). Note that the above derivation of the pseudogap behaviour down to T = 0 relies on the exponential growth of the correlation length and hence would not apply if the ground-state did not have long-range order. For an alternative approach based on thermal fluctuations starting from the T = 0 ordered state, see Ref. [107]. The width of the AF insulating peaks [20,99,100] are of order (gT 2 /ρ s ) 1/2 (for a DMFT analysis see also [61]). It is furthermore interesting to understand the pseudogap regime in terms of lengths instead of the energy ω c [34]. The one-particle thermal de Broglie wavelength v F /(πT ) determines the maximum distance over which an electron wave packet remains coherent despite thermal agitation.
Eq. (21) shows that the imaginary part of the self-energy at zero frequency becomes very large when this length is smaller than the correlation length for spin fluctuations, leading to the decrease of spectral weight at zero frequency. For electrons in the pseudogap regime, the spins look ordered over the length where temperature does not destroy its coherence. This argument holds only in two dimensions. In higher dimensions, the phase space associated with the integration over wave vectors in Eq. (12) diminishes the importance of long wavelength fluctuations [34] for the self-energy. As is well known from DMFT, in infinite dimensions, the self-energy is not influenced by long-wavelength spin fluctuations [9,10].
A similar analysis can be performed for Fermi surface points in the vicinity of a van Hove singularity, i.e. when v F = 0. The expansion of the dispersion relation then is stopped only after the second order term, giving for k AN = (π, 0) after integration over q which yields logarithmic divergences. Even if these di-vergences could be cut off, for ω → 0 so that the imaginary part of Σ is now proportional to ξ 2 instead of ξ (see also [34,44,108], [109,110] for second order perturbation theory and [111] in the context of ferromagnetic fluctuations), i.e. way larger than for momenta away from van Hove points.
To conclude this subsection on the pseudogap insulating regime, we recall some findings and comment on the ability of the description of this regime in self-consistent theories. First we note, that if one uses G (with a finite self-energy) instead of G 0 as fermionic propagator in Eq. (12), one obtains Im Σ(k, 0) ∝ T −2 , hence missing the exponential growth of Eqs. (17,25) [44], see also [99].
On more general grounds, we note that the selfconsistency at the level of vertex corrections should be consistent with the one of the Green function. This remark serves as a warning about imposing self-consistency solely over the Green function. When using approximate schemes which truncate perturbation theory or select a specific class of diagrams without taking vertex corrections into account in a consistent manner, selfconsistency may lead to an incorrect physical description -especially in relation to the opening of insulating gaps. For an insightful discussion of these limitations of selfconsistent perturbation theory, see Sec. 6.1 of Ref. [34]. Observations about the inadequacy of self-consistency in truncated perturbation theory approximations have previously been made also in other contexts, such as the GW approximation and DMFT where it is well-known that self-consistent perturbation theory truncated to secondorder fails in yielding Hubbard bands and the opening of a Mott gap [10,[112][113][114][115]. An additional, however, completely different problem can arise when self-consistent theories (e.g. bold diagrammatic Monte Carlo [116,117]) are applied: due to the multivaluedness of the Luttinger-Ward functional [118][119][120][121][122], intimately related [123] to divergences in vertex functions [124][125][126][127][128][129][130][131], the bold approaches may converge to an unphysical branch.
C. Insights into the metallic regime

Temperature dependence and characteristic scales
Finally, in this section, we provide insights into the physics of the metallic regime. We display in Fig. 17 the quantities Z k and γ k as a function of temperature, at the nodal and antinodal points on the Fermi surface, obtained by a fit to the DiagMC self-energy on Matsubara frequencies. For a precise definition and details about how these parameters are extracted, see App. B. For a Fermi liquid, Z k would correspond to the quasiparticle spectral weight and γ k to the inverse of the quasiparticle lifetime. We are dealing however with a perfectly nested system in which Fermi liquid behaviour does not strictly apply, and this interpretation has to be taken with care, as further discussed in Sec. VII C 3. Indeed, we observe that γ k at the nodal point follows an approximately T -linear temperature dependence for T > ∼ 0.1. The antinodal value of γ k is significantly larger, highlighting momentum differentiation.
The metallic regime 3 was conventionally defined above as the regime where the self-energy has a negative slope at low frequency, thus allowing to define a Z k smaller than unity, as indicated on the upper panel of Fig. 17. From the T -dependence of γ k , we see that region 3 defined in this manner actually consists of two subregimes (see also [132,133]): For 0.1 < ∼ T < ∼ T AN QP (regime 3 a covering most of region 3 ), the inverse quasiparticle lifetime γ k decreases upon reducing T , indicating an increase of quasiparticle coherence upon cooling (as characteristic of a metal). When cooling below T 1/10 a minimum of γ k (and Z k ) is found, below which the lifetime decreases upon cooling as a precursor of the pseudogap regime (regime 3 b).
In order to rationalize these findings, we turn again to spin fluctuation theory, Eq. (12). In contrast to the pseu- dogap regime, it is crucial here to perform the full convolution over all Matsubara frequencies. Using the spectral representation of χ(q.iΩ n ), the single-particle scattering rate at zero frequency is obtained from Eq. (12) as: The key point here is to compare the width in frequency of the two factors entering this expression: 1/ sinh βω and Imχ(q, ω). The width of the former is of order temperature T , while the width of the latter is set by Landau damping, as clear from the Ornstein-Zernike form Eq. (13) applicable close to the antiferromagnetic wavevector which reads for real frequencies: This expression peaks at the characteristic spin fluctuation frequency: Hence, for T < ∼ ω sf , only the low-frequency behaviour of Imχ(q, ω) matters in Eq. (26). This analysis shows that, besides the scale ω c = v F /ξ which controls the pseudogap regime, the spin-fluctuation scale ω sf = γξ −2 is important for the physics of the metallic regime (see [134] and references therein). In Fig. 18, we display these two scales as a function of temperature, with ξ(T ) and γ(T ) determined from a fit to the DΓA susceptibility as described above and in App. A. We compare these two scales to temperature itself, and observe that this indeed defines three distinct regimes, which correspond to a very good approximation to the regimes observed in our numerical results, namely (note that ω sf < ω c ): • T < ∼ ω sf (ξ < ∼ γ/T ): metallic regime 3 a. dogap regime 4 , 5 .
• ω sf < ∼ T < ∼ ω c (v F /T > ∼ ξ > ∼ γ/T ). In this narrow intermediate regime (0.06 < ∼ T < ∼ 0.1, 3 b) the pseudogap is not yet fully opened but the scattering rate no longer behaves as in a metal: it displays a minimum and increases at lower T as the pseudogap regime is entered, as displayed in Fig. 17. We note in passing that the upper boundary of this regime coincides with the temperature at which the double occupancy displays a local maximum (see Fig. 12), corresponding to the onset temperature above which local magnetic moments are formed (see also [56,131]). It also has corresponding signatures in other thermodynamic observables [56].

Which fluctuations dominate the metallic regime ?
Expression (26) also allows to understand an important observation made in Sec. VII A, namely that in the metallic regime the spin fluctuations contributing to the self-energy are not dominated by the vicinity of the antiferromagnetic wave vector. Indeed, if the q integration was limited to q Q in Eq. (26), the correlation length would directly control the temperature dependence of ImΣ. Since ξ(T ) was found to increase approximately exponentially in this regime, this would be inconsistent with metallic behaviour. For the latter to hold, it is crucial that all wave vectors contribute.
In order to further document and validate this point, we display in Fig. 19 a comparison between the momentum dependence of the spin susceptibility as obtained in DΓA and that of its Ornstein-Zernike fit (the latter privileges the vicinity of the antiferromagnetic wave vector). We see (left panel) that at T = 1/8 which lies in the metallic regime, the OZ fit captures very poorly the momentum dependence of the susceptibility far from Q. In particular, the OZ fit misses the increased weight that is visible in the DΓA susceptibility along the diagonals. That weight comes from scattering parallel to the antiferromagnetic zone boundaries, as is apparent already in second-order perturbation theory [109]. In contrast, the OZ form is much more accurate at lower temperatures when ξ has grown to larger values and the system is in the pseudogap regime (right panel): indeed, the physics of the pseudogap is dominated by antiferromagnetic fluctuations (Sec. VII B). This is further documented in Fig. 20 which displays the ratio between the susceptibility integrated over a region of the order of 2ξ −1 around Q to its integral over the whole zone. It is seen that this ratio decreases as T is lowered throughout the metallic regime, while it increases again in the pseudogap regime: hence a large fraction of the spectral weight is missed by focusing on the vicinity of Q in the metallic regime. Fig. 20 also shows that the OZ fit overestimates the relative spectral weight associated with the vicinity of the antiferromagnetic wave vector.
Summarizing, we have shown that a comparison of the two key energy scales ω c = v F /ξ and ω sf = γξ −2 to temperature allows for a determination of the different physical regimes as a function of temperature, in excellent quantitative agreement with our numerical results. The momentum-differentiated metallic regime evolves at low-T into a precursor regime of the pseudogap where the quasiparticle lifetime reaches a maximum. The metallic regime is not dominated by the spin fluctuations associated with the antiferromagnetic wave vector: rather all momentum transfers contribute. This explains, importantly, why quasiparticles whose lifetime increases upon cooling can coexist with antiferromagnetic fluctuations whose correlation length strongly increases (approximately exponentially).

Perfect nesting and non-Fermi liquid behaviour
When Imχ is linear in ω at low frequencies, Fermi liquid behaviour follows from Eq. (26) which yields for T < ∼ ω sf : in which 1/Γ q ≡ Imχ(q, ω)/ω| ω=0 . In the present case however, we expect that strict Fermi liquid behaviour does not hold at low frequency, because of the perfect nesting of the Fermi surface. The whole one-dimensional set of wave vectors with q x = ±q y maps a point on the Fermi surface onto another one. The consequences of perfect nesting in two dimensions and above were explored early on in Ref. [135] and subsequent works. They already manifest themselves at the level of second order perturbation theory, as studied in details in the PhD thesis of Lemay [109] and also recently in Ref. [132]. Perfect nesting leads to singularities in Imχ, which in turn yield a singular behaviour of the low-frequency self-energy for ω < ∼ T . At second order in perturbation theory (denoted by Σ 2PT ) this singular behaviour is of the form [109,132]: with k N = (π/2, π/2), i.e. the nodal point. At finite temperatures, this behaviour is replaced by scaling functions of ω/T , e.g. of the form Im Σ 2PT (k, ω) = T φ(ω/T ) at a generic point on the Fermi surface (with possible logarithmic deviations from ω/T scaling, e.g. Im Σ 2PT (k F , ω = 0) ∼ T ln T [109]). More details about second-order perturbation theory and comparison to our computational results are provided in Appendix C.
Treatments beyond lowest order perturbation theory have been attempted using parquet [136] and fRG [137,138] methods. Symbolic approaches to diagrammatic expansions may also provide a route to fast evaluation of self energies and susceptibilities beyond low orders and can be evaluated on the real frequency axis [139][140][141][142]. While a full solution may not be available, we expect non-Fermi liquid behaviour of the form T α φ(ω/T ), where the exponent α may flow as T is reduced towards and the pseudogap regime is approached. The situation is, to a certain extent, similar to one-dimensional systems.
Resolving the low-frequency behaviour for ω < T in the metallic regime is obviously not accessible to methods which provide data at Matsubara frequencies only. It therefore goes beyond the scope of this work and can be viewed as a challenge for future computational work. However, the dependence on temperature of the parameters Z k and γ k obtained by a fit of the self-energy on Matsubara frequencies is expected to reflect the non-Fermi liquid behaviour through its temperature dependence when ω/T scaling applies, as further detailed in Appendices B and C. Indeed, we can see two hints of this in the data displayed on Fig. 17: (i) Z k does not stabilize to a T-independent value before the pseudogap kicks in, but rather slowly decreases throughout the metallic regime and (ii) the temperature dependence of γ k is better described as T -linear rather that quadratic like in a Fermi liquid. Further details about the proper interpretation of the Matsubara fitted parameters Z k and γ k are given in Appendices B and C.

D. A simple approximation
In the previous sections, we have seen that fluctuations from all momentum transfers q, including the momentum independent background of local fluctuations, are playing a crucial role in the metallic regime. Long-wavelength antiferromagnetic fluctuations associated with q Q become increasingly important as temperature is lowered and are responsible for the formation of the pseudogap. This suggests to try a simple approximation in which one relies on DMFT for the local part of the self-energy, while a weak-coupling approximation Eq. (12) is used to account for the non-local part, namely: In this expression, the last term removes the local contribution from the spin-fluctuation formula, so that the local, k-averaged part of the self-energy is entirely accounted for by DMFT. Please note that, in this subsection, we take the self-energy from paramagnetically restricted DMFT. To evaluate Eq. (31), we use the bare value g 2 = 3U 2 /8 of the coupling to spin fluctuations, and the DΓA spin susceptibility. We recall that the latter is simply given by: χ −1 (q, ω) = χ −1 DMFT (q, ω)+λ with λ the Moriya correction which ensures consistency with the Mermin-Wagner theorem and is determined such that T q,n χ(q, iΩ n ) = χ loc , with χ loc being the local susceptibility as obtained in DMFT from the effective impurity model (for details see App. D 7). Expression (31) is actually a simplified version of the DΓA, in which vertex functions are set to unity (see also App. D 7 to recall of the DΓA equation of motion. For an analysis in the context of superconductivity in cuprates see [143]).
The result of this simple 'DMFT+SF' approximation of Eq. (31) is compared in Fig. 21 (top row) to the Di-agMC benchmark and to DΓA for three temperatures. We see that it performs excellently in the incoherent and in the metallic regime. In the pseudogap regime it is qualitatively reasonable, but overestimates the magnitude of the self-energy at low-frequency (in accordance with the observations made above.) This indicates that, for this rather weak value of the coupling U = 2, the vertex correction terms included in DΓA are not essential in the metallic regime, but become increasingly important as temperature is lowered into the pseudogap regime.
The lower panels in Fig. 21 display the local part of the self-energy as obtained from DMFT together with its approximation from weak-coupling spin fluctuation theory [the local term removed in Eq. (31)], in comparison to DiagMC. This clearly demonstrates that a weak-coupling spin fluctuation approximation does not provide a good estimate of the local component, even at this weak value of U . In contrast, the non-perturbative DMFT provides an excellent description of the k-averaged (local) selfenergy, while non-local contributions can be reasonably accounted for by weak coupling spin fluctuation theory at this value of U .

VIII. CONCLUSION AND OUTLOOK
In this article, we have purposefully focused on an apparently 'simple' regime of the two dimensional Hubbard model (simple square lattice, weak coupling, half filling) with three goals in mind: (i) assessing the ability of state-of-the-art computational methods to address the physics of this model at finite temperature, through the computation of a range of one-particle and two-particle observables (ii) providing an extensive assessment of basically all many-body methods currently available for this purpose, by comparing them to two very different Monte Carlo methods (DQMC and DiagMC) serving as benchmarks and, importantly (iii) investigating the rich physics associated with the different regimes and crossovers found in this model, as it evolves upon cooling from a high-temperature incoherent regime into a momentum-differentiated metal and eventually into an insulating regime with a pseudogap, and elucidating the nature and role of spin fluctuations in these different regimes.
It is satisfying to observe that the two benchmark methods considered in this article are successful at computing the properties of this model throughout the metallic regime and are in excellent agreement with each other. However, we also found that they face rather severe limitations when attempting to reach low temperatures, to the extent that only the onset of the pseudogap regime can be reached. Indeed, the pseudogap is associated with an exponentially growing antiferromagnetic correlation length which requires prohibitively large systems to be simulated with DQMC, while DiagMC, which in contrast works directly in the thermodynamic limit, is faced with convergence issues when summing the perturbative series at low-T . We are hopeful that this can be overcome by using recent improvements of diagrammatic Monte Carlo based on the CDet algorithm [25], together with improved schemes to improve the ability to resum the perturbative series.
Many of the methods considered in this article use the dynamical mean field theory (DMFT) as a starting point and treat spatial fluctuations beyond this starting point to various degrees of approximation. Obviously, the chosen regime of parameters puts DMFT very much out of its 'comfort zone' since the physics at low-T is dominated by long-range spatial fluctuations which in two dimensions also prevent long-range order at any non-zero temperature. Nonetheless we have emphasized and documented that, when properly interpreted, DMFT provides a useful initial description of the different regimes. Furthermore, it provides a highly accurate determination of the local (momentum-averaged) self-energy not only at high-T but also through most of the metallic regime, which is remarkable given that the correlation length reaches values as large as ten lattice spacings there.
Because of the very large correlation length at low-T , cluster extensions of DMFT are not the best route to follow in this parameter regime of the Hubbard model. Although they succeed in describing the momentumdifferentiated metallic regime in satisfactory agreement with the benchmarks, huge cluster sizes would be required to properly capture the pseudogap regime. Cluster extensions of DMFT are much better equipped for addressing strong coupling regimes with shorter correla-tion lengths.
In contrast, we found that extensions of DMFT (ladder-DΓA, TRILEX, ladder-DF and single-shot DB), that make use of response and vertex functions in order to take into account the contributions of fluctuations to the self-energy, perform much better here. We assessed their respective degree of accuracy and found that the DΓA and DF methods are particularly successful at describing all physical regimes of interest in satisfactory agreement with the benchmarks, and are in particular able to properly capture the pseudogap at low temperatures. These two methods therefore provide the arguably best available approximation deep in the insulating regime at low temperature, where currently no reliable benchmarks exist.
We have also considered many-body methods that do not make use of DMFT as a starting point, such as TPSC/TPSC+, fRG and the PA. TPSC played an important role early on in elucidating the physical origin of the weak-coupling pseudogap in relation to longwavelength antiferromagnetic correlations. We found that it captures the different regimes qualitatively, but strongly overestimates the onset of the pseudogap temperature and the correlation length itself. TPSC+, a recent extension of the method, leads to significant improvements. In contrast, the PA underestimates the pseudogap onset temperature, especially at the nodal point, whereas it captures the behavior of double occupancy very accurately. As implemented here, fRG (oneloop Katanin) shows a qualitative agreement with the benchmark until its running coupling constants diverge.
On the physics side, perhaps the most intriguing finding of our study is that, in the metallic regime, singleparticle excitations with a lifetime that increases upon cooling appear to happily coexist with an antiferromagnetic correlation length which increases steeply, reaching about ten lattice spacings at the onset of the pseudogap. We have shown that the explanation lies in the fact that spin fluctuations associated with the vicinity of the antiferromagnetic wave-vector are not the dominant contribution to the self-energy in this regime. In contrast, fluctuations associated with all momentum transfers contribute, including short-range local ones. This also explains why the local (momentum-independent component) of the self-energy is so accurately captured by single-site DMFT, even in the presence of strong antiferromagnetic fluctuations which are responsible for the momentum differentiation on the Fermi surface.
At the lowest temperatures, antiferromagnetic fluctuations eventually dominate and destroy the coherent metal, which gives way to the pseudogap regime (with a narrow precursor regime at which the quasiparticle lifetime reaches a maximum). These observations also provide a rationale for the success of methods such as DΓA and the DF, which incorporate fluctuations of all wave vectors beyond DMFT. In the parameter regime of interest, we have proposed and tested a simplified version of such an approach.
We have documented in detail the rich sequence of crossovers occurring in this model upon reducing temperature, from an incoherent regime at high T all the way down to the pseudogap regime at low-T , through the intermediate temperature momentum-differentiated metal, and assessed the ability of the different computational methods to capture these crossovers. By analysing the contribution of spin fluctuations, we showed that these different regimes are delimited by the comparison of temperature itself to the two characteristic energy scales v F /ξ and γξ −2 , or equivalently of the correlation length ξ to the thermal de Broglie wavelength ∼ v F /T and to the length scale associated with Landau damping ∼ γ/T . Finally, we have also emphasized (Sec. VII C 3 and Appendices B and C) the limitations of imaginary time/Matsubara computational methods in probing the delicate low energy non-Fermi liquid singularities of this perfectly nested model in the metallic regime.
Looking forward, the numerical, theoretical and physical insights obtained in this article should find direct applications in computational studies [143][144][145] of materials for which the interplay between local electronic correlations and long-wavelength spin fluctuations play a crucial role, such as Sr 2 RuO 4 [146][147][148] or the iron-based superconductors [149]. On the model level, our findings will be immediately useful for the interpretation of how physical (density, magnetic, pairing) fluctuations on the two-particle level influence one-particle spectral functions [123,[150][151][152][153][154][155]. Close to a magnetic phase transition, the magnetic susceptibility (bosonic propagator) assumes the following Ornstein-Zernike form [105]: with A being a prefactor, γ the Landau damping constant (the dynamical critical exponent has been set to z = 2 here), and ξ the magnetic correlation length. This form can be used to evaluate the spin-fluctuation diagram given in Eq. (12). In the renormalized classical regime 5 , the sum may be restricted to the lowest Matsubara frequency without great loss of precision, which is, however, not possible in the metallic regime 3 . The estimation of the correlation length from the static magnetic susceptibility in this manuscript has been performed by an empirically more robust formula (especially in the limit of small ξ) of the Ornstein-Zernike form, which incorporates the lattice (tight-binding) structure of the problem [68], neglecting a possible anomalous exponent η (cf. η ≈ 0.037 for the three-dimensional Heisenberg model [157]): Re χ DΓA sp (q, iΩn = 0), 1/T=12.5 and reduces to the original form in the long wavelength limit (small arguments of the sine). For the fits shown we fixed q y = π. Fig. 22 shows data for the static antiferromagnetic susceptibility in DΓA (upper row) and the respective Ornstein-Zernike fits for various temperatures. Please note again, as discussed in the main text, that the Ornstein-Zernike fit is most accurate only around the vicinity of Q = (π, π) and may not capture the physics at other q points, important e.g. in the metallic regime (see Sec. VII C). These fits of the static antiferromagnetic susceptibility are used to determine the prefactor A and the magnetic correlation length ξ. For extracting the Landau damping γ, A and ξ have been fixed and γ has been fitted from the frequency dependence of the propagator. The upper panels of Fig. 23 show these fits of DΓA data for three representative temperatures. The lower panels show the obtained parameters A, ξ and γ as a function of temperature.
Appendix B: Low-frequency expansion of the self-energy In a Fermi liquid, the self-energy on the Fermi surface can be expanded at low (imaginary) frequency as a Taylor series: Inserting this expression into Dyson's equation G −1 = G −1 0 − Σ, we see that the spectral function displays a quasiparticle peak with a spectral weight Z k and width γ k (inverse of the quasiparticle lifetime τ k ) given by: Performing also a Taylor expansion for momenta k close to Fermi surface, the renormalization of the quasiparticle effective mass is obtained as: (B4) In practice, Z k and γ k , as displayed e.g. in our discussion of the metallic regime in Fig. 17, are obtained in this article by a fit of the imaginary part of the self-energy Im Σ(k, iω n ) calculated on the discrete set of Matsubara frequencies, by a fourth-order polynomial. This is common practice in analyzing computational results available in imaginary time at finite temperature. The sign of the slope of the self-energy at low Matsubara frequencies was also used to distinguish between 'metallic' behaviour (Z k < 1) and the pseudogap and incoherent regimes in which the positive slope corresponds formally to Z k > 1.
Care should be applied, however, in interpreting the values of Z k and γ k obtained from such Matsubara fits in terms of quasiparticles properties, when the system does not obey Fermi liquid behaviour at low frequencies or low temperatures. We here illustrate and clarify this point by considering, for simplicity, a linear fit of the imaginary part of the self-energy over the first two Matsubara frequencies ω 0 = πT , ω 1 = 3πT , namely: with: For simplicity, we have dropped the momentum dependence in these expressions: it is understood that the analysis is performed at a given value of k. We now use the spectral representation of the self-energy Σ(iω n ) = dω σ(ω) iωn−ω with σ(ω) ≡ − 1 π ImΣ(ω + i0 + ), yielding: −ImΣ(iω n ) = ω n dω σ(ω) ω 2 n + ω 2 . (B7) Substituting this into Eqs. (B6), after simple algebra, we finally obtain: We now analyze how these two quantities behave in the low-temperature limit. To this aim, we assume that ω/T scaling applies for small ω and T (but an arbitrary ratio ω/T ), and that σ(ω, T ) = T α φ(ω/T ) with φ a scaling function such that φ(x 1) ∼ const. and φ(x 1) ∼ |x| α . The case α = 2 corresponds to a Fermi liquid for which the scaling function is known to be φ(x) = A(π 2 + x 2 ). A value of α < 2 corresponds to non-Fermi liquid behaviour. For an analysis of ω/T scaling and non-Fermi liquid behaviour within second order perturbation theory, see Appendix C.
To obtain the low-T behaviour of −ImΣ(i0 + )| linfit , we observe that the ω/T scaling form can be directly substituted in the above expression while preserving the convergence of the integral. This leads to, for T → 0: .
(B8) This expression demonstrates that a fit to the self-energy on Matsubara frequencies is able to correctly capture the non-Fermi liquid temperature dependence of the scattering rate. Note that for this to hold, it is crucial that ω/T scaling indeed applies. This derivation has been performed for simplicity for a linear fit, but extends to a higher order polynomial fit.
The low-T analysis of Z linfit proceeds along slightly different lines. Indeed, direct substitution of the ω/T scaling form into the equation above would lead to an ultraviolet divergent integral when α > 1. In that case, . The insets show a zoom of the lowest frequency data (the interval range of the y-axis is constant for all k-points shown in order to ensure comparability). Third row: results for k = (3π/4, π/4). Lowest row: Ornstein-Zernike parameters (for the bubble) and energy scales for 2PT (obtained from a Matsubara fit) as functions of (inverse) temperature.
the limit T → 0 can directly be taken to yield: (B9) This is indeed the exact T = 0 expression of the quasiparticle weight resulting from the spectral representation. Note that, for σ(ω) ∼ ω α , at low frequency the above integral converges in the infrared, and that in that case there are coherent quasiparticles, with a finite spectral weight and a scattering rate ∼ T α smaller than their energy. In the opposite case α < 1, the ω/T scaling function must be used which yields the low-T behaviour: At low-T , Z linfit vanishes at T 1−α indicating the breakdown of the quasiparticle concept. For α = 1 (the behaviour obtained by second order perturbation theory away from the nodal point in the fully nested case), a very slow logarithmic vanishing of Z is expected. In conclusion, this analysis demonstrates that when ω/T scaling applies together with a simple power-law behaviour, a fit to the self-energy over Matsubara frequencies is able to pickup the correct T -dependence of the scattering rate and quasiparticle weight in both the Fermi liquid and non-Fermi liquid case, in spite of the fact that the frequency dependence of the self-energy for ω < T is inaccessible from Matsubara frequencies.
Appendix C: Self-energy from second order perturbation theory For further reference, and to facilitate the discussion of the consequences of perfect nesting in Sec. VII C 3, here we show the results of second order perturbation theory (2PT) for the self-energy. For the calculations on the Matsubara axis we use Im Σ(k, iω n ) = U 2 T q,iΩn G 0 (k+q, iω n +iΩ n )χ 0 (q, iΩ n ), (C1) with the (non-interacting) bubble and G 0 being the non-interacting Green function. We calculate this expression with N iω = 200 Matsubara frequencies and a linear momentum mesh of N q = N k = 200. Fig. 24 shows the 2PT data (blue crosses) in comparison to the DiagMC benchmark (black circles and dashed line) for the antinode (uppermost row) and node (second row) as well as for the generic point k = (3π/4, π/4) [third row] for various temperatures. The benchmark result is described quantitatively by 2PT in the incoherent regime and fairly well inside the high temperature metallic regime. In particular, 2PT exhibits a clear momentum differentiation between node and antinode. However it fails in the description of the pseudogap, i.e. the regimes 4 and 5 are absent. This is particularly apparent comparing the emerging energy and length scales obtained from Ornstein-Zernike fits (last row, cf. discussion in Sec. VII C and Fig. 18). Using a Fermi-liquid ansatz (see App. B), one can extract a momentum-dependent quasiparticle weight Z k and inverse quasiparticle lifetime γ k from the Matsubara axis (Fig. 25). As already discussed in the main text, these parameters are only meaningful (in the strict sense of describing Landau quasiparticles) if the low-frequency behavior of the self-energy is the "correct" one, i.e. compatible with the one expected by a Fermi liquid (e.g. in exhibiting a negative slope and a long quasiparticle lifetime).
However, as already discussed in Sec. VII C 3, this may be quite hard to judge by looking at data from the Matsubara axis only. Therefore, we additionally computed the self-energy in 2PT on a continuous imaginary frequency grid (solid orange lines in Fig. 24, see insets for a zoom [158]). One indeed observes that, for frequencies way smaller than the ones accessible by a Matsubara grid |ω| πT , the slope of the imaginary part of the selfenergy changes, at least in this fully nested case, to being positive, yielding clear non-Fermi liquid behavior. This behavior is particularly emphasized at the nodal point (second row).
In order to solidify our analysis, we also computed the self-energy in 2PT on the real frequency axis (see also [109,132]). The upper row of Fig. 26 shows the low-frequency behavior of the imaginary part of the selfenergy on the real axis for a temperature 1/T = 16.0 deep inside the metallic regime for the antinode (red), node (blue) and generic point k = (3π/4, π/4) (light-green). One can observe that, despite their different magnitudes, the functional forms of the antinode and the generic point are very similar. The node, however, exhibits a com- pletely different functional behavior, eventually resulting in |Im Σ(k = N, 0)| > |Im Σ(k = (3π/4, π/4), 0)| for very small frequencies. In order to compare to the Matsubara data, we extract the "offset" −Im Σ(k, ω → 0) from both the real frequency data (crosses in the lowest panel of Fig. 25) and the Matsubara fits. We notice that both methods of extraction are in a very good agreement inside the metallic regime, except at the nodal point, where a severe underestimation by the Matsubara data of ∼ 50% is found. Clear non-Fermi liquid behavior can also be inferred from the real part of the self-energy on the real axis (second row in Fig. 26): at low frequencies its slope is reversed for the node and the generic k-point and non-linear also for the antinode. For comparison we also show (with the dotted-dashed line) the respective low frequency result from Fermi-liquid fit from the Matsubara axis, which are oblivious of the features at |ω| πT . Further insight into the functional form can be gained by careful analytic calculations [109]. They show that indeed at frequencies lower than the temperature |ω| πT (and, therefore, smaller than the first Matsubara frequency), the self-energy takes a clear non-Fermi-liquid form for Fermi-surface vectors k F (see also Sec. VII C): and with k N = (π/2, π/2).
Eventually, with our precise numerical data, we can show that ImΣ(k, ω) indeed exhibits ω/T scaling as discussed in the main text, by plotting ImΣ(k, ω)/T over ω/T (lower panel of Fig. 26), which results in a data collapse over a broad temperature (and frequency) range inside the metallic regime.

Appendix D: Brief description and references of the presented methods
As one of the goals of this manuscript is a synopsis of a variety of different modern many-body techniques, in this Appendix, we will give a short overview on these. Of course, the Appendix cannot serve as an exhaustive review, but will contain the basic idea of each method as well as the necessary references for further information and, if necessary, computational details for obtaining the results presented in this paper.

DiagMC
Diagrammatic Monte Carlo (DiagMC), first developed by B. Svistunov and N. Prokofiev [24], is based on the idea of stochastically sampling Feynman diagrams contributing to a perturbative expansion directly in the thermodynamic limit. In a broad sense, the Diagrammatic Monte Carlo approach consists of the choice of a diagrammatic expansion for a given observable, a numerical algorithm computing the corresponding expansion coefficients and, if necessary, the resummation techniques applied to the resulting perturbative series.
The physically most natural [159], skeleton, expansions are potentially dangerous when applied to the twodimensional fermionic Hubbard model as they can converge to wrong results in the vicinity of half-filling or whenever the local magnetic moment develops [118,[121][122][123]131]. In contrast, bare expansions have a well-defined analytic structure and generally have a nonzero convergence radius in the two-dimensional fermionic Hubbard model. Such expansions can be constructed around an arbitrary starting point [151,[160][161][162]. The default choice [163,164] is to expand around the noninteracting Hamiltonian with the chemical potential shifted by a constant (α = U/2 at half-filling) [165], which effectively eliminates all diagrams with Hartreetype insertions, leading to a much improved variance, which is what we adopt here.
The choice of algorithm used to compute the expansion coefficients determines how large expansion orders can be obtained with reasonable error bars. In this publication we employ the currently most efficient class of algorithms, based on the Connected Determinants Diagrammatic Monte Carlo (CDet) introduced for connected quantities by R. Rossi [25] and generalized to the evaluation of one-particle irreducible self-energy diagrams in Refs. [27,55,166]. The particular self-energy implementation used in this publication is ΣDDMC from Ref. [26] as it is very well-suited for the computation of specific k-points directly in momentum-space. In its original formulation, CDet only works with bare diagrammatic expansions (consisting exclusively of interaction vertices and bare propagators G 0 (r, τ )), although a generalization to semi-bold [160] and bold [117,167] schemes has recently been developed [162]. Contrary to previously developed Diagrammatic Monte Carlo algorithms [116,117,163,164,[167][168][169] which evaluate a single Feynman diagram at each Monte Carlo step, CDet computes a factorial number of relevant diagrams permuted over all internal vertices at exponential cost, scaling with the expansion order m as o(3 m ) for connected quantities [25] and o(m 2 3 m ) for one-particle irreducible quantities [27,55,166] (or, alternatively, scaling as o(m 2 2 m ) and o(m 4 2 m ), respectively [170]). This is achieved by grouping all possible diagram topologies of a particular diagram order into determinants and subsequently subtracting all disconnected, and in the case of the self-energy also all one-particle-reducible, diagrams from the sum by using a recursive formula. As a consequence of the exponential scaling, the total error of a convergent series scales polynomially with computational time [171]. The statistical variance of the algorithm only allows for the computation of a finite number of expansion coefficients, typically m ∼ 10−12 (compared to m ∼ 6−7 for previous Diagrammatic Monte Carlo algorithms).
In challenging cases, such as low temperatures and/or high values of physical coupling, resummation techniques are needed to evaluate the series close to, or beyond the, radius of convergence. This constitutes the only source of systematic error in this approach. A protocol for the resummation of extrapolation, which gives control over this systematic error, has been developed in Ref. [26].
At half-filling, every other order is equal to zero for the double occupancy (even) as well as for the self-energy evaluated at every k-point on the non-interacting Fermisurface (odd). A consequence of particle-hole symmetry at half-filling is that a single calculation of the expansion coefficients at a given temperature yields results for all values of U , provided the series is resummable (away from half-filling this is still the case, however the evaluated density n also changes as a function of interaction strength U ). For all presented observables, the limiting factor to the applicability of this approach at high enough interactions and low enough temperatures is the presence of singularities in the vicinity of the positive real axis in the complex plane of interaction strength U , which render the resummation of the series prohibitively hard. Despite the fact that, in this publication, series have only been evaluated within their radius of convergence, resummation by means of Padé and Dlog-Padé approximants [26,172,173] has been performed in order to accelerate the convergence of series whilst controlling the systematic error of the extrapolations. In the halffilled two-dimensional Hubbard model, this approach has previously allowed for numerically exact results that shed light on the metallic to quasi-antiferromagnet crossover as described by the self-energy [55] as well as various other thermodynamic observables [56,174].

DQMC
The determinantal Quantum Monte Carlo (DQMC) algorithm, formulated by Blankenbecler, Scalapino, and Sugar [23], is a controlled method and has been widely applied to finite-temperature simulations of correlated fermion systems. Its basic methodology is to transform two-body interactions into free fermions coupled with auxiliary fields and then sample the fields to compute fermionic observables. To achieve this, a Trotter decomposition (such as e −∆τĤ = e −∆τĤ0/2 e −∆τĤ I e −∆τĤ0/2 + O[(∆τ ) 3 ] is used in this work, whereĤ 0 andĤ 1 are the tight-binding and interacting parts of the model Hamiltonian), and a Hubbard-Stratonovich (HS) transformation (the discrete spin decomposition [175] for on-site Hubbard interaction is used in this work) are applied within the discretization for the inverse temperature as β = M ∆τ . When convergence ∆τ → 0 can be established, the method is numerically exact.
The systematic error from finite ∆τ is controllable and can be eliminated by extrapolating simulations with several different ∆τ values. In Fig. 27, we present representative results of systematic Trotter errors in DQMC simulations. Two systems, L = 24 with βt = 10 and L = 28 with βt = 15, are studied. We concentrate on the self-energy at k = (π/2, π/2) (upper panel) and the spin susceptibility (lower panel). The results for both parameters explicitly show the convergences to the ∆τ = 0 limit for ∆τ t ≤ 0.06 within our statistical uncertainty. This means that our choice of ∆τ t ≤ 0.02 is reliable.
Further details about DQMC algorithm can be found in several review papers [176,177]. We have also implemented our most recent improvements [178,179] of this method in this work. The minus sign problem is absent for the Hubbard model at half filling studied in this work and DQMC method can reach large system sizes. For the calculation of dynamical quantities, the imaginary-time single-particle Green function and correlation functions are measured in DQMC simulation, and then the imaginary-frequency observables are obtained through the Fourier transformation. The Dyson equation was then applied to compute the self-energy. For the U/t = 2 case, convergence to thermodynamic limit of all the physical observables is observed for up to L = 32 for βt ≤ 10. For lower temperatures, we perform finite-size scalings using second-order polynomials with constraints of monotonic behaviors to reach the thermodynamic limit. We typically use from 3 × 10 5 to 10 4 measurement samples for systems with linear system size from L = 20 to L = 48. In DQMC the number of Matsubara frequencies is determined by the number of imaginary time slices M , which in turn depends on the Trotter step, chosen here to be ∆τ = 0.02 as mentioned earlier.

MFT
The calculation of the static mean-field theory (MFT) susceptibility is analytically equivalent to the (nonselfconsistent) random phase approximation (RPA). The idea of the RPA, originally introduced in [180], is building ladders for the susceptibilities with the bare interaction U in the physically relevant channel, utilizing the non-interacting susceptibility Eq. (C2). The fermionic propagator, thus, is the non-interacting Green function G 0 (k, iω n ). The interacting susceptibility is obtained via: We used momentum grid with a maximum linear mesh sizes of N q = 128 and N k = 128 and the number of fermionic Matsubara frequencies being N iω = 1000.

DMFT
The dynamical mean field theory (DMFT) has become one of the standard techniques for tackling strongly correlated systems over the past decades. Its basic idea consists in mapping the full lattice Hamiltonian Eq. (1) onto a self-consistently determined single site Anderson impurity model. This procedure is exact in infinite dimensions, but represents an approximation in finite dimensions due to its neglection of spatial correlations. Hence, the self-energy in DMFT is purely local. For further reading about the details of the algorithm we recommend the seminal papers [9,10] and the review [28]. The selfenergy and magnetization data presented in this paper have been produced using a state of the art continuous time quantum Monte Carlo impurity solver in its interaction expansion (CT-INT [165,181,182]), which is, like the DMFT self-consistency scheme, entirely implemented in the TRIQS framework [183]. For the calculation of the (magnetic) susceptibility and correlation length we used data from both (i) an exact diagonalization (with four bath sites) impurity solver and the implementation of the Bethe-Salpeter equations presented in [184] and (ii) the CT-INT solver and tprf framework of TRIQS [185], carefully crosschecking that they obtain the same results. As the magnetic correlations (and correlation lengths) grow exponentially reaching low temperatures, we used a very dense momentum grid for these data with a maximum linear mesh sizes of N q = 200 and N k = 200 and the number of fermionic Matsubara frequencies being N iω = 160. In the case of the CT-INT impurity solver, we used N cycles = 1.2 × 10 7 Monte Carlo steps.

DCA
The dynamical cluster approximation (DCA) [29,63,64] is an embedding technique wherein the electron selfenergy is obtained from the solution of an impurity model with N c interacting sites coupled to an infinite bath. As such, the DCA is just one particular generalization to N c > 1 of the 'single site' dynamical mean field method that becomes exact as the number of impurity sites N c → ∞ [10,28,29]. The DCA formulation partitions the Brillouin zone into N c equal area tiles and approximates the self-energy in each tile a, Σ a as a piecewise constant function of momentum as where φ a (k) = 1 if k ∈ a and φ a (k) = 0 for k / ∈ a. Similar to DMFT, DCA invokes a self consistency loop where an initial guess for the impurity model parameters produces a set of Σ a which are then used to update the bath. In this work we obtain results for large values of N c in the paramagnetic phase but note that at weakcoupling and low temperatures studied in this work the DCA expansion in cluster size N c does not appear to be in a scaling regime and therefore one needs much larger clusters in order to estimate the infinite system size limit [55,186].
To solve the N c site impurity problem we use the continuous time auxiliary field method [182,187] with submatrix updates [188]. Similar to other methods that employ Hubbard-Stratonovich decoupling, the primary computational hurdle is finding determinants of an m×m matrix where m is referred to as the expansion order. Unlike Hirsch-Fye solvers that employ a fixed ∆τ the expansion order and time steps in CT-AUX are not fixed, an important improvement that removes the necessity for a ∆τ → 0 extrapolation. In the case of DCA, applying CT-AUX to clusters of size N c results in an expansion order m ∝ βU N c [69]. In this work we have presented data primarily for fixed cluster size N c = 128 on a bipartite cluster. The DCA self-consistency is presumed to be paramagnetic and is iterated ≈ 15 → 20times until the deviation between iterations is much less than the statistical error obtained from ≈ 2 × 10 6 samples of each frequency up to n = 1024. Our DCA and CT-AUX codes are based on the ALPS libraries [189,190].

CDMFT
Cellular dynamical mean-field theory (CDMFT) is the real space cluster extension of DMFT [29,73,77,191]. While proposed around the same time, it has up to today been less spread than the complementary dynamical cluster approximation which is formulated in mo-mentum space. The impurity of the auxiliary Anderson model in CDMFT consists usually, contrary to single site DMFT, of superstructures obtained by upfolding the lattices unit-cell. On the square lattice, e.g., these may consist of N × N quadratic patches. Hence, non-local correlations are included on length scales given by the cluster size as intersite single-particle self energies. Even within the aforementioned patch geometries, which typically retain point group symmetry, translational symmetry is broken (for N > 2 even within the cluster). Re-periodization schemes suffer from ambiguity and may even lead to convergence problems when attempted inside the self-consistency loop. The lattice quantities are therefore approximated from the converged CDMFT solution by restoring the translational invariance for the Green function, the self-energy, or its cumulants [32,66,[191][192][193][194][195][196]. This study uses a cumulant scheme presented in App. B of [67].
A recent study [67] has shown that a so-called centerfocused extrapolation of self energies (and even susceptibilities) to infinite cluster sizes (with a linear regression in 1/N ) yields the best results in comparison to numerically exact diagrammatic quantum Monte Carlo calculations. The CDMFT+CFE self-energies presented in this paper have been obtained by considering quadratic patches of up to N = 8 for the extrapolation scheme. The auxiliary impurity models were solved using a state-of-the-art continuous time quantum Monte Carlo impurity solver in its interaction expansion (CT-INT [165,181,182]) implemented in the TRIQS framework [183]. We used 20 CDMFT self-consistency loops (until convergence). Afterwards we performed 30 simulations with different random seeds starting from this previously converged solution. We take the self-energy as the mean of these 30 simulations. The solver parameters were a cycle length of 300, 10000 warmup cycles and 4 × 10 6 cycles.

Diagrammatic extensions of DMFT
A different route for the inclusion of spatial correlations on top of DMFT is taken by its diagrammatic extensions. Their principle is to extract a higher order correlation function from an auxiliary impurity model and consistently calculate the desired observables from it. For a very recent overview of diagrammatic extension of DMFT see [33]. More detailed information for each method is given below.

DΓA
In the dynamical vertex approximation (DΓA) [43,101] the two-particle analog of the self-energy, the fully irreducible two-fermion scattering vertex Λ, in a DMFT spirit, is assumed to be purely local. In order to obtain the corresponding susceptibilities and self-energies, without the knowledge of leading instabilities, the par- quet equations have to be solved self-consistently (see [197,198] for nanoscopics, [199] for π-tons and [53,200] for the two-dimensional Hubbard model). However, if (like in the present case of the Hubbard model), the leading instability is known, the scheme can be significantly facilitated by considering only the Bethe-Salpeter in the associated scattering channel and the Dyson-Schwinger equation, in order to obtain susceptibilities and selfenergies (for successful applications see [54,58,68,144,145,198,[201][202][203][204][205][206]). In the vicinity of a phase transition, in order to restore the proper sum rules for the susceptibilities (and, thus, the asymptotics of the self-energy and therefore the two-particle self-consistency), within the ladder-version of the DΓA, a Moriyaesque so-called λ correction is used [43,44,101,207]. For the results presented in the main text of the current manuscript, this correction is only done for the spin susceptibility (λ sp , see below). For the calculation of the (magnetic) susceptibility and correlation length we used data from an exact diagonalization (with four bath sites) impurity solver and the implementation of the Bethe-Salpeter equations and Dyson-Schwinger equation presented in [184] (again, like in the DMFT case) were crosschecked with CT-INT. We used momentum grid with a maximum linear mesh sizes of N q = 200 and N k = 200 and the total number of fermionic as well as bosonic Matsubara frequencies being N iω = N iΩ = 160. In the perspective of a future improved version of a Moriyaesque λ-correction [43,44], it is instructive to also compare the results obtained with the correction in both particle-hole channels (i.e. charge and spin, λ sp + λ ch , see Eq. (6) (with r = (ch, sp)) to the ones obtained only with a correction in the spin channel (λ sp , Eq (5) in [44]): (D4) Results for both schemes are presented in Fig. 28. In comparison with DiagMC in Fig. 10, the self-energy obtained with λ sp + λ ch would lead to a slightly better agreement of the imaginary part at low frequencies at the antinode in comparison to the benchmark (left and central panel), but a slightly worse agreement for the susceptibility, the correlation length and T * (right panel). It should be noted that the spin susceptibility and correlation length of DΓA with λ sp agrees almost perfectly with a version of TPSC+ coined TPSC+GG down to low T . A more detailed comparison of different TPSC variants is done in App. D 11. Eventually, in order to be able to compare the simple DΓAesque approximation discussed in Sec. VII D, we recall the Dyson-Schwinger equation of motion for the ladder-version of the DΓA (see [43,101,208], here omitting the Hartree term): where χ refers to the physical, Moriya-corrected susceptibility, F denotes the full two-particle vertex from the self-consistently determined Anderson impurity model in DMFT and G := G DMFT and the last line of the equation accounts for double counting corrections of the local part of the self-energy. γ denotes the electron-boson coupling vertex (already integrated over the internal momentum k). For the simple approximation in Sec. VII D, the electron-boson coupling vertex is set to unity for all momenta and frequencies. A precise definition of the components of the equation of motion can be found at Eq. (4.19) of [208].

TRILEX
The triply irreducible local expansion (TRILEX) method is a relatively recent approach to strong correlations, which utilizes the decoupling of the fermionfermion interaction of the Hubbard Hamiltonian into auxiliary bosons [209,210]. For these bosons, the fermion-boson vertex of a self-consistently determined impurity model is extracted and the lattice polarization and self-energy are calculated via Hedin's equations. The beauty of this approach consists in the facts that (i) the fermion-boson vertex is relatively easy to calculate (compared to the four-point vertex) and (ii) the convergence of the method to the exact solution can be achieved via a cluster extension of TRILEX [211] (we used DCA clusters with two and four sites in this paper). However, the interaction is not unique and comes at the price of introducing an additional (Fierz) parameter α, which in this paper was set in such a way that we decouple only the spin channel (i.e. α = 1/3 in Heisenberg spin decoupling [210]). This choice can be motivated by a socalled fluctuation diagnostics analysis of the self-energy in the pseudogap regime of the Hubbard model [84]. A different route to enhance the results of TRILEX can be made by inserting the fermion-boson vertex on both sites of the Hedin equations (an approach which we called TRILEX Λ 2 , see also [212,213] for a similar idea in the dual theories). The self-energies in the main text using this method show results which turn out to be similar to N c = 4 cluster TRILEX. We used momentum grids with a maximum linear mesh sizes of N q = 128 and N k = 128 and the number of fermionic and bosonic Matsubara frequencies were chosen as N iω = 200 and N iΩ = 20, respectively. For the impurity solver (CT-INT [165,181,182]) we used N cycles = 6.9 × 10 7 Monte Carlo steps.

DF
The dual fermion approach [214][215][216] is a diagrammatic extension of the single site dynamical mean field theory (DMFT), motivated by the idea that non-local corrections to DMFT can be captured by a perturbative expansion around a solution of the dynamical mean field impurity problem.
In the dual-fermion formalism [214,217] one replaces the lattice problem with a lattice of coupled Anderson impurity problems resulting in an action of the form and where g w is the momentum independent Green function of the Anderson impurity problem, and ∆ ω is the hybridization function between the impurity and the bath [28,214]. The dual-fermion action now depends on f and f * which are dual operators obtained via a Hubbard-Stratonovich transformation and in this dual space the interaction terms have become local, collected into the function V i [218][219][220]. Correctly representing V i remains problematic due to the complexity of higher-leg vertex functions [221,222] and so we truncate the vertex at the level of 4-leg operators although higher order contributions may be important [33,223,224]. Once the dual self-energy,Σ(k, ω) is obtained, the lattice self-energy is given by We present DF results from the open-source opendf code [217], starting from a single-site (N c = 1) dynamical mean field solution obtained with a continuoustime auxiliary-field (CT-AUX, [69,182,188]) impurity solver. The method treats all local correlations in a nonperturbative manner via DMFT and then perturbatively includes non-local correlations via a restricted set of selfconsistent ladder diagrams for the non-local ('dual') self energy in the charge and spin channels -this is known as the self-consistent ladder dual fermion method [216]. From this method one can also extract two particle spin susceptibilities which we present. However, there is no self-consistency on the two-particle level and as such the status of two particle susceptibilities from DF is known to be approximate and expected to maintain only qualitative correctness [225]. In all cases the input DMFT solutions (Green functions and Vertices) are obtained using a continuous time auxiliary field method (CT-AUX) [182].
DMFT results for the full vertex F ν,ν ,ω are obtained with a frequency truncation in ν, ν , ω. We present results with truncations |ν c | = |ν c | = 96 in fermionic frequencies and |ω c | = 64 in bosonic frequencies. The DF result is not strongly dependent upon this truncation. More important is the known sensitivity of the result to the momentum resolution [226]. As such we employ a 64×64 k-space grid in the DF dual self-consistency. Further, the result for the lattice self energy is extremely sensitive to the value of the dual self energy due to the inversion in equation (C7). We iterate the DF self-consistency until convergence on the scale of 1 × 10 −8 . We also want to refer to a recent study of the doped Hubbard model in

DB
The Dual Boson (DB) approach [229][230][231] is an extension of the DF approach that additionally to the local Coulomb interaction accounts for the effect of nonlocal interactions in different bosonic channels (ς). Within the DB approach, local electronic correlations are considered exactly in the framework of the (extended) dynamical mean-field theory (EDMFT) [232,233]. Nonlocal collective fluctuations are treated diagrammatically beyond the EDMFT level. For this aim, dual boson fields ϕ ς are introduced in addition to dual fermion variables f ( * ) that are already present in the DF approach. The DB theory is derived analytically using a path integral formalism, so many existing EDMFT-based approaches can be obtained as a certain approximation of the DB method [213,231]. Also, the DB theory fulfills the Mermin-Wagner theorem, which allows to avoid unphysical phase transitions in two dimensions [234]. The action of the DB theory is following (D8) Here, the bare fermionG k,ν,σ and bosonW ς q,ω propagators are given by nonlocal parts of EDMFT Green function and renormalized interaction [231], respectively.
The interaction partF[f * , f, ϕ] of the dual action contains all possible exact local fermion-fermion and fermion-boson vertex functions of the impurity problem. Here, as well as in most of DB approximations, we restrict ourselves to the lowest-order (two-particle) interaction terms that are given by the 4-leg fermion-fermion and 3-leg fermion-boson vertex functions. This truncation of the interaction allows to describe collective charge [235,236] and spin [237,238] degrees of freedom in a conserving way using the ladder DB approximation [229][230][231]239].
In the main part of the text only single-shot ladder Dual Boson results are discussed. These calculations are performed on the basis of the converged DMFT solution of the problem, where the bosonic hybridization function is equal to zero. Importantly, in the latter case the DB theory fully coincides with the DF approach if only the local Coulomb interaction is considered. The corresponding local impurity problem is solved using the open source CT-HYB solver [240,241] based on the ALPS libraries [242]. This solution requires N cycles = 8.1 × 10 7 Monte Carlo steps. After that, we calculate the dual self-energy and polarization operator diagrammatically, and perform only the inner self-consistency loop in order to obtain the dressed Green's function and renormalized interaction using the Dyson equation. For this, we use momentum grid with a maximum linear mesh size of N k = 128, and the number of fermionic and bosonic Matsubara frequencies N iω = 128 and N iΩ = 64, respectively. The expression for the lattice self-energy of the DB approach coincides with the one of the DF theory in Eq. (D7). The lattice polarization function can be found using a similar expression [212]: The fully self-consistent DB calculations can be performed as follows. To obtain the fermionic hybridization of the effective impurity problem, we use the outer selfconsistency condition that equates the local part of the lattice Green function and local impurity Green function k G kωσ = g ωσ . Regarding the bosonic hybridization function, there is no clear receipt how this quantity has to be determined. Here, we investigate two different selfconsistency schemes that fix the bosonic hybridization. For the X-self-consistent (Xsc) result the local part of the lattice susceptibility is equated to the corresponding local susceptibility of the impurity problem q X ς Ω = χ ς Ω . The other self-consistency can be imposed on a renormalized (screened) interaction (Wsc) q W ς Ω = w ς Ω . The renormalized interaction W of the lattice problem can be defined as where U ch/sp = ±U/2 [213]. The EDMFT renormalized interaction can be obtained neglecting the dual contribution to polarization operator in Eq. (D10), so that Π ς (q, iΩ n ) = Π ς EDMFT (iΩ n ). The renormalized interaction of the impurity problem can be found as where Y ς (iΩ n ) is the bosonic hybridization function.
Corresponding results are shown in Figs. 29 and 30. We note as well that the comparisons between self-consistent DB and self-consistent DF schemes are in good agreement, but differ from the exact result. As we point out in the main text, the single-shot DB approach correctly reproduces exact DiagMC results at almost all temperatures. Surprisingly, we observe that the Xsc DB calculations strongly deviate from the exact result presented in both Figures. At the same time, we find that the Wsc DB result for the self-energy agrees with DiagMC calculations even better than the single-shot DB one. However, two particle quantities, such as the lattice susceptibility and double occupancy, get worse when the selfconsistency on the renormalized interaction is utilized. This can be explained by the fact that considered bosonic self-consistencies cannot fix all desired single-and twoparticle quantities at the same time. Therefore, the question of a good self-consistency for the bosonic hybridization function remains open.
Other approximations 11. TPSC and TPSC+ a. TPSC The two-particle self-consistent (TPSC) [89,90,243] approach is a nonperturbative approach that is only valid from weak to intermediate interaction strength; hence, it does not describe the Mott transition. Nevertheless, there are a large number of physical phenomena that it allows to study with an accuracy comparable to other numerically exact methods. In particular, in two dimensions, it describes the entry into the renormalized classical regime of antiferromagnetic fluctuations. There are no adjustable parameters as opposed to, for example, self-consistent renormalized spin-fluctuation theory. It satisfies conservation laws for total spin and total charge, the Mermin-Wagner theorem, the Pauli principle in the form n 2 σ = n σ , as well as the local spin and local charge sum rules. The two local sum rules in addition to the ansatz that the renormalized irreducible spin interaction vertex is given by U sp n ↑ n ↓ = U n ↑ n ↓ suffice to obtain the irreducible vertices and then compute the spin and charge susceptibilities from χ sp (q) = χ (1) (q) 1 − 1 2 U sp χ (1) (q) ; (D13a) χ ch (q) = χ (1) (q) 1 + 1 2 U ch χ (1) (q) .
The local spin and charge sum rules are given by Tr χ sp (q) = n−2 n ↑ n ↓ and Tr χ ch (q) = n+2 n ↑ n ↓ −n 2 , respectively, and the trace is over the spin and momentum space. The above ansatz for U sp is inspired from Ref. 244  At the first-level approximation (hence the superscript "(1)") of TPSC, the irreducible particle-hole bubble diagram is obtained from where the Green function G σ (k) includes a constant self-energy Σ (1) σ and renormalized chemical potential that lead to a non-interacting form G 0 (k). The two irreducible vertices suffice to find an improved self-energy [106] that does not assume Migdal's theorem and that takes into account rotational invariance and crossing symmetry The consistency condition between one-particle and two-particle properties (Galitskii-Migdal formula referred in the main text) 1 2 Tr Σ (2) G (1) = U n ↑ n ↓ (D16) serves as a guide for the domain of validity of TPSC.
(Double occupancy obtained from sum rules on spin and charge equals that obtained from the self-energy and the Green function. When 1 2 Tr Σ (2) G (2) starts to deviate from the above, the method starts to fail.) See [34,247] for detailed comparisons with QMC calculations and other approaches and discussions of the physics. Ref. [31] reviews the work related to the pseudogap and superconductivity up to 2005 including detailed comparisons with quantum cluster approaches in the regime of validity that overlaps with TPSC (intermediate coupling). A pedagogical review of results and derivation appears in Ref. [35].
TPSC numerical results in the present paper have been obtained on a 256 × 256 lattice and 2 × 1024 Matsubara frequencies with tails treated analytically. Due to the self-energy iterations, we use a 128 × 128 lattice and 2×512 Matsubara frequencies in TPSC+ calculations, so that the sufficient precision can be obtained using moderate computational resources, as in TPSC.

b. TPSC+
One of the main limitations of TPSC is that even for weak to intermediate interaction strength, it is not valid deep in the renormalized classical regime [243]. For example, for the half-filled Hubbard model with only nearest-neighbor hopping, one finds that U sp tends to zero as temperature goes to zero, contrary to the fact that U sp and the site double occupancy n ↑ n ↓ = (U sp /U ) n ↑ n ↓ must saturate to a finite value as T → 0 due to virtual states.
To remedy this problem, an improved TPSC has recently been proposed [36]. We refer to this extension as TPSC+. It is based on an extension of Kadanoff and Martin's ideas [248] who treated the normal state of superconductors in such a way that it connects smoothly to the superconducting state described by the BCS equation. The key idea is that the pair susceptibility takes an asymmetric form χ(q) = − Tr[G(k)G 0 (q − k)]. It is called the pairing approximation or GG 0 theory in this context and it has been extensively used by Levin's group to study pairing pseudogap and related phenomena [249,250].
We apply that idea to the repulsive Hubbard model and use the asymmetric form in the particle-hole bubble where the self-energy Σ (2) entering G (2) has the same form as Eq. (D15): Contrary to the original Kadanoff-Martin approach, the irreducible vertices are computed from the same sum rules and ansatz but the susceptibilities in Eq. (D13) are obtained from χ (2) instead of from χ (1) . The tilde symbol ( ) indicates that Σ (2) , G (2) , and χ (2) self-consistently depend on each other. Note that the trace also includes the spin. The advantages of this approach are as follows. (a) The generalized Stoner criterion for the phase transition temperature T N becomes identical to the mean-field gap equation in the antiferromagnetic state with the interaction vertex reduced from the bare U to U sp [36].
(b) The Mermin-Wagner theorem and the Pauli prin-ciple are satisfied. Analytical arguments that demonstrate these results proceed in a manner analogous to those in TPSC. (c) One-and two-particle properties are consistent in the sense that 1 2 Tr Σ (2) G (2) = U n ↑ n ↓ is satisfied exactly with double-occupancy equal to that obtained from the local spin and charge sum rules. (d) Analytical arguments analogous to those in TPSC show that there is a pseudogap in two dimensions that is a precursor to the antiferromagnetically ordered state at zero temperature. Because the susceptibility dressed by self-energy remains finite at zero temperature, renormalized U sp remains finite while pushing T N to zero. Furthermore, at zero temperature, the size of the pseudogap ∆ pseudogap (0) becomes equal to the finite magnetic gap ∆(0).
On the down side, the susceptibilities at zero wave vector do not vanish at finite frequency, as conservation laws say they should. However, the values of the zero-wave vector susceptibilities at finite frequency are generally small. Finally, the f -sum rule is also slightly violated.
Heuristically, another explanation of the approach comes from the expected cancellations between vertices and self-energy. Consider the Bethe-Salpeter equation for the four-point function (susceptibilities) where we have used the identity GG −1 = δ(1 − 2). Suppose a quasiparticle picture applies, namely G = ZG 0 where Z is the quasiparticle weight. In some approximation, the Z −1 from the vertices δG −1 δφ cancels Z from one of the Green functions. In some sense, U sp and U ch account for the average momentumspace screening of the interaction vertices while Z −1 accounts for their frequency dependence.
Numerical results for TPSC+ in the present paper have been obtained on a 128 × 128 lattice and 2 × 512 Matsubara frequencies. Fast Fourier transforms with a uniform grid for imaginary time τ are used in the self-consistent calculations of Σ 2PT and χ 2PT .
Eventually, we would like to comment on result from a TPSC+ variety, which utilizes the full Green function G instead of G 0 and is therefore coined TPSC+GG. In Fig. 28 we already saw that there is a very good agreement of the TPSC+GG with DΓA for the spin susceptibility and correlation length. Fig. 31 now also shows that at low T the double occupancy agrees well with the benchmark, although it does not open the pseudogap. n is the number of positive fermionic frequencies that determined the parametrization of the two-particle vertex (see Ref. [258] for the definitions). The rest function contains (4n + 1) × (2n) × (2n), the K2-function and the fermion-boson vertex (4n + 1) × (2n), the K1-function (128n + 1), and the self-energy (8n) frequencies. The fermionic momentum dependence of the vertices and response functions is accounted for by a form factor expansion, where we considered only the local s-wave contribution since at half filling the physics is dominated by antiferromagnetic fluctuations. The remaining momentum dependence of the vertices, the response functions and the self-energy are calculated on (kx × kx) equally spaced momentum patches, while the Green's functions and their summation in the particle-hole and particle-particle excitation are calculated on a (px × px) grid. The grid for the vertices and response functions is refined around q = (π, π) adding (k refine x × k refine x ) patches from the (px × px)-grid. This accounts for 216 additional patches with respect to the ones already included in the (kx × kx) grid.
(see also [254][255][256][257] for the context of the pseudogap). Its starting point is an exact functional flow equation, which yields the gradual evolution from a microscopic model action to the final effective action as a function of a continuously decreasing energy scale. Expanding in powers of the fields one obtains an exact hierarchy of flow equations for the n-particle irreducible vertex functions, which in practical implementations is truncated at the two-particle level. Neglecting the renormalization of three-and higher order particle vertices yields approximate 1-loop flow equations for the self-energy and twoparticle vertex [259]. The underlying approximations are devised for the weak to moderate coupling regime, where forefront algorithmic advancements brought the fRG for interacting fermions on 2D lattices to a quantitatively reliable level [49,258]. The substantial improvement with respect to previous fRG-based computation schemes relies on an efficient parametrization of the two-particle vertex, where we combine the so-called "truncated unity" fRG [260][261][262] using the channel decomposition in conjunction with a form factor expansion for the momentum dependence with the full frequency treatment [263] including the high-frequency asymptotics [102,221]. We here use the Katanin replacement [264] in the flow equation for the two-particle vertex as a first step towards the multiloop extension of the fRG which would allow to sum up all the diagrams of the parquet approximation with their exact weight [48,258,265]. In order to account for the form-factor truncation, the self-energy flow is determined by the scale derivative of the Schwinger-Dyson equation [49,50] replacing the conventional one-loop flow equation (in Fig. 32). We refer to Refs. [49,50,258] for the details of the algorithmic implementation and in particular also for the computation of the susceptibilities from the flow, the technical parameters are reported in Tab. II.

PA
The parquet approach is a diagrammatic scheme first introduced by DeDominicis and Martin in 1964 [38]. It relies on the classification of diagrams contributing to the full 1-particle irreducible 2-particle vertex in terms of their 2-particle reducibility. These diagrams can be either fully 2-particle irreducible or 2-particle reducible in one of three channels [39,266,267]. The method is unbiased with respect to the channels and exact if using the exact fully 2-particle irreducible vertex as an input.
The parquet approximation (PA) applied here, consists in setting this fully 2-particle irreducible vertex to the bare Hubbard U i.e. its lowest order contribution. This is correct only up to fourth order in perturbation theory. The PA has previously been used to study several condensed matter problems including the Hubbard model [199,200,[268][269][270][271][272]. The here presented results have been obtained with the TUPS [53] implementation which relies on the method of Truncated Unities [52,262] and vertex asymptotics [102,273] in order to represent vertex functions as well as a previously developed parallelization scheme in order to perform computations efficiently [269,270]. All results are converged in the number of discrete lattice momenta (N q ) and positive fermionic Matsubara frequencies (N f + ) used as well as in the number of basis functions (N F F ) inherent to the Truncated Unity method [53]. Concretely we use a linear mesh size of N q = 48 and N f + = 40 at β = 16/t and fewer frequencies and momenta for higher temperatures as well as N F F = 9 for all calculations. The results for β = 20/t have been extrapolated to infinite momentumgrid size. The double occupancy was calculated using Eq. (8). Please note, that in the PA the sum rule relating the double occupancy obtained from one-particle and two-particle quantities is fulfilled by construction. The PA also fulfills [51,53] the Mermin-Wagner theorem [21,22].