Numerical approaches for calculating the low-field dc Hall coefficient of the doped Hubbard model

Using determinant Quantum Monte Carlo, we compare three methods of evaluating the dc Hall coefficient $R_H$ of the Hubbard model: the direct measurement of the off-diagonal current-current correlator $\chi_{xy}$ in a system coupled to a finite magnetic field (FF), $\chi_{xy}^{\text{FF}}$; the three-current linear response to an infinitesimal field as measured in the zero-field (ZF) Hubbard Hamiltonian, $\chi_{xy}^{\text{ZF}}$; and the leading order of the recurrent expansion $R_H^{(0)}$ in terms of thermodynamic susceptibilities. The two quantities $\chi_{xy}^{\text{FF}}$ and $\chi_{xy}^{\text{ZF}}$ can be compared directly in imaginary time. Proxies for $R_H$ constructed from the three-current correlator $\chi_{xy}^{\text{ZF}}$ can be determined under different simplifying assumptions and compared with $R_H^{(0)}$. We find these different quantities to be consistent with one another, validating previous conclusions about the close correspondence between Fermi surface topology and the sign of $R_H$, even for strongly correlated systems. These various quantities also provide a useful set of numerical tools for testing theoretical predictions about the full behavior of the Hall conductivity for strong correlations.


I. INTRODUCTION
Transport measurements are among the most common and accessible experimental probes, and are often among the first to be performed following the discovery of new materials. Yet, the theoretical investigation of normal state transport properties of quantum materials presents a number of unique challenges. Fermi liquid theory and the associated Boltzmann transport theory, which provide the theoretical framework for the understanding of ordinary metals, are known to break down in certain regimes. In the case of the high-T c cuprates, this is evidenced by the linear-in-T longitudinal resistivity, which violates the Mott-Ioffe-Regel (MIR) limit and has been synonymous with what has been called "strange metallicity" [1][2][3][4][5][6][7]. In addition, the Hall coefficient R H for cuprates has a strong temperature dependence [8,9], in contrast to predictions of Fermi liquid theory.
The Hubbard model, despite its simple form, successfully captures some of the non-Fermi liquid signatures in the normal state of cuprates. Strange metallic resistivity, without a signature of saturation at the MIR limit, has been successfully observed in the Hubbard model, both numerically [10] via determinant quantum Monte Carlo (DQMC) [11,12] simulations and in cold atoms experiments [13,14]. R H investigated in another recent DQMC work [15] also shows strong temperature dependence and a non-trivial peak at temperature T ∼ t, the kinetic energy scale, which may be connected to the rise of R H in cuprates like LSCO as temperature decreases from the ultra high temperature limit [16]. The close relationship between the experimental observations of cuprates and the theoretical results from the Hubbard model motivate us to continue investigating transport properties of the Hubbard model, specifically R H . We seek to better understand novel transport phenomena in materials without quasiparticles, and the relationship to the evolution of the electronic structure for understanding intertwined phases in the cuprates [17].
A calculation of the Hall coefficient is more complex compared to the longitudinal resistivity. Previous works have investigated R H in the t − J and Hubbard models for both high frequency and dc limits [18][19][20][21][22][23][24]; however, much of the work for strongly correlated models has involved certain simplifying assumptions, approximations, and limiting cases [22,25]. A faithful comparison of methods that can be used to calculate the dc Hall coefficient in a strongly correlated framework has been lacking.
Measuring R H in a numerically exact way poses challenges, limited by both the speed and efficiency of numerical techniques. For example, exact diagonalization is restricted to small lattice sizes [26]; density-matrix renormalization group (DMRG) is limited to a small number of exited states and may not be stable for calculating transport properties of 2D systems, especially in the metallic phase [27]. Numerical simulations for larger lattice sizes can be performed using quantum Monte Carlo (QMC) simulations in imaginary time for temperatures where the fermion sign problem is not too severe [28]. In the work presented here we use DQMC, a particular flavor of QMC.
One approach for obtaining R H = −B −1 ρ xy via QMC simulations is to explicitly couple the Hubbard model to a finite magnetic field. Current-current correlators χ αβ (τ ) (α,β = x or y direction) measured in imaginary time are then analytically continued to real frequency to obtain all components of the conductivity tensor σ αβ (ω). In this approach, explicitly adding a magnetic field B raises the computational complexity by requiring complex (as opposed to real) calculations. Apart from the inherent difficulty in properly incorporating the magnetic field B due to considerations of gauge invariance, this procedure also suffers from the need to analytically continue both the diagonal and off-diagonal components of σ αβ concurrently [29,30].
In an alternative approach, one could consider the zerofield limit by expanding the off-diagonal part of χ αβ up to linear terms in B. This method still requires analytic continuation, but avoids measurements in a finite field. However, in this approach one must evaluate a correlation function of higher order fermion operators (six fermions in the Hubbard model), which can increase error propagation. In addition, by introducing an extra imaginary time and an extra space index, the simulation becomes computationally more expensive. We provide additional detail in Appendix A.
In Ref. [25], another route was laid for studying R H numerically. As an application of the recursion method [31,32], this technique expands the Kubo formula of dc Hall conductivity in a Liouvillian representation into terms determined by magnetization matrix elements and Liouvillian matrix elements (or recurrents) [33] in a Krylov basis. By expanding H consist of thermodynamic susceptibilities, the expansion avoids the need for analytic continuation [25,34]. We refer to this method as the recurrent expansion. One drawback of this method is that the expansion is only conditionally convergent and its truncation error can be hard to estimate for strongly interacting systems.
In previous work [15], we investigated the dc Hall coefficient R H of the Hubbard model using DQMC to evaluate the leading order of the recurrent expansion R showing a strong temperature dependence -increasing with decreasing temperature -mimicking the behavior seen in cuprates [16]. Despite the strong temperature dependence deviating from Fermi liquid behavior, the sign of R H displayed a surprisingly close relationship with the Fermi surface topology, which has been usually understood as a feature of free electrons. As the interaction increases or the doping decreases towards half filling, R H changes sign concomitant to changes in Fermi surface topology. We argued that the "Hall coefficient sign -Fermi surface topology" correspondence may apply even for very strong interactions and low doping, in close proximity to a Mott insulator.
Higher-order corrections in the recurrent expansion of R H [25,34] could be large enough to qualitatively change this behavior. In the expansion described in Refs. [25,34], the higher order k th magnetization matrix elements and recurrents are constructed from correlators containing operators proportional to d k J α /dt k , where J α is the current operator along the α = x or y direction. d k J α /dt k may produce terms that include a number of fermion operators, which makes these higher order recurrents more computationally expensive to measure in comparison to R (0) H . Since the convergence rate of the recurrent expansion is hard to determine away from weak coupling, and higher order corrections are expensive to calculate, we consider the two approaches mentioned previously, which focus directly on the field response of χ αβ . As exact expressions measured in a well-controlled algorithm, they can be compared with our result for R (0) H [15]. The imaginary time dependence of χ αβ also contains real-time dynamic information about the conductivities.
In this work, we use numerically exact DQMC simulations to evaluate the dc Hall coefficient R H in the weak-field limit using multiple methods: • Recurrent Expansion -leading order R (0) H in terms of thermodynamic susceptibilities [15,25,34], • Zero Field (ZF) -the three-current linear-response of the off-diagonal part of the correlator χ αβ to first order in the magnetic field, χ ZF αβ . • Finite Field (FF) -directly evaluating χ αβ for a gauge invariant Hamiltonian in weak finite-fields on a finite-size lattice [35], χ FF αβ . We compare results from the latter two methods directly in imaginary time, finding a high degree of consistency, demonstrating that the DQMC algorithm is wellequipped to handle orbital effects of magnetic fields. To avoid the caveats of analytic continuation, we estimate various proxies for R H from the three-current correlation function. We find reasonable consistency with previous results for R (0) H [15]. These findings reaffirm the correspondence between the sign structure exhibited by R (0) H and the topology of the underlying Fermi surface, even in the limit of strong correlations that lack well-formed quasiparticles. In addition, we find that R H varies more slowly in Matsubara frequency than the individual longitudinal or transverse conductivities. We speculate that the cancellation of strong Matsubara frequency dependence of the individual conductivities also may be related to the observed correspondence between R H and the Fermi surface topology [15].
The remainder of this paper is organized as follows. In section II, we discuss the inclusion of orbital magnetic fields into the DQMC algorithm, and provide an expression of the zero-field linear response χ ZF xy and show the comparisons between the ZF and FF results in imaginary time. In section III we construct proxies for esti-mating the Hall coefficient from χ ZF xy and χ xx (taken as the ZF longitudinal response) and discuss the comparisons between them. We close with a discussion of our results and the challenges that remain for an evaluation of the full frequency dependence of the conductivities in the Hubbard model in a magnetic field.

