Thermodynamics of correlated electrons in a magnetic field

The Hofstadter-Hubbard model captures the physics of strongly correlated electrons in an applied magnetic field, which is relevant to many recent experiments on Moir\'e materials. Few large-scale, numerically exact simulations exists for this model. In this work, we simulate the Hubbard-Hofstadter model using the determinant quantum Monte Carlo (DQMC) algorithm. We report the field and Hubbard interaction strength dependence of charge compressibility, fermion sign, local moment, magnetic structure factor, and specific heat. The gross structure of magnetic Bloch bands and band gaps determined by the non-interacting Hofstadter spectrum is preserved in the presence of $U$. Incompressible regions of the phase diagram have improved fermion sign. At half filling and intermediate and larger couplings, a strong orbital magnetic field delocalizes electrons and reduces the effect of Hubbard $U$ on thermodynamic properties of the system.


I. INTRODUCTION
Strong magnetic fields allow us to probe the phase diagram of strongly correlated materials and uncover novel phases. For example, a magnetic field induces charge/pair density wave order in cuprate superconductors [1][2][3][4]; magnetic field is a convenient tuning parameter for accessing quantum critical points [5,6]; and field-induced reentrant superconductivity also has been reported in uranium compounds [7,8]. With the recent proliferation of experimental evidence for fractional quantum Hall effect, superconductivity, and other correlated electron phases [9][10][11][12][13][14] in graphene Moiré superlattices, there is renewed interest in studying the behavior of strongly correlated electronic systems in strong magnetic fields.
Properties of non-interacting electrons in a twodimensional periodic lattice under the influence of a strong magnetic field are fairly well-understood. In this system, the competition between lattice and magnetic length scales leads to the fractal Hofstadter butterfly spectrum with recursive magnetic subband structure [15,16], which generalizes the idea of Landau levels in a free electron gas. The Chern numbers associated with these magnetic subbands provide an elegant explanation of the integer quantum Hall effect [17]. Experimentally, the Hofstadter Hamiltonian has been realized in ultra-cold atoms loaded on optical lattices [18], and direct observation of the Hofstadter spectrum has been reported in Moiré superlattices in graphene with high resolution [19,20].
The most natural framework for understanding the simultaneous influence of magnetic field and Coulomb interaction on electrons in a periodic lattice is to take * jxding@stanford.edu † tpd@stanford.edu the Hofstadter Hamiltonian and add to it a Hubbard interaction term. In the literature, this is sometimes called the Hofstadter-Hubbard or Hubbard-Hofstadter model. This model has been investigated using Hartree-Fock mean-field theory [21,22], exact diagonalization [23,24], dynamical mean-field theory [25,26], and in the large U limit via renormalized mean-field theory [27]. Aside from exact diagonalization, which is limited to small system sizes, all methods used to study the Hubbard-Hofstadter model have been approximate and don't capture the full extent of quantum fluctuations. It is not conclusive, for example, whether interactions change or preserve the gap structure of the Hofstadter butterfly [22,24,26]. Determinant quantum Monte Carlo (DQMC) [28][29][30] is an unbiased and numerically exact algorithm for studying quantum systems at finite temperature. It employs a discrete Hubbard-Stratonovich transformation to reduce the quartic Hubbard interaction term to quadratic at the cost of introducing a fluctuating auxiliary field. This auxiliary field is then sampled using the Metropolis-Hastings algorithm. DQMC has been employed successfully in the (zero-field) Hubbard model to study spin and charge excitations [31,32] and superconducting fluctuations [33], as well as find evidence for fluctuating stripes [34] and T -linear resistivity [35]. The DQMC method is especially powerful at half-filling in the absence of kinetic frustration [36], where the fermion sign problem is absent due to particle-hole symmetry, even in the presence of a magnetic field. This allows simulations to be performed at much lower temperatures, providing access to properties more reflective of the ground state.
At half filling, the Fermi surface of the non-interacting Hofstadter model consists of a finite number of Dirac points at even-denominator rational fractions of magnetic flux per plaquette [37]. The ground state of the Hubbard-Hofstadter model is thus expected to remain a Dirac semi-metal up to some finite coupling strength U c . At half a magnetic flux quantum per plaquette, the model also is known as the π-flux model [38], and has been studied extensively numerically. As interactions are turned on, the π-flux model exhibits a quantum phase transition of the chiral Heisenberg Gross-Neveu universality class at U c ≈ 5.6t into an antiferromagnetic Mott insulator (AFMI) [39][40][41][42][43]. Since the π-flux model corresponds to the Hubbard-Hofstadter Hamiltonian threaded with maximum possible flux, we may think of the zerofield Hubbard model on a half-filled square lattice as the "0-flux model", which exibits a metal−AFMI transition with U c = 0 [29,30]. Our simulations address intermediate field strengths between the 0-flux and π-flux Hubbard model, which, to the best of our knowledge, has not been studied via DQMC.
In this work, we study the Hubbard-Hofstadter model using DQMC and present the evolution of thermodynamic properties of correlated electrons in an orbital magnetic field B. We demonstrate that the gross structure of magnetic Bloch bands and band gaps determined by the non-interacting Hofstadter spectrum is preserved in the presence of U . Moreover, we determine that the many-body fermion sign is directly connected to electronic charge compressibility. Finally, focusing on the half-filled AFMI, we find that an orbital magnetic field tends to delocalize electrons and thus effectively lower the influence of Hubbard U .

