Real-space cluster dynamical mean-field theory: Center focused extrapolation on the one- and two particle level

We revisit the cellular dynamical mean-field theory (CDMFT) for the single band Hubbard model on the square lattice at half filling, reaching real-space cluster sizes of up to 9 x 9 sites. Using benchmarks against direct lattice diagrammatic Monte Carlo at high temperature, we show that the self-energy obtained from a cluster center focused extrapolation converges faster with the cluster size than the periodization schemes previously introduced in the literature. The same benchmark also shows that the cluster spin susceptibility can be extrapolated to the exact result at large cluster size, even though its spatial extension is larger than the cluster size.


I. INTRODUCTION
Even after decades of intense research, the single band Hubbard model in finite dimensions larger than one remains an unsolved cornerstone paradigm in theoretical solid state physics. As a non-perturbative technique dynamical mean-field theory (DMFT) 1-3 led to a significant leap forward as it provided an approximation which is nowadays widely used to treat material realistic lattice models 4 . In this approximation, the lattice problem is mapped to a self-consistent auxiliary quantum impurity model, leading to a local approximation of the self-energy.
Single site DMFT has led to many successes, e.g. a description of the Mott-Hubbard metal-to-insulator transition (MIT) from a weak coupling metallic over to a strongly correlated Fermi liquid and eventually Mott insulating state 3,5 within a single theoretical framework. The locality of the self-energy is however a major limitation for some of the most interesting problems as for example cuprate high-T c superconductivity [6][7][8][9] where the pseudo-gap phase is characterized by a strong nodeantinode differentiation. Note, however, that in materialrealistic multi-orbital models even local DMFT selfenergies can lead to k dependent effects in the spectral function A(k, ω) due to orbital mixing as, e.g., shown in 10 . Quantum critical behavior [11][12][13] in low dimensions and/or at low temperatures is another example where the physics is dominated by non-local correlations beyond the mean-field description of spatial fluctuations within DMFT.
In this paper, we focus on the real space CDMFT method 19 . While large cluster sizes were studied in DCA up to convergence and compared to exact methods like diagrammatic Monte Carlo in some parameter regimes, e.g. 57,58 , large CDMFT clusters have not received the same attention. The central issue is that CDMFT breaks translational invariance by definition for any finite cluster. One therefore needs to use reperiodization schemes to restore the translation symmetry of the lattice self-energy when approximating it from the cluster self-energy 19 . Several ways to do this have been discussed in the literature 24,59,60 .
In this work, we solve large CDMFT clusters of up to 9 × 9 sites for the half-filled Hubbard model, and benchmark their convergence against exact results obtained with diagrammatic Monte Carlo (DiagMC) [61][62][63] in its connected determinant formulation (CDet 64 for one particle reducible quantities and ΣDDMC 65-67 for one particle irreducible quantities).
We show that if we use a center-focused extrapolation (CFE) to approximate the lattice self-energy, instead of averaging over the whole cluster as done in previous works, we obtain a much better approximation. This CFE converges faster and is in excellent agreement with the DiagMC result at moderate temperature. We further present results at lower temperature T and larger U than the range currently accessible with diagrammatic Monte Carlo methods.
The paper is organized as follows. In Section II we introduce the 2D Hubbard model and present the CDMFT formalism. We then present in Section III the center focused extrapolation technique, and benchmark it against DiagMC. In Section IV we discuss the spin-spin correlation function and finally summarize our results in Section V.

