The anisotropic Harper-Hofstadter-Mott model: competition between condensation and magnetic fields

We derive the reciprocal cluster mean-field method to study the strongly-interacting bosonic Harper-Hofstadter-Mott model. The system exhibits a rich phase diagram featuring band insulating, striped superfluid, and supersolid phases. Furthermore, for finite hopping anisotropy we observe gapless uncondensed liquid phases at integer fillings, which are analyzed by exact diagonalization. The liquid phases at fillings 1 and 3 exhibit the same band fillings as the fermionic integer quantum Hall effect, while the phase at filling 2 is CT-symmetric with zero charge response. We discuss how these phases become gapped on a quasi-one-dimensional cylinder, leading to a quantized Hall response, which we characterize by introducing a suitable measure for non-trivial many-body topological properties. Incompressible metastable states at fractional filling are also observed, indicating competing fractional quantum Hall phases. The combination of reciprocal cluster mean-field and exact diagonalization yields a promising method to analyze the properties of bosonic lattice systems with non-trivial unit cells in the thermodynamic limit.


I. INTRODUCTION
Since the discovery of the quantum Hall effect [1][2][3], the lattice geometry's influence on charged particles in magnetic fields has been the subject of extensive research. Prototypical models such as the non-interacting Harper-Hofstadter model (HHm) [4,5] exhibit fractionalization of the Bloch bands with non-trivial topology (see Fig.  1), manifesting in quantum (spin) Hall phases [6][7][8]. Ultracold atomic gases with artificial magnetic fields [9][10][11][12][13] enabled the experimental study of the non-interacting model [8,[14][15][16], while the effect of strong interactions on the band properties remains an open problem. While heating processes in the regime of strong interactions still represent a problem for cold atom experiments with artificial magnetic fields, recent experimental progress gives hope that this can be controlled in the near future [17,18].
For bosons in the HHm with local interaction, i.e. the Harper-Hofstadter-Mott model (HHMm), previous theoretical studies found fractional quantum Hall (fQH) phases, which have no counterpart in the continuum for strong fields [19], using exact diagonalization (ED) [19][20][21], composite fermion theory [19], and the density matrix renormalization group (DMRG) on a cylinder [22]. Composite fermion studies also found evidence of a bosonic integer quantum Hall phase in bands with Chern number two [23], also observed with ED [24] in the presence of next-neighbor hopping. In a recent DMRG study [22] a bosonic integer quantum Hall groundstate was also observed in the standard HHMm at filling ν = 2. However, the composite fermion approach is biased by the choice of the wavefunction [19,23], while ED on small finite systems suffers from strong finite-size effects [20,21,24]. where the arrows indicate the direction of the corresponding hopping processes. The 4 × 4 cluster employed in the RCMF approach is also shown (gray shaded area). (b) Single-particle dispersion for Φ = π/2 and tx = ty = 1. The precession of thê h k,q vector [Eq. (6)] is shown for three states (red, blue and green) when varying k. The vector-colors indicate the values of k (see colorbar).
The issue becomes especially challenging when going to strong fluxes, where extrapolation to the thermodynamic limit is impossible [25] and the overlap of the finite-size groundstate with Laughlin [20,21] or composite fermion [19] wavefunctions quickly decreases. DMRG, on the other hand, is restricted to cylinder geometries [22,26,27] and encounters huge numerical difficulties in the case of critical phases. Variational Gutzwiller mean field studies also found evidence of fQH phases [28][29][30], as well as striped vortex-lattice phases [28], but the variational basis is restricted by construction. The results of a recent cluster Gutzwiller mean field (CGMF) study [31] are likewise hard to interpret since the method breaks the translational invariance and the topology of the system.
To overcome these problems, we develop a reciprocal cluster mean field (RCMF) method, directly defined in the thermodynamic limit, which preserves the topology of the lattice, and yields excellent agreement with numerically exact results for the Bose-Hubbard model. Further, we introduce an observable for the measure of topological properties in the presence of interactions.
We systematically map out the phase diagram of the strongly interacting HHMm as a function of the chemical potential and the hopping anisotropy. The phase diagram features band insulating, striped superfluid, and supersolid phases. At integer fillings we further observe highly anisotropic gapless uncondensed liquid phases, which are analyzed using exact diagonalization. For fractional filling, we find incompressible metastable states, indicating competing fQH phases. We define the respective order parameters, and present spatially resolved density, condensate-density, and current patterns. Finally, we discuss how on an infinite cylinder with a single unitcell in the y-direction the liquid phases become gapped and show a quantized Hall response to the adiabatic insertion of a magnetic flux. This paper is organized as follows. In Sec. II the HHm and HHMm are discussed, and the method for measuring non-trivial topological properties is introduced. The RCMF method is derived and discussed in Sec. III, while the results for the HHMm are presented and discussed in Sec. IV. Finally, Sec. V is devoted to conclusions.