II. METHODS
We study the single-band Hubbard-Hofstadter model on a two dimensional square lattice where t is the hopping integral between the nearest neighbor sites ij , µ is chemical potential, and U is the on-site Coulomb interaction strength. c † iσ (c iσ ) is the creation (annihilation) operator for an electron on site i with spin σ =↑, ↓ and n iσ = c † iσ c iσ measures the number of electrons of spin σ on site i. As this model only has nearestneighbor hopping, it preserves particle-hole symmetry at half-filling with µ = 0. A uniform, orbital magnetic field is introduced by the Peierls substitution via the phase where the integral is taken over the shortest straight line path, Φ 0 = h/e is the magnetic flux quantum, and R i is the position of site i. We choose the symmetric gauge A = (−yx + xŷ)B/2 and do not include any Zeeman coupling terms.
We simulate the Hamiltonian in Eq. (1) on a finite cluster with lattice constant a = 1, and N x and N y sites in the x and y directions, respectively. N = N x N y denotes the total number of sites. We implement modified periodic boundary conditions consistent with magnetic translation symmetry [44]. Requiring that the wave function be singlevalued on the torus gives the flux quantization condition Φ/Φ 0 = n f /N , where Φ = Ba 2 is the flux through a plaquette and n f is an integer.
Allowing the hopping integral to carry a complex phase requires us to modify the standard DQMC algorithm to use complex numbers, which increases the run-time of our algorithm ∼ 3-fold. The complexified DQMC algorithm retains the same O(M 3 L) scaling as the real DQMC algorithm, where M = N x = N y is the linear size of the lattice, and L is the number of imaginary time discretization steps. Unless otherwise specified, all Hubbard-Hofstadter DQMC simulations are performed on a N x = N y = 8 square cluster. Error bars in DQMC results, when shown, denote ±1 standard error of the mean, estimated by jackknife resampling. Detailed simulation parameters are listed in Supplementary Note 1. Non-interacting results, where shown, are obtained from diagonalizing the Hofstadter model on a 40 × 40 cluster in order to minimize finite size effects.

