Ultracold atoms in optical lattices are the subject of intense research as these systems are excellent proxies to study phenomena of quantum many-body physics. The strongly correlated regime is of particular interest as parameter domains were earlier inaccessible in condensed matter systems. For bosonic atoms, the physics of the system is described by the Bose–Hubbard model (BHM) [1, 2]. The interplay between the hopping strengths of the bosons onto neighboring lattice sites and the on-site repulsive interaction leads to the Mott insulator (MI) and superfluid (SF) phase in various parameter domains [3, 4]. Such systems allow unprecedented control over the system parameters for easy tunability in the experiments. Optical lattices are crystals of light and are generated by superimposing counter-propagating laser beams. The lattice parameters, like the depth of the potential wells, can be controlled by changing the intensity of the laser beams, and the lattice periodicity can be varied by choosing different wavelength lasers. Additionally, different lattice geometries like square, triangular, honeycomb, and kagome are synthesized for various choices of the angles between the lasers, and the lattice dimensionality can be increased by employing more sets of lasers in other orthogonal directions. Furthermore, these systems are clean and free of any defects. However, if needed, the defects can be introduced in a controlled manner, thus allowing a systematic study of the associated effects. This has allowed the study of various novel quantum phases and the experimental observation of quantum phase transitions with good control on the system parameters [5, 6]. They have been used extensively to study the quantum phases which have eluded experiments. Examples include the much-studied supersolid phase with dipolar bosons [7–9] and in various lattice geometries [10, 11]. The supersolid phase has been experimentally observed recently for the system of dipolar gases [12, 13] and with cavity-mediated interactions [14]. Other studies include the fermionic and bosonic mixtures and the physics of phase separation [15–17], collective excitations [18, 19], the effect of the synthetic magnetic field on the supersolid phase and the Bose glass phase in a disordered system [20, 21], Kibble–Zurek mechanism, and non-equilibrium dynamics [22–24].

Atoms being charged neutral, no Lorentz force is generated when an external magnetic field is applied. However, there is a host of phenomena in condensed matter physics associated with external magnetic fields. A prime example is the physics of quantum Hall states. This is now resolved with the recent experimental realization of synthetic gauge fields using lasers [25–27]. Thus, the system of neutral bosons in optical lattices is an ideal setup to study the physics of quantum Hall effect. This is due to the possibility of realizing high flux density per plaquette of the lattice, which is otherwise difficult in condensed matter systems [28]. This has opened the door to understand the strongly correlated physics of topological states like the integer quantum Hall and fractional quantum Hall states, proposals for their realization [29], and various theoretical studies on FQH states [30–33]. Recently, FQH states have been realized experimentally using ultracold atoms in optical lattices [34].

The integer quantum Hall effect (IQHE) is observed in electronic systems under the influence of strong magnetic fields at very low temperatures, where the resistivity of the system does not assume continuous values but is quantized in integer multiples of e ²/h [35]. As a function of the magnetic field, the resistivity has plateaus corresponding to the integer values. Later, the discovery of the fractional quantum Hall effect (FQHE) in 1982 observed the resistivity quantized as the fractional multiple of e ²/h [36]. It is understood to be arising from the electronic correlations, which play no role in the IQHE [37]. The FQH states are characterized by the fractional values of the filling factor ν, which quantifies the number of the flux quanta attached per electron. Earliest identified FQH states had ν = 1/3 and 2/3. Subsequently, several states were identified with ν = 1/5, 2/5, 3/7, 4/9, …, 4/3, 5/3, 7/5, 5/2, … [38]. The Jain sequence predicts FQH states for ν = p p + 1 , where p = 1, ± 2, ± 3, … [39] and the Read–Rezayi sequence appears at fillings ν = p/2, where p = 1, 2, 3, … [40, 41]. Various studies predicting the FQH states in the lattice systems have led to immense interest in such systems. These are described by the tight-binding model for the electronic motion in the magnetic field originally studied by Harper and Hofstadter [42, 43].

In this work, we numerically study the FQH states exhibited by the ultracold bosons in a square optical lattice with a synthetic magnetic field described by the bosonic Harper–Hofstadter Hamiltonian. We use the cluster Gutzwiller mean-field (CGMF) and the exact diagonalization (ED) methods to probe the FQH states near the vacuum state. We have also studied the FQH states in the bosonic Harper–Hofstadter Hamiltonian extended by incorporating dipolar interactions. We focus, in particular, on the ν = 1/2 FQH state. This study adds to the previous literature on FQHE with dipolar interactions in the identification of the ground state among competing FQH and superfluid states obtained using the CGMF method. The superfluid state, obtained from the exact diagonalization method, has number density that is an integer multiple of the inverse of the system size and hence may not represent the ground state for the system in the thermodynamic limit. Thus, with ED, the metastability of the state is not determined. However, in order to identify the topological order of the state, we ultimately use exact diagonalization for a given number density. The FQH states are identified by studying their characteristic decaying trend of the two-point correlation function. These states are topological in nature and have fractional values of topological charge. The identification of the topological nature for these states is carried out by calculating the many-body Chern number (MBCN) [44].