II. MODEL AND METHOD
We study the 2D single-orbital Hubbard model on the square lattice where c † i,σ (c i,σ ) denotes the creation (annihilation) operator for an electron with spin σ on site i, with the density operatorn i,σ = c † i,σ c i,σ and the chemical potential µ. t and U are nearest-neighbor hopping amplitude and onsite Hubbard repulsion respectively. We study the half-filled case which corresponds to µ = U/2. For a single-orbital unit cell the electronic dispersion relation in momentum space reads ε(k) = −2t(cos k x + cos k y ).
In the following, we compute propagators of single-and two-particle excitations in the CDMFT approximation 19 . CDMFT is a DMFT approximation on the super-lattice made of real space clusters on the original Bravais lattice 19 , where the cluster sites within the super-cell play the role of orbitals. Hence, the Green's functions and selfenergy contains inter-sites elements. They are computed using a self-consistent quantum impurity model, which is solved by a continuous time Quantum Monte-Carlo solver using an interaction expansion 68 , implemented with the TRIQS library 69 .
For the present study we consider square clusters of up to N × N = 9 × 9 atoms (see Fig. 1 for a schematic illustration). Technically, the equations for the Green's functions and self-energies are identical to those of a multiorbital (single-site) DMFT with a density-density interaction. The local lattice Green's function in Matsubara frequencies is computed by integrating over the reduced Brillouin zone (RBZ) of the super-lattice : (2) with i, j indexing cluster sites, dispersion relationε ij (k) (with k ∈ RBZ), and self-energy Σ ij (iω n ). The Weiss field is given by the standard equation 19 The self-energy acquires intersite components. Hence it can capture correlation effects up to distances of about the linear size of the cluster. The cluster self-energies break translational invariance for any finite cluster size. Transformation to lattice quantities in the Brillouin zone of the original lattice therefore requires restoring translational invariance for the Green's function, the self-energy, or its cumulants 18,19,23,24,59,[70][71][72] (see also Appendix B).

III. CENTER FOCUSED EXTRAPOLATION
In a very large cluster, we expect on general grounds the Green function at the center of the cluster to converge faster with cluster size than at the boundary. In this section, we use this insight to improve on the reperiodization within the CDMFT algorithm. We begin with the study of single-particle cluster quantities, and analyze the convergence of the local and the nearest neighbor self-energy w.r.t. cluster size. In a second step we use these converged self-energies to approximate the lattice self-energy, and benchmark it against exact DiagMC results.

A. Cluster quantities
Let us first concentrate on the cluster Green's and spectral functions. In Fig. 2 we show site-resolved single particle spectral functions obtained from an N ×N = 8×8 CDMFT calculation for different values of U . Exploiting cluster symmetry, we show only the upper right quadrant of the 64 cluster sites, where equivalent spectra are plotted in the same color. The violation of the translational invariance is most pronounced for U/t = 4 (solid lines) as can be seen by comparing the spectra in the center at R = (0, 0) with those in the corners at R = (3, 3). For U/t = 2 all spectra exhibit a finite quasiparticle weight at F , while for U/t = 12 they are all gaped.
In Fig. 3 we show the corresponding site-dependent value of G(R, τ = β/2). The spectra in the cluster center are generally more correlated/insulating than those at the edges/corners of the cluster. Given the fact that the exact solution has a critical value of U c = 0 at T = 0, the single particle spectrum of the central site appears to be a better approximation than the average over the whole cluster (dashed line in Fig. 3).

B. Convergence of cluster quantities
We now study systematically the convergence of the onsite and first neighbor self-energy with cluster size, and compare it to the exact lattice DiagMC result in the regime where it is available. In the left panels of Fig. 4 we plot the (purely imaginary) onsite self-energies of the four inner-core sites of the even clusters as a function of Matsubara frequency for U/t = 4 and two different temperatures. In the right panels we plot the (purely real) nearest neighbor self-energy corresponding to the inter-site correlation between two of the innermost cluster sites. The overall cluster-size dependence is pronounced in all plots and stronger for the lower temperature. In the insets of At high temperature (βt = 5, upper panels of Fig. 4) we compare CDMFT results to DiagMC benchmark data for Matsubara frequencies iω n up to n = 10 (red color). For both the onsite and nearest neighbor self-energies, the extrapolated CDMFT results are compatible with the DiagMC data within error bars. For the onsite self-energy (upper left panel), the extrapolated (and Di-agMC) self-energy is quite different from the largest considered cluster N × N = 8 × 8, especially at the first n = 0 Matsubara frequency. For the nearest-neighbor self-energy, the extrapolation effect extends to Matsubara frequencies as high as n = 10.
At lower temperature (βt = 12.5, lower panels of Fig. 4), which is out of reach of the DiagMC algorithm, the effects of the extrapolation are even stronger, and enhance the impact of the correlations.
FIG. 6. Imaginary (left) and real (right) part of the lattice self-energy Σ(k) at the first and second Matsubara frequency n = 0, 1 along a path in the first Brillouin zone, for an inverse temperature βt = 5 and different values of the interaction U/t. We plot the lattice self-energy obtained by conventional periodization of the 8 × 8 cluster (red), Fourier transformation of the real-space correlations relative to the central site for the 8 × 8 cluster (blue), Fourier transformation using the extrapolated real-space self-energy (black), and the numerically exact DiagMC data (green).

