Weak topological insulating phases of hard-core-bosons on the honeycomb lattice

We study the phases of hard-core-bosons on a two-dimensional periodic honeycomb lattice in the presence of an on-site potential with alternating sign along the different y-layers of the lattice. Using quantum Monte Carlo simulations supported by analytical calculations, we identify a weak topological insulator, characterized by a zero Chern number but non-zero Berry phase, which is manifested at either density 1/4 or 3/4, as determined by the potential pattern. Additionally, a charge-density-wave insulator is observed at 1/2-filling, whereas the phase diagram at intermediate densities is occupied by a superfluid phase. The weak topological insulator is further shown to be robust against any amount of nearest-neighbor repulsion, as well as weak next-nearest-neighbor repulsion. The experimental realization of our model is feasible in an optical lattice setup.


I. INTRODUCTION
Since the discovery of topological insulators (TIs), tremendous effort has gone into the understanding of this novel phase of matter, both theoretically and experimentally [1]. The hallmark of a topological insulator is the existence of a bulk band gap, similar to an ordinary insulator, along with protected gapless surface states. While conducting surface states can also be observed in normal band insulators, the signature that makes the TIs unique is the topological protection of the surface states by time reversal symmetry. The application prospects of TIs in spintronic devices and quantum information technology have projected them as worthy candidates for frontier research in condensed matter physics.
A system with a non-trivial value of ν 0 is known as a strong topological insulator (STI), where gapless surface states are manifested on each two-dimensional (2D) surface of the system. On the other hand, if we try to form a three-dimensional structure by stacking layers of 2D STIs along some particular direction, we end up with a system that exhibits gapless states on some of its 2D surfaces (depending on the stacking orientation), while the other surfaces remain gapped. This system is referred to as a weak topological insulator (WTI) for which the strong topological index ν 0 measures to be zero, but some of the weak topological indices ν attain non-trivial values. In general, a d-dimensional WTI can be visualized as a system constructed by stacking (d − 1)-dimensional STIs.
According to the tenfold periodic table of topological phases [6], in the case of BDI symmetry class, a STI phase is manifested in one-dimension (1D) with a Z topological invariant, whereas in 2D the strong topological index is zero. If we now think of a system obtained by stacking L y 1D BDI chains with topological index ν, along the ydirection, the resulting 2D system will behave as a WTI. This system will manifest conducting edge states only along the edges localized at the two ends of the lattice in the x-direction. The weak topological index ν x in such a case can be measured by averaging the strong topological index over L y layers.
The main difference between a strong and a weak topological insulator lies in the robustness of their edge states. While for STIs a symmetry-preserving disorder can never create a gap at the gapless edges, this is not the situation for WTIs. Apparently, the protection of the edge states in WTIs requires lattice-translational symmetry. Therefore, it is only natural to think that even a small disorder in the system can destroy the topological phase. However, in Ref. [7], the authors have shown that the conducting edge state can actually persist in the presence of disorder, as long as the time-reversal symmetry and the bulk gap are preserved. While the experimental verification of STIs have been performed in diverse classes of materials [8][9][10] since its theoretical prediction, there are only a few examples of materials exhibiting WTI phases [11][12][13]. Further venues are needed for the study of WTI phases [14], their stability properties [15,16], and the effects of interactions [17].
In recent years, the study of topological phases in bosonic systems has been a center of attraction in condensed matter physics. Due to the condensation property of bosons, the realization of topological phases requires interaction among the particles. This could, in fact, enhance the richness of the various topological phases observed in a bosonic system. Moreover, recent advancements in optical lattice experiments have created a promising platform, where different phases of interacting and non-interacting bosonic systems can be realized in a controlled manner.
Inspired by these current advances, in this paper we study hard-core-bosons (HCBs) on a 2D periodic honeycomb lattice, subjected to an on-site potential with alternating signs along the different y-layers. Using quantum Monte Carlo (QMC) technique supported by analytical calculations, we find that the phase diagram exhibits three insulating phases at densities 1/4, 1/2 and 3/4, separated from each other by a superfluid region. Depending on the choice of the on-site potential form, either the insulator at 1/4 or 3/4 filling is found to be a WTI, which manifests a nontrivial Berry phase and the existence of edge states along the x-edges of the lattice. We show that the topological phase is robust against any amount of nearest-neighbor (NN) repulsion, as well as The y-layers are indexed by ℓ = 1, 2, . . . and are highlighted by blue (red) backgrounds indicating the presence of an on-site potential W0 (−W0) along the layer. The dashed black lines represent the bonds connecting the lattice across the boundaries. Each unit cell, delineated by a green rectangle, contains four sites labeled by the numbers 1 to 4. weak next-nearest-neighbor (NNN) repulsion among the HCBs.
The paper is organized as follows. In Section II we present the model, the numerical techniques and the relevant order parameters. In Section III we discuss the phase diagram of the model and analyze the different phases. The edge states of the insulating phases of the model are analyzed using QMC methods in Section IV. In Section V we calculate the topological invariants for the insulating phases. Next, in Section VI the effect of NN and NNN repulsion on the WTI is presented. Lastly, in Section VII we conclude. In Appendix A we analyze the band structure of the model.

