The Weak, the Strong and the Long Correlation Regimes of the Two-Dimensional Hubbard Model at Finite Temperature

We investigate the momentum-resolved spin and charge susceptibilities, as well as the chemical potential and double occupancy in the two-dimensional Hubbard model as functions of doping, temperature and interaction strength. Through these quantities, we identify a weak-coupling regime, a strong-coupling regime with short-range correlations and an intermediate-coupling regime with long magnetic correlation lengths. In the spin channel, we observe an additional crossover from commensurate to incommensurate correlations. In contrast, we find charge correlations to be only short ranged for all studied temperatures, which suggests that the spin and charge responses are decoupled. These findings were obtained by a novel connected determinant diagrammatic Monte Carlo algorithm for the computation of double expansions, which we introduce in this paper. This permits us to obtain numerically exact results at unprecedentedly low temperatures $T\geq 0.067$ for interactions up to $U\leq 8$, while working on arbitrarily large lattices. Our method also allows us to gain physical insights from investigating the analytic structure of perturbative series. We connect to previous work by studying smaller lattice geometries and report substantial finite-size effects.


I. INTRODUCTION
The two-dimensional Fermi-Hubbard model [1][2][3][4][5] plays the role of the fruit fly within strongly-correlated electron models. On the one hand, the model is a rich platform to investigate fundamental questions about the properties of interacting quantum systems. On the other hand, it is also believed to be directly relevant to the study of correlated electronic materials. It is for example conjectured to capture essential physical properties of hightemperature superconductors such as cuprates [6,7] even though the exact connection is subject to ongoing debate [8][9][10]. Many relevant experiments find intricately convoluted orders of spin and charge stripes, charge and pair density waves, and unconventional superconductivity which undeniably adds to the difficulty of disentangling the roles of individual phenomena [11][12][13][14][15]. While the Hubbard model has mainly been the subject of theoretical studies, it can now also be directly simulated at moderately-high temperatures by means of cold atoms trapped in optical lattices [16][17][18][19][20][21][22]. More recently, efforts in quantum computing technology have been aimed at the Hubbard model with the hope of finding effective algorithms for noisy near-term quantum hardware [23][24][25][26][27][28]. Large collaborative projects have been formed with the aim of comparing the results of numerous state-of-the-art algorithms in order to better understand the model as well as provide unbiased consensus benchmarks for future experimental and theoretical studies [29][30][31][32]. Despite the unquestionable and persistent interest of multiple scientific communities, a consensus phase diagram for the Hubbard model in the most intriguing parameter regimes is still missing [4,5].
Without doubt, a number of additional important unsolved questions remain to be answered before victory over the Hubbard model can be declared. To name but a few: Does the model host high-temperature superconductivity? If yes, how is it affected by competing magnetic instabilities? What is the relation between the onsets of strong spin and charge correlations? How does entering the strongly correlated regime manifest itself? Is a long correlation length a necessary ingredient? What is the relationship between the pseudogap and the coupling strength or magnetic correlation length?
Many previous works have addressed some of these questions. One of the best understood regimes of the model is half-filling (one electron per lattice site on average), where the ground state is an antiferromagnetic insulator and the Mermin-Wagner theorem prohibits longrange order at any finite temperature. As a result, a series of crossovers form from a metallic regime to a quasi-antiferromagnetic insulator with exponentially large correlation lengths [32][33][34][35], which becomes maximal at a particular finite interaction for any non-zero temperature. The nature of the insulating gap also changes as function of interaction, from Slater-like [36] to Mott-like [37], and reaches the Heisenberg regime in the infinite-interaction limit [38].
At finite temperature, one finds magnetic correlations in the spin and charge channels instead (together with superconducting ones). These have been documented by means of low-order diagrammatic Monte Carlo [54], determinant Quantum Monte Carlo (DQMC) [41], dynamical mean field theory (DMFT) and its cluster extensions [55][56][57][58][59], and by minimally entangled typical thermal states (METTS) [60]. These methods have also been able to capture the commensurate to incommensurate crossover, particularly in the spin channel, and as a function of temperature, doping and next-nearestneighbor hopping. Despite these advancements, the relation between the appearance of strong correlations in the spin and charge channels is not yet fully understood. Another open question is the connection between the onset of the strongly correlated regime and the occurrence of long magnetic correlation lengths, which is particularly difficult to answer by simulating small finite-size clusters.
The biggest challenge in this regime of the Hubbard model is the conciliation of the finite-temperature results with the ground-state calculations [4]. For many finite-temperature algorithms it is notoriously difficult to reach low enough temperatures to capture the relevant physics in order to make the connection, mainly due to the appearance of the fermionic sign problem. Almost all algorithms, whether ground state or finitetemperature, have to resolve to some sort of approximation. Notably, most of these cannot be formulated in the thermodynamic limit and thus suffer from finite-size effects [31,41,44,60], finite-momentum resolution [59], and/or effects coming from the choice of boundary conditions. In particular, there is a number of methods restricted to the analysis of the Hubbard model on widthfour cylinder geometries [31,41,60], whose correspondence to the model in the thermodynamic limit has not yet been thoroughly studied, mainly due to the lack of benchmarks from the latter.
In this work, we study the doped Hubbard model at finite temperature on large systems that do not suffer from finite-size effects This allows us to compute the spin and charge susceptibilities with a high momentum resolution, as functions of density, temperature and interaction strength. We observe three distinct regimes as we move in parameter space, as sketched in Fig. 1. We denote the first regime weak coupling (blue). It is found predominantly at small interactions and has short correlation lengths. Secondly, we observe a regime of what we call long correlations (red) centered around half-filling and at finite interactions. The system is almost longrange ordered with its magnetic correlation length growing rapidly with decreasing temperature. Notably, we find the correlations to be most pronounced at significantly lower interaction strengths than the values usually cited by literature [20,21]. Finally, the strong coupling (purple) regime is found predominantly at higher values of interactions and in a broader range of dopings. The system is more localized, marked by low values of the double-occupancy and short-range correlations, with the correlation length being only mildly affected by decreasing temperature. We also report the appearance of a Mott gap at half-filling in this regime. By evaluating the momentum-resolved spin susceptibility, we also investigate the nature of the magnetic ordering tendencies when the correlation length is sufficiently large, uncovering a commensurate to incommensurate crossover as a function of temperature, doping and interactions. Additionally, we find that spin correlations have much higher onset temperatures as compared to charge correlations that never become long-ranged in the range of temperatures we studied. This suggests that the spin and charge responses are decoupled over an extended region of parameter space.
The findings presented in this manuscript were made possible by a novel connected determinant diagrammatic Monte Carlo algorithm (CDet) for the computation of perturbative double expansions, which, in particular, allows us to construct series at fixed density away from half-filling. This permits us to compute numericallyexact results for large lattice sizes and at unprecedentedly low finite temperatures down to T ≥ 0.067 for interaction strengths U ≤ 8. As our algorithm allows us to study arbitrary lattice sizes and geometries we are able to quantify the finite-size effects which are present in var-ious other methods. Additionally, we show that much of the physics we have observed has a direct correspondence with the analytic properties of spin and charge susceptibilities as functions of complex interaction strength.
The paper is structured as follows: In Sec. II we introduce the Hubbard model and outline the main ingredients of our theoretical formalism. In Sec. III we present our key findings, focusing on mapping out the three distinct correlation regimes in Sec. III A, the commensurate to incommensurate crossover in Sec. III B, the decoupling of responses found in the spin and charge susceptibilities in Sec. III C, we discuss finite size lattice effects in Sec. III D, and we elaborate physical insights from the analytic structure of perturbative series in Sec. III E. We provide further details of our novel Diagrammatic Monte Carlo algorithm in Sec. IV and finish with a cumulative discussion of our results together with an outlook for future studies in Sec. V.