C. Center focused extrapolation (CFE) of the self-energy
After extrapolating the self-energy at the center of the cluster, we approximate the lattice self-energy in real space by a center focused extrapolation (CFE) defined in the following. For each lattice displacement R with vector elements greater than or equal to zero, we take the large N extrapolation of the cluster self-energy between the central site R c (c.f. Fig. 1) and R + R c , i.e.
For negative displacements we infer the value by using the rotational symmetry of the cluster. This procedure differs substantially from the conventional periodization, where averages over the whole cluster are performed 19 .
In practice, given that we solve CDMFT only up to N × N = 9 × 9, our strategy is two-fold. For the short distances R, we use the extrapolated self-energy values computed in the previous subsection. For the largest R however, where no 1/N extrapolation is possible given the lack of data points, we use instead the value of the largest cluster. This approximation is justified because the self-energy in real space Σ(iω n , R) decays on length scales smaller than our largest cluster. Indeed in Fig. 5 we show Σ(iω 1 , R) for βt = 5 and different values of the interaction U/t. Even for the largest interaction value U/t = 4 do we observe a decay of the absolute value of the self-energy at the first Matsubara frequency to zero within less than 5 lattice constants. In Fig. 6 we compare the Fourier transform of this estimator of the lattice self-energy (black) with the exact DiagMC computation (green) for the first two Matsubara frequencies (ω n={0,1} ), at the temperature where it is available. We further compare these to the selfenergy obtained by using the same procedure directly on the N ×N = 8×8 results, i.e. without the large N extrapolation (blue dashed), and to the "standard" periodization of the N × N = 8 × 8 self-energy (red) 19 . We find again an excellent agreement, within error bars, between the CDMFT+CFE and the exact DiagMC results. In particular, for the largest available U , we see a substantial improvement over the conventional self-energy periodization, and a sizable correction of the self-energy data obtained without extrapolation. We see in the second to lowest panel on the right side (βt = 5, U/t = 4, n = 0) that the CFE procedure significantly improves even details of the self-energy along the path Γ-X-M -Γ (i.e. nonmonotonous behavior between X and M points).

D. Lattice spectral function of CDMFT+CFE
The single particle spectra corresponding to the selfenergy data of Fig. 6 are presented in Fig. 7. In addition we present here also the spectrum for the cumulant periodization scheme applied to the N × N = 8 × 8 cluster (see Appendix B). The qualitative differences between the CFE results and the conventional schemes using data without extrapolation are sizable and most visible around the Fermi level. The gaps around X = (0, π) and between M and Γ at (π/2, π/2) are more pronounced due to the larger value of the CFE self-energies. This shows that the CFE has a qualitative impact on computed observables, even at high temperature.