II. CURRENT-CURRENT CORRELATION FUNCTIONS IN THE PRESENCE OF A MAGNETIC FIELD
In this section, we first discuss the inclusion of magnetic fields into the Hubbard model and derive an expression for the off-diagonal component of the current-current correlation function, to linear order in the magnetic field. We compare this directly to the current-current correlation measured under the lowest nonzero allowed field in imaginary time.
Here and throughout the paper, we have neglected Zeeman coupling of applied magnetic fields to spins and focus solely on the orbital contributions relevant for the Hall conductivity. The Hamiltonian of the Hubbard model in the presence of an orbital magnetic field is where t is nearest-neighbor hopping energy, µ is the chemical potential, U is the on-site repulsive interaction, c † r,σ (c r,σ ) is the creation (annihilation) operator for an electron at position r with spin σ, and n r,σ ≡ c † r,σ c r,σ is the number operator, with real-space lattice position r given by r = xe x + ye y , where e x and e y are unit vectors and the lattice constant is set to 1. The model is placed on a square lattice and we use periodic boundary conditions such that c r+Lxex ≡ c r+Lyey ≡ c r , unless otherwise specified, where L x and L y are the linear size of the system in the x and y directions, respectively. Here, θ r1,r2 = r2 r1 eA(r) · dr is the Peierls phase and A is the vector potential, and this Peierls phase integral is calculated using the shortest straight line path. Finally, the current density operator [36] is given by To circumvent the problem that the magnitude of a uniform field is limited to integer multiples of the flux quantum Φ 0 = 2π/e on a torus, we take the magnetic field to have a finite wavevector B(r) = B cos(qy) with q = 2π/L y in finding the zero-field linear response expression. With this choice, the field magnitude B can be arbitrarily small while maintaining periodic boundary conditions. We have verified that for the systems sizes investigated here, the value of q does not affect the overall results (See Appendix C).
Expanding the current-current correlation function to first order in B yields, where T τ is the imaginary time ordering operator, β = 1/k B T , T is the temperature, and we use the gauge where A(r) = − (B sin(qy)/q) e x . Here, · · · B denotes the expectation value taken with the full Hamiltonian, while · · · is taken with the B = 0 Hamiltonian. [37] Next, by making use of translation and reflection symmetries of the unperturbed Hamiltonian, we find where is the zero-field linear response, which we evaluate using DQMC simulations. We note in passing that a similar expression was derived for the case of a continuum model in the limit q → 0 [38][39][40].
In addition to χ ZF xy , we also consider the current-current correlation function in a uniform magnetic field. As discussed above, the smallest uniform magnetic field that can be applied corresponds to a single flux quantum through the system B = Φ 0 /V , where V = L x L y is the area of the system. We use the gauge A = B(−ye x + xe y )/2 and the corresponding modified periodic boundary conditions c r+Lxex ≡ c r e −i eBLxy/2 , c r+Lyey ≡ c r e i eBLyx/2 , following Ref. [41]. We then define the finite-field currentcurrent correlation function as We evaluate χ FF xy by performing a separate set of simulations in the presence of a magnetic field, and compare with χ ZF xy . The two expressions are expected to agree in the thermodynamic limit. Technical details of simulations can be found in Appendix A. While = 1 for convenience, in natural units the unit of conductivity σ is e 2 / and the unit of magnetic field strength B is /(ea 2 ) = /e for lattice constant a = 1. Therefore, our unit for R H is e −1 .
Results for χ ZF xy (τ ) and χ FF xy (τ ) (both purely imaginary) plotted against imaginary time are shown in Fig. 1. Generally, the transverse conductivity is reduced for large Comparisons between zero-field (ZF, solid lines) Im χ ZF xy (Eq. 3) and finite-field (FF, dashed lines) Im χ FF xy (τ ) (Eq. 6) in imaginary time for calculations on a 6 × 6 lattice. a A comparison at fixed filling (ρ = 0.9) and temperature (T /t = 0.5) for interaction strengths U/t = 4 − 16. b A comparison for fixed U/t = 12 and T /t = 0.5 for fillings ρ = 0.85 − 0.95. c A comparison for fixed U/t = 16 and ρ = 0.95 for temperatures T /t = 0.5 − 2. The ZF and FF values have the same units -e 3 t 2 . Error bars represent ±1 standard error of the mean, determined by jackknife resampling [42]. U , Fig. 1a, and as half filling is approached, Fig. 1b, as expected when charge fluctuations are suppressed at large U and particle-hole symmetry is restored at halffilling. The zero field (ZF) result χ ZF xy (τ ) also qualitatively matches the finite field (FF) result χ FF xy . The significant features, including doping dependence, temperature dependence, and imaginary time dependence, agree quite well between ZF and FF, which implies that when converted to σ xy (ω)/B, the two methods should also produce similar frequency dependent features. We verified that the small discrepancies between ZF and FF results are reduced for larger lattice sizes, as shown in Fig. 6 of Appendix C. We henceforth suppress the labels FF and ZF from χ ZF/FF xy unless it is needed. We observe that ∂ τ Im χ xy at τ = β/2 tends to increase with increasing U or decreasing doping, and may even change sign for certain parameters (see Fig. 5 in the Appendix for more details). As the off-diagonal conductivity is related to χ xy by (Appendix B) by symmetry, χ xy (τ = β/2) = 0. For U = 0, in the thermodynamic limit, ∂ τ χ xy does not depend on τ . In this case, Eq. 7 leads to ω Im σ xy (ω) ∝ δ(ω), as expected for infinite lifetime quasiparticles in the non-interacting limit. In Fig. 1a-b, we see that χ xy measured under the weakest interaction strength (U/t = 4) and lowest filling (ρ = 0.85) shows similar τ behavior to that expected for U = 0, and shows strong deviations with increasing U or as the system approaches a Mott insulator at half filling. The large curvature of Im χ xy (τ ) at small τ reflects the features of σ xy (ω) at ω ∼ U due to transitions to the upper Hubbard band. This becomes more pronounced at low temperatures and higher U as evident in Fig.1c.