A. Harper-Hofstadter model
To facilitate the discussion for the strongly-interacting system, we first review the non-interacting HHm on the square lattice. The Hamiltonian is given by with hopping amplitudes t x/y and annihilation (creation) x,y . Each plaquette is pierced by a flux such that a phase Φ is picked up when going around it, as illustrated in Fig. 1a. For Φ = 2π/N Φ the unit cell can be chosen as N Φ sites in the y-direction.
Eq. (1) is diagonalized by the transform b l (k, q) = (2π) −1 x,j e −i(kx+q(l+jNΦ)) b x,l+jNΦ , where l ∈ [0, N Φ − 1] and k(q) are the momenta in x(y)direction. For even N Φ the Hamiltonian reduces to and For Φ = π/2, used below, the system has three isolated topologically non-trivial bands, see Fig 1b. Here we use the notation of Ref. 16 where the central (super)band contains twice the number of states as compared to the two other bands. For a discussion of how the hopping anisotropy t x /t y = 1 affects the bandstructure, see Appendix A.
The Hamiltonian H Φ and Eq. (2) can be rewritten in the compact notation where v k,q is a vector of scalars and h k,q is a vector of operators The operator h k,q fully determines the momentum dependence of the non-interacting system, and we can apply the concept of parallel transport [32]. The local Berry curvature at the point (k, q) is proportional to the rotation of the unit-vector under an infinitesimal momentum shift. In fact, ifĥ k,q shows a non-trivial winding under transport on a closed path through the Brillouin zone, the Berry-curvature cannot be continuously deformed to a trivial one and the system is topologically non-trivial. The Chern number of the nth band is given by the number and direction of closed loops ofĥ k,q , i.e.
where γ n is the solid angle subtended byĥ k,q when taking the expectation value with respect to the single-particle eigenstates of the nth band and sweeping the momenta through the Brillouin zone. This is shown in Fig. 1b. For the lowest band A 0 (k, 0) and A 1 (k, 0) are shown while B(k, 0) varies only slightly:ĥ k,q performs one anti-clockwise loop, corresponding to a Chern number of c 0 = −1. Equivalently, for the central bandĥ k,q performs a double clockwise loop (c 1 = 2), while the highest band again has c 2 = −1.
The connection between the winding ofĥ k,q and the Hall conductivity can be seen from the example of the integer quantum Hall effect, i.e. the lowest band being completely filled with non-interacting fermions. Adding a magnetic flux Φ y piercing a torus of size L x × L y in y-direction can be achieved by transforming the hopping amplitudes as t x → t x e iΦy/Lx for hopping processes in +x, while taking the complex-conjugate in the opposite direction. The effect of this transform on the Hamiltonian (5) amounts to v k,q → v k−Φy/Lx,q , which is manifested in a translation of the vector h k,q with respect to the case without flux at each momentum (k, q), i.e.
where |Ψ(Φ y ) is the many-body groundstate under the flux Φ y . The sole effect of Φ y is therefore a transform of the many-body groundstate such that at each momentum (k, q) Eq. (7) is fulfilled, resulting in a rotation ofĥ k,q . Inserting a flux of Φ y = 2πL x /4 yields the transform A 0 (k, q) → A 1 (k, q), A 1 (k, q) → −A 0 (k, q), and therefore n l (k, q) → n l+1 (k, q). Adding a magnetic flux of Φ y = 2πL x /4 is equivalent to translating the manybody groundstate by one site in the y-direction.
If the lowest band is completely filled, the total number of particles on the torus is L x L y /4. Therefore adiabatically inserting a flux of Φ y = 2πL x /4 results in L x L y /4 particles being translated by one site in y-direction, or equivalently a total number of L x /4 particles being transported once around the periodic boundary in the ydirection. Consequently, adiabatically inserting a flux of Φ y = 2π results in a quantized total transverse transport of a single particle around the periodic boundary.

B. Harper-Hofstadter-Mott model
We proceed with the study of the HHMm with interaction U , chemical potential µ, and magnetic flux Φ = π/2, in the hard-core limit U → ∞.
In contrast to the non-interacting case, for a finite interacting system the Berry curvature is defined with respect to boundary twisting angles [25], i.e., where A j (θ x , θ y ) = i Ψ (θ x , θ y ) |∂ θj |Ψ (θ x , θ y ) is the Berry connection, Ψ is the many-body groundstate, and θ x , and θ y are twisting angles of the boundary conditions in x-and y-direction, respectively (i.e. T x/y Ψ(θ x , θ y ) = e iθ x/y Ψ(θ x , θ y ), where T x/y is a translation by the system size L x/y in x-, and y-direction, respectively). The twisted boundary conditions can be implemented in the same way as the magnetic flux discussed in Sec. II A by transforming the hopping as t x → t x e iθx/Lx and t y → t y e iθy/Ly for hopping processes in +x, and +ŷ-direction, respectively, while taking the complexconjugate in the opposite directions. The interaction and chemical potential terms in Eq. (8) remain unchanged. The only effect of adding the twisting angles (θ x , θ y ) to the infinite system is, as in Sec. II A, v k,q → v k−θx/Lx,q−θy/Ly .
In other words, if T θx,θy is the momentum-space translation operator which transforms each momentum as k → k + θ x /L x , and q → q + θ y /L y , we have For the Berry-curvature The Berry curvature is therefore fully determined by the response of the periodic-boundary many-body groundstate Ψ(0, 0) to a translation in momentum.
If we define P h.c. as the projector onto the Hilbert space of hard-core bosons (where multiple occupancy in position space is forbidden), the interacting many-body Hamiltonian [Eq. (8)] can be written as with particle-number operator N . The full momentum dependence of the hard-core bosons is therefore contained in the term P h.c. h k,q P h.c. . Furthermore, for any hardcore boson many-body eigenstate Ψ we have P h.c. |Ψ = |Ψ . As in the non-interacting case therefore a nontrivial winding of Ψ(0, 0) h k,q Ψ(0, 0) in momentum space indicates a non-trivial topology of the many-body groundstate. It should be emphasized that this measure is different from summing over the individual singleparticle Chern numbers of the occupied bands, since no projection onto non-interacting bands is involved.
For a further discussion of the measurement of topological properties with theĥ k,q -vector, see Appendix D.