II. MODEL AND FORMULATION
We consider HCBs in a 2D periodic honeycomb lattice, as depicted in Fig. 1, governed by the Hamiltonian creates (annihilates) a HCB at site i, i, j represent NN pairs of sites, t is the amplitude of NN hopping, W i is the on-site potential at site i and µ denotes the chemical potential. We take the NN hopping as the unit of energy and set t = 1 for our numerical calculations. In our study W i forms a periodic potential along the y-direction with a period of two lattice sites, i.e., we take W i = W 0 (−W 0 ) for layers (along the y-direction) labeled by odd (even) values of ℓ. We shall assume that the lattice constant is a = 1 throughout.
To study the various phases of the Hamiltonian in Eq. (1), we use the Stochastic-Series-Expansion (SSE) technique [18,19], a quantum Monte Carlo method, employing directed loop updates [20,21]. To capture the ground state-properties of a L × L honeycomb lattice using SSE, all simulations have been done at low enough temperatures such that the inverse temperature β ∼ L [22].
To construct the phase diagram using SSE we use four order parameters: average density ρ, superfluid density ρ s , structure factor S(Q) and dimer structure factor S D (Q).
The average density of a system containing N s sites isρ = in i /N s , wheren i =d † id i gives the number of HCBs (either 0 or 1) at site i. To calculate the superfluid density using SSE, we employ the following expression in terms of the winding numbers Ω x and Ω y along x and y-directions [19], where · · · represents ensemble average. For example, the winding number Ω x can be calculated by counting the total number of operators N + x (N − x ) transporting particles in the positive (negative) x-direction, according to the formula Ω where L x is the length of the lattice along the x-direction.
Next, the structure factor per site is expressed as, where r i = (x i , y i ) is the position of site i. To calculate the structure factor for particles in a L × L honeycomb lattice, we can always use a transformation on the lattice to straighten the bonds along the y-direction, such that the resulting lattice looks like the one depicted in Fig. 1 b. With the use of the position vectors r of this new transformed lattice, the allowed values of the wavevector Q coincide with those of an L × L square lattice, i.e., Q = (2πp/L, 2πq/L), where p = 0, 1, · · · , L − 1 and q = 0, 1, · · · , L − 1. To detect the presence of diagonal long-range orders in the system, we have calculated S(Q) for all possible values of Q and identify the ones at which the structure factor displays peaks. Lastly, we define the dimer structure factor as whereD α =d † αLdαR +d † αRdαL is the dimer operator defined on the α-th NN bond aligned along x-axis with α L , α R being the two lattice sites attached to this bond. In Eq. (4), the summation runs over N b NN bonds oriented along x-axis and the vectors R represent the position coordinate corresponding to the midpoints of these bonds in the transformed lattice in Fig. 1 b. The dimer operator is chosen in a way such that it will give a nonzero expectation value only when a dimer is formed, i.e., when the constituent particle hops back and forth along the NN bond.