IV. SPIN-SPIN CORRELATION FUNCTION
Let us now consider the spin susceptibility of the cluster The computation of the full lattice susceptibility, which uses the impurity two-particle vertex function as an input to the lattice Bethe-Salpeter equation 3 , is currently computationally too demanding for large clusters. However, in the limit of infinite cluster size, we expect the two quantities to coincide. In Fig. 8 (upper panel), we show the imaginary time onsite and nearest-neighbor spin-spin correlation functions for U/t = 4 and βt = 7.5, for different cluster sizes. For these parameters both components are i) sizable, ii) fully dynamic (i.e. beyond static meanfield approaches), and iii) antiferromagnetic (as expected for the present model). We see a quick convergence of the onsite susceptibility with N . The nearest-neighbor component for R = 1 has a stronger cluster dependence.
In the central panel of Fig. 8 we show the spatially resolved dynamical susceptibility χ(R, τ ) for the 8×8 cluster. Remarkably, the dynamic part of the susceptibility seems to decay much faster with distance than the static part (i.e. already χ(R = (1, 1), τ ) is barely dependent on τ ). This could indicate the presence of two length scales, one for coherent (i.e. dynamic) and one for non-local correlations that can be captured by static mean fields. The real-space decay of the static (i.e. τintegrated) spin susceptibility χ(R, Ω = 0) is presented on the bottom panel of Fig. 8 and compared to that of the corresponding self-energy at the smallest Matsubara frequency. The susceptibility data is well approximated by an Ornstein-Zernike form A exp(−R/ξ) (ξ/R) with a correlation length of ξ = 6.89. The fact that we manage to capture the exponential decay of spin-spin correlations reaching far beyond the scale of the clusters under consideration is indeed quite remarkable. A similar analysis for the self-energy proves to be less feasible, given the fact that it exhibits a much stronger radial dependence. We observe, however, that the self-energy decays overall on smaller length scales than the susceptibility.
In order to assess the validity of these susceptibility results we compare it for different cluster sizes N × N = 4 × 4, 6 × 6, 8 × 8 to the DiagMC benchmark in Fig. 9, where plot χ(R, iΩ n ) as a function of Matsubara frequency for the two displacements R = (0, 0) (left) and R = (1, 0) (right).
As for the self-energy, we perform an extrapolation to infinite cluster size (c.f. insets of Fig. 9) and obtain excellent agreement, within error bars, with the DiagMC data. We observe a strong dependence on N only for the displacement R = (1, 0) at the first Matsubara frequency.
Our analysis of the spin-spin susceptibility suggests that real-space cluster approaches such as CDMFT are valid even at temperatures and/or interaction values for which the magnetic correlation length exceeds the cluster size, as long as the self-energy decay in real-space is sufficiently captured. However, we expect that the feedback of this long range magnetic mode onto the electronic self-energy is not correctly described at this cluster size.