III. RECIPROCAL CLUSTER MEAN FIELD
The previously employed CGMF method [31,33] breaks translational invariance by applying the meanfield decoupling approximation only to the hopping-terms at the boundary of the cluster, while the hopping terms within the cluster are treated exactly. The simplest case where this can be observed is when the symmetrybreaking field is zero, reducing the lattice to a set of decoupled clusters with open boundaries. This violation of translational invariance breaks the symmetries of the dispersion and thereby its topological properties. In order to mitigate such artifacts we develop a mean-field decoupling based on the concept of momentum coarse-graining, introduced in the context of the dynamical cluster approximation [34].
We term this method as "reciprocal cluster mean field" (RCMF). It crucially preserves both the translational invariance and the topology of the system. For topologically trivial translationally-invariant systems it yields more accurate results than previous mean-field methods (see Appendix F). It is well-suited for cases where the underlying symmetries of the dispersion are indispensable to understand the physical properties, such as, e.g., topological insulators. For benchmarks of the method on the hard-core Bose-Hubbard model and the chiral ladder with artificial magnetic fields [35,36], see Appendix F.
To illustrate our procedure let us first start from a general non-interacting hopping Hamiltonian with hopping amplitudes t (x , y ),(x, y) in position-space and dispersion k,q in reciprocal space. As in the dynamical cluster approximation [34], the main idea of RCMF consists in projecting the N ×M lattice system onto a lattice of N c × M c clusters (later we will take N, M → ∞, but the method is also well-defined for finite systems). Each cluster is spanned by the internal cluster coordinates X and Y , such that we can decompose the position coordinates x and y on the lattice into wherex andỹ are inter-cluster coordinates. In the same way the momenta in x and y-direction -k and q, respec-tively -are decomposed as where K and Q are the cluster momenta in reciprocal space. Through a partial Fourier transform, the creation and anihilation operators in reciprocal space can be written in the mixed representation where b K,Q (x,ỹ) annihilates a boson with clustermomenta K and Q on the cluster located at (x,ỹ) [34]. The central idea of the momentum coarse-graining consists of projecting the dispersion of the lattice k,q onto the clusters in reciprocal space. This can be done by a partial Fourier transform of the dispersion onto the subspace of cluster-local hopping processes, giving the intra-cluster dispersion¯ K,Q as representing hopping processes within the cluster, while the remainder δ K,k,Q,q = K +k,Q +q −¯ K,Q represents all other hopping processes between different clusters [34]. Now we can decompose the Hamiltonian of Eq. (10) into where, using Eq. (11), the part H c is cluster-local, while ∆H contains the coupling between different clusters where in the second line we introduced the mixed representation of δ K,k,Q,q , Our goal is to derive an effective Hamiltonian which is cluster local through a mean-field decoupling approximation of ∆H. To this end we decompose the creation/annihilation operators into their static expectation values and fluctuations, i.e.
The standard procedure of the mean-field decoupling approximation consists of neglecting quadratic fluctuations. Furthermore, we assume translational invariance between the different clusters, i.e. that the condensate φ K,Q is independent of the cluster location As derived in detail in Appendix E, this approach reduces a general system with local interactions and Hamiltonian local systems with effective mean-field Hamiltonian, The symmetry breaking field F X, Y is given by and the effective hopping amplitudes are defined as Within the RCMF approach the effective free energy of the lattice system is given by where Ω is the free energy of the cluster local Hamiltonian of Eq. (16). Note that Eq. (19) is consistent with the standard lattice free-energy within the single-site meanfield approximation [37]. In fact, requiring stationarity in the symmetry breaking field F X, Y , taking into account Eq. (17), reproduces the standard mean-field self-consistency condition Here, . means taking the expectation value with respect to the mean-field Hamiltonian [Eq. (16)].
We note in passing that the treatment of the symmetry-breaking field F is identical to the way it should be implemented in a dynamical cluster approximation extension of bosonic dynamical mean-field theory [38][39][40].
Groundstate phase diagram of the HHMm in two dimensions with hard-core bosons and flux Φ = π/2 in terms of µ/tmax and (tx − ty) /tmax. The observed phases are band insulating (BI, light blue), supersolid (SS, dark blue), striped superfluid (SF, white), gapless uncondensed liquid (Liquid, pink), and fractional quantum Hall (fQH, dark gray). The dashed regions indicate where the RCMF groundstate has a non-zero condensate order parameter but is very close in energy (< 3%) to metastable uncondensed states. At zero anisotropy the striped superfluid undergoes phase separation between vertically (for tx > ty) and horizontally (for tx < ty) striped order (black vertical line), while for µ = 0 the density is homogeneous and fixed to n = 1/2 in all phases (green dashed line).