III. PHASE DIAGRAM
To construct the phase diagram we study the Hamiltonian in Eq. (1) at various values of W 0 by varying the chemical potential µ. Fig. 2 depicts the variations of HCB density ρ, superfluid density ρ s , structure factor S(0, π) and dimer structure factor S D (0, π) as a function of µ, where the value of W 0 is fixed at 6.0. The three plateaus in the ρ − µ curve (apart from the trivial ones at ρ = 0 and ρ = 1) clearly indicate the presence of three incompressible insulators at densities 1/4, 1/2 and 3/4. At these plateaus the superfluid density becomes zero, whereas in the intermediate regions it attains some non-zero value, thus separating the three insulating regions by a superfluid phase. To understand the nature of the insulators we have calculated the structure factor S(Q) and dimer structure factor S D (Q) for all possible values of Q. We find that S(Q) peaks only at wavevector Q = (0, π), whereas S D (Q) displays peaks for both Q = (0, π) and Q = (π, 0) with the same peak value (in Fig. 2 We note that the results in Fig. 2 are independent of the sign of W 0 , i.e., whether we choose W i to be positive (negative) for odd (even) values of ℓ in Fig. 1 or the reverse scenario, Fig. 2 remains unaltered. In the following, we assume W 0 > 0 to analyze our results, but a similar analysis can be extended to the reverse scenario as well. In Section IV we will see that the sign of W 0 nevertheless plays a role in the characterization of the different phases.
Since the on-site potential for the layers labeled by even values of ℓ is −W 0 , upto half-filling the particles will prefer to occupy these layers only, keeping the odd ℓ layers completely empty. Due to the presence of NN hopping, at 1/4-filling, it is energetically favorable for the system to fill the upper two sites of each unit cell (i.e., site 1 and 4) by one particle only, so that this particle can hop back and forth between sites 1 and 4 of two adjacent unit cells to further lower the energy of the system. As a result of this hopping process dimers are formed between two sites belonging to the upper half of two different unit cells. Due to the formation of these dimers there is no net flow of HCBs in the x or y-directions, which makes the phase insulating in nature. The structure of this dimer insulator at 1/4-filling is depicted in Fig. 3 a. We note that at each even ℓ level we have one dimer which involves two boundary sites when the open boundary condition is applied along x-direction with zigzag edges. Now, let us analyze the structure factor, Eq. (3), for Q = (0, π), Although the summation in Eq. (5) is over all possible pairs of sites in the lattice, only those pairs for which both sites are occupied will have a non-zero contribution.
Since for the dimer insulator at ρ = 1/4 (see Fig. 3 a), all particles reside on the even layers only, for all contributing pairs y i − y j is even, so Eq. (5) reduces to At 1/4-filling there are N s /4 particles in the system and each of them participates in N s /4 pairs in the summation (including the case where i = j) with a +1 contribution to the structure factor. Therefore, for the dimer insulator at ρ = 1/4, S(0, π) attains the value This result matches well with the result in Fig. 2 b. It is clear from the discussion above that at 1/4-filling, as long as the particles are constrained to reside on alternate layers, the value of S(0, π) will be 0.0625. This value is independent, e.g., of whether the particles form a dimer insulator or arrange themselves in a charge-density-wave (CDW) pattern.
To manifest the formation of dimer insulator at density ρ = 1/4, we next calculate the dimer structure factor S D (Q) as prescribed in Eq. (4). Since in our system the dimers are formed along the NN bonds oriented along x-direction of the lattice, we have defined the dimer structure factor such that it will detect dimers along these bonds only. Now, if we think about the dimer- FIG. 4. Pictorial description of hopping processes in the superfluid region in-between the dimer insulator at ρ = 1/4 and CDW structure at ρ = 1/2.
insulator structure corresponding to ρ = 1/4 (as depicted in Fig. 3 a) for the transformed lattice in Fig. 1 b, it is easy to see that for any two dimers with midpoints As a result, the dimer structure factors for Q = (0, π) and Q = (π, 0) reduce to the exact same expression, In the dimer insulator phase, with N b being the total number of NN bonds along x-direction, there are N b /2 dimers in the system. Each of these dimers will participate in N b /2 pairs (of dimers) in the summation, with a +1 contribution towards the dimer structure factor. Therefore, the dimer structure factor reduces to, which is indeed attained in Fig. 2 b at filling 1/4. Next, at 1/2-filling, the layers labeled by even ℓ values are completely filled, such that the upper two sites of each unit cell are occupied by two HCBs. Therefore, at this density the dimers of 1/4-filling disappear completely and we have a CDW, similar to the one depicted in Fig. 3 b, which is insulating in nature. This can further be verified from Fig. 2 b, where we can see that at ρ = 1/2 the dimer structure factor [S D (0, π)] vanishes and the structure factor [S(0, π)] shows a peak with value 0.25. The half-filled system contains N s /2 particles in total and each of them participates in N s /2 pairs of sites, which has a non-zero contribution towards the structure factor. So, Eq. (6) in this case simply becomes which is the maximum value attained by S(0, π). Finally, the structure of the insulator at 3/4-filling is depicted in Fig. 3  cells this means that the upper two sites (sites 1 and 4) of each cell are fully occupied and the lower two sites (sites 2 and 3) share a single HCB. Again by virtue of NN hopping, the particles at odd ℓ-levels can further lower their energy by hopping back and forth between sites 2 and 3 of each unit cell. As a result, dimers are formed in the lower half of each unit cell. The main difference between the dimers formed at ρ = 1/4 and ρ = 3/4 is that, the dimers at 1/4-filling are formed between two sites belonging to two different unit cells, whereas the sites involved in 3/4-filling are residents of the same unit cell. Since the number of dimers formed in this case coincides with the one for ρ = 1/4, the dimer structure factor attains the same peak value 0.25 in this situation as well. The value of the corresponding structure factor can be extracted by realizing that out of the 3N s /4 particles in the system, N s /4 reside on the odd ℓ layers, whereas N s /2 particles are located at even ℓ layers. So, in total there are 2(N s /4)(N s /2) pairs for which the separation between the particles along the y-axis (i.e., y i − y j in Eq. (5)) is odd. Clearly each of these pairs will contribute −1 to the structure factor (as e iπ(yi−yj) = −1 for these cases). On the other hand, for (N s /4) 2 +(N s /2) 2 number of pairs, y i −y j is an even multiple of the lattice constant, which gives rise to a positive contribution to the structure factor. Therefore, the structure factor for this insulator attains the value, which coincides with the one for 1/4-filling, Eq. (7). Next we turn to discuss the superfluid phase. In the intermediate regions between the insulating phases, where the superfluid density is finite, both the structure factor S(0, π) and dimer structure factor S D (0, π) admit nonzero values, see Fig. 2. Interestingly, in these intermediate regions, we observed anisotropy in the superfluid density where ρ y s , the superfluid density along y-direction of the lattice, is much larger than the one along the xdirection, ρ x s . So while the superfluid retains some additional structure from the two neighboring insulators in the transition region, it can still superflow in both directions, as we now argue. The mechanism is described in Fig. 4, where we take for example the transition between the insulators at fillings 1/4 and 1/2. Consider a situation where we add one HCB to the dimer insulator at ρ = 1/4. At this range of fillings the particles will naturally prefer to occupy the even layers, having the lower on-site potential, so the extra particle chooses to reside on one of the sites attached to bond b 1 in Fig. 4 a, effectively generating a doubly occupied bond. Now, this extra particle can hop through the lattice giving rise to superfluidity in two different ways. Firstly, the extra particle can follow a two-step hopping process similar to the one depicted by blue arrows in Fig. 4 a. By virtue of this process, effectively the doubly occupied bond b 1 has hopped to bond b 2 (Fig. 4 b) giving rise to superfluidity along y-direction. Since the difference of energies of the particle at the initial (or final) and intermediate step of this process is 2W 0 , the energy gained by the particle during this process is ∼ t 2 /(2W 0 ). Secondly, the particle can also follow a three-step hopping process depicted by the green arrows in Fig. 4 a. This process results in a configuration as shown in Fig. 4 c, where the doubly occupied bond b 1 has effectively hopped along x-direction of the lattice to the bond b 3 , contributing to a non-zero superfluid density ρ x s . On account of the fact that both the intermediate sites involved in this hopping process have energies higher than the initial or final sites, by an amount of 2W 0 , one can see that the energy gain in this process is ∼ t 3 /(4W 2 0 ). Therefore, in comparison the doubly occupied bond can always gain more energy by hopping in the y-direction of the lattice than in the xdirection. Consequently, anisotropy is developed in the superfluid density with ρ y s > ρ x s . Similar arguments hold for the superfluid regions between any two insulators.
The complete phase diagram of the model in the (µ, W 0 ) plane is depicted in Fig. 5 in the atomic limit, i.e., when the hopping t is turned off (Fig. 5 a) and for t = 1.0 (Fig. 5 b). The phase boundaries are obtained from QMC (solid lines with points), performed for a 20 × 20 periodic honeycomb lattice with inverse temperature β = 120.
In the presence of finite NN hopping t = 1.0 (Fig. 5 b), for W 0 = 0, the superfluid phase fills the range between the Mott lobes at densities ρ = 0 and ρ = 1. Beyond some critical value of W 0 , additional insulating lobes start to appear at densities 1/4, 1/2 and 3/4 separated by superfluid regions. The insulator at half-filling is a CDW, whereas the other two are dimer insulators. One can see that the phase boundaries obtained from QMC are more or less consistent with the calculated band edges from Appendix A (dashed lines in Fig. 5 b), except in the neighborhood of W 0 = 1, where they slightly deviate. Since the analytical results demonstrate the phase boundaries for a lattice in the thermodynamic limit, we expect the deviations to be smaller for larger system sizes. Within the error bars of our QMC calculations the critical value of W 0 beyond which the insulating phases appear, is 1.4 for the CDW at half-filling, whereas for the dimer insulators it appears to be 1.5. However, the calculated band edges predict that the tips of all three insulating lobes lie at W 0 = 1.
Next, in the atomic limit (Fig. 5 a), we see that only the CDW insulator at half-filling survives and both the dimer insulators vanish completely. Indeed, the dimers in the dimer insulators are formed by virtue of the NN hopping in the presence of a finite W 0 . On the other hand, in the atomic limit the CDW insulator appears as soon as we have a non-zero W 0 . In fact, the presence of NN hopping can destroy this structure by transforming it into a superfluid, as is indeed observed in Fig. 5 b. Beyond some critical value of W 0 the CDW phase sets in because at this stage W 0 dominates over hopping t and it becomes energetically favorable for the particles to be frozen in this structure instead of moving around in a superfluid phase. The boundaries of the CDW phase in Fig. 5 a can be determined by considering the change in the total energy of the system when we introduce an additional HCB in the half-filled system manifesting the CDW phase. At half-filling all the particles occupy the sites with on-site potential −W 0 . As a result, in this situation the total energy of the system is simply E[N s /2] = −µN s /2 − W 0 N s /2. Now, if we try to add another HCB to the system, this additional particle will have to occupy a site with on-site potential W 0 . So, in this scenario the total energy of the system will be E[N s /2 + 1] = −µ(N s /2 + 1) − W 0 N s /2 + W 0 . Thus, the change in the total energy of the system to add an additional HCB is ∆E = −µ + W 0 . As long as ∆E > 0 the phase remains stable against the addition of an extra particle; therefore, the upper phase boundary of the CDW is given by the line µ = W 0 . Similarly the lower phase boundary can be determined by following the same procedure for the case when we reduce one particle from the half-filled system, for which the boundary will be given by the line µ = −W 0 .
The calculated band edges (see Appendix A) appear as dashed lines in Fig. 5 b. The nature of the different phases is further elucidated by the calculated spectrum, presented in Fig. 6. For W 0 < t a Dirac semi-metal is observed at ρ = 1/2, while for W 0 = t it is replaced by a nodal line semi-metal. Finally, for W 0 > t, the phases at ρ = 1/4 and ρ = 3/4 develop a full gap, and the dimer insulators are formed. Since the dimer insulators have no corresponding atomic limits, they appear to be more interesting to investigate. In the next section, we shed some light on the nature of these insulators by exploring their edge structure.

IV. EDGE STATES
To further explore the nature of the dimer insulators, we measure the shift in the average density ρ when we switch to open boundary conditions along the xdirection, i.e., by turning off the horizontal dashed bonds in Fig. 1 b [ 23,24]. We observe that as a function of the chemical potential µ the average density ρ remains unaltered except for ρ = 1/4, which splits into two plateaus corresponding to densities ρ 1 = 0.225 and ρ 2 = 0.275, see Fig. 7 a. Further investigation reveals that the values of ρ 1 and ρ 2 depend on the size of the system: for an 8 × 8 system ρ 1 = 0.1875 and ρ 2 = 0.3125, whereas for a 10 × 10 system ρ 1 = 0.2 and ρ 2 = 0.3. Generally speaking, we find that for a honeycomb lattice with N e edge sites (depicted by red circles in Fig. 8 c) and N s total number of sites, the plateau at ρ = 1/4 splits into ρ 1 = ρ − N e /(2N s ) and ρ 2 = ρ + N e /(2N s ).
We now argue that the splitting of the plateau under open boundary conditions corresponds to the number of in gap edge states. At 1/4-filling the dimer structure can be thought of as a series of 1D Su-Schrieffer-Heeger (SSH)-like chains (depicted by blue dashed rectangle in Fig. 3 a), where dimers of strength ∼ t are formed. The dimers in each chain are weakly coupled to each other via third-order hopping through the intermediate sites with higher on-site potential, t x ∼ t 3 /(2W 2 0 ). Now, each of these SSH-like chains will give rise to a pair of degenerate in-gap edge-states under open boundary condition, localized at the two ends of the chain. As is clear from Fig. 3 a, for a L x ×L y lattice, there are L y /2 different SSH chains weakly connected to each other via second order hopping, t y ∼ t 2 /(2W 0 ). Due to this inter-chain interaction, the degeneracy of the edge-states will be lifted and the resulting L y in-gap states will now form bands on the two edges, each with a bandwidth ∼ 4t y . The two plateaus therefore correspond to the situation when: (1), none of the edge sites are occupied; and, (2), all of the edge sites are completely occupied, beyond some critical chemical potential determined by the edge bandwidth. In our notations this reduces to ρ 1 = ρ − 1/(2L x ) and To help visualize this, in Fig. 8 a, we plot the local density profile of a 20 × 20 honeycomb lattice under open boundary condition along the x-direction, by choosing a value for chemical potential, µ = −6.5, in the lower plateau ρ 1 . We can clearly see that the lower plateau corresponds to a situation where bulk sites of the even ℓ-layers have density 0.5 each, while the edge sites have density close to zero. Contrarily, Fig. 8 b depicts the density profile corresponding to µ = −5.8, a point in the upper plateau ρ 2 . The density of the edge sites now becomes close to 1, while the density of the bulk sites remains 0.5 as before. So, at the ρ 1 (ρ 2 ) plateau, each of the edge sites have 1/2 HCB less (more) compared to the sites in the bulk. The difference of these two local density profiles is shown in Fig. 8 c, which clearly demonstrates that the transition between the two plateaus, in Fig. 7 a, corresponds to the occupation of the in-gap edge states. All this indicates towards the possible topological nature of the dimer insulator at 1/4 filling-fraction.
Next, we study the effect of the reversal of the sign of W 0 . In this case, at 1/4-filling, the dimers are formed at odd ℓ-layers and hence no dimers are split by the opening of the boundary conditions [see Fig. 3 c with black (filled) sites replaced by white (empty) sites]. On the other hand, at ρ = 3/4, all the sites residing on the odd ℓ-layers are completely occupied and the dimers are formed in the even ℓ-layers [see Fig. 3 a with white (empty) sites replaced by black (filled) sites]. Therefore, in this case, it is the dimer insulator at 3/4 filling which involves the formation of edge states. As a result, under open boundary condition along x-direction, the plateau corresponding to ρ = 3/4 splits into two plateaus (as shown in Fig.7 b) corresponding to densities ρ 1 = 0.725 and ρ 2 = 0.775, while the other parts of the ρ − µ curve remain unchanged.
Finally, we note that by opening the boundaries of the lattice along the y-direction, we end up with a honeycomb lattice with armchair edges and no dimers are split by this change. Therefore, in this case the ρ − µ curve remains unaffected by the open boundary condition for both positive and negative W 0 . Thus, edge states are manifested only along the zigzag edges in the x-direction.

V. THE TOPOLOGICAL INVARIANT
In order to confirm the topological nature of the dimer insulators from a bulk perspective, one must calculate the topological invariant for the system under periodic boundary conditions. We now turn in this direction. For a L x × L y Honeycomb lattice, the lattice points on the discrete Brillouin zone (k x , k y ) can be expressed as, Let ψ ±,± denote the normalized column eigenvectors corresponding to the energy bands E ±,± of Eq. (A1) (see Ap- and where det + (det − ) corresponds to the permanent (determinant) of the matrix. It is important to note that in the definition of the link variable, the permanent is applicable when the particles under consideration are HCBs, whereas for fermions U involves determinant. Finally, the Chern number can be calculated as, where F xy is the lattice field strength defined as, with −π < 1 i F xy (p, q) ≤ π. For the three gapped phases at densities 1/4, 1/2 and 3/4, we calculate the Chern number [25] for W 0 > t to determine the topological nature of these phases. Despite of the presence of the edge states for the 1/4 (or 3/4) dimer insulator, the Chern number turns out to be zero for all of the three insulators.
As mentioned at the end of Section IV, the edge states are observed only along the zigzag edges of the lattice under open boundary condition along the x-direction. This is related to the fact that the dimer insulator structure can be effectively described as 1D SSH-like chains stacked in a 2D lattice, connected via weak hopping. Therefore, in order to probe the topological nature of the dimer insulators, we calculate Berry phase for each k y value separately according to the formula, It turns out that the Berry phase (γ) is independent of k y (or q). Fig. 9 depicts the variation of the Berry phase as a function of W 0 /t calculated in the gapped region for three different densities ρ = 1/4, 1/2 and 3/4. We can see that for positive values of W 0 (W 0 > t), the Berry phase is quantized at π for the dimer insulator at ρ = 1/4 and remains zero for the one at density 3/4. The situation is reversed when we reverse the sign of W 0 (W 0 < −t). On the other hand, in both of these cases γ remains zero for the insulator at ρ = 1/2. This identifies the dimer insulators at ρ = 1/4 and 3/4 as weak topological insulators (WTIs) for W 0 > t and W 0 < −t respectively. The physical situation for positive W 0 (W 0 > t) is illustrated in Fig. 3, which depicts the underlying SSH chains for densities ρ = 1/4 (Fig. 3 a) and ρ = 3/4 (Fig. 3 c). Our model could be viewed as a series of 1D SSH chains, weakly-coupled via t y , each with alternating tunnelings t, t x . The Hamiltonian for a single 1D SSH chain belongs to the BDI symmetry class [6] and admits a topological invariant. Depending on the sign of W 0 , the chains for either ρ = 1/4 or ρ = 3/4 manifest their topological phase, which leads to a WTI at the corresponding density, as long as t y is weak enough.

VI. EFFECT OF INTERACTIONS
In this section, we discuss the effect of interactions between HCBs on the phase diagram and the WTI phase. First we consider NN repulsion between HCBs by adding to the Hamiltonian, Eq. (1). Fig. 10 compares the variations of the average density ρ as a function of the chemical potential µ for three different values of V 1 with periodic and open boundary conditions. The blue curves in Fig. 10 depict the average density of a 20 × 20 periodic honeycomb lattice, whereas the brown curves in the insets show the alteration of the plateau at ρ = 1/4 under open boundary condition. One can see that as we increase the NN repulsion the width of the plateaus corresponding to ρ = 1/4, 1/2 and 3/4 increases. Since the band gap in any insulating phase is determined by the width of the plateau in the ρ − µ curve, the gaps corresponding to the three above-mentioned insulating phases simply gets larger for larger NN repulsion. In other words, the insulating phases become more stable in the presence of NN repulsion. This can be understood by realizing that the introduction of NN repulsion effectively increases the energy cost of adding another particle to the insulating structures. Therefore, energy minimization forces the system to be in the insulating phases for wider ranges of the chemical potential, thus increasing the band gap of these phases. Under open boundary condition along the x-direction, for each value of V 1 , the plateau at ρ = 1/4 further splits into two plateaus corresponding to densities 0.225 and 0.275 similar to the case with zero NN repulsion. This means that the topological nature of the dimer insulator at ρ = 1/4 remains unaffected by the presence of NN repulsion.
We now argue that the WTI phase is in fact robust against any amount of NN repulsion, as elucidated by considering the spatial structure of the insulating phase. As discussed in Section III, the WTI is a dimer insulator, where each dimer is formed by a particle hopping back and forth between the two sites belonging to a NN bond aligned along x-direction. Since the dimers are situated at alternate y-levels, the particles in two neighboring dimers do not feel any repulsion, as they never reside on two NN sites. Consequently, NN repulsion cannot restrict the hopping process involved in the formation of a dimer and thus the WTI remains uninfluenced.
This argument is also valid for the WTI phase at 3/4filling, which is the particle-hole conjugate of Fig. 3 a. While NN repulsion is felt between the particle in each dimer and the particles frozen in the in-between layers (depicted by white circles in Fig. 3 a and reflecting filled sites in its particle-hole conjugate), this does not disrupt the formation of dimers. In fact this configuration is the minimum-energy configuration at density ρ = 3/4, even in presence of NN repulsion. While at this filling, the particle in each dimer encounters 2V 1 repulsion from the two occupied NN sites, this would not depend on which of the two sites of the dimer it occupies. The particle will therefore prefer to hop back and forth between these sites, rather than choosing a particular site to reside in, thereby lowering the energy of the system.
The above discussion motivates the consideration of the effect of NNN interactions, described by an additional term As depicted in Fig. 11 a, in presence of weak NNN repulsion the variations of the order parameters with the chemical potential remain almost unchanged for a periodic honeycomb lattice. With open boundary conditions the plateau corresponding to the dimer insulator at ρ = 1/4 splits into two parts (inset of Fig. 11 a) similarly to the non-interacting scenario. This demonstrates that the WTI phase is robust against small NNN repulsion values. Now, with the increase of the ratio V 2 /V 1 the hopping process of the constituent particle of each dimer gets disrupted. Consequently, beyond some critical value of this ratio, it is energetically favorable for the particle to be localized at one of the two sites of the dimer instead of hopping back and forth. This way the particles can avoid NNN repulsion felt between two neighboring dimers along the y-direction. Such a configuration is depicted in Fig. 12 b. In Fig. 11 b the dimer structure factor S D (0, π) and the structure factor S(0, π) are plotted as a function of V 2 /V 1 . We can see that with increasing value of V 2 /V 1 , the dimer structure factor decreases from 0.25 to a value close to zero, whereas the structure factor remains constant at 0.0625. Since the dimers are destroyed for larger values of V 2 , the dimer structure factor obviously decreases. Nevertheless as the particle number in each y-layer is fixed the structure factor remains unaltered. Hence, it can be concluded that the WTI transforms into a normal insulator (similar to the one in Fig. 12 b) for larger values of NNN repulsion.
As can be seen from Fig. 11 c, the NNN interactions are observed to stabilize additional insulating plateaus at fillings 1/8, 3/8, 5/8 and 7/8 for a periodic honeycomb lattice. At ρ = 1/8 only half of the dimers in Fig. 3 a are formed so that no NN or NNN repulsion is felt between two HCBs. The structure corresponding to this insulator is not unique. One of its possible structures is demonstrated in Fig. 12 a for an 8 × 8 lattice. On the other hand, at filling-fraction 3/8, half of the dimers of Fig. 3 a are occupied by an extra particle each, thus transforming a dimer into a pair of particles. Fig. 12 c depicts one of the many possible structures corresponding to this insulator. The dimers (pair of HCBs) in this insulating phase are distributed in a way such that no two dimers (pair of HCBs) are NNN of each other. Despite of the repul-sion felt by the HCBs, this is in fact the minimum energy configuration of the system at ρ = 3/8. One should note that the insulator at density 5/8 (7/8) is exactly the same as the one at ρ = 3/8 (1/8) once the even and odd layers, as well as particles and holes, are interchanged. The detailed characterization of these additional insulating plateaus would be interesting to pursue in future.

VII. CONCLUSIONS
To summarize, we have studied HCBs in a periodic honeycomb lattice with NN hopping (t) and alternating positive and negative on-site potential (W 0 ) along different y-layers, using SSE QMC supported by analytical calculations. The system reveals the existence of three insulating phases for W 0 > t: a CDW at 1/2-filling and two dimer insulators at densities 1/4 and 3/4. Depending on the on-site potential pattern, one of the dimer insulators turns out to be a WTI, with a zero Chern number but a non-trivial Berry phase. The underlying 1D STI can be effectively thought of as weakly coupled SSH chains, where the intermediate layers being either completely empty (for ρ = 1/4) or fully occupied (for ρ = 3/4). Although our study involves HCBs, it is important to note that the WTI phase persists in case of fermions as well. Since the energy bands are well-separated in the regime W 0 > t, the topological phase becomes oblivious to the exchange statistics of the constituent particles.
With recent advancements, optical lattice with ultracold atoms would be a perfect tool to actualize our model experimentally. Experiments on hexagonal lattices in this framework have already been around for a while [26][27][28][29]. The on-site potential of the lattice sites can also be tuned in these experiments, making it possible to achieve the pattern required by our model. Additionally, the measurements of Berry phase [30] as well as Chern number [31] in an optical lattice setup have also been performed successfully. Thus, we believe that our model is a promising and interesting candidate to realize weak topological insulating phase in optical lattice experiments. In addition, certain features of our model, including the band spectrum and the presence of edge states, could also be probed using driven-dissipative exciton-polariton microcavity lattices [32].
Besides the WTI phase, our model exhibits a rich phase diagram, which includes intriguing phases such as bosonic Dirac semi-metal and nodal line semi-metal among others. Since the main focus of our current work is the WTI phase, a detailed study of the other novel phases is outside the scope of this paper. It would be interesting to investigate these phases in more detail and examine how these phases are affected by the presence of off-site interaction in the system. Furthermore, it would be worthwhile to analyze the dependence of the phase diagram on the on-site repulsion U , when the HCBs are replaced by soft-core-bosons, as well as the persistence of the WTI in the U → 0 limit. In this Appendix we present details of the band structure of the Hamiltonian in Eq. (1). The absence of interaction between two different HCBs makes it possible to express the Hamiltonian in single-particle basis. Consequently, the energy spectrum obtained from this Hamiltonian is the same for a HCB and a fermion. For a unit cell containing four sites, as depicted in Fig. 1, the Hamiltonian in Eq. (1) can be expressed in the momentum space as, The four branches of the energy spectrum corresponding to this Hamiltonian are, where ǫ(k) = W 2 0 + t 2 + 4t 2 cos 2 √ 3 2 k y , η(k) = 2t W 2 0 + 4t 2 cos 2 3 2 k x cos 2 √ 3 2 k y . (A3) The Brillouin zone can be chosen to lie between − π 3 < k x < π 3 and 0 < k y < 2π √ 3 . In order to study the energy spectrum, we divide the parameter space into three regions: W 0 < t, W 0 = t and W 0 > t. For the halffilled system, we expect the lowest two bands E −,+ and E −,− to be occupied. In the region W 0 < t, a pair of Dirac nodes appear in the energy bands E −,− and E +,− at k x = 0 and two different k y values (Fig. 6 a). The positions of these nodes shift along the k y -axis as a function of W 0 and t values. This can be understood by realizing that to have nodes in the E ±,− energy bands at k x = 0 we must have ǫ(0, k y ) = η(0, k y ), which can be simplified to the solution, Therefore, for W 0 < t the half-filled system behaves as a bosonic Dirac semi-metal due to the presence of the pair of Dirac nodes at momenta (0, k ± y ). Next, we note that at k y = π √ 3 , for all possible values of k x , we have and At W 0 = t, the two point nodes thus coalesce and disappear at k x = 0, k y = π √ 3 and a line node is formed at k y = π √ 3 , as shown in Fig. 6 b. As a result, at half-filling we have a nodal line semi-metallic phase for W 0 = t.
To understand the situation for W 0 > t of the halffilled system, it is important to note that, at k y = π √ 3 the lower band E −,− attains its maximum value (−|W 0 − t|), while the upper band E +,− reaches its minimum (|W 0 − t|). The band gap at half-filling is thus given by, Hence, for W 0 > t, the half-filled system develops a gap of width 2|W 0 − t| and becomes insulating (see Fig. 6 c). The analysis of a 3/4-filled system for different values of W 0 and t can be done by comparing the energy bands E +,− and E +,+ . The top of the lower band E +,− lies at k x = π 3 , k y = 0, with the energy value whereas the bottom of the upper band E +,+ occurs at a different momentum k x = 0, k y = π √ 3 with energy, The indirect band gap is therefore expressed as, So, the system at 3/4 filling-fraction will behave as an insulator as long as ∆ 3/4 > 0, which simplifies to W 0 > t.
For W 0 = t we have E +,− ( π 3 , 0) = E +,+ (0, π √ 3 ) = 2t and the band gap becomes zero, while for W 0 < t the band gap ∆ 3/4 becomes negative. Consequently in the region W 0 ≤ t, both of the energy bands E +,− and E +,+ become partially filled and the system behaves as a bosonic semimetal. A similar analysis of the system at 1/4 fillingfraction leads to the same conclusions as density 3/4.