Fractional quantum Hall states with variational Projected Entangled-Pair States: a study of the bosonic Harper-Hofstadter model

An important class of model Hamiltonians for investigation of topological phases of matter consists of mobile, interacting particles on a lattice subject to a semi-classical gauge field, as exemplified by the bosonic Harper-Hofstadter model. A unique method for investigations of two-dimensional quantum systems are the infinite projected-entangled pair states (iPEPS), as they avoid spurious finite size effects that can alter the phase structure. However, due to no-go theorems in related cases this was often conjectured to be impossible in the past. In this letter, we show that upon variational optimization the infinite projected-entangled pair states can be used to this end, by identifying fractional Hall states in the bosonic Harper-Hofstadter model. The obtained states are characterized by showing exponential decay of bulk correlations, as dictated by a bulk gap, as well as chiral edge modes via the entanglement spectrum.

Introduction.−The physics of interacting topological phases of matter, which includes phenomena like fractionally charged quasiparticles, exotic exchange statistics and long-range (topological) entanglement has led to the formulation of many new concepts in physics, most notably topological order [1].Additionally, using these phenomena offers a promising route for robust quantum information technology [2,3].The paradigmatic family of quantum states for these phases of matter are the fractional quantum Hall states [4,5] and their lattice generalization as fractional Chern insulators [6].Given this, an enormous amount of theoretical effort, driven by the fundamental interest in the fractional quantum Hall states has been put forward since their discovery and continues to this day.Moreover, within the last years advances, in capabilities of experiments on cold atoms, like the possibility to include artificial gauge fields have put the extremely rich physics of the fractional Hall states into the grasp of current technologies [7][8][9][10][11].Recently, a crucial first experimental step was made with the realization of a Laughlin state of two atoms in such an experiment [12].A major part of the quest for the realization of the fractional quantum Hall physics in well controlled synthetic quantum matter experiments is the identification of specific Hamiltonians where these phases arise.As these states generically emerge from an intricate interplay of magnetic fields and interactions between the particles, the candidate Hamiltonians are challenging for state of the art numerical methods to treat.Within the class of tensor network methods for such numerical studies, the infinite Projected Entangled-Pair State (iPEPS) [13] framework has unique advantages for this task.As it * corresponding author: weerda@thp.uni-koeln.deoperates directly in the thermodynamic limit in two dimensions it can eliminate finite size effects that can significantly bias other related methods, e.g.cylinder-MPS, leading to possibly altered phase boundaries or even artificial phases that are non-existent in the two dimensional limit.Thus, the iPEPS ansatz constitutes an important tool for the determination of the phase structure of manybody Hamiltonians.However, doubts had been raised about the general applicability of the PEPS ansatz for the numerical simulation of gapped chiral states because of a no-go theorem for the case for free fermions [14].In a recent work this issue was studied more closely from a numerical perspective [15].It was shown that for a truncated parent Hamiltonian of a chiral spin liquid [16,17] the iPEPS ansatz was indeed successful in finding and representing this ground state numerically.In this work, we take a further step in utility of the iPEPS framework, by treating for the first time mobile bosons in an external gauge field, showing that chiral gapped topological phases arise.We are able to treat flux values that are in the range of today's experimental technology.This setting necessitates the use of large unit cells in the iPEPS ansatz due to the magnetic unit cell over which the Hamiltonian terms vary due to the external gauge field.Crucial for this investigation is the fact that we make use of a fully variational ground state search, as attempts with simple imaginary time evolution turn out to yield inconclusive results.By scanning through different parameter values of the Harper-Hofstadter Hamiltonian we are able to find incompressible plateaus at the filling fraction expected for bosonic fractional Hall states.We focus on the paradigmic case of the Laughlin state in the case of hard-core bosons, but show that also for soft-core bosons and higher fillings we can find incompressible plateaus corresponding to more exotic fractional Hall states.The chiral nature of these states is verified with the use of the entanglement spectrum.We thus for the first time represent bosonic fractional quantum Hall states with PEPS.Model and Methods.− We study a two dimensional system of bosonic particles on a lattice, which are interacting on-site and are subject to a magnetic flux ϕ through each plaquette of the lattice.This situation is described by a interacting Harper-Hofstadter Hamiltonian in which the hopping amplitude t is dressed by the semiclassical gauge field.The gauge field is chosen such that it produces the desired magnetic flux per plaquette, ϕ = □ A ij .This Hamiltonian can be realized in cold atom experiments [8,9,11], and was the basis for the recent preparation of a two-body fractional quantum Hall state [12].In the limit of large fillings it also describes the physics of Josephson Junction Arrays [18].For low densities and small values of ϕ, where the magnetic length ℓ B = 1/ √ ϕ is much larger than the lattice spacing, we approach the continuum setting, in which the presence of bosonic fractional Hall states is well established [19].However, at large flux values at which current cold atom experiments can operate and lattice effects play a central role, presence of truly lattice versions of the fractional quantum Hall states in the phase diagram remains a question of ongoing research.Studies based on exact diagonalization [20][21][22][23][24][25][26] as well as (infinite) density matrix renormalization group ((i)DMRG) [27][28][29][30] and tree-tensor networks (TTN) [31,32] have suggested a multitude of fractional quantum Hall states in such systems.Even though these methods provide inside into the possible phases of the Hamiltonian in the two dimensional thermodynamic limit, it is widely acknowledged that they suffer from strong finite size effects.In particular for elongated cylinder geometries, as usually employed by the (i)DMRG methods of treating two dimensional Hamiltonians, these strong finite size effects lead to an overestimation of the stability of gapped phases as well as the emergence of charge density wave orders, breaking translation invariance, as discussed in [29,30].Generically, it is not clear what the influence of finite size effects is on all the possible low-energy phases close to the true ground state is.Thus if these phases are in close competition, the phase with energetically most favorable finite size effect can be accidentally identified as the ground state.The circumstances where these considerations can become crucial include quite important models like the doped Hubbard model or models relevant for the search for spin liquids in materials just to name a few.Here, we thus aim to provide important clarification on the situation by showing that the iPEPS ansatz, using variational optimization, is capable of tackling such problems directly in the thermodynamic limit, and hence without spurious finite size effects.The iPEPS ansatz [13,33] for the many-body wavefunction used in this investigation is parameterized by a set of tensors (A 1,1 , . . .A m,n ) that form a m × n unit cell (illustrated with the grey line) which is periodically repeated on the infinite lattice: The individual tensors A i,j are rank 5 tensors with 4 virtual indices of bulk (bond-)dimension χ B as well as a physical index whose dimension, d, is that of the local Hilbert space of the physical system.In order to contract this infinite tensor network for the calculation of norms and expectation values, we make use of the corner transfer matrix renormalization group (CTMRG) method [34][35][36][37][38][39][40][41].The CTMRG is a power method that generates effective environments approximating the contraction of subsections of the infitite tensor network.The accuracy of this approximation is controlled by the bond dimension of the environment tensors, which we call the environment bond dimension χ E .Beyond the calculation of expectation values, the environment tensors generated with the CTMRG method can also be utilized in other ways, e.g. in the evaluation of the entanglement spectrum [42,43], cf.supplementary material [44].In order to allow us to for the first time treat more complicated Hamiltonians with a semi-classical gauge field and (large) magnetic unit cell and hosting chiral topological ground states with iPEPS we need to perform a variational ground state search.Within the ansatz-class of iPEPS of a fixed unit cell size and bond dimension χ B , the variational state-optimization [45,46] using automatic differentiation [47], implemented as explained in [48], allows us to find the optimal approximation of the ground state.This variational framework for iPEPS has been advantageous in treating e.g.frustrated spin system including recently spin liquids [15,45,49,50].For Hamiltonians, as in eq. ( 1) a treatment using the alternative methods based on imaginary-time evolution (simple update) did not yield useful results.This further illustrates the importance of the framework of variational optimization of the iPEPS ansatz.
Results.− We initially consider the U − → ∞ limit of hardcore bosons with a flux of ϕ = π/2 per plaquette.Note, that for this choice of flux, the noninteracting Hamiltonian has remarkably flat lowest and highest bands, making it a promising choice for exhibiting fractional quantum Hall physics.The flatness of a band can be quantified by the ratio of its gap to the next band and its bandwidth [51][52][53].Further, this large choice of flux is within the reach of current cold atom experiments [7][8][9][10][11][12]54].To search for fractional Hall states in this model we first map out the average occupation per site ⟨n⟩ of the ground state as a function of the chemical potential µ, see fig. 1.The calculations were performed with bulk bond dimensions χ B = 4 and χ B = 5 and values of environment bond dimension up to χ E = 120, which showed converged behaviour for the observables, as shown in figure 1.We used the Landau gauge for the vector potential in eq. ( 1) and multiples of the 4 × 1 magnetic unit cell for the unit cell of the iPEPS ansatz, cf.supplementary material for more details on gauge choices and unit cells used in this letter [44].We find two incompressible plateaus in the occupation, one at low and one at high filling.The states on these plateaus show a homogeneous density.The average occupation of the lower plateau corresponds exactly to the occupation for a ⟨n⟩ 2πϕ = ν = 1/2 Laughlin state of bosons.We observe a symmetry of occupation around µ = 0, cf fig. 1.This can be understood by relating the Hamiltonian from eq. ( 1) for a particular choice of flux and chemical potential Ĥ(ϕ, µ) to the same Hamiltonian with reversed sign for the chemical potential Ĥ(ϕ, −µ) using a particle-hole transformation as well as time reversal Θ. Applying these transformations results in with N the number of sites in the system.This also relates the eigenstates of Ĥ(ϕ, µ) and Ĥ(ϕ, −µ) and hence can be used to understand precise nature of the plateau at high filling of fig. 1 as a Laughlin state of holes at hole-filling ν = ⟨n h ⟩ 2πϕ = 1/2, with ⟨n h ⟩ = 1 − ⟨n⟩.We note further, that the incompressible plateaus in fig. 1 make up more than 20% of the parameter-space between the empty and filled states.We attribute the fact that these plateaus make up such a sizable fraction to the flatness of the lower and upper bands of the Hamiltonian at ϕ = π/2.
We now characterize the properties of the states at the ν = 1/2 filling in order to verify that this is indeed a chiral Laughlin state.We do so by first inspecting the entanglement spectrum, which should correspond to the gapless spectrum of the chiral edge modes, by following Li and Haldane [42].The entanglement spectrum is defined as the spectrum of the entanglement Hamiltonian H ent that, in turn, can identified from the reduced density matrix ρ l of a bipartition of the system ρ l = e −Hent .We can access the entanglement spectrum for our system by utilizing the bulk-boundary correspondence for PEPS [55] and do so directly in the thermodynamic limit [43] as well as on a cylinder of finite circumference.A more detailed description the methods used for the computation and details on the approximations made can be found in the supplementary material [44].
In fig.2a we can see that we find a almost perfectly chiral entanglement spectrum for the states on the incompressible plateau at the ν = 1/2 filling, confirming the chiral nature of the edge modes of the states on the incompressible plateaus in fig. 1.For the the states on the incompressible plateau at high filling, which we expect to be a Laughlin state of holes, we find that they have a chiral entanglement spectrum with reversed chirality, see fig.2a, as expected for holes.The fact that the entanglement spectra shown here are not perfectly chiral at momentum zero stems from the numerical approximations made, but the chirality can systematically be improved by increasing the approximation parameter ϑ as illustrated in fig.2b.This parameter corresponds to the bond dimension of the boundary Hamiltonian from the PEPS bulk-boundary correspondence, cf.suppl.material for more information.We also perform the entanglement spectrum on a finite cylinder of circumference 10.This discretizes the edge spectrum and allows us to identify the degeneracy of the lowest energy edge excitations, cf.fig.2c.These degeneracies (1, 1, 2, 3, 5, ...) are the ones expected for a chiral boson modes, the edge mode of the Laughlin state.
Additionally we also calculate bulk correlation functions for the incompressible states at ν = 1/2.We find a exponential decay of correlations at short distances, consistent with the bulk gap of the Laughlin states, with a correlation length of ξ/ℓ B = 1.32 in units of the magnetic length ℓ B = 1/ √ ϕ.This is compatible with the result obtained with tree tensor networks for finite systems at lower flux [31].Furthermore, we observe a finite tail in the correlation function at long distances.This tail is a numerical artifact and can be identified as consequence of the fact that we are representing a gapped chiral state with PEPS: it has been observed as well in the context of related spin systems [15].Increasing the environment bond dimension χ E , for a fixed value of χ B changes the details of the artificial long range tail, while leaving the short range exponentially decaying correlations invariant.The value at which these tails appear depends on χ B .
In order to show that we are not limited to hard-core bosons we relax the U − → ∞ condition and move to large but finite interactions U/t = 50.In our calculations we allow up to two bosons per site in this case, given At short distances the correlations decay exponentially with distance, consistent with a bulk gap of the state.We extract a correlation length of ξ/ℓB = 1.32 which is consistent with previous results [31].At long distances we find a small aritficial tail of in the correlation functions, that is due to the fact that we represent a chiral topological gapped state with PEPS.that even higher fillings are energetically strongly suppressed.In this situation we observe, for a magnetic flux of ϕ = π/2, a similar average occupation per site as in the case of hardcore bosons, with a slightly higher occupation at a given µ stemming from the fact that states including two bosons on a single site are now possible.We again find incompressible plateaus at ν = 1/2 as well as at its hole analog, cf.supplementary material [44].If we instead reduce the magnetic flux to ϕ = π/3, we additionally find a plateau at filling factor ν = 2.This plateau is expected to be associated with exotic states that belong to the Jain sequence obtained in the composite fermion picture [56].The ν = 2 state from this sequence is a bosonic integer quantum hall state, which is not intrinsically topological but rather a symmetry protected topological state [25,26,[57][58][59].
Conclustion & Outlook.− In this letter we demonstrated for the first time that the variational iPEPS framework allows for the investigation of Hamiltonians of mobile, interacting particles with semi-classical gauge fields that host chiral topological states as ground states, such as fractional Hall states and Chern insulators.The iPEPS method is a crucial complement to other well established methods for this task (such as exact diagonalization or other tensor network ansatzes) as it avoids finite size effects that can alter the correct phase diagram.Further, as it can give genuine information about the bulk physics, it can be used to benchmark whether finitesize experiments are actually displaying bulk behaviour or are dominated by their boundary.The class of Hamiltonians describing mobile, interacting particles in a semiclassical gauge field, of which we have treated a paradigmatic example in this letter, comes up in the description of several physical systems that have attracted enormous amounts of interest in recent years [8,12,[60][61][62][63].Furthermore key characteristics, such as the importance of flat bands which aid the emergence of interesting correlated phases, are often shared among these different systems.In the bosonic Harper-Hofstadter model, the (relatively) flat bands of the Hamiltonian can be understood as leading to the prominence of interaction driven physics and hence to fractional Hall states [64].Correspondingly, in twisted bilayers systems the isolated moiré flat bands are the reason for many candidate ground states such as spin-liquid states, quantum anomalous Hall insulators, and chiral d-wave superconductors [65].An important illustrative example is twisted WSe 2 which can be described by a moiré Hubbard model [66,67], in which the hopping on the triangular lattice is enhanced by a spin-dependent phase, tunable in experiments.Given the similarity of this Hamiltonian to the one studied in this letter, this puts the investigation of the physics of twisted bilayer systems well within reach of the variational iPEPS framework, promising an unbiased numerical perspective on the physics of these systems.Furthermore, recent advances drastically improved the ability to calculate excitation spectra and structure factors in the variational iPEPS framework [68][69][70].This allows direct connection to experimental results [71] and thus should establish the variational iPEPS framework clearly as a mayor tool of investigation into the above mentioned systems.
For treating the Harper-Hofstadter model using hardcore bosons at a flux per plaquette of ϕ = π/2 we used the conventional Landau gauge choice.This leads to a 4 × 1 magnetic unit cell.For the ground state search we applied both a 4 × 1 and an enlarged 4 × 4 unit cell.While both unit cell choices lead to the emergence of the Laughlin plateau, their use turns out to be convenient for different purposes.The enlarged unit cell ensures a more stable and less biased optimization, albeit slower, and is therefore the choice for scanning the different filling regimes.The original magnetic unit cell offers instead the advantage of a full periodicity of the iPEPS in one direction, which is instrumental for the calculation of the entanglement spectra (see next section), and can be safely used within the stable plateau.The same holds for the calculations of the soft-core bosons (with large interactions) at flux per plaquette ϕ = π/2.
For the case of soft-core bosons at flux per plaquette of ϕ = π/3 we choose a again can choose a Landau gauge and choose an enlarged unit cell of the iPEPS of 6 × 6.Here we observe the emergence of the incompressible plateaus shown in the main text.We additionally use a custom gauge choice with a magnetic unit cell of 3 × 2 corresponding to the same value of the flux per plaquette: by using this as iPEPS unit cell, we are able to confirm the emergence of all plateaus with faster convergence.The custom gauge choice is illustrated in in Fig. 5.We note in addition that we did not to enforce local symmetries on the tensors in our calculation as it was done in previous works [15].

II. CALCULATION OF THE LOW-LYING ENTANGLEMENT SPECTRUM IN THE THERMODYNAMIC LIMIT
In order to characterize the edge modes of a PEPS state we got from our variational optimization procedure, we can calculate the entanglement spectrum (ES).Following Li and Haldane [42], this is corresponding to the * corresponding author: weerda@thp.uni-koeln.despectrum of the edge modes of the system.The entanglement spectrum is defined as the spectrum of the entanglement Hamiltonian H ent , that in turn can identified from the reduced density matrix ρ l of a bipartition of the system ρ l = e −Hent .We can access the ES for our system by utilizing the bulk-boundary correspondence for PEPS [55], and can do so directly in the thermodynamic limit [43].
Given a PEPS network, we bipartition it into a left and right side network.These two networks have open physical indices, as well as open virtual indices at the boundary defined by the bipartition.Tracing out all physical indices of this left (right) part of the PEPS network we obtain σ l (σ r ).We can think of σ l (σ r ) as the reduced density matrices of the virtual indices at the boundary for the left and right half of the original PEPS.The PEPS bulk-boundary correspondence allows us to express the reduced density matrix ρ l of the left-side physical indices in terms of these objects: where U is an isometry that preserves the spectrum.
With the CTMRG procedure, we construct environment tensors that can be used to approximate σ l and σ r .Specifically, the CTMRG we produces environment tensors with bond dimension χ E , that approximate semiinfinite stripes of a double layer PEPS network.Thus a product of these environment tensors can be used gives an bond dimension χ E MPO-approximation of σ l or σ r , as follows: For practical purposes we note that the operator σ T l σ r σ T l has the same eigenvalues as the operator σ l σ r : see Ref. [73] for mathematical details.Note that therefore the approximation parameter from the main text is related to the MPO-bond dimension ϑ = χ 2 E .Therefore, to calculate the entanglement spectrum, we only need to find the eigenvalues of the product of the matrix product operators (MPOs) for σ l and σ r , which is an MPO itself.We obtain the lowest part of the ES by first using the VUMPS algorithm to approximate the leading eigenvector of σ l σ r [74].Starting from this leading eigenvector we use the quasi-particle excitation ansatz to find the energies of the low lying excitations at a given momentum [75,76].The bond dimension χ ev we choose for this leading eigenvector needs to be viewed as an additional approximation to be controlled.To this end, we note that a change in the bond dimension of this leading eigenvalue χ ev does not qualitatively change the chiral character of the entanglement spectrum of the fractional Hall states treated in this work, but can lead to slightly more fluctuations in the spectral values, as illustrated in fig.6.In the main text we show the result of these ES calculation and highlight especially the bond dimension of the MPO, that we use for the approximation of ρ l and 2 Laughlin state and the corresponding hole analog but no additional plateaus at higher filling.
hence H ent , serves as a crucial approximation parameter, as it can be used to increase the chirality of the spectrum.This can be understood from the fact that the MPO bond dimension can be seen as controlling the length scale of interactions of the boundary Hamiltonian following the PEPS bulk-boundary correspondence.For a truly chiral entanglement spectrum, and hence boundary modes, it can be argued by analogy to [79], that a truly nonlocal boundary Hamiltonian is needed [15] which is only archived at infinite MPO bond dimension.We further note in this context that for different Hamiltonians that permit additional symmetries such that σ l = σ r , one only needs to calculate the spectrum of one of these MPOs to get the spectrum of ρ l .

III. ENTANGLEMENT SPECTRUM ON A CYLINDER OF FINITE CIRCUMFERENCE
In addition to the calculation of the entanglement spectrum in the thermodynamic limit, we can gain information about the nature of the boundary modes by considering a system on a cylinder with finite circumference.The discretization of the available momenta makes the counting of spectral degeneracies feasible: indeed, depending on what exactly is the nature of the boundary mode, we expect different counting patterns, e.g., (1, 1, 2, 3, 5, . . . ) for a chiral boson -corresponding to the integer partitions of the added momentum in discrete units [78].
Technically, this can be done by using Lanczos methods to look for the largest eigenvalues of the ρ l on a cylinder.Alternatively, one can also use again the calculation of the leading eigenvector in combination with the calculation of the excitations of periodic systems using open boundary MPS [77].Once the leading eigenvalues and eigenvectors are determined with one of the methods above, one can use the periodic translation operator to identify the momentum associated with the different states.In the main text we showed the lowest branch of the entanglement spectrum.Here we show a larger window, where we can see an additional branch at negative momenta.This branch is actually comprised of two quasi-degenerate branches, for one additional or one missing particle in the system: Thus we observe a doubled counting, i.e., (2, 2, 4, 6, 10, ...).Higher branches exhibit more complicated countings and are affected by numerical spurious effects.
Lastly, it should be mentioned the even for the calculation of the entanglement spectrum on a cylinder, we use the environment tensors from the CTMRG calculated on a plane, which is technically an additional approximation.

IV. OCCUPATION PLATEAUS FOR SOFTCORE-BOSONS AT ϕ = π/2
When investigating the bosonic Harper-Hofstadter model we investigated the case of soft-core bosons in the case of strong interactions U/t = 50, allowing for at most 2 bosons per site.The average occupation per site is shown in fig.8.We find a situation close to the one of hard-core bosons, which is not enirely unexpected at very strong local two-body interactions.Again, a incompressible plateau emerges at low and high filling corresponding to the Laughlin states of particles and holes respectively.We notice that, while the exact particle hole transformation presented in our discussion of the hard-core case does no longer hold here, it is reasonable to expect that the arguments extend qualitatively to the case of large finite interactions.

V. NOTES ON THE CONVERGENCE WITH CHIRAL GAPPED STATES
In order to reach convergence during the optimization, one may encounter instabilities when the ground state is a gapped chiral topological state.In this scenario one can usually get back to a stable optimization by increasing the environment bond dimension before the instability occurs.It is also helpful to use the most accurate projectors possible in this case.The numerical problems can be attributed to the artificial tail in the correlation functions that develops, when treating such gapped chiral topological states.In order to include these artificial tails -which are crucial for representing such states with PEPS -in the calculation of the energy, one needs a large environment bond dimensions.This also leads to a situation in which the CTMRG procedure need a large number of steps to converge.This is the main reason why the optimization of these states -especially in the case of soft-core bosons -can be challenging in practice as oftentimes a large number of optimization steps (on the order of hundreds to thousands) are necessary at quite large environment bond dimensions.

Figure 1 :
Figure 1: Average occupation per site ⟨n⟩ as a function of µ/t for hard-core bosons with bond dimensions χB = 4(green circles) and χB = 5 (yellow diamonds).Note that the two incompressible plateaus are found at the filling expected for the ν = 1/2 Laughlin state and its hole analog.The regime given by the fractional Hall plateaus makes up more than 20% of the parameter-space in between the empty and completely filled state.

Figure 2 :
Figure 2: (a): The orange dots show the low lying part off the entanglement spectrum of a state on the incompressible plateau at low filling in fig. 1.We clearly find a chiral entanglement spectrum indicating that this is a fractional quantum Hall state.The blue dots show the low lying part off the entanglement spectrum of a state on the incompressible plateau at high filling with reversed chirality.(b): The low lying part of the entanglement spectrum is shown for three values of the approximation parameter ϑ, illustrating that we can systematically increase the chirality of the entanglement spectrum.(c): lowest branch of the entanglement spectrum on a 10 site cylinder.The counting associated with the chiral boson on the edge of the Laughlin state is observed.

Figure 3 :
Figure 3: Correlation function of the plateau state at low filling, cf.fig.1.At short distances the correlations decay exponentially with distance, consistent with a bulk gap of the state.We extract a correlation length of ξ/ℓB = 1.32 which is consistent with previous results[31].At long distances we find a small aritficial tail of in the correlation functions, that is due to the fact that we represent a chiral topological gapped state with PEPS.

Figure 4 :
Figure 4: Softcore boson occupation per site for U t = 50 and flux ϕ = π/3.We see that an additional incompressible plateau at ν = 2.The plateau at ν = 2 corresponds to bosonic integer quantum Hall states.

Figure 5 :
Figure 5: Illustration of the custom gauge choice.The blue numbers illustrate a choice for the the value of gauge fields on the edges.If we consider a single plaquette, and add up all the values for the gauge fields they add to π/3 mod 2π.Note that we must take a minus sign for those gauge fields corresponding to arrows that point clockwise for the plaquette in question.

Figure 6 :
Figure 6: Illustration showing that the bond dimension χev used for the MPS approximation of the leading eigenvector in the entanglement spectrum calculation does not qualitatively change the spectrum, while value of the bond dimension used for the MPO approximation of the boundary Hamiltonian indeed increases the chirality of the spectrum.

Figure 7 :
Figure 7: Soft-core boson occupation between the Mott lobes for U t = 50.We see again the plateau corresponding to the ν = 12 Laughlin state and the corresponding hole analog but no additional plateaus at higher filling.

Figure 8 :
Figure 8: Soft-core boson occupation between the Mott lobes for U t = 50.We see again the plateau corresponding to the ν = 12 Laughlin state and the corresponding hole analog but no additional plateaus at higher filling.