IV. RESULTS
In our RCMF approach the HHMm [Eq. (8)] is reduced to an effective 4×4 cluster Hamiltonian, for details see Appendix G. For a comparison of our RCMF results with exact diagonalization at zero hopping anisotropy, see Appendix C. In Fig. 2 we present the groundstate phase diagram in terms of the chemical potential µ/t max and the hopping-anisotropy (t The phases at densities n and 1 − n are related by a CT -transformation consisting of a particlehole transform combined with complex-conjugation (see Appendix B). The symmetry around the (t x − t y ) = 0 axis corresponds to gauge invariance, since t x and t y can be exchanged in combination with a lattice-rotation of π/2. At n = 0 and n = 1 we find topologically trivial band insulators (BI). Below we discuss the other resulting phases in more detail.

A. Condensed phases
At moderate values of µ we observe superfluid phases with striped density and condensate density modulation. For t x > t y this is a vertically striped superfluid (VS-SF), with vertically striped density distribution ρ(x, y) and condensate-density distribution ρ c (x, y) = |φ x,y | 2 , as shown in Fig. 3 together with the spatially resolved particle current J(x, y). The net current is zero, as expected for an infinite system. Locally, however, there are chiral currents around two plaquettes in the horizontal direction. We therefore introduce the striped-superfluid order parameter J str = For t x < t y the superfluid phase is horizontally striped (HS-SF), with the patterns of Fig. 3 rotated by π/2 compared to the VS-SF. Since at t x = t y the system is invariant under a π/2-rotation, for |µ|/t max 2 the VS-SF and HS-SF undergo phase separation.
At |µ|/t max 2 and for low anisotropy we find a supersolid phase (SS) with lower free energy than the striped phases. The density distributions ρ and ρ c spontaneously break translational invariance, having a period larger than the unit cell (see Fig. 3). A similar spontaneous breaking of translational invariance has already been observed in the staggered-flux bosonic Harper-Mott model [41] and the bosonic Hofstadter model on a dice lattice [42], and has recently been measured experimentally in spin-orbit coupled Bose-Einstein condensates [43]. The SS exhibits chiral currents around single plaquettes, with position-dependent amplitudes, as captured by the order parameter J ss = x,y cos π 2 (x + y) (J x (x, y) − J y (x, y)). In all phases, at µ = 0 the density distribution is homogeneous, ρ(x, y) = 1/2, while ρ c (x, y) remains modulated.
The phase transition between the striped superfluids and the SS phase is characterized by a kink in the average condensate density n c , see Figures 4a and 4b. For n c > 0, the striped superfluid order parameter J str is only zero at |t x − t y |/t max = 1 (where the lattice is a set of trivial one-dimensional chains), exhibiting a kink at the phase transition to the SS, where also J ss becomes non-zero.

B. Uncondensed phases
At density n = 1/2 ( Fig. 4a) and stronger anisotropy we find a phase with zero condensate density (n c = 0). In Figs. 4c and 4d we show a sweep in µ for (t x − t y )/t max = −0.8, where we observe plateaus with zero n c , zero current, and homogeneous density distribtution ρ(x, y) = ν/4, with fillings ν = 1, 2, 3. In these phases, since F x,y = 0, the RCMF Hamiltonian of Eq. (16) reduces to a finite 4 × 4 torus without any external variational parameter. In order to further analyze these phases we therefore turn to ED using twisted boundary conditions in order to analyze finite-size effects (see Sec. II B and Appendix C). If the phases are gapped, one expects the manybody gap to stay finite for all twisting angles (θ x , θ y ), while in gapless phases the groundstate mixes with excited states.
As can be seen in Fig. 5 for t x < t y the groundstate remains gapped with respect to boundary twisting in the x-direction with θ y = 0, while it mixes with the excited states for twisting in the y-direction. For t y < t x the behavior is reversed. This is also consistent with the correlations | b † x,y b 0,y | decreasing exponentially to zero as a  function of x for t x < t y , while staying finite throughout the system for t x > t y (see Appendix C). In contrast to the two-dimensional Bose-Hubbard model without magnetic flux, which in the superfluid groundstate always shows condensation as long as both hopping amplitudes are finite [44], the HHMm therefore shows a transition at finite t x /t y to a highly anisotropic uncondensed gapless liquid. The fact that these phases are adiabatically connected to the one-dimensional limit (t x = 0 or t y = 0), where hard-core bosons are in a superfluid phase, as well as the highly anisotropic correlations | b † x,y b 0,y |, possibly point to an unconventional one-dimensional superfluid order.
As a function of the hopping anisotropy, the liquid phases occur where the lowest band is particularily flat either in k-or q-direction, suppressing condensation in the minima of the dispersion (see Fig. 6a and Appendix A). While the system has zero current everywhere (due to the periodic boundaries), a signature of the response of the liquid to the magnetic field is found by analyzing current-current correlations (see Appendix G 2 for details), resulting in two counter-propagating currents which cancel each other, shown in Fig. 6b.
In Fig. 6c we show the projection of the groundstate onto the three non-interacting bands n 0 , n 1 , and n 2 , for the same parameters as in Figs. 4c and 4d. At ν = 1 the lowest band shows unit filling. As shown in Fig.  6a, this phase appears in the same regions of µ as the integer quantum Hall plateau of non-interacting spinless fermions (see Appendix A)). At ν = 3 the holes show unit filling in the lowest band, due to the CT transform (see Appendix B).
As can be seen in Fig. 6d, the vectorĥ k,q shows the same behavior as for the lowest non-interacting band in all three liquid phases (shown for ν = 2), in contrast to the trivial BI at n = 0, 1 and the one-dimensional superfluid at t x = 0 or t y = 0. For ν = 1, as in the case of non-interacting fermions discussed in Sec. II A, this winding indicates the transverse transport of a single particle if a magnetic flux of Φ y = 2π is inserted. At ν = 3, the transverse transport consists of a single hole. This is consistent with the band fillings in Fig. 6c and the CT transform discussed in Appendix B, i.e. the reversal of the Hall conductivity σ xy (ν = 3) = −σ xy (ν = 1). As these phases are gapless in the two-dimensional ther- modynamic limit, the quantization of the Hall conductivity is not topologically protected by edge modes and therefore sensitive to disorder, as is the case for metallic Fermi-liquid-like phases of hard-core bosons [24,45,46]. By contrast, in the case of a cylindrical geometry, i.e. L y = 4 and L x → ∞, the response to the twisted boundaries in x-direction while θ y = 0, shown in Figs. 5a and 5c, indicates that all three plateaus are gapped. As a function of θ x the vectorĥ k,q shows a complete loop and appears to be robust against local perturbations (see Appendix D). What the nature of the phases at ν = 1, 3 is in such a quasi-one-dimensional setup remains to be investigated: the non-trivial winding indicates a gapped phase with odd Hall conductivity, which is expected to show intrinsic topological order and fractional quasiparticle excitations for bosons in two dimensions [47,48]. The non-degenerate groundstate we observe (which for bosons is only expected at even Hall conductivities) apparently is at odds with this prediction. However, the argument of Ref. 47 relies on the fact that the quasiparticle excitations need to behave as fermions under exchange in such an odd Hall conductivity phase. In the cylinder, however, we approach the one-dimensional limit, where hard-core bosons naturally behave as free fermions also in the absence of fractionalization. This possibly explains why we observe an ED groundstate which remains gapped and non-degenerate for all accessible system sizes L x × 4 (see Appendix C).
At ν = 2 (n = 1/2), the Hamiltonian of Eq. (8) is CT -symmetric. It directly follows that σ xy (ν = 2) = 0. This is consistent with the bands being equally filled with particles and holes (see Fig. 6b), resulting in a zero net Hall conductivity (see Appendix B). In two-dimensional systems such a CT -symmetric phase is expected to be topologically trivial [49][50][51], in line with the gaplessness observed in Figs. 5c and 5d.
On a cylinder, however, this phase is gapped as the one-dimensional limit is approached, where CTsymmetric phases can have a non-zero topological Z invariant for non-interacting fermions [49]. The nontrivial winding implies the quantized transport of a single particle-hole pair under the adiabatic insertion of a flux of 2π, resulting in a total zero Hall conductivity. Wether this particle-hole transport is a consequence of topologically protected edge modes on the cylinder remains to be investigated on larger system sizes, possibly using DMRG [22,26,27].
Whereas away from integer fillings the groundstate is always symmetry-broken, it is always possible within RCMF to find (metastable) stationary solutions with zero symmetry-breaking field (F = 0) and therefore n c = 0, as shown in Fig. 4d. While at large hopping anisotropy fractional fillings are largely supressed, at low anisotropy the F = 0 solution shows plateaus at any filling commensurate with the 4 × 4 cluster, i.e. ν = m/4 with integer m, as shown in Fig. 7c.
As mean-field approaches such as RCMF tend to overestimate the stability of symmetry-broken phases, we critically analyze the difference in free energy between the groundstate and the metastable plateaus in Fig. 7. As can be seen, this quantity shows local minima at integer (ν = 1, 2, 3) and half-integer (ν = 1/2, 3/2, 5/2, 7/2) fillings, while quarter fillings correspond to local maxima. This is consistent with the argument that without longrange interactions it costs a negligible energy to compress the ν = 1/4 to the ν = 1/2 Laughlin liquid [21]. The energy difference is particularily low in the vicinity of the liquid phases indicating that these plateaus might extend to lower values of hopping anisotropy. Furthermore, at low anisotropy (where the lowest band is particlularily flat, see Fig. 6c and Appendix A) the metastable plateau at ν = 1/2 (and ν = 7/2) is very close in free energy to the groundstate. This plateau has been shown to correspond to a fQH phase in ED [19][20][21], variational Gutzwiller mean field [30], and DMRG [22] studies. As shown in Fig. 7c, at zero anisotropy, the condensate fraction of the groundstate shows a local minimum at both ν = 1/2 and ν = 2, further indicating that it might converge to zero with increasing cluster size. These are also the fillings where at zero anisotropy there is the largest discrepancy between RCMF and ED results on small finite systems (see Appendix C). It should however be noted that when assuming a cylinder geometry, RCMF observes both a ν = 1/2 and ν = 2 plateau, in agreement with ED. To conclude, there are regions of the phase diagram, where the symmetry-broken groundstate and the metastable plateaus are too close in free energy to dismiss finite size effects. We denote these regions (identified by the condition |Ω GS − Ω F =0 | < 3% of the groundstate energy) as dashed areas in the phase diagram of Fig. 2.