A. Interacting Gap Structure
In Fig. 1, we show the electron density n vs. chemical potential µ at different field strengths for U/t = 0 − 8. In Fig. 1(a) we observe an electron density plateau where there are energy gaps between magnetic Bloch bands, n = 2Bν/Φ 0 , ν ∈ Z. As U increases in Fig. 1(b)-(e), these plateaus are weakened and pushed outward in chemical potential, but inflections of the n vs. µ curves are still visible at the same density values, indicating that degeneracy of Landau levels is not modified by U . Since we are at relatively high temperature, only the most prominent band gaps (with Chern number C = ±1) n = 2B/Φ 0 remain visible at larger U values. Additionally, at half filling, as U increases, a Mott gap appears and widens for all values of magnetic field. In Fig. 1(d)-(e), for U/t ≥ 6, when the Mott gap is well-defined, it decreases monotonically as the magnetic field increases, consistent with previous exact diagonalization results [24]. The same trend can be seen in a "correlated Hofstadter butterfly" plot, as shown in Supplementary Fig. S1 and described in Supplementary Note 2. We will discuss later, and in more detail, the behavior of the Mott gap.
It is instructive to plot our data as Wannier diagrams [45], i.e. color intensity plots of charge compressibility χ = ∂ n /∂µ as a function of electron density n and magnetic field strength B. Charge compressibility, or thermodynamic density of states, is directly measurable in experiments [9,13]. In a non-interacting system, at Curves with the same color have the same magnetic field strength Φ/Φ0 across all panels.Each curve is plotted with an offset 2B/Φ0 in order to improve visibility of inflection points. Error bars, corresponding to ±1 standard error of the mean, estimated by jackknife resampling, are smaller than the size of data points.All plots have inverse temperature β = 4/t.
zero temperature, charge compressibility is equivalent to the single-particle density of states. We measure charge compressibility in DQMC simulations as where n i = n i↑ +n i↓ . In Fig. 2, we show Wannier diagrams for U/t = 0 − 6. For all values of U , we observe local minima of χ (indicating incompresssible states) along straight lines satisfying the Diophantine equation where r and s are integers, and n 0 = 2 is the electron density of the completely filled system. This is consistent with what we expect from the Hofstadter spectrum in the non-interacting system [45]. The most prominent incompressible state with r = 1, s = 0 remains clearly visible up to U/t = 8. Less prominent incompressible states with r = 2 and r = 3 persist to U/t = 4 and U/t = 2, respectively. These results show that the integer quantum Hall states for r ≤ 3, s = 0 have no weak coupling instabilities with respect to Hubbard repulsion; the r = 1, s = 0 state remains stable up to large U . At half filling, the vertical compressibility minima indicative of the Mott gap becomes visible for U/t 4. Thus, we argue that the gross structure of magnetic Bloch bands and band gaps determined by the non-interacting Hofstadter spectrum is preserved in the presence of U , with the Mott gap at half-filling when U/t 4 superimposed as an additional feature. Due to the sign problem, our DQMC simulations are restricted to relatively high temperature βt ≤ 5, so we cannot resolve conclusively how much U changes the fine structure of the Hofstadter spectrum. All plots have inverse temperature β = 5/t. Magnetic field strength is displayed as Φ/Φ0, and electron density is shown as n /n0, where n0 is the electron density of a completely filled system, which in our case is 2. The system is particlehole symmetric, so only the range n /n0 ∈ [0, 0.5] is shown. Dotted cyan lines indicate where we expect local minima of charge compressibility to occur from the non-interacting model.

