Various Topological Mott insulators and topological bulk charge pumping in strongly-interacting boson system in one-dimensional superlattice
In this paper, we study a one-dimensional boson system in a superlattice potential. This system is experimentally feasible by using ultracold atomic gases, and attracts much attention these days. It is expected that the system has a topological phase called topological Mott insulator (TMI). We show that in strongly-interacting cases, the competition between the superlattice potential and the on-site interaction leads to various TMIs with non-vanishing integer Chern number. Compared to hard-core case, the soft-core boson system exhibits rich phase diagrams including various non-trivial TMIs. By using the exact diagonalization, we obtain detailed bulk-global phase diagrams including the TMIs with high Chern numbers and also various non-topological phases. We also show that in adiabatic experimental setups, the strongly-interacting bosonic TMIs exhibit the topological particle transfer, i.e., topological charge pumping phenomenon, similarly to weakly-interacting systems. The various TMIs are characterized by topological charge pumping as it is closely related to the Chern number, and therefor the Chern number is to be observed in feasible experiments.
pacs:67.85.-d, 03.75.Lm, 05.30.Jp, 73.21.Cd
Keywords: Bose-Hubburd model, Topological phase, Superlattice potential, Topological charge pumping.
Recent years, topological phase is one the most interesting subjects in condensed matter physics. It is generally defined as a state characterized by a nontrivial topological number even though it has no local order parameters. Topological phase is expected to form even in one dimensional (1D) system as suggested in the celebrated work by Thouless . The topological phase in 1D system originates from geometrical similarity of the (1+1)D spacetime to 2D space, in which well-known topological states of matter, e.g., the integer quantum Hall (IQH) [2, 3] and fractional quantum Hall (FQH) states [4, 5] form. Inspired by the important observation by Thouless , certain 1D models have been studied from the view point of topological phase [6, 7, 8, 9, 11, 12, 13, 14]. Also the experiments on cold atomic gases in an optical lattice have started to “quantum simulate” such 1D systems . As one of the recent remarkable successes in the experiments, we notice the realization of the topological Thouless pumping [16, 17]. The topological Thouless pumping is a phenomenon in which a spontaneous atomic transportation takes place by changing a Floquet parameter characterizing topological properties of the Hamiltonian. The experimental successes stimulate theoretical study of the topological phase of 1D cold atomic system in an optical lattice.
Motivated by these theoretical observations and experimental successes, we shall study the systems of strongly-interacting Bose gases on 1D superlattice in this work. In experiments, interactions between cold atoms can be controlled by using optical experimental techniques, e.g., Feshbach resonance [18, 19] and orbital Feshbach resonance . The target system is described by the Bose-Hubbard model (BHM) with an applied modulate chemical potential term. Interestingly, the BHM with a modular chemical potential is expected to have non-trivial topological states [6, 7, 8, 9, 11, 12] and it can be quantum simulated by the recent experiments .
Most of the previous studies have focused on the existence of non-trivial topological phases. Some numerical studies by using the density-matrix renormalization group method (DMRG) confirmed the existence of the non-trivial topological phase called topological Mott insulator (TMI) [6, 7, 9]. TMI is classified by a topological number such as the Chern number and it has a gap in the bulk but a gapless excitation in the (spacial as well as phase diagram) boundaries. A quantum Monte-Carlo simulation was also carried out to detect the topological phase . However, the existence was verified only in a limited parameter regime, and detailed global phase diagrams are still lacking. In particular, strongly-correlated bosonic topological states with high Chern number in the bulk have not been clarified yet in a global parameter regime, nor it is understood well how the competition between the superlattice potential and the on-site repulsion determines the ground state of the system. Also, there is one important question, i.e., how the topological phases in the obtained phase diagrams are related with the topological charge pumping as bulk topological properties.
In this paper, we shall study the above problems in the strongly-interacting boson system by using the exact diagonalization [22, 23], and show explicitly relation between the equilibrium topological phases and the topological charge pumping in the adiabatic process by following Ref.. From the relation, the global phase diagram plays a role of a guide for detecting various topological charge pumping in the experiments.
This paper is organized as follows. The target boson model in the 1D optical superlattice is explained in Sec. 2. Section 3 studies the ground states of the system for the hard-core boson limit. We explain the single particle (SP) equation, which is related to the famous Harper equation in 2D electron lattice model in uniform magnetic fields, and we present the SP spectrum. We observe the ground states and their topological order by using the exact diagonalization. In Sec. 4, we clarify the ground state properties of the soft-core boson system, where the on-site interaction plays an important role to determine the ground states and their topological properties. In particular, we shall show global phase diagrams of the Chern number. The phase diagrams includes rich topological phases. In Sec. 5, we discuss a relationship between the ground states obtained by the exact diagonalization and the dynamical topological charge pumping. We shall show the Chern number of the many-body interacting ground-state is directly connected to the particle transfer performed in an adiabatic pumping cycle of superlattice. Section 6 is devoted for conclusion. Detailed calculations concerning to the Chern number and the topological charge pumping are given in appendices.
2 Bose-Hubbard model on 1D superlattice
We consider a dilute Bose gas system in a superlattice, whose Hamiltonian is given as follows,
where is the boson annihilation (creation) operator at superlattice sites , the density operator , is the hopping amplitude, is the on-site repulsion. The superlattice is created in cold atom experiments by using two different standing-wave lasers as shown in Fig. 1. The parameter in [Eq.(1)] is related to the amplitude of the superlattice potential. The parameter is the superlattice period, which is a tunable parameter in experiments, and we call modulate parameter. On the other hand, is the phase shift between the two standing-wave lasers. is a twisted phase coming from the twisted boundary condition [25, 26] and is the system size. In this paper, we consider the 1D system with the periodic boundary condition. From the view point of the topological state, adiabatic Floquet parameters of the present system are and , i.e., the Hamiltonian [Eq.(1)] is invariant under the transformations of the adiabatic parameters such as and , independently. From this fact, the adiabatic parameters span a 2D periodic parameter space, i.e., a torus . In Sec. 5, we shall study the topological charge pumping, which takes place by varying the parameter as a function of time such as .
3 Phase diagrams of hard-core Bose-Hubbard model
In this section, we shall consider the hard-core boson limit , i.e., multiply occupied states are prohibited, and then we drop the on-site interaction term in . In this limit, the system approaches to a non-interacting fermionic system . Under this condition, the SP spectrum of the model determines the ground state of the system. The previous studies [6, 7, 8, 11] showed that the SP spectrum of the present system is given by the solutions of the Harper equation, i.e., the Hofstadter butterfly , which describes the 2D lattice electron system in uniform magnetic fields. Interestingly, this correspondence can be understood by the consideration of the dimensional extension from spacial 1D system to spacial 2D system . Actually, by substituting the SP wave function into with the BHM Hamiltonian of Eq.(1), the SP equation is obtained as follows,
Then, let us compare the above SP equation (LABEL:SingleE) with the Harper equation ,
where is the hopping in the direction and is the magnitude of the applied magnetic flux per plaquette. As shown in Ref. , the modulate parameter and the phase shift in correspond to and the -component wave number , respectively. Furthermore, there exists correspondence such as , . As a future work, it is interesting to study the above correspondence from a relativistic field-theoretical view point .
From the SP equation of the BHM of Eq.(LABEL:SingleE), it is naturally expected that there exist a band insulator at specific filling factors, [filling factor is an average particle number per site]. As an example, Fig. 2 (a) shows the energy spectrum obtained by solving Eq.(LABEL:SingleE) for the case of , with various . We observe that as increasing , the SP spectrum splits into three bands, i.e., the spectrum has two band gaps, and each band has the -filling. Then in the hard-core boson system with , the first band is fully occupied, and a band insulating forms. Similarly for the system with . [See later discussion and Figs. 2 (c) and (d).] It is known that the band insulating ground state has non-trivial topological nature indicated by a non-vanishing integer Chern number .
where is the non-degenerate ground state depending on the adiabatic parameters and .
In what follows, we focus on classifying the phase at vanishing temperature (i.e., the ground state) of the system. In the bosonic system described by the Hamiltonian [Eq.(1)], three kinds of the ground states [6, 7, 8, 11] are expected to appear, i.e., superfluid (SF), trivial Mott insulator (MI), and TMI. In particular, the TMI is identified as the state with non-vanishing integer of Eq.(3) and without SF order. The TMI is regarded as an analogous state of the IQH states [3, 2].
We identify the TMIs by using the exact diagonalization [22, 31, 23], which is an efficient method to study the bulk properties of the system. Figure 2 (b) exhibits our numerical results, i.e., the phase diagram in the - -plane for the case of . For calculating the Chern number define by Eq.(3), we employed the methods proposed in Ref.  and used the discretized adiabatic parameter space with mesh. (We took as suggested in Refs. [6, 33].) As Fig. 2 (b) shows, the states at specific fillings have a non-vanishing Chern number while the others exhibit vanishing Chern number. The Chern number is quantized as at in the finite- region. Here we note that even for fairly small values of , a finite energy gap exists as seen from Fig. 2 (a), and this gap accompanies the non-vanishing . For example, the TMI at with forms for at which the first-excited energy band appears. We also exhibit density snapshots for with and in Figs. 2 (c) and (d). The ground states are non-degenerate and the (discrete) translational invariance is apparently broken there. This feature is different from that of the ordinary IQH state, in which the continuous translational invariance is preserved.
4 Phase diagrams of soft-core Bose-Hubbard model
In this section, we shall study the soft-core BHM on the superlattice. In the soft-core case, multiply occupied states are allowed at each site and therefore the on-site interaction plays an important role to determine the ground state of the system. In fact, it is expected that the competition between the on-site interaction, , and the superlattice, , leads to rich phase diagrams. That is, the TMIs with various Chern numbers form in the phase diagram in the -plane. For the numerical methods, we employ the exact diagonalization as in the study on the hard-core case. In this section, we focus on the global phase diagrams for the cases and .
To begin with, we obtained the phase diagrams in the --plane with . In Figs. 3 (a) -(c), the phase diagrams of three cases are shown, where the particle filling is denoted by as before.
Figure 3 (a) for the case shows that there are two phases, i.e., the SF, which is characterized by a finite value of the single particle density matrix (SPDM) , ( is the floor function), and the TMI phase with . We determined the ground-state phase diagram as follows; if the ground state has a non-vanishing , we regard the ground state as a TMI even though the SPDM has a finite value. That is, we regard the finite SPDM as a finite system-size effect in this case. (See later discussion.) On the other hand, if the ground state has a finite SPDM with , we regard the ground state as a SF state. Hereafter, we denote the TMI with the () as TMI ().
As seen in Fig. 3 (a), the SF phase forms for large . The typical behaviors of the SPDM along the line in Fig. 3 (a) are shown in Fig. 4. The result exhibits the system size-dependence. On the other hand for the small regime, the system exists in the TMI(+1). Figure 5 (a) shows the system size dependence of along the line in Fig. 3 (a). Interestingly, its system size dependence is much smaller than that of the SPDM. We also plot the “finite-size scaling” of the critical point obtained by in Fig. 5 (b). The result indicates that the critical point is almost independent of the system size. From the data, we can conclude that the Chern number is a good order parameter for identifying phase boundaries of the system. Similar numerical result concerning to the system size dependence of Chern number was reported in Ref. , though the target model is a spin model. In addition, we measured the energy gap between the ground state and the first excited state  along the line in Fig. 3 (a). The result is shown in Fig. 6 (a). There, gapless excitation in SF cannot be clearly seen due to the finite-size effect in the small system.
Next, we consider the case for . As shown in Fig. 3 (b), we obtain the rich phase diagram. As far as we know, this phase diagram is one of new findings in this work. There are four phases, i.e., the SF phase exists for large , and interestingly the TMIs with and appear in the small regime. In particular, the TMI with is an interesting phase, where the system permits the double occupancy at each site as shown in Fig. 7 (a) since , and then the system essentially behaves as a two-species boson system. Each species of boson fully-occupies the lowest band in the SP spectrum in Fig. 2 (a), therefore, .
On the other hand for the case , since the on-site interaction is dominant, the multiple occupancy does not appear. Bosons behave like fermions as shown in Fig. 2 (b). This situation leads to full-occupancy up to the second-lowest band of the SP spectrum. Thus, the TMI state exhibits . Similar discussion was given in Ref. . Furthermore, we found that the TMI phase with forms as an intermediate state between the two TMIs with and . We consider that in this state with , the multiply-occupied sites appear in addition to the fully-occupied lowest band of the SP spectrum.
Here, we mention the possibility of the coexistence of the SF and TMI (non-vanishing ). In our simulation, the SPDM certainly exhibits a small but finite value in the TMIs near the SF regime. For example, as one of the possible regime of the coexistence, we indicate the area such as and in Fig. 3 (b) by the dotted ellipse. There, a part of bosonic atoms form the TMI(+1) and the others may Bose condensate i.e., SF. However, due to the finite-size effect, the present calculation of the small system cannot reveal the precise behavior of the SPDM. This problem will be studied in more detail in the near future.
It is interesting to measure the energy gap along the lines in the parameter space, on which a phase transition between two different TMIs takes place. As a typical example, we plot the energy gap as a function of along the line in the phase diagram in Fig. 3 (b). The result is shown in Fig. 6 (b). For calculation of , we chose the minimums of in the adiabatic-parameter space . We find that the gap apparently closes at two transition points between the TMI(-1) and TMI(+1) and between the TMI(+1) and TMI(+2). The result was independent of the size of the tridiagonal matrix of the Lanczos algorithm. Our numerical study obviously captures the level-crossing of the lowest and first-excited states. This level-crossing induces the change of the topological number of the ground state .
For the line in the phase diagram in Fig. 3 (b), we show typical distributions of the Berry curvature, i.e., the integrand of Eq.(3). See Fig. 6 (c). Its integral over gives the non-vanishing integer Chern number. From the data, the Berry curvature generated from the gauge field () has no topological defect (quantized vortex), but it measures the density of the magnetic flux penetrating the surface of the 2D torus . Then, non-vanishing Chern number indicates the existence of a hiding monopole that exists in the interior of the torus . In other words, one can define a 3D gauge field induced by in the interior of the torus, which corresponds to the magnetic monopole.
Let us see how the phase diagram changes as the filling factor increases further. Figure 3 (c) shows the phase diagram for the filling and . Here, we have a richer phase diagram than the lower filling cases in Figs. 3 (a) and (b). In the phase diagram Fig. 3 (c), the TMIs with the larger integer form for the large and small regime. In Fig. 3 (c), the highest value of is . This value corresponds to the total particle number in the unit cell, i.e., three particle. In general, for and where and are co-prime integers, the possible highest value of is expected to equal the total particle number in the unit cell.
To get the physical picture of the state with large , we apply the same discussion of the TMI (+2) in Fig. 3 (b) to the present case, i.e., the system permits higher occupancy as increasing the value of with keeping fixed.
Here we remark on the phase boundary of the SF in the small and small regime in Fig. 3 (c). In that parameter region, we cannot clearly identify the phase boundary between the non-topological MI and SF state since both phases have . Furthermore, the typical behavior of the SPDM shown in Fig. 3 (d) indicates that we can obtain the phase boundary only approximately. Therefore, we show the approximate phase boundary between the MI and SF with the dotted line in Fig. 3 (c).
We shall show the obtained phase diagrams for other values of and with various fillings. Figures 8 (a)-(c) show the global phase diagrams for . Similarly to the case shown in Fig. 3, phase diagrams are getting complicated as the filling increases. We note that as increasing , the TMIs with the larger values of appear for large as the on-site energy of the multiply-occupied states is smaller than there. Similarly, we obtained the phase diagrams for the case of , and the results are shown in Figs. 9 (a)-(d). Generally speaking, the phase structures and their dependence on and exhibit similar behaviors with those of the previous cases and . Unfortunately, the present study using the exact diagonalization of the small system is not good enough to get more detailed phase diagrams. Further numerical studies with large system size, e.g., by using the DMRG, is desired.
Finally, we consider the filling-factor dependence of the system, i.e., the phase diagram in the --plane. Due to the on-site interaction, the SP picture is no longer meaningful, i.e., the energy spectrum of the ground state and the energy gap cannot be described by the SP equation of Eq. (LABEL:SingleE), even though TMIs may form in certain parameter regions.
We obtained the global phase diagram in the --plane, which is shown in Fig. 10 (a). The system size , , and . The results indicates the existence of two types of non-trivial TMI. We call the first one stable TMI (STMI) phase, in which is robust for the change of the value of as shown in Fig. 10 (b). On the other hand, we call the second one random Chern number TMI (RTMI) phase, in which takes various integer values as shown in Fig. 10 (c). As seen from Fig. 10 (b), the STMI phase forms when the filling-factor and the modulate parameter are tuned to produce the band-insulator regime, and the Chern number is determined as in the previous discussion based on the intuitive picture shown in Fig. 7. On the other hand, the RTMI forms when the parameters and are not located at the band insulator of the SP spectrum, and the interplay of the on-site repulsion and superlattice plays an essential role there. As far as we know, the phase diagram in Fig. 10 (a) is one of new findings in this work.
To study detailed properties of the RTMI, we measured the energy gap for the cases and . The obtained results are shown in Fig. 10 (d). It is obvious that the energy gaps are small for most of the values of . This indicates that the RTMI ground states are very close to the first-excited state. Then, the level-crossing frequently occurs as varies and states with various Chern number appear as the ground state.
5 Topological charge pumping and Chern number
In Sec. 4, we clarified the ground-state phase diagram of . Here we expect that each TMI with different exhibits different transport properties when one varies the cyclic pumping parameter adiabatically, e.g., as (), where is the period of one pumping cycle. Such a charge pumping phenomenon was studied both theoretically [1, 36, 38, 37, 24] and experimentally [16, 17, 39] for similar models to the present one. In this section, we shall show a connection between the Chern number [Eq.(3)] and the particle transfer.
To begin with, we introduce a current operator,
From the current operator , the current density at time can be directly calculated . By using the genuine ground state denoted by , which satisfies the Schrödinger equation and an adiabatic instantaneous ground state denoted by , the current density is expressed as
Here, we assume that the many-body interacting system has a finite energy gap and introduce energy spectrum , where the corresponding eigenstate is the instantaneous normalized eigenfunctions satisfying . For , the eigenfunction is nothing but the adiabatic instantaneous ground state. Then, the genuine ground state is approximately expanded in terms of the states as 
where we have assumed
We are interested in the current density averaged over ,
Then, we obtain the total particle transfer for one pumping cycle as follows,
By introducing the Berry connection () by , is expressed as
Finally, by changing variables from to , is expressed as
The above relation can be regarded as the interacting many-body version of the similar relation in the SP picture obtained by using a single Wannier state and the Bloch wave function , which is applicable for a non-interacting SP system.
Here, we note that the data shown in Fig. 6 directly show the time evolution of because . corresponds to the number of the transfer particle to a nearest-neighbor unit cell. The time evolution of originates from the form of the superlattice potential depending on the parameter . As one of typical example of the behavior of , we focus on the TMI (-1) state in the left panel in Fig. 6 (c). We consider that the superlattice potential can be described as [16, 17], where the first term is the standard optical lattice potential with amplitude and the second term is another lattice potential with amplitude as shown in Fig. 1. We set a typical ratio . In Fig.11, we show with and the corresponding particle density at denoted by in two nearest-neighbor unit cells. From the data, transfer property (hopping tendency) can be intuitively understood. In Fig. 11 (a), the particle density of the two lattice sites between the nearest-neighbor unit cells are almost unity. Then, particle is hard to hop between the nearest-neighbor unit cells due to strong and the particle current is suppressed. The result directly appears as the small value of the Berry curvature at in the left panel in Fig. 6 (c). On the other hand in Fig. 11 (b), the particle density of the two lattice sites between the nearest-neighbor unit cells are much less than unity and particle is easy to hop between the nearest-neighbor unit cells. The condition leads to significant particle current. The tendency clearly appears in Fig. 6 (c). The value of the Berry curvature near is negatively large, i.e., the current to a left nearest-neighbor unit cell is large. The total amount of the negative current corresponds to the Chern number ().
In fact, the current density is determined by the wave functions of the genuine and instantaneous ground states, but we think that the above observation sheds light on the intuitive understanding of the relation between the charge pumping and the Berry connection for the interacting Bose particle cases. Here, we should mention the study on the system with edges in Ref. . There, focusing on the bulk edge-correspondence in topological charge pumping, (for related experiments on cold atoms, see Refs. [41, 42]) it was shown that the Berry connection in the temporal gauge is directly related to the shift of the center of mass of the system, and the charge pumping is derived by that Berry connection. For the present system without edges, we think that there exists a direct relation between the Berry connection and the charge pumping. This problem is under study, and we hope that the results will be published in the near future.
Equation (12) shows that TMIs exhibit various amounts of the charge transfer in the topological charge pumping depending on the value of the Chern number . Therefore, the obtained phase diagrams of in Sec. IV can be a guide for detecting new properties of topological charge pumping in the strongly-correlated bosonic systems in experiments on ultracold atoms.
In this paper we studied the strongly-interacting boson systems in superlattice potential, which are feasible in 1D optical superlattice system of ultracold atoms. The SP properties of the system are described by the Harper equation, and the system exhibits the band insulator in the hard-core boson limit. The band insulator has a non-trivial Chern number, and it was calculated by obtaining many-body wave functions with varying two adiabatic parameters in . For the hard-core boson case with using the exact diagonalization, we first verified the locations of the band insulator corresponding to the TMI state, which is characterized by a non-vanishing Chern number. Then, we studied the soft-core boson case. There, the competition (interplay) between the superlattice amplitude and on-site interaction plays an important role. As a result, we observed the various TMIs with different integer Chern numbers. For the modulate parameter and , we clarified the global phase diagrams of the Chern number by using the exact diagonalization. The numerical results obtained for different system sizes show that the behaviors of the Chern number are almost independent of the system size. Therefore, the obtained phase diagrams in the present study captures the essential properties of the system. Interestingly, these obtained phase diagrams exhibit a very rich phase structure including various TMIs with large Chern numbers. From the obtained results, we see that the TMI with a high Chern number tends to form when the superlattice amplitude is getting large. We conclude that this behavior originates from the higher occupancy of bosons per lattice site. Furthermore, we clarified particle-filling dependence of the ground state for by studying the --phase diagram. We found that there are two types of non-trivial TMI, i.e., the RTMI and the STMI. The results concerning to the STMI are in fairly good agreement with the DMRG result for a similar model in Ref. , while the existence of the RTMI with the interesting behavior of the Chern number is one of new findings of the present work. We also found that level-crossing frequently occurs in the parameter region of the RTMI.
In Sec. 5, we studied the relationship between the particle transfer and the Chern number of many-body wave functions. We conclude that in order to detect the various TMIs in real experiments, the charge pumping measurement by adiabatic time-evolution of the parameter [16, 17] is useful. To measure the TMIs in real experiments, the Streda formula  is expected to be useful as pointed out in Ref. .
We acknowledge Y. Takahashi and S. Nakajima for helpful discussions. Y. K. acknowledges the support of a Grant-in-Aid for JSPS Fellows (No.17J00486).
Appendix A. Derivation of Eq. (6)
where and are time-dependent expansion coefficients. Since the genuine ground state satisfies the Shrödinger equation and adiabatic setup allows us to approximate , and () along time evolution, we can show that the -th coefficient satisfies the following time differential equation,
Appendix B. Derivation of Eq. (7)
where the term was canceled. Here we drop the last term in Eq. (A.16) since we assume the energy gap between the ground state and the excited states, , is large enough. Then we notice the following equations,
The above equations are obtained from and . Then substituting the above equations in Eq. (A.16) and using the complete relation of the state bases , the current is expressed as
-  Thouless D J 1983 Phys. Rev. B 27 6083
-  Thouless D J, Kohmoto M, Nightingale P and den Nijs M 1982 Phys. Rev. Lett. 49 405
-  Kohmoto M 1985 Ann. Phys. (N. Y). 160 343
-  Laughlin R B 1983 Phys. Rev. Lett. 50 1395
-  Wen X -G 2004 Quantum Field Theory of Many-body Systems: From the Origin of Sound to an Origin of Light and Electrons, Oxford Graduate Texts (OUP Premium, New York)
-  Zhu S L, Wang Z D, Chan Y H and Duan L M 2013 Phys. Rev. Lett. 110 075303
-  Deng X and Santos L 2014 Phys. Rev. A. 89 033632
-  Lang L J, Cai X and Chen S 2012 Phys. Rev. Lett. 108 220401
-  Matsuda F, Tezuka M and Kawakami N 2014 J. Phys. Soc. Japan 83 083707
-  Li T, Guo H, Chen S and Shen S Q 2015 Phys. Rev. B 91 134101
-  Xu Z, Li L and Chen S, 2013 Phys. Rev. Lett. 110 21530
-  Ganeshan S, Sun K and Das Sarma S 2013 Phys. Rev. Lett. 110 180403
-  Hu H, Cheng C, Xu Z, Luo H G and Chen S 2014 Phys. Rev. B 90 035150
-  Roscilde T, Faulkner M F, Bramwell S T and Holdsworth P C W, New J. Phys. 18, 075003 (2016).
-  Lewenstein M, Sanpera A and Ahufinger V 2012 Ultracold Atoms in Optical Lattices: Simulating Quantum Many-body Systems (Oxford University Press)
-  Nakajima S, Tomita T, Taie S, Ichinose T, Ozawa H, Wang L, Troyer M and Takahashi Y 2016 Nat. Phys. 12 296
-  Lohse M, Schweizer C, Zilberberg O, Aidelsburger M and Bloch I 2016 Nat. Phys. 12 350
-  Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys.80 885
-  Inouye S, Andrews M R, Stenger J, Miesner H -J, Stamper-Kurn D M and Ketterle W 1998 Nature 392 151
-  Zhang R, Cheng Y, Zhai H and Zhang P 2015 Phys. Rev. Lett. 115 135301
-  Chen Y -A, Huber S D, Trotzky S, Bloch I and Altman E 2010 Nat. Phys. 7 13
-  Zhang M and Dong R X 2010 Eur. J. Phys. 31 591
-  Raventos D, Gras T, Lewenstein M and Julia-Diaz B 2017 J. Phys. B 50 113001
-  Hatsugai Y and Fukui T 2016 Phys. Rev. B 94 041102(R)
-  Rousseau V G, Arovas D P, Rigol M, Hebert F, Batrouni G G and Scalettar R T 2006 Phys.Rev.B 73 174516
-  Niu Q, Thouless D J and Wu Y S 1985 Phys. Rev. B 31 3372
-  Rigorously, the connection is achieved though the Jordan-Wigner transformation.
-  Hofstadter D R 1976 Phys. Rev. B 14 2239
-  Kraus Y E and Zilberberg O 2012 Phys. Rev. Lett 109 116404
-  Kuno Y and Ichinose I, in preparation.
-  Noack M and Manmana S R 2005 AIP Conf. Proc. 789 93-163
-  Fukui T, Hatsugai Y and Suzuki H 2005 J. Phys. Soc. Japan 74 1674
-  Zeng T -S, Zhu W and Sheng D N 2016 Phys. Rev. B 94 235139
-  In our exact diagonalization, the energy gap can be obtained from a tridiagonal matrix with the optimal size of the dimension from the Lanczos algorithm.
-  A similar results are plotted in Ref..
-  Wei R and Mueller E J 2015 Phys. Rev. A 92 013609
-  Qian Y, Gong M and Zhang C 2011 Phys. Rev. A 84 013608
-  Wang L, Troyer M and Dai X 2013 Phys. Rev. Lett. 111 026802
-  Kraus Y E, Lahini Y, Ringel Z, Verbin M and Zilberberg O 2012 Phys. Rev. Lett. 109 106402
-  Asbóth J K, Oroszlány L, Pályi A 2016 A Short Course on Topological Insulators: Band Structure and Edge States in One and Two Dimensions (Berlin: Springer)
-  Gaunt A L, Schmidutz T F, Gotlibovych I, Smith R P and Hadzibabic Z 2013 Phys. Rev. Lett. 110 200406
-  Mukherjee B, Yan Z, Patel P B, Hadzibabic Z, Yefsah T, Struck J and Zwierlein M W 2017 Phys. Rev. Lett. 118 123401
-  Streda P 1982 J. Phys. C Solid State Phys. 15 L717