V. CONCLUSION
We derived the reciprocal cluster mean field method and applied it to the groundstate phase diagram of hardcore bosons in the Harper-Hofstadter-Mott model at flux Φ = π/2. The bosons exhibit band insulating, striped superfluid, and supersolid phases. At finite anisotropy and integer filling we further found anisotropic gapless uncondensed liquid groundstates characterized by a non-trivial winding of the newly introduced vectorĥ k,q . We further analyzed the properties of these phases using exact diagonalization. At fillings ν = 1 (3) this corresponds to integer particle (hole) filling of the lowest band, while the ν = 2 phase is CT symmetric with zero Hall response. We also observed metastable fractional quantum Hall phases predicted by other methods [19][20][21][22]30], which do not correspond to the groundstate (most likely due to finitesize effects), but are very close in free energy. Finally, we discussed how the liquid phases at integer fillings become gapped on a cylinder with just one unit-cell in the y-direction and show a quantized Hall response to the adiabatic insertion of a magnetic flux. These properties, which are not expected for the full two-dimensional system, seem inherent to the quasi-one-dimensional nature of the cylinder geometry and need to be further investigated on larger system sizes. The combination of reciprocal cluster mean-field and exact diagonalization provides a promising venue for the numerical simulation of bosonic lattice systems with larger unit cells in the thermodynamic limit. The HHm can be solved by diagonalizing the Hamiltonian of Eq. (2), yielding three topologically non-trivial bands (see Fig 1b). For the gauge used in this work, the non-trivial topology arises in k-direction, while in qdirection the dispersion has a trivial cosine-shape. Both the topology and the four minima of the dispersion are independent of the anisotropy between the hopping amplitudes t x and t y . The bandwidths of the three bands, on the other hand, are affected by the ratio t x /t y .
In order to analyze this, we introduce the quantities ∆E k and ∆E q for the lowest band, where ∆E k is the bandwidth in k-direction, i.e. q) is the dispersion of the lowest band, and max k/q corresponds to taking the maximum with respect to k and q, respectively. The bandwidth in q-direction, ∆E q is defined analogously. As shown in the upper panel of Fig. 6c, for t x t y the lowest band is particularily flat in the k-direction (∆E k 1), while for t y t x it is particularily flat in the q-direction (∆E q 1). Note that this is not to be confused with the "flatness" of the bands that typically supports fractional quantum Hall effects, which would consist in max [∆E k , fact this quantity is low in the region of low anisotropy). Instead, having only ∆E k 1 or ∆E q 1 will result simply in supressing the condensation of bosons in the minima of the dispersion.
Another quantity affected by the anisotropy is the gap between the lowest and the central band. The simplest many-body problem where this plays a role is the case of spinless non-interacting fermions, which exhibit an integer quantum Hall phase for integer filling of the lowest band, i.e. if the chemical potential µ lies within the (anisotropy-dependent) gap, see Fig. 5c.