A. Hubbard model on the square lattice
The grand-canonical Hamiltonian for the Fermi-Hubbard model [1][2][3][4][5] reads: whereĉ † andĉ denote the fermionic creation and annihilation operators, k = (k x , k y ) the reciprocal lattice momentum, σ ∈ {↑, ↓} the fermionic spin, r = (x, y) labels lattice sites, U denotes the onsite repulsion strength, µ the chemical potential, and the square lattice dispersion relation is given by k = −2 t (cos k x + cos k y ), where t is the nearest-neighbor hopping amplitude (we set t = 1 in our units), andn rσ counts the number of particles with spin σ at site r.
In the majority of this work (with the exception of Sec. III D) we will present results for the model on a periodic square lattice of fixed linear system size L = 64 as this is found to be sufficient to eliminate finite-size effects for the range of temperatures we consider. This allows to easily gather statistics for all the L 2 values needed to construct the momentum-resolved susceptibilities, but there exists no technical limitation for our algorithm in considering larger systems or even directly the thermodynamic limit. We study thermal equilibrium properties at nonzero temperature T . In the following, the thermal average of an operatorÔ is denoted by Ô = 1 Z Tr e −Ĥ/TÔ where the grand-canonical partition function is defined as Z = Tr e −Ĥ/T , and we also use the imaginary-time Heisenberg picture for operators,Ô(τ ) = e τĤÔ e −τĤ .

B. Probes of spin and charge correlations
The main probes of spin and charge correlations we use in this work are the momentum-resolved spin susceptibility (in real and reciprocal space), which provides a quantitative measure of magnetic ordering tendencies, and the momentum-resolved charge susceptibility, which characterizes the response of the system to an inhomogeneous density perturbation. The spin susceptibility in real space is: whereŜ z (r) = 1 2 (n r↑ −n r↓ ) is the z-component of the spin operator. The Fourier transform of χ sp (r), χ sp (q), diverges at the onset of long-range magnetic order of wave vector q. As antiferromagnetic correlations are relevant near half-filling, we also consider the staggered spin susceptibility, defined in real space for r = (x, y) as: The charge susceptibility is defined by: where δn(r) = σn rσ − n, and n is the average number of particles per site n = σ n rσ . The doping is defined as 1 − n, the deviation from half-filling. The Fourier transform of χ ch (r), χ ch (q), diverges at the onset of longrange charge order of wave vector q.
To probe local correlations in our model we use the double occupancy: which quantifies the formation of local moments. The entrance into the Mott-insulating regime can be characterized by a plateau in the density n as a function of the chemical potential µ. Alternative quantities to probe spin and charge correlations are the equal-time structure factors: as well as their reciprocal space counterparts S sp (q) and S ch (q). These have been used for this goal in previous studies [30,41,59,60], along with the dynamic structure factor. However, the structure factors, while being experimentally relevant, have the disadvantage of not being directly related to ordering tendencies at finite temperature, where many phases are competing. Consequently, in this work we limit ourselves to the study of susceptibilities.
The correlation length is extracted by means of fitting momentum-space cuts of the susceptibilities using double Lorentzians with constant offsets. This procedure was also used in Ref. [41,59] for the structure factors. Unlike these previous works, we use a simple linear scale with cutoffs for our colormaps in both real and reciprocal space in order to give a fair account of the magnitude of correlations.

C. Theoretical formalism
Traditional finite temperature perturbation theory is defined in the grand-canonical ensemble. More specifically, one can write a so-called bare expansion for any quantity O in the form of a power series of the interaction strength U at fixed chemical potential µ as: To achieve better series convergence properties, it is advantageous to take into account the change in density at the mean-field level by introducing a Hartree shift to the chemical potential [61,62], which is diagrammatically equivalent to eliminating diagrams with tadpole insertions: where n 0 (µ 0 ) is the average number of particles per site at chemical potential µ 0 and zero interaction strength.
When studying spin and charge correlations as functions of an increasing on-site repulsion U , it is particularly important to disentangle variations due to rather trivial changes in particle density from proper correlation effects. For this reason, we introduce in this work the fixed-density expansion: the chemical potential µ = µ(n, U ) is renormalized as a function of U in such a way that the average number of particles per site, n, does not change with U : We refer to this expansion as the "fixed-density" scheme, and it is the method of choice in this work, while we still sparingly use the Hartree-shifted scheme in the vicinity of half-filling and at large values of U , as well as for cross-checking our results. The perturbative series are numerically computed to 8 − 10 expansion orders using a Monte Carlo procedure. We give more details about the formalism and the numerical methods in Sec. IV.

III. RESULTS
A. The three intermediate-temperature regimes of the Hubbard model We commence by inspecting the spin susceptibility χ sp (q) for all q wave-vectors in the first quadrant of the Brillouin zone of a 64 × 64 square lattice with periodic boundary conditions. A typical result for χ sp (q) is shown in Fig. 2 for U = 5, T = 0.1 and density n = 0.8. We would like to note that no momentum interpolation procedure has been used and each point corresponds to a uniquely calculated momentum value. For the parameters that we have investigated, the spin susceptibility displays one or several peaks close to q = (π, π). In Section III B we discuss how the location of the maxima and the associated commensurate to incommensurate crossover depend on the coupling U , the tempera- ture T and the density n. However, let us first focus our attention on the intensity of these peaks and the associated correlation length ξ, extracted by fitting the peaks with a double Lorentzian (plus constant offset) Ornstein-Zernike form (see right panels of Fig. 2).
[63] The correlation length and the maximum value of the spin susceptibility across the Brillouin zone are shown in Fig. 3 and are complemented with double occupancy and chemical potential data in Fig. 4. Together we use these quantities to identify three distinct intermediate-temperature regimes of the Hubbard model: a weak-coupling regime with short magnetic correlation length at small values of U , a regime with long magnetic correlation length at intermediate values of U and close to half-filling, and a strong-coupling regime at large U with strong shortrange spin correlations (see a sketch in Fig. 1). We discuss these three regimes in detail below.

