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é 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. The Hofstadter–Hubbard model on 2D square lattices is a paradigmatic model to study the interplay of electron correlations and external magnetic field. The authors use quantum Monte Carlo to study the thermodynamic properties of the Hofstadter Hamiltonian at intermediate to strong coupling, finding that a strong orbital magnetic field delocalizes electrons and reduces the effective Hubbard interaction.

S trong 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 two-dimensional 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 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 evendenominator 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.

Results and discussion
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. 1a we observe an electron density plateau where there are energy gaps between magnetic Bloch bands, hni ¼ 2Bν=Φ 0 ; ν 2 Z. As U increases in Fig. 1b-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 Fig. 1 Electron density vs. chemical potential. Electron density 〈n〉 vs. chemical potential μ a in the non-interacting system, and b-e with Hubbard U/t = 2-8. 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. half filling, as U increases, a Mott gap appears and widens for all values of magnetic field. In Fig. 1d, 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 44 , 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 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 44 . 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.
Fermion sign. An important quantity in QMC simulations of interacting fermions is the fermion sign. The Hubbard-Hofstadter model is sign-problem-free at half filling on a bipartite lattice. But the fermion sign problem 45 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 46 , 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 47 .
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 Wannier diagrams a in the non-interacting system and b-d with Hubbard U/t = 2-6. Gray regions in b-d are parameter regions where we don't have simulation data. All plots have inverse temperature β = 5/t. Magnetic field strength is displayed as Φ/Φ 0 , and electron density is shown as 〈n〉/n 0 , where n 0 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〉/n 0 ∈ [0, 0.5] is shown. Dotted cyan lines indicate where we expect local minima of charge compressibility to occur from the non-interacting model. Fig. 3 Fermion sign-charge compressibility correspondence. a Charge compressibility χ = ∂〈n〉/∂μ and b 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.
value at other electron densities and that of the standard zerofield 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 signcompressibility correspondence has been reported 30,[48][49][50] , 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 [51][52][53] suggesting that the fermion sign is not merely a coincidental barrier to accessing lowtemperature physics, but may be reflective of intrinsic physics of model Hamiltonians.
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 54 .
In the presence of an orbital magnetic field, the non-interacting density of states at half filling is modified significantly. As can be seen in Fig. 2a, 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 hm 2 z i, 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 measures the degree of spin localization. It is 0.5 in the noninteracting system and approaches 1 in the U/t → ∞ limit. In the zero-field Hubbard model, hm 2 z i has features at T~U associated with the formation of local moments, and at T~J, associated with short-or long-range ordering of local moments 36,55 . 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 e iQÁðR i ÀR j Þ hðn i" À n i# Þðn j" À n j# Þi: ð4Þ 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,56 .
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. 4 Evolution of local moment and magnetic structure factor in the half-filled Hubbard-Hofstadter model. Temperature and field dependence of local moment hm 2 z i 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. Figure 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 55,57 . When a magnetic field is turned on, as shown in Fig. 5b, 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. 5a, 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-5, along with the decreasing width of the Mott gap in Fig. 1d, 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.

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-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 58,59 . 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 60,61 .

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 y iσ (c iσ ) is the creation (annihilation) operator for an electron on site i with spin σ = ↑, ↓ and n iσ ¼ c y iσ c iσ measures the number of electrons of spin σ on site i. As this model only has nearest-neighbor 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. (5) 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 62 . Requiring that the wave function be single-valued 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 runtime 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.

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.