Appendix B: Charge conjugation relations of hard-core bosons
For hard-core bosons a particle-hole transform (i.e. simultaneous b † → b and b → b † ) is equivalent to an inversion of the flux in the Hamiltonian of Eq. (8), i.e. Φ → −Φ. A direct consequence of this, is that the manybody groundstates at densities n and 1 − n are related by the CT operation, where C is the particle-hole transform, and T the complex conjugation operator. This implies that the Hall conductivity σ xy is anti-symmetric under the transform n → 1 − n [52], i.e. σ xy (n) = −σ xy (1 − n). (B1) This effect is known as the charge conjugation symmetry of hard-core bosons [52]. Further, the uncondensed phase at ν = 2 (n = 1/2) discussed in Sec. IV B is by definition CT -symmetric with zero Hall conductivity. This implies that in the chiral current patterns of Fig. 6b, the "charge" transport of the "particle" and "hole" channels will always cancel each other.

Appendix C: Comparison with exact diagonalisation
In Fig. 8a we compare our RCMF results with ED results on finite systems. The 4 × 4 ED system differs from the F = 0 solution of RCMF by a renormalization of the hopping according to Eq. (G1), resulting in a shift in chemical potential of the plateaus. We present a sweep of the density in chemical potential without hopping anisotropy (i.e. t x = t y ). As can be seen the only regions where we see a large discrepancy with respect to ED are around fillings ν = 1/2 and ν = 2. These are the fillings where the metastable plateaus are particularly close in energy to the symmetry-broken groundstate (see Fig. 7).
We further compare the ED results with RCMF results on a cylinder with just 4 sites and periodic boundary conditions in y-direction. This can easily been done by modyfying the coarse-graining procedure of Eq. (12), which is now only integrated over k. This results in a new cluster-hoppinḡ As can be seen in Fig. 8a, in this case also RCMF shows a fractional plateau at ν = 1/2 and a plateau at ν = 2, indicating that at zero anisotropy these phases are much more robust in the cylinder geometry than they are on the infinite square lattice. Apart from the response to twisted boundaries discussed in Sec. IV B, the anisotropic gapless nature of the two-dimensional uncondensed phases can also be observed in the scaling of the many-body gap as a function of L x while keeping L y = 4 fixed. As shown in Fig.  8b the manybody gaps remain essentially constant if t x is (sufficiently) smaller than t y , while it decreases in a non-monotonous way if t x is larger than t y . If the same scaling is done in y-direction the situation is reversed. This also implies that on the cylinder these phases are gapped for t x < t y . The same behavior can also be observed in the correlations | b † x,y b 0,y | in a system with L y = 4 and L x = 8, shown in Fig. 8c, which quickly drop to zero as a function of x for t x < t y , indicating a gapped phase on the cylinder. For t x > t y on the other hand, the correlations stay finite throughout the whole system hinting at the anisotropic gapless nature of the phase in two dimensions (and on the cylinder if t x is sufficiently large). the absence of U (1) symmetry-breaking. Computingĥ k,q in the phases with F = 0 by taking expectation values for the discrete momentum values of the cluster (K and Q) then is equivalent to taking the same expectation values with respect to the infinite lattice. By looking at the values ofĥ K,Q at these discrete momenta and extrapolating its rotation on the infinite lattice, we are thereby able to measure the topology of the infinite lattice in a way that is not limited by finite-size effects. This is shown in Fig. 9a, where theĥ k,q vector is compared in a 4 × 4 and 8 × 4 periodic system, respectively, for filling ν = 2 and (t x − t y )/t max = −0.5, yielding excellent agreement. In Figs. 9b and 9c we show the precession ofĥ k,q for filling ν = 1/2 on a 4×4 system, and ν = 3/5 on a 5×4 system. At ν = 1/2 the system is in a fQH phase showing a topological winding. At ν = 3/5, where for bosons no fQH phase is possible, the vector does not show any closed loop and has a net geometric angle of zero.
By using twisted boundaries in x-direction (i.e. varying θ x ), while keeping θ y = 0, we measure the response of the vectorĥ k,q on a cylinder to a magnetic flux piercing the system in y-direction in Fig. 9d. As discussed in Sec. II A, under the insertion of a flux of Φ y = θ x = 2π, the winding of the vector indicates an adiabatic translation of the manybody groundstate by one site in y-direction. For ν = 1 (shown in Fig. 9d) this translates into a quantized transverse transport of a single particle around the periodic boundaries in y-direction, while for ν = 3 a sin-gle hole is being transported. In the case of ν = 2 the total charge transport is zero, as a particle-hole pair is transported.
We further analyze the stability of the winding number against local perturbations. To that end we introduce a local shift in chemical potential ∆, such that the chemical potential is shifted on site (X, Y ) = (0, 0) on the cluster, i.e.
In Fig. 10 we show the response of the uncondensed phases on a 4 × 4 system at ν = 2 to this local perturbation. As can be seen in Fig. 10a, while the individual momentum values of the h k,q -vector change with the strength of the perturbation ∆, the total winding of the vector around the origin remains stable even at such large values as ∆ ≈ t max , while the single-particle gap (i.e. the size of the plateau) decreases (see Figs. 10b and 10c). At larger values of ∆ the phase has a non-zero condensate order parameter. A similar behavior can also be observed in the ν = 1, 3 phases. In the full two-dimensional system we expect this robustness to vanish as the system-size is increased and the manybody gap goes to zero. In the cylindrical geometry however this points to a stability of the winding against disorder, indicative of topological protection.
which, after a Fourier transform to position space, and dropping the (x,ỹ)-notation, yields the effective meanfield Hamiltonian where the symmetry breaking field F X, Y is given by and If instead of a pure hopping Hamiltonian, the Hamiltonian also includes local (interaction) terms, e.g.
the local part H int is already inherently cluster-local and can be absorbed into H c in Eq. (13), such that the effective Hamiltonian becomes Taking into account the constant shift of Eq. (E3), the free-energy of the full lattice system under the mean-field decoupling approximation can now be expressed as where Ω is the free energy of the cluster with the Hamiltonian of Eq. (E6). With Eqs. (E4)-(E7) we have everything in place in order to formulate the full RCMF approach, see Eqs. (16)- (20) in Sec. III.