The weak-coupling regime
For U 3 and T 0.1, the system is in a weakly correlated regime. The double occupancy gradually decreases with increasing U and slowly decreases with decreasing density, see Fig. 4 for T = 0.2. At the same time, the maximal value of the spin susceptibility as a function of q and the associated correlation length slowly increase with U as shown in Fig. 3. In this regime, the double occupancy decreases with increasing temperature (see inset in Fig. 4) because of the Pomeranchuk effect: an increase of the interaction strength U leads to a higher degree of localization and a favorable gain in entropy [33,[64][65][66][67]. This behaviour can be understood from the Maxwell relation ∂S/∂U T = −∂D/∂T U . We observe this Pomeranchuk effect for all densities at least down to n = 0.775. Close to half-filling, we can observe a small change of slope in the density versus chemical potential curve, indicating that the compressibility at half-filling is smaller than in the doped system, see right panel of Fig. 4.

The long correlation length regime
As the coupling U is further increased, the maximal value of the spin susceptibility as a function of q increases rapidly, even more so as one approaches half-filling, see Fig. 3. The position of the maximum for T = 0.2 is around U 4 − 5 and is only weakly dependent on the density. This maximum is in good agreement with the prediction of Ref. [54], which found the onset of magnetic spin correlations to occur at highest temperatures for U ∼ 4. Note that this maximum is significantly lower than the single-particle bandwidth 8t [20,21] and seemingly slowly decreases with decreasing temperature.
In this intermediate coupling regime, the correlation length becomes large (see Fig. 3), especially close to half-filling and as the temperature is lowered. It then decreases rapidly when the density is reduced and for T ≥ 0.1 it is only a couple of lattice sites long beyond 10% doping. As a result of this more pronounced difference in correlation length between different density levels, the double occupancy decreases more rapidly for occupancies closer to half-filling (see Fig. 4). This is natural as longer correlation lengths tend to favour more localized spins.
The curves of the density versus chemical potential seem to indicate that there is a wider region close to halffilling with a plateau where the compressibility would be small (Fig. 4). It is however difficult to obtain results in this regime, mainly due to the very long correlation lengths.

The strong local correlations regime
Finally, for U 6 we enter a regime of strong correlations (at temperatures T ∼ 0.2). The double occupancy that has decreased to a fraction (about 25%) of its noninteracting value and has become less dependent on the density. At the same time, these (quasi) antiferromagnetic correlations are very local and the spin susceptibil- ity and magnetic correlation length now decrease with increasing U , see Fig. 3, and display a slower increase as temperature is decreased as compared to the previous regime. Also, the temperature dependence of the maximum of the spin susceptibility as a function of q becomes weaker.
Another clear indication that this regime has strong correlations comes from the behavior of the density versus chemical potential. In this regime there is a clear plateau at half-filling showing a charge gap of the order ∆ ∼ U/2. This is compatible with a Mott insulating state at half-filling [68][69][70][71][72].

Crossover phase diagram
Our results for the magnetic correlation length are summarized in Fig. 5 with the help of an intensity map in the plane of doping and interaction strength as obtained from fitting (Q, π) as well as (Q, Q) momentum cuts. Both plots clearly display a dome centered around half-filling and U 5 where the correlation length is largest. The dome separates a regime of long correlation length for U ≤ 4 − 5 where spin-correlation theory may provide a reasonable description of the physics with a regime at larger U ≥ 6 where the physics is qualitatively different and dominated by short-range (quasi) antiferromagnetic correlations. It is interesting to note that in a recent CDet study [73] performed in the same range of parameters (T = 0.2 and U 4) a maximum in the entropy which develops away from half-filling was observed in exactly this regime of long correlations giving evidence for a crossover from metallic to non-Fermi-liquid behaviour.
From the right panel of Fig 4 we observe no evidence of phase separation occurring at temperature T = 0.2 (or at T = 0.1 which we have also investigated). The only region where we cannot make any conclusive statements is in the immediate vicinity of half-filling and at interaction strengths U ∼ 4 − 6 (see also Fig. 19), which has also been conjectured to be the most likely place for the occurrence of spatial phase separation by previous studies [12,43,44,[50][51][52][53]73].

B. Commensurate to incommensurate spin correlation crossover
In this section, we discuss the precise nature of the spin correlations in regimes where the correlation length is sizeable. We systematically investigate the structure of the spin susceptibility in momentum space and the staggered spin susceptibility in real space (see Fig. 6 for a graphical comparison to the non-staggered spin susceptibility) on large 64 × 64 lattices to study the commensurate-incommensurate crossover as a function of temperature, density and interaction strength.
Let us commence by showing the density dependence of the spin susceptibility χ sp at temperature T = 0.2. We set the interaction to U = 5 [74], where we have shown that the magnetic correlations extend furthest in real space (see Fig. 3). The results are shown in Fig. 7. Close to half-filling, the spin susceptibility in momentum space is strongly peaked at (π, π). Correspondingly, in real space, the staggered susceptibility shows extended uniform antiferromagnetic correlations. As the density is reduced, the peak first remains at (π, π) and below a critical density n 0.95 gradually splits into four separate peaks at (π ± δ s , π) and (π, π ± δ s ), compatible with the square lattice symmetry. The peaks become weaker and wider as the density is further reduced. The resulting real-space staggered spin susceptibility displays incommensurate spin correlations for densities below n 0.9 with domain walls separating π-shifted regions with antiferromagnetically correlated spins. The correlation length quickly becomes shorter as the density is reduced. Let us readily mention here that the onset of incommensurate spin correlations is never accompanied by a significant redistribution of the charge in the range of parameters that we have studied, see Section III C.
In Fig. 8 we additionally study the dependence of the crossover on interaction strength 3 ≤ U ≤ 7 for T = 0.1 and n = 0.875. We find that the four peaks present at weak-coupling gradually split further away from (π, π) as the interaction is increased. This corresponds to a gradually shrinking length of domain walls in real space. The increased splitting with increased interaction suggests that the leading wave-vector in the infinite interaction limit is likely different from (π, π) for densities away from half-filling, as was also found in Ref. [75]. We also observe four sub-leading momentum-space peaks forming at vectors (Q, Q) and hinting at possible diagonally-striped phases in the ground state which are commonly found in mean-field based approaches [75]. We have equally found a commensurate to incommensurate crossover as a function of temperature T , which is described in more detail in Appendix A (in particular see Fig. 20).
The commensurate to incommensurate crossover has previously been addressed by other numerical methods: DQMC has been used at finite temperature and predominantly on the 16×4 cylinder geometry, first for the threeorbital Hubbard model [76] and then for the single-band Hubbard model (that we focus on in this publication) in Ref. [41]. The authors studied temperatures T ≥ 0.22 for interactions in the range of 5 ≤ U ≤ 7. They found commensurate spin correlations in the vicinity of halffilling (n = 1) and in the incommensurate correlations in the doped regime without next-nearest-neighbor hopping t = 0 [77]. These findings have been confirmed by a recent DCA study on 8 × 8-sized embedded clusters in Ref. [59], which was performed for U = 6 and down to T = 0.167. Both of these studies have only found weak hints of correlations in the charge susceptibility. Another recent study using METTS on 16 × 4 cylinders did also find incommensurate correlations at U = 10 and 1/16th doping (n = 0.938), but at much lower temperatures T < 0.05. The authors also reported sizeable corresponding charge correlations picking up below these temperatures. It is also worth mentioning that these magnetic crossovers become a actual phase transition in the three-dimensional Hubbard model, where they have been studied by means of DMFT and its diagrammatic extensions [78].
Our results are summarized in Fig. 9, which shows the dominant wave-vector for different densities and values of U at temperature T = 0.2 and along the q = (Q, π) and q = (Q, Q) cuts. We find that the leading wave-vector is commensurate, q = (π, π), close to half-filling and becomes incommensurate either when the doping level is increased or when the interaction strength U is larger. We also observe that the incommensurability appears earlier along the (Q, π) than along the (Q, Q) cut. In real space, we generally find that, while the leading vector is always of (Q, π) nature, domain walls tend to form diagonally. This observation has also been made in previous studies Ref. [41,54,59] and was explained as the result of a superposition of a horizontal and a vertical stripe-pattern.

C. Incommensurate spin correlations with no charge redistribution
In this section, we want to clarify whether there is a particular charge response connected to the onset of incommensurate spin correlations discussed above.
Previous works [41,59,60] suggest that the formation of incommensurate spin correlations with a wave-vector (π ± δ s , π) is accompanied by incommensurate charge correlations with wave-vector (±δ c , 0) where δ c 2δ s . In Ref. [41] no conclusive evidence of charge correlations was found at finite temperature beyond what the authors identified as boundary effects. In [60] such a (±δ c , 0) peak became apparent below T = 0.05 for n = 15/16 and U = 10 and for a cylindrical width-four geometry. In Ref. [59] the authors found a broadened maximum around the (0, 0) wave-vector which they fitted with a double-Lorentzian that revealed two distinct maxima at (±δ c , 0). It is apparent that more conclusive evidence for such incommensurate peaks in the charge response is needed.
In Fig. 10, we fit the zero-frequency charge susceptibility χ ch (q) around (0, 0) along the cuts (Q, Q) (top) and (0, Q) (bottom) using double-Lorentzians with constant offsets, much like we have done for the spin counterpart in the previous sections. In the left panel, we use the same density n = 0.8 as Ref. [59] but at lower temperature T = 0.1 and slightly lower U = 5 (chosen to be close to the maximum of the spin correlations). Similarly to Ref. [59], we find a very flat peak around (0, 0) which, when fitted with our procedure, reveals weak maxima at incommensurate wave-vectors (±δ c , 0). In the right panel, we regard a higher density n = 0.892(6) at the same U = 5 and even lower temperature T = 0.067, where we can more clearly identify a broad peak at (±δ c , 0). The overall shape of the curve is reminiscent of what was observed in Ref. [60]. We proceed to investigate how the charge response depends on the parameters of our model. Fig. 11 shows the charge susceptibility χ ch in real and reciprocal space obtained from a Hartree-shifted (non-fixed density) series expanded around the noninteracting density n 0 = 0.8 and at temperature T = 0.067. The series is evaluated at interactions U = {0, 4, 5, 6, 7}, which corresponds to actual densities n = {0.8, 0.860(2), 0.892(6), 0.932(11), 0.978(20)}. We can identify four maxima at incommensurate wave-vectors (π ± δ, π) along with curved broad ridge features connecting them to wave-vectors (±δ c , 0). In real space we observe the formation of diagonal features which exhibit grouping of positive and negative values. It is, however, impossible to identify particular domain walls within these results.
Let us emphasize that these ridges in χ ch are characteristic of the Lindhard function and are already present in the non-interacting system [79]. They become less prominent at larger values of U and wash out as the temperature in increased. For temperatures down to T 0.067 and interactions U 5 at densities 0.8 n 0.875 we found a 20 − 25% increase in the maximum of the charge susceptibility in reciprocal space from our calculations as compared to low-order RPA. At the level of RPA these ridges are determined by nesting lines in the Brillouin zone which connect Fermi momenta with collinear Fermi velocities [56,79].
This absence of a significant charge response is further confirmed by studying the maximum value of the charge susceptibility χ ch (q max ) as a function of U and T as present for density n = 0.875 in the left panel Fig. 12. The behavior of χ ch (q max ) is in striking contrast to that of χ sp (q max ) (see Fig. 3). It does not display any maximum as a function of U and in general has a very weak temperature dependence.
We complete this analysis by studying the position of the leading wave-vectors given by δ s and δ c , along the (Q, π) and (Q, 0) directions, respectively, as well as their Left: Maximum of the charge susceptibility χ ch (qmax) for n = 0.875 as a function of the interaction strength U and temperature T . Right: Leading wave-vector computed from fitting with a double-Lorentzian (with constant offset) from the spin (δs) and charge (δc) susceptibility as well as the ratio between the two (RQ = δs/δc). Data is shown as a function of U for n = 0.8 and T = 0.1.
ratio R Q = δ s /δ c at n = 0.8 and T = 0.1 (see right panel of Fig. 12). Again, the quantities related to the spin and charge response behave qualitatively differently. Indeed, we find that δ c gradually increases with U while δ s stays roughly constant up until U ∼ 4 − 5, which is where the strong correlation regime starts, and then swiftly shifts further away from the (π, π) wave-vector. As a result, the ratio R Q starts decreasing and then increases beyond U = 4 − 5 until it reaches a value close to 0.5 at U = 7. While this is the value expected for a stripeordered phase, we see no evidence for significant charge correlations. An interesting question is whether the ratio remains fixed at 0.5 at larger values of U and whether stronger charge correlations eventually appear.
To summarize, our results support that, for the range of parameters that we have studied, the charge correla-tions remain uniform over the lattice even when longerrange spin correlations develop. From our observations we conclude that this absence of charge stripe correlations persists in the phase diagram at least down to T ∼ 0.10 and for values of U < 7.

D. System size dependence
In this section, we discuss the dependence of spin and charge correlations on the lattice size and shape. Such an analysis is important because multiple ground-state as well as finite-temperature methods can only study relatively small lattices thus making controlled extrapolations to the thermodynamic limit difficult. In particular, ground-state DMRG calculations are predominantly carried out on elongated geometries, mostly on width-four cylinders [31,[40][41][42], even though slightly larger widthfive and width-six cylinders have also been recently investigated [31].
In contrast, a recent VAFQMC study of the ground state [44] has analyzed system sizes of up to 16 × 16 as well as performed an extrapolation with system size. The authors found that stripe ordered phases compete closely with superconducting ones whilst spatial phase separation also occurs for a large range of parameter space.
At finite temperatures, DQMC calculations have also been mainly performed on the 16 × 4 lattice and benchmarked against an 8 × 8 lattice geometry [41]. The vertical stripy patterns found for 16 × 4 has not been reproduced by the 8 × 8 lattice, which lead the authors to conclude that a square lattice should be treated as a superposition of two stripes, one horizontal and one vertical, as the rotational symmetry is not broken. Another finite-temperature method METTS [60] has studied 32×4 sized cylindrical geometries with open boundary  conditions in the long direction and periodic boundary conditions in the short one. While the method is very effective in tackling any given temperature, it is extremely hard to treat cylinders with width larger than four. It is therefore of great interest to quantify the systematics with respect to an infinite lattice. Finally, a study using DCA [59] has been able to inspect up to 8×8 clusters embedded in a bath, which the authors suggest effectively reached the thermodynamic limit for the temperatures they computed. It should be stated that despite this major technological success the momentum-resolution of such results is still limited by the size of the cluster.
In Fig. 13, we study elongated geometries, including 16 × 4 and 32 × 4, but also study 32 × 8 and 32 × 16. We use the 64 × 64 lattice as a benchmark. At density n = 0.875 and for temperature T = 0.1 we chose U = 4 which is around the maximum of correlation length. We find that both 16 × 4 and 32 × 4 qualitatively reproduce the right picture in the sense that they reveal incommensurate correlations. However, we find that the maximum in reciprocal space is underestimated by roughly 40% and the length of the horizontal domain is roughly 9 sites as opposed to 13 for the thermodynamic limit result. We further see that the relative strength of highly non-local correlations gets enhanced on these cylindrical geometries. From the 32 × 8 and 32 × 16 geometries we see that even for sufficiently large geometries with broken rotational symmetry the horizontal striped patterns yield to a more complex two-dimensional pattern, which can be best described by the broadened reciprocal space peaks at (Q, π) vectors. The fact that the real-space pattern for 32 × 16 shows the circular character found in the 64 × 64 lattice rather than well defined vertical domain walls is indicating that the system does not have a strong nematic response.
In Fig.14, we additionally study square geometries of 8 × 8 and 16 × 16 for the same parameters. We observe that for too small lattices, such as 8 × 8 the incommensurability is effectively hidden as can be easily understood by regarding the finite momentum resolution in reciprocal space. Already for a 16 × 16 lattice the incommensurability is more apparent in real-space, while the reciprocal space plot is more inconclusive and only shows a broadened peak around the antiferromagnetic (π, π) vector. All computations in this work have been done on lattices with periodic boundary conditions. We do not attempt to study the effects of the choice of different boundary conditions in this work, however such a study has been recently performed in [59] and consider- able effects due to the choice of these have been reported.
We have also investigated the system size effects in the spin channel at higher temperatures (see Figs. 22,23 of Appendix B) and system size effect in the charge channel (see Figs. 24,25 of Appendix B) and found them to be a lot less severe in both cases. We conclude that finite size effects clearly worsen as magnetic correlations increase, which is the case when approaching half-filling and lowering the temperature.

E. Physical insights from complex plane structure
In this section, we will show that the magnetic and charge properties of the Hubbard model can to a large extent be inferred from the analytic structure of χ sp and χ ch , seen as a function of the complex interaction strength U , which one can attempt to reconstruct from the knowledge of perturbative series obtained with diagrammatic Monte Carlo (see Sec. II C, Sec. IV and Appendix C for more details on the method and formalism).
Using Padé approximants on the series coefficients, we have found the positions of the closest poles of χ sp in the complex U plane as obtained from fixed-density series within CDet. This is shown in the top plots of Fig. 15 at T = 0.2 for χ sp (q max ). For a given density n, the complex plane has a pole on the negative real U axis, most likely connected to a finite-temperature Kosterlitz-Thouless phase transition to the superconducting state of the negative-U (attractive) Hubbard model [80]. The re- maining two poles are on the positive U side and symmetrically placed about the real axis. Quite independently from the doping the real part of the poles is always close to U 5. As a result, one expects that the spin susceptibility will have a maximum around this value, which is indeed what we observe, see Fig. 3. The poles get closer to the real axis as the density approaches half-filling, thus yielding a larger spin susceptibility. Despite our ability to identify up to five poles in the complex plane (due to our choice of Padé approximants) the two additional poles appear significantly further away from the origin and are seemingly unstable with respect to slight changes in the series coefficients coming from a stochastic sampling of their error bars.
The right panel of Fig. 15 shows the same poles for the charge susceptibility. There, the positive U poles remain further from the real axis and are much less sensitive to the choice of density n, a different behavior leading to a smaller charge response of the system. The negative pole, on the other hand is found at roughly the same values of U as compared to its spin counterpart, which further corroborates the suggestion that this pole corresponds to an actual phase transition.
If we instead investigate the complex U plane of the Hartree-shifted series (bottom plots of Fig. 15) we find the positive poles for the spin susceptibility to be a lot closer to the real-axis as compared to the fixed-density series poles. This is logical as the density changes with U and the series are thus sensible to regimes with very high values of the spin susceptibility in the vicinity of half-filling.
The same analysis can be performed to investigate the the momentum dependence of the dominant peak in the spin susceptibility. Fig. 16 shows the closest poles of χ sp (q = (Q, π)) and χ sp (q = (Q, Q)) for two densities, n = 0.975 (top) and n = 0.725 (bottom). For densities in the vicinity of half-filling the closest poles are those that have a wave-vector q = (π, π) in accordance with the system developing commensurate spin correlations. As the density is reduced, the poles move further away from the real axis but also reorganize in a way that the closest poles to the origin are associated with an incommensurate wave-vector (Q, π). This makes sense since at low densities we find incommensurate correlations for practically all attainable values of U .

A. Diagrammatic Monte Carlo
In this work, we introduce and use a novel version of the diagrammatic Monte Carlo (DiagMC) method [81][82][83][84][85][86][87][88][89]. DiagMC has been used for lattice and continuum models with short and long-range interactions [90][91][92][93], for interacting topological models [94], for real-time propagation [95][96][97][98][99][100][101][102], and in combination with extensions of DMFT [103][104][105]. The main idea of DiagMC is to write a diagrammatic expansion for intensive physical quantities and to sample the diagrams of this expansion with a Monte Carlo procedure. As the expansion of any physical quantity can be written down directly for the thermodynamic limit, there is no overhead in considering arbitrary system sizes, thus circumventing the associated fermionic sign problem [106]. The Connected Determinant Diagrammatic Monte Carlo (CDet) algorithm [62] goes one step further in this direction by allowing to efficiently sum all Feynman diagrams topologies at given space-time positions of the interaction vertices, therefore reducing the computational effort.

B. Efficient evaluation of chemical-potential diagrammatic insertions
The technique we introduce in this work allows us to efficiently sum all Feynman diagrams for fixed vertex positions and arbitrary chemical-potential diagrammatic insertions, allowing in particular to fix the density perturbatively. More specifically, using the notation of Sec. II C, it is possible to write a double series in the interaction strength U and the chemical-potential shift α U We are able to efficiently compute the coefficients of this expansion, O kj , numerically, to high order. Diagrammatically, the double α-U expansion is given by bare < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 A N S h k S M X x R h a c K 3 Y s x i 7 / 2 0 7 j A = " > A A A C U H i c b V F B a x N B F H 6 b W o 1 p 1 a h H L 6 O h U K G E X Q l W P B W 8 e I x g 2 k I 2 h L e T t 8 n Q 2 d l l 5 m 1 L W P Y n e u n N 3 9 G L B 6 W d p C u 2 q Q + G + f i + b 9 6 8 + S Y p t H I c h j + D 1 t a j 7 c d P 2 k 8 7 O 7 v P n r / o v n x 1 7 P L S S h r J X O f 2 N E F H W h k a s W J N p 4 U l z B J N J 8 n Z l 5 V + c k 7 W q d x 8 5 2 V B k w z n R q V K I n t q 2 p 3 H G s 1 c k 4 h n p B l F v E C u T D 2 t 4 r J A a / O L e j / O k B d J W t n 6 I G Y s 3 z f W f 8 5 Z f m E 2 v G F 9 E H q j b X q / / T z t 9 s J + u C 7 x E E Q N 6 E F T w 2 n 3 0 v e V Z U a G p U b n x l F Y 8 K R C y 0 p q q j t x 6 a h A e Y Z z G n t o M C M 3 q d a B 1 G L P M z O R 5 t Y v w 2 L N 3 j 1 R Y e b c M k u 8 c z W w 2 9 R W 5 P + 0 c c n p p 0 m l T F E y G X l 7 U V p q w b l Y p S t m y p J k v f Q A p V V + V i E X a F G y / 4 O O D y H a f P J D c P y h H 3 3 s D 7 4 N e k e D J o 4 2 v I F 3 s A 8 R H M I R f I U h j E D C D 7 i C 3 / A n u A x + B d e t 4 N b 6 d 4 f X c K 9 a n R u 5 R b X K < / l a t e x i t > h n " (r, ⌧ ) n # (0, 0)i :

< l a t e x i t s h a 1 _ b a s e 6 4 = " 9 u p e k c Z t 0 4 M 8 V 4 k d z u Q j w m o O p W 0 = " > A A A C T n i c f V F N S y N B E O 3 J q h v j V 3 Y 9 e m k N g o K E G Q n u 4 k n w 4 t E F o 0 I m h J p O T d L Y 0 z N 0 1 + w S h v m F X s T b / o y 9 e F D E 7 c Q 5 + I U P G h 6 v X l V X v 4 4 y J S 3 5 / l + v 9 m V u f u F r f b G x t L y y u t b 8 9 v 3 c p r k R 2 B W p S s 1 l B B a V 1 N g l S Q o v M 4 O Q R A o v o q v j a f 3 i N x o r U 3 1 G k w z 7 C Y y 0 j K U A c t K g i a E C P V L I w y E q A h 6 O g Q p d D o o w z 8 C Y 9 E + 5 E y Z A 4 y g u T L k X E u S 7 l f U T p 1 / u + c 5 m q s m b h 4 N m y 2 / 7 M / D 3 J K h I i 1 U 4 H T R v w 2 E q 8 g Q 1 C Q X W 9 g I / o 3 4 B h q R Q W D b C 3 G I G 4 g p G 2 H N U Q 4 K 2 X 8 z i K P m 2 U 4 Y 8 T o 0 7 m v h M f d l R Q G L t J I m c c 7 q w f V u b i h / V e j n F P / u F 1 F l O q M X z R X G u O K V 8 m i 0 f S o O C 1 M Q R E E a 6 X b k Y g w F B 7 g c a L o T g 7 Z P f k / P 9 d n D Q 7 v z q t I 4 6
V R x 1 t s G 2 2 A 4 L 2 A 9 2 x E 7 Y K e s y w a 7 Z P 3 b P H r w b 7 8 5 7 9 J 6 e r T W v 6 l l n r 1 C r / w f 9 m 7 X i < / l a t e x i t > h n " (r, ⌧ ) n " (0, 0)i : Feynman diagrams with chemical-potential insertions, as shown in Fig. 17 for two building blocks of the perturbative expansions for the spin and charge susceptibilities (see Appendix C). The algorithm we introduce to efficiently compute these expansions is a generalization of CDet, similar in spirit to the method of Ref. [87]. We provide more details of our algorithm in Appendix C. We can use the chemical-potential insertions to work at fixed particle number in the grand-canonical ensemble by writing the expansion for the density: and by solving, order by order in U , the following equation for α(U ): where n(µ 0 , 0) ≡ n is a constant value and α(U ) = ∞ k=0 U k α k . The first term of the expansion for α(U ) is given by the mean-field value α 0 = n/2, but in general there is no analytical solution for higher-order corrections, apart from half-filling where α k =0 = 0. The computed α(U ) is then substituted into Eq. (10) to obtain a perturbative series in only one variable, U , for an arbitrary quantity O: The freedom in choosing the chemical-potential insertions α(U ) is used in this work to cross-check the results obtained by the fixed-density expansion. FIG. 18. Partial sums ( kmax k=0 U k O k ) for the spin susceptibility (left) χsp(q) and charge susceptibility χ ch (q) (right) at q = (π, π) for n = 0.875, U = 4.3, T = 0.2 obtained from the fixed-density expansion (green) and the Hartree-shifted expansion (blue). Extrapolated results using [4,5] Padé approximants are shown as horizontal bands.

C. High-order expansions and resummation
The series of Eq. (13) has a non-zero radius of convergence, therefore it provides an explicit definition of O(U ) for complex values of U . The radius of convergence is determined by the position of the closest singularity in the complex U plane: Fig. 16 and the discussion of Sec. III E provide a physical interpretation of this otherwise purely-mathematical fact [107]. For U larger than the convergence radius, we can apply analytic continuation as the singularities of the susceptibilities are not on the positive U axis. If that would not be the case, we would have a finite-temperature phase transition instead of a crossover, and the latter scenario is supported by the numerical data of Fig. 16 and physical expectations. In this work, analytic continuation is performed with the Padé approximants method [108], which constitutes the only potential source of systematic errors.
The convergence and the uncertainty of our results is checked by comparing different Padé schemes [109] to high orders (typically 8−10), taking into account the statistical noise on our numerical estimation of O k , and comparing the results of the Hartree-shifted and the fixeddensity expansions for the same chemical potentials. See Fig. 18 for representative examples of partial sums for perturbative series that were obtained in our calculations.
We expect, on physical grounds, the region of extremely-long correlation lengths to pose a currently insurmountable obstacle to our resummation schemes, and the numerical data of Fig. 19 confirms this picture. This issue is common to all numerical techniques, as shown by the multi-method half-filling benchmark study [32]. On the other hand, Fig. 19 also explains why previous CDet studies [62,109] were able to resum Hartree-shifted expansions far into the strong-coupling regime (up to U = 7 for n = 0.95 at T = 0.2 in Ref. where r l is a lattice site and τ l ∈ 0, 1 T is the imaginarytime. The integrals of Eq. (14) are evaluated with the Many-Configuration Markov-chain Monte Carlo method introduced in Ref. [110] and determinants are evaluated using a fast principal minor algorithm [111,112]. As briefly mentioned before, if the quantity O is well-defined in the thermodynamic limit, its series coefficients O kj are also well-defined, and there is no problem to extend the integration range of Eq. 14 to larger system sizes, or even directly to the thermodynamic limit.

V. CONCLUSIONS
In this paper, we have established numerically-exact results for the momentum-resolved charge and spin susceptibilities of the two-dimensional Fermi-Hubbard model. We have systematically studied a range of couplings 0 ≤ U ≤ 8, densities 0.7 ≤ n ≤ 1.0 and temperatures down to T = 0.067. We have shown that there are three distinct regimes separated by rather sharp crossovers: A weak-coupling regime with short magnetic correlation length, an intermediate coupling regime where the correlation length rapidly increases when the temperature is lowered or when the density approaches half-filling, and a strong coupling regime with short-range magnetic correlations. These regimes are also characterized by different behaviors of the double occupancy. In the weak coupling regime, the double occupancy de-creases with increasing temperature, displaying a Pomeranchuk effect. At strong coupling instead, the electrons are more localized and a Mott gap appears at half-filling.
In the regime with long magnetic correlation length, we established the existence of a crossover from q = (π, π) antiferromagnetic correlations to incommensurate correlations. This crossover is mainly seen as a function of increasing doping, but also as temperature is lowered or as the interaction is increased. The most common incommensurate vectors are (π, π ± δq), (π ± δq, π) but we found hints of the existence of peaks at (π ± δq, π ± δq), which we expect to become stronger in parts of the phase diagram at lower temperatures. We have shown that the different phenomena we observed can, to a large extent, be read from the analytic structure of χ sp and χ ch as function of a complex interaction strength.
In the parameter range that we investigated, we found that the charge susceptibility behaves very differently from the spin susceptibility and saw no evidence of a reinforcement of the charge response in regimes where the (incommensurate) magnetic correlation length increases rapidly. We attribute the existence of weak peaks at (±δq c , 0) to the presence of similar peaks at the noninteracting level in the Lindhard function. This absence of charge stripe correlations above T 0.067 raises interesting questions about the conciliation of the finite temperature properties of the system and its ground-state properties [30,31]. A natural scenario would be that, as temperature is lowered further, charge-stripe correlations eventually appear, as this is what happens in the width-four cylinder geometry [60], but further work is needed to confirm this picture. An additional question is whether such an ordered stripe phase would appear at finite temperature or only exists in the ground state.
Our findings from benchmarking different lattice geometries show that, especially as temperature is lowered and correlations increase, finite-size effects do lead to sizeable quantitative, and in the case of very small lattice sizes even quantitative, differences in results. It would be interesting to understand how these effects compare to the small energy differences that are observed between different competing ordered ground-state phases that have been studied on width-four cylinders and the 16 × 16 lattice [40][41][42]44].
In this work, we have limited ourselves to the discussion of the Hubbard model without next-nearestneighbor hopping t , which is however believed to play a crucial role in being able to connect the model to experimental findings on cuprate superconductors with effective t in the range of −0.4 ≤ t ≤ −0.1 [113]. Indeed, there has been an extensive amount of neutron scattering experiments performed on cuprates studying commensurate and incommensurate magnetic correlations [114][115][116][117][118][119][120][121][122]. In particular, it was found that electron doped cuprates show commensurate peaks [115,116], while experiments on hole-doped cuprates found the peaks to be incommensurate [119][120][121]. Another interesting difference was found between La-based cuprates and LBCO, where for the former both charge and spin ordering vectors increase with doping, while for the latter the charge vector decreases and the spin one increases. For the t = 0 case studied here we find behaviours consistent with La-based cuprates, similarly to Ref. [41]. It is, of course, impossible to talk about cuprates without mentioning superconductivity and its interplay with magnetic order. In that respect it would be of great interest to study superconducting d-wave correlations on equal footing with spin and charge ones, much like what has been recently carried out in Ref. [59]. It would be equally important to relate our findings with single-particle properties and to explore a possible connection to pseudogap physics. We, however, leave these investigations for future work.
The results obtained in this work were, to a great degree, enabled by algorithmic progress, in particular the development of a connected determinant diagrammatic Monte Carlo algorithm for double expansions, in U and an arbitrary chemical potential shift α(U ). Importantly, this allowed us to make the choice of the function α(U ) a posteriori and without any additional computational effort. Here, we have limited ourselves to only two such choices, the Hartree-series and the fixed-density series. However, one can generate an arbitrary number of alternative expansions from different choices of α(U ). This, on one hand, provides an additional degree of crosscontrol over the systematics coming from the resummations of perturbative series. On the other hand, this freedom of choice can potentially be exploited through optimisation schemes such as machine learning based algorithms.  diagonal domain walls across the crossover. The length of the domain around the center stays relatively unaffected by temperature at 5 lattice sites.
We also study the charge susceptibility in real and reciprocal space for the same set of parameters, where the spin correlations clearly undergo a commensurate to incommensurate crossover (see Fig. 21). The reciprocal space charge susceptibility χ ch (q) shows little structure except for a clear peak around (π, π). At the lowest available temperature T = 0.067 we observe the start of the formation of ridges as discussed in relation with Fig. 11, as well as the corresponding weak maxima at wavevectors (δ c , 0) and (0, δ c ). The real space charge susceptibility χ ch (r) always takes positive value at r = (0, 0), which is expected because because the system is holedoped n < 1 and the double-occupancy must be positive D > 0. There is no indication of domain walls forming, where the correlator would change sign as seen for the spin susceptibility for these parameters. We observe a similar behavior for all the temperatures that we have investigated here T ∈ {0.20, 0.167, 0.125, 0.10, 0.067}. It does therefore seem that the charge response is only very weakly affected by the formation of incommensurate spin correlations in this regime.
Appendix B: System size dependence of the spin susceptibility at higher temperatures In Fig. 22 and Fig. 23 we perform an additional study of the dependence of the spin susceptibility χ sp in real and reciprocal space on variations in system size and lattice geometry (rectangular vs. square). We have picked the same density n = 0.875 as in Figs  however, we have chosen to investigate here the higher temperature T = 0.22 (in accordance with Ref. [41]) and a somewhat higher interaction U = 5 as it is closer to the maximum in correlation length as a function of interaction strength. We observe that all geometries manage to reproduce the correct value and position of the maximum in reciprocal space, which is centered around (π, π), corresponding to commensurate correlations. This is consistent with the fact that lattice-size and boundary effect play an increasing role as correlations increase, which in what we universally observe when lowering the temperature.
In Fig. 24 and Fig. 25 we equally analyse the system size effects of the charge susceptibility χ ch , at density n = 0.875, temperature T = 0.1 and interaction strength U = 4. We observe that the overall real space picture is qualitatively similar for most lattices, even though the small geometries 16 × 4 and 8 × 8 show stronger nonlocal correlations as well as an enhanced peak around wave-vector (0, 0) as compared to the result on the large 64 × 64 lattice.
Appendix C: Chemical-potential insertions as double expansions

Connected Determinant Diagrammatic Monte Carlo
In this section we briefly reproduce the formalism of Ref. [62] in a form that is most convenient for its generalization. With the notation of Sec. II C, for any quantity O and for any real α, we define where hereafter a dependence on the chemical potential at U = 0, by µ 0 , is implied in all expressions. It is possible to express O α;k as the integral of the sum of all connected Feynman diagrams with {(r 1 , τ 1 ), . . . , (r k , τ k )} as space-time positions for the vertices where V = {(r 1 , τ 1 ), . . . , (r k , τ k )}, r l is a lattice site and τ l ∈ 0, 1 T the imaginary-time. The integrand O α ({(r 1 , τ 1 ), . . . , (r k , τ k )}) can be computed numerically  exactly for any quantity with exponentially scaling computational cost by using the CDet technique. As an example, we present here the computation of the density, O = n. The CDet algorithm allows to obtain the sum of all connected Feynman diagrams, n α (V ), by removing disconnected diagrams from the sum of all, connected and disconnected, Feynman diagrams, which we denote here by a α (V ), and which has the advantage to be easily expressed in form of determinants thanks to the Wick's theorem. Indeed, one has the explicit form: where V ≡ {(r 1 , τ 1 ), . . . , (r k , τ k )}, |V | is the cardinality of the set V corresponding to the expansion order, A α (V ) is a (|V | + 1) × (|V | + 1) matrix given by: We further set r |V |+1 ≡ 0, τ |V |+1 ≡ 0, and the non-interacting Green's function is defined by: where T is the time ordering operator and the average is computed at µ = µ 0 and U = 0. Then, n α (V ) is obtained from the recursive elimination of disconnected diagrams from the sum of all diagrams: where we have introduced The determinants of the matrices A α (S) and Z α (S) for all subsets S ⊆ V , needed for a α (S) and z α (S), can be computed in a computational time proportional to O(2 k ), where k ≡ |V |, by using a fast principal minor algorithm [111,112]. This leaves the recursive step, Eq. (C6), as the main computational bottleneck with a complexity of O(3 k ) (or, alternatively, O(k 2 2 k ) using fast subset convolutions [123]).
For completeness, we give the explicit expressions for the spin and charge susceptibilities at spin balance in terms of the F σσ functions: where F σσ (r, τ ) = δn σ (r, τ ) δn σ (0, 0) .
In order to obtain F σσ α (V ), analogously to what was done before, we introduce where B α is a (|V | + 2) × (|V | + 2) matrix defined by (B α (V )) uv = G (0) ((r u , τ u ), (r v , τ v )) − α δ u,vδu,|V |+1δu,|V |+2 , where r |V |+1 = r, τ |V |+1 = τ , r |V |+2 = 0, τ |V |+2 = 0, and A α is a (|V | + 1) × (|V | + 1) matrix such that (Ã) uv = (B) uv . We also define: We can then obtain F σσ from the recursive relation: which has the diagrammatic interpretation of the elimination from the sum of all, connected and disconnected, symmetrized Feynman diagrams for n σ (r, τ )n σ (0, 0) , A σσ α (V ), of two classes of diagrams: the non-necessarily connected diagrams contributing to n σ (r, τ ) n σ (0, 0) (second line in the previous equation), and the remaining disconnected diagrams (third line in the previous equation). The equivalent CDet recursion formula for single Hartree-type expansions of the spin and charge susceptibilities has been first introduced in Ref. [35]. and Eq. (C6) shows that n α (V ) is a polynomial of α of degree, at most 2 |V |, and from the Feynman-diagram representation it is easy to understand that the actual degree is only |V |: The same reasoning applies to any quantity O: O α (V ) = |V | j=0 α j O j (V ), where O j (V ) is the j-th chemical potential insertion in the sum of all connected Feynman diagrams for fixed space-time vertex positions V , from which we can determine the contribution to the physical quantity O after integration over space-time coordinates. In the following, we limit our discussion to O = n.
Our strategy to compute n j (V ), for j ≤ |V |, is the following: we first determine a j (S) and z j (S) for all S ⊆ V and j ≤ |V |; then, we apply Eq. (C6) seen as a polynomial equation to recursively obtain n j (V ). We use the following identities: where V = {(r 1 , τ 1 ), . . . , (r k , τ k )}. As briefly mentioned in the main text, one of the features of the algorithmic advancement we discuss here is that it allows to crosscheck the results obtained by different choices of α. This freedom has been proven to be extremely useful [61,90], even if in these previous works it was impossible to use the full functional freedom of α(U ) as a function of U , as there was no technical means to achieve this goal efficiently. In this work, we have limited ourselves to comparing the Hartree-shift choice, α = n 0 /2, and the fixed-density choice, i.e. α(U ) such that the density does not change as a function of U . See Fig. 26 for a representative comparison of the two set of series for different wave-vectors of the upper quarter of the Brillouin zone. These have two fairly different starting points, yet are converging to the same results within error bars.
It is also possible to generalize this formalism to arbitrary symmetry-breaking field insertions. A constant symmetry-breaking field shift has already been shown to be an efficient way to build a convergent perturbative series in the superfluid phase [124]. We further remark that a similar technique to the one we introduce in this work can be applied to a more general shift in the singleparticle propagator, which would allow even more flexibility in the control of the convergence of the series. For completeness we note that alternative exponentially scaling algorithms which alter perturbative expansions whilst still grouping Feynman diagrams into determinants have been introduced in Refs. [87,88,125].