III. PROXIES
We wish to obtain the zero frequency limit of the transverse conductivity from the imaginary time result shown in Fig. 1. Normally, this procedure for the longitudinal conductivity σ xx involves inverting [29,30] may be employed to obtain σ xx (ω) from χ xx (τ ). However this is not the case for the imaginary part of the transverse frequency dependent conductivity Im σ xy , and as a result MEM techniques encounter difficulties. This can be seen in Fig. 1 wherein Im χ xy (τ ) can change sign in the range [0, β/2), while the kernel K(τ, ω) in Eq. 7 does not change sign.
In addition, the six-fermion correlator in χ ZF xy is computationally expensive to measure and suffers from large numerical errors. Therefore, in this section, to compare our χ ZF xy result obtained from Eq. 5, with R H obtained in Ref. [15], we construct proxies for dc R H using χ ZF xy , and compare these proxies to R (0) H . We consider two types of proxies that we derive and discuss below: • D type -stemming from an analogy to Drude theory: expressed as which is obtained by inserting the Drude formulas into Eq. 7 and 8 and taking the limit γ → 0, where γ is the scattering rate. Another candidate proxy D γ is constructed by assuming γ to be non-zero and fitting χ xx and χ ZF xy using Eqs. 7, 8, 10, and 11. Results of proxies D and D γ are shown in Fig. 2a-c. • M type -determined by extracting the zero Matsubara frequency limit of: where χ αβ (i ω n ) in Matsubara frequency is defined as the Fourier transform of the imaginary time data, given by χ αβ (i ω n ) = β 0 dτ e iωnτ χ αβ (τ ). Here we utilize a cubic spline interpolant. Another candidate proxy M2 is defined as where ω 1 = 2π/β, the smallest nonzero Matsubara frequency. Results of proxies M1 and M2 are shown in Fig. 2d-f. H from previous work [15], at least within a factor of order unity. Proxy M2, on the other hand, consistently produces results that are much smaller (about 1/4 of M1) in magnitude. Calculation details of these proxies are in Appendix D.
For a Fermi liquid with a momentum-independent scattering rate γ, theoretical work [8,[43][44][45][46] has demonstrated a Drude-like ω dependence of the conductivities in Eq. 10 and 11 for ω γ. Therefore, we expect proxy D to be in good agreement with the true R H for a Fermi liquid where the scattering is isotropic and weak. Details of the derivation of Eq. 9 can be found in Appendix D. Figures 2a-c  H . Proxy D γ assumes γ to be non-zero and finite. As shown in Fig. 2a-c, the results for proxy D γ are largely the same as those of proxy D, implying that under the assumptions of Drude theory, our estimation of R H is not sensitive to changes in relaxation rate. As we can see from Eq. 10 and Eq. 11, any direct effect of γ cancels out in R H .
There are some limitations for D-type proxies. Conductivities can deviate significantly from Drude theory for strongly interacting systems, leaving the approximate R H far from the true results. In addition, from solely Eq. 7 and 10, one can never obtain a Im χ xy (τ ) that changes sign as a function of τ in the range τ ∈ [0, β/2), which is an important feature in our Im χ xy (τ ) data due to interaction effects. Now we switch to the M type proxies. Proxy M1 using Eq. 12 is exact for dc R H in the zero-temperature limit T Λ (ω n Λ), where Λ is defined as the scale on which χ αβ begins to deviate from its low-frequency behavior. This is because (see Appendix B and [21]). Proxy M2 approximates σ xx (ω = 0) with and uses ω = ω 1 to estimate Eq. 14 at ω = 0. Equation 16 is also exact at T Λ, and is often used as a proxy for σ xx in other work [10].
Figures 2d-f show proxies M1 and M2. Values of M1 are overall close to those for both proxy D and D γ , even for strong interactions up to U = 16 t. Indeed generally, as long as the rapidly varying components in σ xy (iω) and σ xx (iω) cancel out in R M1 H (iω), proxy M1 is very accurate regardless of the explicit form of frequency dependence of the conductivities. In other words, proxy M1 only requires the ratio of the conductivities, R M1 H (iω n ), to vary slowly with ω n , as we show in Fig. 3, where we plot R M1 H (i ω n ) against n (similar to Ref. [21]). Proxy M2 is not able to make use of cancellation between the transverse and longitudinal conductivities, so it still requires T Λ for both conductivity components. The difference between values of M1 and M2 in our simulations H . M2 has been multiplied by a factor of 4 to demonstrate that it has similar temperature and doping dependence, as well as sign change structure, as other proxies, but with a reduced magnitude. Panels in the same row share the same legend. For all proxies, χ ZF xy (τ ) is measured on a 6 × 6 lattice, and χxx(τ ) is measured on a 8 × 8 lattice. Pairs of a dotted and solid lines in panels a-f share the same parameters. a-c Proxy D (solid) and proxy Dγ (dotted). d-f Proxy M1 (solid) and proxy M2 (×4, dotted). g-i Dashed lines are R (0) H , defined by the leading order term in the expression for RH in Refs. [25] and [34], constructed from thermodynamic susceptibilities. These are evaluated on 8 × 8 lattices. Error bars represent ±1 standard error determined by jackknife resampling.
indicates that Λ T , so that σ xy (iω) and σ xx (iω) vary significantly within the scale set by the smallest non-zero Matsubara frequency ω 1 .
Considering that D and M type proxies have different assumptions and approach R H from quite different aspects, it is remarkable that D, D γ , and M1 are overall comparable to each other. Even more remarkably, all the proxies compare particularly well to the sign changing structure to R (0) H . Also considering the previous comparison of ZF results with those of FF in Sec. II, we conclude that the previous method we used to calculate R (0) H [15] as an approximation for R H is reliable for large U , de-spite our neglect of higher order corrections. As we have shown in this section, the proxies constructed using χ ZF xy and χ xx are also a useful approximation when direct analytic continuation to find σ xy (ω) is challenging.
We remark that within the temperatures accessible for DQMC simulation, our R H results do not show a strong doping dependence away from half-filling, specifically in the region around 20% hole doping as referenced in cuprate experiments [47][48][49].  Fig. 2. χ ZF xy was calculated on a 6 × 6 lattice, while χxx was calculated on an 8 × 8 lattice. Error bars for n > 0 are smaller than the points.

IV. DISCUSSION
Through a consideration of various methods to evaluate the Hall coefficient, we have demonstrated agreement between the sign change structures of χ ZF xy , χ FF xy , and R (0) H , supporting prior claims that the sign of the dc Hall coefficient has a close relationship to the Fermi surface topology in the strongly correlated zero-field Hubbard model [15]. The rough agreement between proxy D, D γ , M1 and R (0) H , as well as the significantly dif-ferent result of proxy M2, suggests that making use of cancellation between transport quantities may simplify evaluations of difficult multifermion correlation functions such as the Hall coefficient and allow us to construct a good description of transport without quasiparticles. The following facts lend support to our idea. While longitudinal resistivity in the Hubbard model shows typical non-Fermi liquid behavior [10], R H shows relatively flat ω n -dependence in Fig. 3. In this work, proxy D and D γ provide similar results. This likely results from our assumption that γ is the same for the two conductivities, leading indirectly to simplifications that allow D γ to mimic D due to apparent cancellation of lifetime effects. In constructing proxy M2, there are no such cancellations, and the proxy is fragile and easily fails. Finally, Refs. [25] and [34] showed how ratios of conductivities like the Hall coefficient or thermal Hall coefficient reduce to expressions constructed from simple thermodynamic susceptibilities. It is an open and intriguing question whether a Fermi-liquid like correspondence between Fermi surface topology and R H in Ref. [15] is due to such cancellations between conductivities, rather than necessarily Fermi-liquid-like ω-dependence of each conductivity.
It remains a challenge to perform analytic continuation directly from Matsubara frequencies or imaginary time data. A promising approach would be to use techniques designed to treat non-positive-definite spectra [50,51], or other methods of analytic continuation [52][53][54][55]. If a reliable method of analytic continuation can be found, then our evaluation of the exact three-current linear-response χ ZF xy through the numerically exact and unbiased DQMC algorithm will allow us to find the exact σ xy (ω) spectra for all frequencies for the Hubbard model. Even in the absence of reliable analytic continuation methods, our χ ZF xy results will still be a benchmark for any theory that proposes a frequency dependence of σ xy for the strongly correlated Hubbard model or similar models.

Appendix A: simulation details
We use DQMC simulations to measure χ ZF xy (τ ) for the Hubbard model in Eq. 1 without B and χ FF xy (τ ) under eB = 2π/V , respectively, on a 2D square lattice with the corresponding periodic boundary conditions [41]. Strategies for evaluating Green's functions in this work are the same as those of Ref. [10].
We can write down an alternate expression for χ ZF xy : Because of C 4 symmetry, we confirm that χ ZF xy given in Eqs. 5 and A1 give identical results and report the average of the two quantities to reduce sampling errors and improve statistical measurements. Similarly, since χ FF xy (τ ) = −χ FF yx (τ ), we measure and report (χ FF xy − χ FF yx )/2. In our simulations for χ ZF xy , we measure T τ J(τ )J(τ )J(0) at discretized τ values. The integral over τ is computed by finding two cubic spline interpolants using data points in the ranges τ ∈ [0, τ ] and τ ∈ [τ, β] respectively (here T τ J(τ )J(τ )J(0) is discontinuous at τ = τ due to time ordering) and then integrating the interpolants.
To tune the chemical potential µ for a specific target filling level ρ tar at a specific temperature, we use DQMC to calculate ρ for a range of chemical potentials µ (at 0.05t intervals for ZF case and 0.1t intervals for FF case) and obtain the best µ by linear interpolation of the ρ versus chemical potential curve. In the ZF case, µ is tuned on a 8 × 8 lattice and the same µ is used for both 8 × 8 and 6 × 6 lattice sizes. In the FF case, µ tuning for 6 × 6 and 8 × 8 lattices are done separately, although in practice, the optimal µ is almost identical for the two lattice sizes. For our parameters, we can obtain ρ to within a tolerance of O(10 −3 ) of the target density ρ tar (written as ρ throughout this work). An example of this tolerance in the ZF case is shown in Fig. 4. Including the effects of the Trotter error, DQMC statistical error, and the density shift between lattice sizes for specific values of µ, our density is accurate to O(10 −3 ).
Regarding the Trotter error, we define dτ as the interval between imaginary time data points. In previous work [15], we set a minimum partition of imaginary time L = β/dτ = 20 and a maximum dτ = 0.1/t for all interactions and temperatures. In this work, for U/t = 4 − 8, all ZF calculations (χ ZF xy , χ xx , R H , and ρ ) use the same dτ as in the previous work [15], while FF calculations (χ FF xy and ρ ) also have maximum dτ = 0.1/t, but have minimum L = 10. So the imaginary time spacing of χ FF xy is larger than χ ZF xy at the highest temperatures, as shown in Fig. 5a-d for U/t = 4 and Fig. 6 for U/t = 8. For U/t = 12 − 16, χ ZF xy , χ xx , χ FF xy , and R H are all obtained with a maximum dτ = 0.05/t and minimum L = 20, in order to reduce the Trotter error. Measurement of ρ in tuning of µ is done using a minimum L = 20, and a maximum dτ = 0.01/t (ZF) and dτ = 0.02/t (FF). In summary, U t(dτ ) 2 ≤ 0.08 in this work, comfortably below the conventionally adopted limit U t(dτ ) 2 ≤ 1/8. We consider Trotter error to be negligible. The accuracy of our results is affected primarily by the limitations of individual proxies, as discussed in the main text. In addition to Trotter error discussed here, our only other source of systematic error is finite-size effects discussed in Appendix C.
We run up to approximately 600 independently seeded Markov chains for χ ZF xy , about 40 Markov chains for χ xx and R  Fig. 1 with a direct comparison between the ZF Im χ ZF xy (τ ) and FF Im χ FF xy (τ ). Calculations were performed on 6 × 6 lattices using the DQMC algorithm. Error bars are ± 1 standard error determined by jackknife resampling. The same color (marker) within each panel represents the same parameters, as referenced in the legends in the final column. Rows have fixed interaction strength U and columns have the fixed temperature T .
tions to conductivities, we assume the thermodynamic limit. In that limit, χ ZF xy is considered to be the linear response to an infinitesimal uniform magnetic field B, i.e. χ ZF xy = lim B→0 χ xy . On the other hand, in defining χ ZF xy in Eq. 5, we used a finite q and a non-uniform magnetic field to ensure that the expression is well-defined on a finite lattice. This approach is reasonable, as the conductivities obtained in this way converge to their values in the thermodynamic limit as q → 0.
These are Eq. 7 and 8 in the main text. Comparing Eq. B2 and Eq. B4 and considering Eq. B10, we also conclude that Appendix C: Finite size effects A comparison of ZF and FF data obtained using different lattice sizes for U/t = 8 and U/t = 12 is shown in Fig. 6 and 7, respectively. Fig. 6 reveals that the FF results have smaller finite-size effects than the ZF results. As we have discussed in the main text, we expect the ZF and FF results to be identical in the thermodynamic limit, and that trend is indeed shown in Fig. 6 and Fig. 7.
For proxy D γ , We fit a few values of χ ZF xy and χ xx near τ = β/2 to Eq. 10 and 11 through Eq. 7 and 8. In doing so, we assume that the γs are equal in σ xx and σ xy . We also make use of the fact that values of χ xx (τ ) and χ xy (τ ) near τ = β/2 are determined more predominately by the low-frequency Drude-like behavior of conductivities, compared with other τ values. By fitting χ xx we extract γ and Ω xx ; and by fitting χ ZF xy we extract Ω xy /B. We choose to use σ xx rather than σ xy to find γ for several reasons: DQMC measurements of χ xx (τ ) have smaller numerical errors than χ ZF xy (τ ), and we could use the second derivative of χ xx to estimate γ, but we need to use the third derivative of χ xy to estimate γ. Correspondingly, using χ ZF xy would have required us to fit more data points around β/2 to calculate γ. In general, we need to use at least two data points for the χ xx fit, and at least one data point for the χ ZF xy fit. When numerical errors are large, we choose to include more data points in order to obtain a accurate fit. The downside of including points further away from τ = β/2 is that we must assume Drude-type behavior holds on a wider frequency range for σ xx and σ xy . Note that χ xx (τ ) and χ ZF xy (τ ) are symmetric and antisymmetric about τ = β/2, respectively, so we only use data points on one side of τ = β/2.
For the M-type proxy, when we calculate χ αβ (i ω n ), we use a cubic spline to fit χ αβ (τ ) and insert 10000 sampling points on the imaginary time axis, and integrate the oscillatory function as a function of τ using the composite trapezoidal rule.
For proxy M1, the errors from χ xx (i ω n ) are neglected and the error bar is ±1 standard error determined by jackknife resampling of χ ZF xy (i ω n ). The cubic spline extrapolation to ω n = 0 utilizes data on the first three nonzero ω n . The error bars for proxy M2 are constructed by error propagation of the standard errors in χ ZF xy (i ω 1 ) and [χ xx (τ = β/2)] 2 , which are themselves determined by jackknife resampling.