Appendix F: Benchmarking RCMF
In order to benchmark RCMF we turn to the Bose-Hubbard model with hard-core bosons on a twodimensional square lattice using a 4 × 4 cluster Hamiltonian. In Fig. 11a we show RCMF results for the condensate density ρ c = X,Y |φ X, Y | 2 as a function of chemical potential for t x = t y = 1 and compare with standard single-site mean field, CGMF [33] on a 4 × 4 cluster, and numerically exact path integral quantum Monte Carlo (QMC) [44,53] results. As expected, RCMF shows better agreement with QMC than the two other mean-field methods. In contrast to CGMF, which due to the breaking of translational invariance converges towards a weakly position-dependent (unphysical) condensate φ X, Y , the condensate in RCMF is completely homogeneous.
We also compare RCMF results with QMC for anisotropic systems in Fig. 11b, observing stronger deviations with increasing anisotropy |t x −t y |. This is related to the use of a square symmetric 4 × 4 cluster, while the bandwidths in k-and q-direction are no longer equal. As the one-dimensional limit (t x = 0) is approached, meanfield methods are always expected to behave worse, since quantum fluctuations play a bigger role. However, the results are still qualitatively correct, and we conclude that RCMF works reasonably well also for anisotropic systems.
In order to ensure that RCMF can properly treat artificial gauge fields, we simulate the two-leg ladder of Refs. 35 and 36 with a magnetic flux of Φ = π/2 per plaquette and hard-core bosons using a 2 × 8 cluster. This ladder corresponds to the Harper-Hofstadter-Mott model where the x-direction is restricted to just two sites. It shows Mott phases at density n = 0.5 and superfluid phases otherwise, with both phases exhibiting Meissner and vortex current-patterns depending on the anisotropy [36]. The Meissner phases can be found for anisotropies where for the gauge of Ref. 35 the non-interacting groundstate momenta -i.e. the momenta where the dispersion has (degenerate) global minima -are k gs = ±π/4. These momenta are fully captured by the 2 × 8 cluster with cluster-momenta K = mπ/4, where m = 0, 1, 2, ..., 7. On the other hand, in the anisotropy-region where the vortex phases appear, k gs varies as a function of the hoppinganisotropy [35] and can no longer be represented within a 2 × 8 cluster. This is shown in Fig. 11c, where the RCMF chiral current is compared to DMRG results [36] both in the Mott (n = 0.5) and superfluid (n = 0.25) regime. Here, J y (l, y) is the current in y-direction on the yth site of the ladder-leg l. The RCMF results agree very well in the Meissner phases, while they cannot capture the vortex phases. This is a good example of what RCMF can do and what not: for RCMF to work it is indispensable that the cluster is both an integer multiple of the unit cell and that the groundstate momenta of the non-interacting model can be reproduced exactly by the grid of cluster momenta spanned by K and Q. This is a direct consequence of Eq. (E2), which in reciprocal space implies φ K+k,Q+k = φ K,Q δk ,0 δq ,0 . This condition, which was used to decouple the clusters in RCMF, can only be fulfilled if the momenta of the minima of the dispersion (which are the momenta where condensation will occur for local single-site interactions) can be reproduced by the momenta of the cluster (K, Q). If this is the case, as seen in Fig. 11c, the deviation from the DMRG results on the chiral current [36] is below 1%.