B. Fermion Sign
An important quantity in QMC simulations of interacting fermions is the fermion sign. The Hubbard-Hofstadter fermion sign s plotted against electron density n at magnetic field strength Φ/Φ0 = 12/64, inverse temperature β = 6/t and Hubbard interaction U/t = 4. This system is particle-hole symmetric, so we only plot the density range n ∈ [0, 1]. Dashed lines indicate electron densities at which χ reaches local minina and s reaches local maxima. Error bars, corresponding to ±1 standard error of the mean, estimated by jackknife resampling, are smaller than the size of data points. model is sign-problem-free at half filling on a bipartite lattice. But the fermion sign problem [46] fundamentally prevents us from obtaining high quality simulation data at low temperatures and away from half filling. Thus, any insight into factors affecting the severity of the sign problem is valuable. Since the fermion sign problem is NP-hard [47], we do not expect a general solution to the fermion sign problem to exist. Nevertheless, as the sign problem is representation-dependent, it is possible to reduce or completely remove the sign problem for specific classes of non-generic Hamiltonians [48].
In this work, we find a correlation between the fermion sign and the charge compressibility. In Fig. 3, we show the fermion sign s and charge compressibility χ, both plotted against n , for one representative set of parameters. Local minima of charge compressibility in this interacting system exactly correspond to local maxima of the fermion sign. At these local maxima, the fermion sign may be an order of magnitude improved over its value at other electron densities and that of the standard zero-field Hubbard model. For an extended figure demonstrating that this correspondence is general across our parameter space and not a finite size artifact, see Supplementary  Fig. S2 and Supplementary Note 3. Our results may mean that although the Hubbard model in general suffers from a sign problem, it is possible to obtain good results when we are precisely located on an integer quantum Hall plateau. Since similar sign-compressibility corre-spondence has been reported [30,[49][50][51], it appears that the improvement of fermion sign in insulating phases is quite general, consistent with our intuition that fermionic statistics become less important in localized states. Our results also relate to recent work [52][53][54] suggesting that the fermion sign is not merely a coincidental barrier to accessing low-temperature physics, but may be reflective of intrinsic physics of model Hamiltonians.

C. Half-Filling
Finally, we focus on half filling, where we believe interesting interplay between Hofstadter physics and Hubbard physics occurs. In the absence of a magnetic field, the ground state of the half-filled Hubbard model is an AFMI at any nonzero value of U , i.e. (U c = 0) [29,30], due to perfect nesting of the Fermi surface and a logarithmically divergent single-particle density of states. In the limit of strong interactions U t, the half-filled Hubbard model maps to the Heisenberg model with antiferromagnetic nearest neighbor spin exchange energy J = 4t 2 /U [55].
In the presence of an orbital magnetic field, the noninteracting density of states at half filling is modified significantly. As can be seen in Fig. 2(a), the density of states/charge compressibility does not change monotonically with field, but instead shows prominent minima at Φ/Φ 0 = p/q, where p and q are co-prime and q is even, corresponding to a non-interacting ground state with q inequivalent Dirac cones [37]. The large-field limit corresponds to the π-flux model, in which a semimetal−AFMI transition occurs at U c ≈ 5.6t. As the orbital magnetic field significantly changes the non-interacting density of states at half filling, we expect that critical U c should exhibit B field dependence. It would be interesting to investigate if U c changes monotonically with B, or if it exhibits non-monotonicity commensurate with the oscillatory behavior of the density of states. We defer the mapping of this U -B phase diagram to future work.
For the remainder of this section, we focus on the parameter region U/t ∈ [6,10]. Here, the system is safely an AFMI at all field strengths. We examine the evolution of local magnetic moment m 2 z , antiferromagnetic structure factor S(π, π) and specific heat c v = ∂ E /∂T with magnetic field strength, and see that these thermodynamic quantities all consistently show that a strong orbital magnetic field tends to modify the AFMI by delocalizing electrons and thereby reducing the effect of U on the low-energy properties of the Mott insulating phase. For finite-size analysis of thermodynamic observables at half filling, see Supplementary Fig. S3 and Supplementary Note 3.
In Fig. 4, we show the temperature and field dependence of the local moment. The local moment or sublattice magnetization z at half filling. Insets display the field dependence of magnetic structure factor S(π, π) at βt = 16 for magnetic field strength Φ/Φ0 ∈ [0, 0.5]. (a)-(c) Correspond to Hubbard interaction strength U/t = 6 − 10, respectively. Curves with the same color and marker type have the same magnetic field strength across all panels. Error bars, corresponding to ±1 standard error of the mean, estimated by jackknife resampling, are smaller than the size of data points.
measures the degree of spin localization. It is 0.5 in the non-interacting system and approaches 1 in the U/t → ∞ limit. In the zero-field Hubbard model, m 2 z has features at T ∼ U associated with the formation of local moments, and at T ∼ J, associated with short-or longrange ordering of local moments [36,56]. We see that at fixed U , increasing magnetic field strength reduces the local moment monotonically at all temperatures, with the effect largest below temperatures T ∼ J. The zero-field and π-flux limit of local moment data are consistent with previous work [43].
The magnetic structure factor is the Fourier transform of the real-space spin-spin correlation function When the system has long-range antiferromagnetic order, the structure factor is strongly peaked at the ordering wave vector Q = (π, π), with peak height scaling linearly with lattice size [36,57]. In our simulations, we find that at all B and U values, the magnetic structure factor is sharply peaked at (π, π), consistent with the system being in the AFMI phase. Insets to Fig. 4 show that at all U , the magnetic field monotonically reduces S(π, π), indicating that the magnetic field reduces the AFMI ordering tendencies. Fig. 5 shows the evolution of the low temperature peak of c v with magnetic field and Hubbard U . We calculate c v numerically by measuring energy as a function of temperature E(T ) and taking the finite difference ∆ E /∆T . In the zero-field half-filled Hubbard model, at large U , the specific heat has a "two peak" structure, with a broad high temperature peak at T ∼ U associated with charge fluctuations and and a narrow low temperature peak at T ∼ J associated with spin fluctuations [56,58]. When a magnetic field is turned on, as shown in Fig. 5(b)-(c), U is large enough that the two peaks remain well-separated. The high temperature peak doesn't move, while the low temperature peak shifts to higher temperatures. In Fig. 5(a), the two peaks are initially well-separated at low fields, but at Φ/Φ 0 1/4, the low temperature peak shifts upwards and merges with the high-temperature peak, complicating our interpretation. This is likely due to U/t = 6 being low enough for the system to not simply map to the Heisenberg model, and for the system to be close to the AFMI phase transition at π-flux.
Since the low temperature peak in specific heat is associated with the spin exchange energy J, we are tempted to say that the orbital magnetic field increases J. However, this interpretation may be overly naive. We believe the more accurate statement is that the orbital magnetic field tends to delocalize electrons, and thus, effectively lower the influence of U on low energy properties of the system. Insets to Fig. 5 show that a magnetic field increases (in magnitude) kinetic energy in the insulating phase, which supports this interpretation. This runs contrary to our usual intuition that an orbital magnetic field localizes electrons by winding them up into Landau orbits. However, here our starting point is a correlated insulator, rather than free electrons (or a Fermi liquid). Our results in Figs. 4 and 5, along with the decreasing width of the Mott gap in Fig. 1(d)-(e), suggest that in the AFMI phase, the orbital magnetic field tends to delocalize electrons, increase kinetic energy, and lower the effective influence of U . We observe that the influence of magnetic field is suppressed as U increases. As U/t → ∞, the influence from the B field will diminish and become negligible, since in the atomic limit, no hopping exists, and the orbital magnetic field cannot have an influence on the system.

IV. CONCLUSIONS
In this work, we implemented DQMC to simulate the Hubbard-Hofstadter model and directly investigate field dependent thermodynamic properties of correlated electrons, specifically focusing on the charge compressibility, local moment, magnetic structure factor, and specific heat. By examining charge compressibility, we find that magnetic Bloch bands are smeared out by the Hubbard interaction, but the non-interacting band gaps away from half-filling persist at temperatures accessible to DQMC in the presence of local correlations. At half filling, we find that the orbital magnetic field and Hubbard potential act antagonistically. At intermediate to strong coupling U/t ∈ [6, 10], strong orbital magnetic fields reduce the apparent width of the Mott gap, reduce the magnitude of the local moment and magnetic structure factor, increase kinetic energy, and shift the low-T peak of specific heat to higher temperatures. Together, these phenomena indicate that an orbital magnetic field tends to delocalize electrons and reduce the effect of U . From the algorithmic perspective, we find that the fermion sign in DQMC simulations is improved significantly when the physical system is incompressible.
For future work, we are interested in mapping out the U -B phase diagram for U/t ∈ [0, 6] and in studying the evolution of U c with magnetic field strength. At smaller values of U , we can achieve lower temperatures in DQMC simulations, but finite size effects also become more significant. Pinning down the precise location of the phase transition will require simulations on much larger lattices, as well as careful finite size scaling analysis.
As this work is partially motivated by understanding the large variety of exotic states in Moiré materials [9][10][11][12][13][14], another potential direction for future work is to use methods developed here to study the Hubbard-Hofstadter model on the honeycomb or triangular lattice. It has been suggested that the single-band Hubbard model on a triangular lattice is directly applicable to transition metal dichalcogenide heterobilayers [59,60]. However, due to the linear dispersion relation of Dirac fermions, the Coulomb interaction is poorly screened in graphene. It is likely that some multi-orbital, extended Hubbard model may be required to capture the full effect of electron correlations, for example, in twisted bilayer graphene [61,62].

V. DATA AVAILABILITY
Aggregated numerical data and analysis routines required to reproduce the figures can be found at https: //doi.org/10.5281/zenodo.6383764. Raw simulation data that support the findings of this study are stored on the Sherlock cluster at Stanford University and are available from the corresponding author upon reasonable request.

VI. CODE AVAILABILITY
The most up-to-date version of our DQMC simulation code can be accessed at https://github.com/edwnh/ dqmc.

VII. AUTHOR CONTRIBUTIONS
EWH wrote the DQMC simulation code. JKD performed simulations and data analysis. TPD conceptual-ized the work. All authors (JKD, WOW, BM, YS, EWH, TPD) participated in discussions and manuscript writing.

VIII. COMPETING INTERESTS
The authors declare no competing interests. Half-filling data shown in main text Figs. 4 and 5 are obtained from simulations performed with 2 × 10 4 to 4×10 4 warm-up sweeps and 2×10 5 to 5×10 5 measurement sweeps. We run 20 to 200 independently seeded Markov chains for each set of parameters. At half filling, we enforce a smaller imaginary time discretization interval in order to reduce effects from Trotter error. Here the imaginary time discretization interval ∆τ ≤ 0.05/t, and the number of imaginary time slices L = β/∆τ ≥ 10. We check that ∆τ = 0.05/t is sufficiently small that Trotter error is negligible in thermodynamic observables and do not affect the conclusions presented in this paper.
In all simulations, multiple equal-time measurements are taken in each full measurement sweep through the auxillary field. Each Markov chain with M measurement sweeps collects M L/5 equal-time measurements. The mean and standard error of observables are estimated via jackknife resampling of independent Markov chains.

SUPPLEMENTARY NOTE 2: CORRELATED HOFSTADTER BUTTERFLY
This section shows the "correlated Hofstadter butterfly" plot, which is a color intensity plot of charge compressibility χ = ∂ n /∂µ as a function of chemical potential and magnetic field strength B. Charge compressibility is treated as a proxy for local density of states in the correlated system, since we don't have enough data to perform Maximum Entropy analytic continuation to obtain the "true" local density of states. In Fig. S1(c)-(d), we see the width of the apparent Mott gap shrinks as field strength increases. This is consistent with main text Fig.1(d) gap, determined by performing cubic spline interpolation on the χ(µ) data, and solving for the chemical potential at which χ = 0.05/t. Black dashed vertical lines mark the zero-field Mott gap boundary, demonstrating that the Mott gap shrinks as orbital magnetic field increases. All subplots have inverse temperature β = 4/t. Magnetic field strength is displayed as Φ/Φ0, and electron density is shown as n /n0, where n0 is the electron density of a completely filled system, which in our case is 2. The system is particle-hole symmetric, so only the range n /n0 ∈ [0, 0.5] is shown. Results are obtained on 10 × 10 lattice, with inverse temperature β = 4/t, and Hubbard interaction strengths U/t = 4 − 6.
This system is particle-hole symmetric, so we only plot the density range n ∈ [0, 1]. Curves with the same color have the same magnetic field strength Φ/Φ0 across all panels. Vertical dashed lines are color-coded according to magnetic field strength, and indicate electron densities at which the charge compressibility reaches its first local minina, and where the fermion sign may have a corresponding local maxima or inflection point. Magnetic field strength is displayed as Φ/Φ0. Error bars (±1 standard error of the mean, estimated by jackknife resampling) are smaller than the size of data points.

SUPPLEMENTARY NOTE 3: FINITE SIZE ANALYSIS
Since our DQMC simulations are limited to relatively small clusters, it's important to check if the conclusions of this paper correctly extrapolate to the thermodynamic limit. Extended data plot Fig. S2, obtained on a 10 × 10 lattice, shows that charge compressibility minima and fermion sign maxima occur at electron filling factors n = 2Bν/Φ 0 , ν ∈ Z, regardless of lattice size. The sign-compressibility correspondence is general across field strengths and Hubbard interactions and is not a finite-size artifact.
In Fig. S3, we show the lattice-size dependence of local moment and specific heat for the half-filled Hubbard-Hofstadter model with U/t = 6. The conclusions that 1) local moment decreases with field and 2) the lowtemperature specific heat peak moves to higher temperature are qualitatively true, independent of the lattice size used in DQMC simulations. In general, we expect finite size artifacts to be most severe at small U , when our sys-tem is least localized, with most extended wavefunctions. As U increases, we expect finite size effects to be reduced, since U localizes electrons, and the system becomes insensitive to boundary conditions. We also expect finite size effects to be most important at low temperatures, when the system is closest to its ground state. Thus, the finite-size effects shown in Fig. S3 is the worst case for the parameters reported in the main text at half filling, when the system is an AFMI. Even in this worst possible case, on all lattice sizes we simulate, the local moment and specific heat display field dependence trends as reported in the main text. Thus, we believe our finding that an orbital magnetic field delocalizes electrons and weakens the effect of U at half filling is valid in the thermodynamic limit. * jxding@stanford.edu † tpd@stanford.edu Solid, dashed, dot-dash, and dotted lines denote data obtained on 8 × 8, 10 × 10, 12 × 12, and 16 × 16 lattices, respectively. Magnetic field strength is displayed as Φ/Φ0. Error bars (±1 standard error of the mean, estimated by jackknife resampling) are smaller than the size of data points. On the 10 × 10 lattice, Φ/Φ0 = 1/8 and Φ/Φ0 = 3/8 can't be achieved exactly, so these values are represented by data obtained with field strengths Φ/Φ0 = 13/100 and Φ/Φ0 = 38/100, respectively.