Consider a system of bosonic atoms loaded in a 2D square optical lattice in the presence of a synthetic magnetic field. The system is described by the bosonic Harper–Hofstadter Hamiltonian, which is similar to the BHM except that the hopping term acquires a lattice-site-dependent Peierls phase. In the Landau gauge, the Hamiltonian is given by H ̂ BHM = ∑ p , q − J e i 2 π α q b ̂ p + 1 , q † b ̂ p , q + J b ̂ p , q + 1 † b ̂ p , q + H . c . + U 2 n ̂ p , q n ̂ p , q − 1 − μ n ̂ p , q , (1) where p (q) is the lattice site index along the x (y) directions, b ̂ p , q ( b ̂ p , q † ) is the annihilation (creation) operator at the lattice site (p, q), n ̂ p , q is the number operator, J is the hopping strength, U is the on-site interaction strength, μ is the chemical potential, and α is the flux quanta per plaquette or unit cell of the lattice. Thus, Φ = 2πα is the phase acquired by an atom while traversing around a plaquette of the lattice.

The ground state is obtained by diagonalizing the Hamiltonian using the CGMF method with the appropriate cluster size. In the CGMF method, the lattice is tiled with clusters, and neighboring clusters are coupled through the mean field of the lattice sites at the boundary of clusters. Thus, the hopping term in the Hamiltonian in Eq. 1 is evaluated exactly for the intra-cluster hopping, while it is modified using mean-field prescription for the inter-cluster hopping [45]. In this approximation, the bosonic creation and annihilation operators of the inter-cluster hopping terms are written in terms of a mean-field plus a fluctuation. The bosonic operator corresponding to lattice sites surrounding the edges of a cluster is replaced by b ̂ p , q = ϕ p , q + δ b ̂ p , q , (2) where ϕ p , q = ⟨ b ̂ p , q ⟩ , and the expectation value is taken with respect to the ground state wavefunction. A similar expression is obtained for the creation operator b ̂ p , q † . This approximation simplifies the Hamiltonian of the system to a direct sum of the cluster Hamiltonians. The Hamiltonian of the cluster is given by H ̂ C = − ∑ p , q ∈ C ′ J e i 2 π α q b ̂ p + 1 , q † b ̂ p , q + J b ̂ p , q + 1 † b ̂ p , q + H . c . + ∑ p , q ∈ C U 2 n ̂ p , q n ̂ p , q − 1 − μ n ̂ p , q − ∑ p , q ∈ δ C J e i 2 π α q ϕ p + 1 , q * b ̂ p , q + J ϕ p , q + 1 * b ̂ p , q + H . c . , (3) where the first summation is restricted over lattice sites (p, q) such that their neighboring lattice sites (p+1, q) or (p, q+1) also lie within the cluster C. The last summation is over lattice sites (p, q) that lie at the boundary of the cluster C. Using the Gutzwiller ansatz, the ground state wavefunction of the system is defined as the tensor product of cluster wavefunctions expressed in terms of the coupled basis as follows: | Ψ  GW 〉 = ∏ C | ψ c 〉 = ∏ C ∑ n 1 , n 2 … , n K L C n 1 , n 2 . . , n K L C | n 1 , n 2 … , n K L 〉 C , (4) where the cluster is chosen to be of dimension K × L. Then, for a M × N system, the number of clusters required to tile the system is M × N/(K × L). To obtain the ground state numerically, we start with an initial guess for the mean field ϕ, and the corresponding Hamiltonian matrix is calculated. The ground state wavefunction of each cluster is obtained by diagonalization of the cluster Hamiltonian. New values of ϕ are calculated for the lattice sites at the edge of the cluster, and these are used to update the cluster Hamiltonian as an improvement over the initial guess. This process is repeated for all the clusters and continued till self-consistency is achieved. For the present computations, we consider self-consistency is reached when the mean field converges upto an accuracy of ∼ 1 0 − 9 . It is to be noted that in the limit of tiling the lattice by a single cluster, that is, there is no mean-field, the CGMF method is equivalent to the exact diagonalization method. More details about the numerical procedures used in this research can be found in [46, 47].