Appendix G: RCMF approach for the Harper-Hofstadter-Mott model
The HHm has groundstate momenta k gs = 0, ±π/2, π and q gs = 0. Since the momenta of the groundstate are independent of the anisotropy, we do not encounter the difficulties described in Appendix F for the vortex phases of the chiral ladder when using finite clusters. In order to fulfill Eq. (E2) by reproducing the minima of the dispersion (see Appendix F) a multiple of 4 sites in X direction is needed, since for 4 sites K is a multiple of π/2. We also need a multiple of 4 sites in Y direction in order to fully capture the 1 × 4 unit cell. In this work we restrict ourselves to the minimal cluster, i.e. 4 × 4.
Since the mean-field decoupling is performed in the thermodynamic limit, the sum overk andq in (12) can be replaced by an integral and computed analytically. In this configuration the coarse-graining described in Sec. III leads to the cluster-hoppinḡ where t (X , Y ),(X, Y ) is the original hopping of the HHm [Eq. (1)] with periodic boundary conditions, which plugged into (16)(17)(18) yields the effective RCMF Hamiltonian for the Harper-Hofstadter-Mott model on a 4 × 4 cluster.

Observables
The free energy Ω of Eq. (19) represents the free energy of the lattice system in the thermodynamic limit within the RCMF approximation. Using functional derivatives of Eq. (19) we can compute expectation values with respect to the full lattice system. According to the selfconsistency condition [Eq. (20)], this is trivial for the condensate φ X, Y = b X, Y .
Accordingly, we get for the condensate density ρ c (X, Y ) = φ X, Y 2 , and the total condensate density per site n c = 1 N c M c X,Y ρ c (X, Y ).
Also for the particle density we get an equivalence between the full lattice and the 4 × 4 cluster, since ρ(X, Y ) = − δΩ δµ X,Y = n X, Y , and, accordingly, for the total particle density per site n = 1 N c M c X,Y ρ(X, Y ).

Current-current correlations
In the uncondensed phases described in Sec. IV the current within the system is zero. However, as the system is not band insulating and the kinetic energy is non-zero this must result from two counter-propagating modes which cancel each other out. We can analyze these modes by measuring current-current correlations between neighboring bonds. By Eq. (G3) we can write the currents on the full lattice J x/y as expectation values of lattice-operatorsĴ x/y with respect to the RCMF Hamiltonian, i.e. J x/y (x, y) = Ĵ x/y (x, y) withĴ x/y (x, y) defined according to Eqs. (G2)-(G4). We now can look at the following current-current correlations: ∆J x,x (x, y) = Ĵ x (x, y)Ĵ x (x + 1, y) ∆J x,y (x, y) = Ĵ x (x, y)Ĵ y (x, y) ∆J x,y−1 (x, y) = Ĵ x (x, y)Ĵ y (x, y − 1) These three quantities are enough to describe the current patterns depicted in Fig. 6d: if ∆J x,x (x, y) is positive the currents J x (x, y) and J x (x+1, y) point in the same direction, if it is negative they point in opposite directions. The same is also true for ∆J x,y (x, y) and ∆J x,y−1 (x, y). Assuming a finite current in +x direction on a given site (x, y) and extracting the sign of the currents on the neighboring bonds through the correlations introduced above, one can therefore easily draw one of the two counterpropagating current patterns. The other pattern results from simply inverting the direction of all currents.