V. CONCLUSION
In this work, we have revisited the CDMFT calculations for the single band Hubbard model on the 2D square lattice at half filling. We performed a detailed momentum and real-space analysis of the spectral properties for different cluster sizes of up to 9 × 9 sites. Using a systematic benchmark with the exact DiagMC result at the lowest temperature for which it is obtainable, we have shown that an approximation scheme of the lattice selfenergy based on the center of the cluster is superior to the conventional periodization approaches based on averages over the cluster. In fact, the broken translational symmetry of the cluster Green's function and self-energy is simply a manifestation of the bulk-like nature of the central cluster sites, making them the proper basis for the approximation of the respective lattice quantities. We have further shown that the exponential decay of spinspin correlations is very well captured by CDMFT calculations, even when correlations extend far beyond the size of the cluster. Finally, CDMFT calculations and the CFE extrapolation can be carried out to temperatures currently inaccessible to exact diagrammatic Monte Carlo techniques, making the CDMFT+CFE procedure a powerful computational tool to access the physics of non-local correlations beyond dynamical mean-field theory.  In this appendix, we consider antiferromagnetic (AFM) ordering in our CDMFT computations.
In Fig. 10 we provide the phase boundary T N (U ) as determined by the CDMFT with the criterion: AFM or-der for m AFM ≥ 0.1 and PM one for m AFM ≤ 0.1. We consider cluster sizes up to 8 × 8. Due to commensurability with the AFM checkerboard symmetry, we restrict our analysis to even N (even though odd clusters could also be easily accessed in self-consistent cycles analogous to AFM single site DMFT). We observe a monotonous reduction of the Néel temperature with increasing cluster size, with a pronounced U -dependence: the reduction of T N in the strong coupling regime is significantly stronger compared to the small U region. This reduction is expected, and compatible with the exact Néel temperature of the 2D square lattice, which is known to be T N = 0 from the Mermin-Wagner theorem 74 . In this supplemental section we present directly postprocessed CDMFT data without extrapolation for various cluster sizes. In order to compute momentum dependent spectral functions in the Brillouin zone of the 2D square lattice we follow the cumulant periodization scheme of Ref. 24. Instead of periodizing the self-energy directly, its cumulant M (iω n ) = (iω n + µ − Σ(iω n )) −1 is periodized to obtain the lattice Green's function G(k, iω n ) by with R i , R j being the real-space positions of the sites i and j in the cluster, and the origin in R = (0, 0) corresponding to the central site for N odd and to one of the innermost sites for N even. We note that also the application of this scheme turned out to be problematic in some situations -which actually motivates the present analysis to account for the real-space anisotropies of CDMFT -but is used here as a reference since it was shown to perform better than re-periodizations of the self-energy or the Green's function 24 .
FIG. 11. Momentum resolved spectral functions along a closed path inside the first Brillouin zone (left) and full spectral function of the local Green's function (right), for U/t = 4, inverse temperature βt = 12.5, and different square clusters, as obtained by using the cumulant periodization. The dotted line corresponds to ω = 0. While the single-site DMFT (N = 1) is clearly in the metallic phase, the quasiparticle peak vanishes with increasing N , and N × N = 9 × 9 is already close to the insulating phase characterized by a gap in the local spectral function.
In Fig. 11 we show results for the single-particle spectral function as obtained from CDMFT for U/t = 4 and an inverse temperature βt = 12.5 by employing (B1) and subsequent analytical continuation with the Maximum Entropy method 75,76 . Each of the nine horizontal panels contains a momentum resolved intensity map for  A(ω, k) (left) on the path Γ-X-M -Γ in the cubic Brillouin zone and the k-integrated local spectral function A(ω) (right). The evolution from top ("DMFT") to bottom ("9×9") illustrates the cluster size dependence of the spectrum at fixed interaction and temperature. As function a of N , we observe an MIT transition which reflects in A(ω) as a spectral weight transfer from sharp quasiparticle excitations around the Fermi-level ε F to incoherent Hubbard bands, in agreement with previous studies on smaller cluster sizes 22,25 . However, the opening of the gap in A(ω) is a gradual process crossing a pseudogaplike regime between N × N = 4 × 4 and N × N = 7 × 7. This effect originates from a momentum selective opening of the gap (see also Ref. 77). Closer inspection of A(ω, k) for the intermediate regime of N × N = 5 × 5 and 6 × 6 clearly reveals the absence of a Fermi surface at the antinodal point X = (π, 0), while the quasiparticle band crossing ε F in the nodal direction around (π/2, π/2) remains sharp.
To complement our analysis of the MIT we also study the U and β dependencies of A(ω). Figure 12 shows A(ω) for N × N = 7 × 7 and U/t = 4. The temperature dependence, especially of the spectral weight at ε F , indicates the presence of a fully developed gap (in agreement with the corresponding N × N = 7 × 7 momentum resolved spectrum) which is barely closed by incoherent thermal smearing. Data without analytical continuation are reported in the lower panel of Fig. 13, with G(τ = β/2) = A(ω = 0).
Here our results for the MIT as a function of N , U , and temperature are summarized. Remarkably, the slope of the T c vs. U c transition line is inverted with respect to the single-site DMFT already for the smallest N × N = 2 × 2 cluster, in agreement with previous works 22 . Increasing N leads to a simultaneous increase of T c at fixed U and decrease of U c at fixed T which is consistent with previous calculations directly in the thermodynamic limit, e.g. DΓA 78,79 and DiagMC 67,80 , where the critical value shows the behavior U c → 0.