For specific values of the filling factor ν = ρ/α, where ρ is the number density, FQH states may be favored as the ground state. The quantum Hall states are identified by the plateaus in the plot of bosonic density ρ as a function of μ, which can also be referred to as the compressibility plot. Thus, the plateaus are the indicators of incompressible ground states and the signature of quantum Hall states. Among the states, those with the integer (fractional) value of ν can be the integer (fractional) quantum Hall states. The compressibility plot is obtained by the CGMF method with the appropriate cluster size. For the present work, we consider α = 1/4 and choose the Landau gauge. It should be mentioned that we choose experimentally achievable high values of α [26], as the lattice structure becomes more prominent at high flux strengths compared to the continuum limit at low values of α. For various choices of the flux strength α, FQHE has been studied extensively and has been reported in various works [31, 46–49]. For this choice of the magnetic field and the gauge employed, the magnetic unit cell consists of four lattice sites along the y direction. So the cluster size along the y direction should be in multiples of four lattice sites in order to respect the magnetic translational symmetry. To further verify that the plateaus represent quantum Hall states, we use the exact diagonalization method to compute properties of the quantum Hall states. In this method, all the terms of Hamiltonian are calculated exactly for a system with a fixed number of particles, that is, we work in the canonical ensemble. The quantum Hall states have a characteristic property related to the decaying behavior of the two-point correlation function defined by ⟨ b ̂ † ( x ′ , y ′ ) b ̂ ( x , y ) ⟩ , where the expectation is calculated with respect to the ground state wavefunction. It has a power-law decay at the edge, while in bulk, it initially decays as an exponential, followed by a non-monotonic trend, and finally with a power-law tail [50]. ⟨ b ̂ † x , y b ̂ 0 , y ⟩ ∝ 1 / x α , y ∼ 0 , e − x / ξ , y ∼ M / 2 , x < M , 1 / x + 2 y α , y ∼ M / 2 , x ∼ M . (5) This arises because the quantum Hall state is gapped in bulk but has a gapless edge at the boundary of the system.

To identify the topological order of the state, we calculate the MBCN using the method of [51]. On a torus geometry, the generalized boundary conditions are given by the single-particle translation along the lattice resulting in the original wavefunction but upto a phase or twist angles at the boundary. Depending on the direction along which the twist angle is applied, it could either represent magnetic flux along the axis of the torus or through the center of the torus. The ground state is a function of the twist angles, and the knowledge of the ground-state manifold across the grid of twist angles (θ _x , θ _y ) is necessary to calculate MBCN. At the boundaries of the lattice, twist angles θ _x (θ _y ) are applied in the hopping amplitude along x (y) directions, and then exact hopping is allowed at the boundary using periodic boundary conditions (pbc). The geometry of the system is thus that of a torus. With the implementation of the twist angles, the Hamiltonian is given by H ̂ = ∑ p , q [ − J e i 2 π α q e − i 2 π δ x L θ x b ̂ p + 1 , q † b ̂ p , q + J e − i 2 π δ y L θ y b ̂ p , q + 1 † b ̂ p , q + H . c . + U 2 n ̂ p , q n ̂ p , q − 1 − μ n ̂ p , q ] , (6) where δ _xL is the Kronecker delta function. The ground-state manifold of the ν = 1/2 quantum Hall state is doubly degenerated on the torus geometry, so we calculate two states from the degenerated ground-state manifold Ψ₀(θ _x , θ _y ) and Ψ₁(θ _x , θ _y ) at various values of twist angles in the range [0,1] and calculate the ground-state manifold projector at (θ _x , θ _y ): P θ x , θ y = Ψ 0 θ x , θ y Ψ 0 θ x , θ y + Ψ 1 θ x , θ y Ψ 1 θ x , θ y . (7) To fix the gauge for the ground-state manifold, we choose two reference multiplets from the ground-state manifold at two different values of the twist angles Φ j ( θ x 1 , θ y 1 ) and Φ k ′ ( θ x 2 , θ y 2 ) , where j, k can take values 1 or 2. This serves to define two different gauge references with associated scalar fields Λ_Φ = det ⟨Φ _j |P(θ _x , θ _y )|Φ _k ⟩ and Λ Φ ′ = det ⟨ Φ j ′ | P ( θ x , θ y ) | Φ k ′ ⟩ , respectively. The gauge fields corresponding to Φ and Φ′ are not defined on the entire grid spanned by (θ _x , θ _y ) but form two complementary regions, where one of the fields is well defined but not the other. The MBCN can then be calculated by counting the number of branch vortices in the argument field Ω defined by Ω = arg det ⟨ Φ j ′ | P θ x , θ y | Φ k ⟩ , (8) in either of the two complementary regions (identified by the regions where either Λ_Φ or Λ_Φ′ field is vanished). While counting the number of branch vortices, the vortices and anti-vortices should be summed with appropriate signs according to the vorticity [31, 48].

We have studied the FQH states of the bosonic Harper–Hofstadter Hamiltonian, which can be experimentally realized with ultracold atoms in optical lattices. We find FQH states as the ground state of the system for various parameter regimes. To calculate the ground state, we have used the CGMF method, which describes the quantum correlations well, and we find that incompressible states occur as the ground states for different regions in the parameter space. The filling factors of these states are similar to those of the FQH states. To verify that the states are indeed FQH states, we identify the topological order of the ν = 1/2 state by calculating the MBCN. With the introduction of dipolar interactions truncated to the nearest neighbor, we show that the extended Hamiltonian also supports the FQH states as the ground state. Thus, FQH states are robust against the dipolar interactions; in fact, we find that the dipolar interactions stabilize these states against the SF state. We have verified the robustness of the ν = 1/2 FQH state against the finite size effects and its stability against the tail of the dipolar interaction. This suggests the possibility of experimentally observing the ν = 1/2 FQH state with dipolar condensates in the optical lattice. Recent experimental realization of FQH states for ultracold gas of ⁸⁷Rb atoms in the optical lattice [34] is a first step and can set a platform to further investigate the FQHE with dipolar condensates in the near future.