Winding number order in the Haldane model with interactions

We study the Haldane model with nearest-neighbor interactions. This model is physically motivated by the associated ultracold atoms implementation. We show that the topological phase of the interacting model can be characterized by a physically observable winding number. The robustness of this number extends well beyond the topological insulator phase towards attractive and repulsive interactions that are comparable to the kinetic energy scale of the model. We identify and characterize the relevant phases of the model.


Introduction
The integer quantum Hall effect initially [1] and the more extended set of topological insulators more recently [2], represent a family of many-body systems with exotic transport properties. These properties originate from the symmetries of the system and a robust and hidden topology that appears in its bands [3]. More precisely, if such two-dimensional non-interacting models are embedded in a torus and diagonalized in momentum space, they give rise to a collection of bands, k n ( ) e , and eigenstates, k n ( ) y for integer n, which are characterized by the Chern numbers When the Fermi energy of the system lays in the gap between two bands and the total Chern number of the occupied bands is non-zero, then its quantum phase corresponds to a topological insulator. Its key physical property is that as an insulator it is non-conducting in the bulk, but it has robust conducting states that live on the boundary of the lattice due to the non-trivial band topology. These edge states are well known for their robustness under perturbations: the non-local nature of the topological phase that supports them, allows edge states to survive under extreme conditions of distortion, impurities and other external perturbations. This resilience of topological edge states is highly favorable when considering their realization and detection in the laboratory. An important open question in topological insulators is what happens to their topological properties in the presence of interactions. One special situation is when we force the energy bands of the insulator to be extremely narrow, having k n ( ) e almost independent of k, and introduce a strong repulsion or attraction. A paradigmatic case is the fractional quantum Hall effect [4], where states with fractional filling support anyonic excitations that could be used for quantum computing [5,6].
A complementary situation is when there are obvious dominant energy scales and the interaction term is comparable to the kinetic term that originates the topological phase. This regime, which has been addressed in earlier theoretical works [7,8], is relevant to what could be experimentally achieved in present quantum simulations of topological insulator with ultracold atoms. Such simulations have been recently realized in the laboratory for the Harper-Hofstadter Hamiltonian [9][10][11] and the Haldane model [12].
In this work we employ the Haldane model to study the role of interactions in the survival or destruction of the topological phases. This model [13] is a particular instance of a topological insulator or anomalous quantum Hall effect. It is realized on a honeycomb lattice with complex hopping amplitudes, as shown in figure 1. We rely on a variant of the model introduced in [14], which is amenable to being implemented using ultracold atoms in optical lattices. It was also shown in [15] that it can introduce new regimes, such as a topological semimetal. Its physical implementation naturally gives rise to nearest-neighbor interactions [16]. This is exactly the setting we simulate here.
Our study relies on an established tool for characterizing non-interacting topological phases, the so called winding number. This is a mathematical object which can be constructed when the unit cell of the lattice contains a pseudospin degree of freedom (see section 2) and the eigenstates of the model are represented by unique Bloch vectors where 0 y is the ground state of a momentum-dependent Hamiltonian H k 2 2  Î´. For the non-interacting case the winding number is identical to the Chern number [17]. It has been shown that ν is an observable that can be directly measured in ultracold atoms experiments [14]. This approach has been generalized to topological superconductors [18] and to composite problems (i.e. unit cell dimensionality 2 n ), where it still signals the existence of topological phase [17]. In this respect, we regard the winding number as a global order parameter in its own right that can be measured experimentally.
The main result of this work is that the topological phase that originates in the non-interacting topological insulator regime extends towards both attractive and repulsive interactions, and is at all times detectable through the winding number (2). For very repulsive interactions a charge density wave (CDW) order develops, as already indicated by exact diagonalizations in finite lattices with a fixed number of particles per site [7,8]. For very attractive interactions the lattice fills up, as we work at zero chemical potential regime. This phenomenology can be observed with a very straightforward mean-field theory that works with the vector field (3) as an order parameter. This theory shows both the survival of the winding number as well as the phase transitions into the topological phase and into the CDW. These mean-field predictions are confirmed by matrix-product-states (MPSs) ansatz which is optimized using DMRG for ground state computations with up to 50 lattice sites, a size that doubles what is currently achieved through exact diagonalizations. Finally, the destruction of the  topological insulator that takes place for large attractive interactions is modeled using an extension of the meanfield theory by Poletti et al [19], which shows the appearance of a superconducting phase in such regimes.
The structure of this paper is as follows. In section 2 we introduce the form of the Haldane model used here. Section 3 introduces a mean-field theory that works with the spin field (3) as order parameter and thus gives at all times a proper value of the winding number. Section 4 introduces MPSs, the variational ansatz that we use to get numerical estimates of the ground state properties in finite-size problems. Since both our mean field and our DMRG fail for very attractive interactions, section 5 introduces a BCS ansatz that we can use to describe that regime, generalizing the ideas in [19]. Section 6 introduces the main results, describing the phase transitions that take place in this model and how they are characterized through the winding number, the mean field and DMRG, or the BCS theory. Finally, section 7 discusses further implications of this work and possible outlooks.

The model
We work with a variant of the Haldane model that introduces flux as a phase on the nearest neighbor hopping of a honeycomb lattice, as explained in [14]. The model is extended to consider also nearest-neighbor interactions that take place between the different sublattices of the Haldane model. Overall h.c.
and a nearest-neighbor interaction [14]; the intralattice hopping parameter in the A B , sublattice is t a b , , and the nearest-neighbor interaction strength is U. Finally, we include a possible lattice imbalance, ò, that will be made zero in this work.
The hopping part of the Hamiltonian, t and t a b , , implements a Dirac-type Hamiltonian with two distinct singularities in the Brillouin zone (BZ). The t term gives rise to the kinetic terms of the effective model, while the t a , t b terms implement a momentum dependent mass. The combination of both terms can be rewritten as a pseudospin model are the energies of the two bands at the given quasimomentum k. At half filling, that is one atom per unit cell, the ground state is obtained by placing one particle in the lowest band, which is a single particle state with an orientation of the pseudospin a b k , As widely explained in [14], this model could be implemented using two separate optical lattices that store atoms in two different internal states, one for sublattice A and a different state for sublattice B. The nearest neighbor hopping between sublattices would then be implemented by laser assisted tunneling between internal states. As a consequence, if the photon-assisted tunneling is implemented by a Raman process, the tunneling amplitude t xy will carry a phase that is proportional to the phases of the laser in the two lattice centers, x and y. If we use a single laser propagating along the direction and momentum p, In our setup we have chosen one fixed direction which is perpendicular to one of the lattice edges, as sketched in figure 1(b). The complex hoppings may then be written -. Note also that by using two overlapping sublattices we also allow for the presence of strong nearest-neighbor interaction between sublattices, which is proportional to the square of the overlap between lattice wavepackets (see [16]) and can be rather large.
More recently the group of Prof Esslinger has demonstrated an alternative implementation of the Haldane model using periodically shaken optical superlattices. There, a square lattice is first distorted into a brick lattice by the introduction of an additional laser, creating a model that behaves like graphene. This lattice is then periodically shaken so that certain hopping amplitudes are enhanced and others are reduced, implementing the Haldane model [12]. A clever combination of the shaking directions and phases allows pretty much the same flexibility as our proposal above. The only drawback is that nearest neighbor interactions are weaker in the shaken lattice than in the overlapping lattice setup that we introduced. This impedes further enhancements or a significant reduction of the hopping rates. Hence, the translation of our results to that particular experiment is only useful in the weakly interacting regime, U t | |  .
3. Mean field theory 3.1. General methodology Our approach towards a mean-field study of the interacting Haldane model builds on our knowledge of the noninteracting solution: in the non-interacting case the ground state wavefunction is exactly parameterized by the pseudospin orientation S k ( ) of the fermion that populates the k quasimomentum in the BZ. Consequently, we may write the ground state wavefunction at half filling as In the interacting case we expect the nearest-neighbor repulsion to have a moderate effect on the fermion distribution, more or less preserving the Fermi sea and the topological phase associated to the vector field S. In other words, we will estimate variationally the ground state properties of the interacting model using the wavefunction (8) and the orientations of the pseudospins as variational parameters. This mean-field approach is rather unique. The variational ansatz for GS | ñ is a product state over the Hilbert spaces of each quasimomentum mode, described by a quasicontinuous field S k ( ) that is defined on N equispaced points of the BZ. Thus, even though our Hilbert space contains 2 N parameters, our semiclassical description only employs N 2 values associated to the orientations of the S on the N points in momentum space. Because the complexity of this description is linear in the problem size, we can afford to treat rather large lattices with it, exceeding those of the MPS simulations by an order of magnitude. Nevertheless, the outcome is numerically the same for these larger lattices, as we will show below, and the optimization becomes slower.
It is important to remark that this method becomes exact in the U=0 limit. Moreover, the mean-field variational wavefunction (8) continuously interpolates between the topological phase and other phases that are to be found, such as the CDW that we expect in the limit U→+¥, where atoms localize in either of the sublattices. This variational ansatz, however, does not capture other orders, such as a superconducting phase that should arise for large negative U t 2| | -. To analyze this phase we employ below a different ansatz. with three displacements v 0,1,2 that connect one site to its neighboring cells. The Fourier transform of this Hamiltonian becomes

Spin model
since the central exponential, which has been summed over s, has led to a (·) d that enforces the conservation of momentum. Finally, we have the weight function The second term induces CDW order, forcing all atoms to be on one sublattice or the other. The third term scatters particles to different momenta. Both terms counteract the effect of the magnetic field p S p ( ) ( ) e that tries to enforce the topological phase.

MPSs ground state 4.1. General methodology
To confirm the mean-field predictions one could perform exact diagonalization for computing the ground state wavefunction and evaluating the winding number. However, this becomes unfeasible for two combined reasons.
The first one is that we need a minimum size of the lattice to detect the winding number. If we perform a simulation with L L 1 2 unit cells in position space, using periodic boundary conditions, it would amount to sampling L L 1 2 points in the BZ. A quick numerical calculation revealed that this sampling is insufficient to capture the winding number accurately if L 3  . This means that the smallest simulation that could conceivable reproduce the non-interacting result would be a 4×4 cells lattice with 32 sites.
The second reason is that given that baseline, performing simulations at half-filling becomes unfeasible using an ordinary computer and algorithms. Already storing the ground-state wavefunction for the 4×4 lattice demands several gigabytes and we were not able to do ground state computations without resorting to state-ofthe-art diagonalization software and a supercomputer with parallelized Lanczos. Note indeed that already the largest clusters that are found in the literature [8] at half filling have sizes around 24 sites or 3×4 cells, which are insufficient for the winding number computation.
A powerful alternative to exact diagonalizations are tensor-network state methods. In particular, we have used MPSs to represent the ground state and a two-site DMRG-type algorithm [20,21] to compute the minimum energy state within that variational ansatz. While the MPS is a variational form of a many-body wavefunction that is built by contracting tensors (matrices) in a 1D scheme, it has been used also for efficiently approximating 2D structures [21][22][23][24]. Using the ordering of lattice sites and tensors in figure 2 our ground state is written in the form is the occupation number of the mth lattice site. Using a rather straightforward optimization procedure based on local updates of the tensors, it is possible to minimize the energy, computing the optimal tensors áY Yñ áY Yñ

Winding number from MPS
After computing the ground state we still have to recover the winding number as measured in a time-of-flight experiment [14]. In such experiments, the expansion of the atoms in the lattice makes them adopt a Fourier transform of their original wavefunctions [25],   already for very small samplings [14]. Moreover, for our theoretical studies with small lattices, it suffices to have M of the order of the lattice size itself. Indeed, when the winding number converges, it does so for M N ( )  points.

MPS technical remarks
There are several important technical remarks regarding the MPS simulations. The first one concerns the use of conserved quantities in the numerical simulation. For instance, we are free to impose a condition on the total number of particles in the lattice. We have chosen not to do this, looking for the ground state of the full Hamiltonian, which is equivalent to minimizing the free energy at zero chemical potential. This means that in these simulations an increase or decrease of the filling is possible. As we will see below, such a change does not affect the identification of a topological phase through the observable winding number-a signature of the robustness of this quantity, as shown already in [14].
The second remark regards the size and type of lattice. MPS simulations were done for lattices with 4×4, 4×5 and 5×5 cells (32, 40 and 50 sites), the latter showing the clearest signal. We explicitly use open boundary conditions for the MPS wavefunction for several reasons. First of all, open boundaries are the ones that best reflect a potential experiment with ultracold atoms. Second, we wish to show that the proposed measurement of the winding number is a robust and powerful technique that does not restrict to abstract geometries, such as tori and spheres. Third, since we already have periodic boundaries in the mean field theory, agreement between this theory and the open boundary conditions MPS represents a strong signature that our approach is correct.
Finally, we have to comment on the size of the MPS, the so called bond dimension. This dimension is the size of the matrices in (16) and it relates to the maximum amount of entanglement that is available in a bipartition of the state. The simulations that we show here are done with bond dimension 100 c = . This restriction is based on the need to scan in detail the whole parameter space and the use of long-range interactions induced by the 2D-to-1D mapping. However, in this particular study, for the observables that we computed, including density correlations and the full winding number, we have verified that convergence starts already at very low bond dimensions, such as 30 for the 5×5 lattice. This is one further evidence of the utility of DMRG or MPS snake ansätze for small two-dimensional systems, a fact well known in the literature [27,28]

BCS ansatz
We can further analyze our system, in the case of attractive (U 0 < ) interactions, by means of a BCS ansatz. In order to do so, and to facilitate a comparison, we will follow the procedure used in [19], which deals with a similar tight-binding Hamiltonian. Up to trivial factors, the Hamiltonian presented there corresponds to the 0 F = case in equation (7). The BCS ansatz is based on a Wick expansion of the interaction term in momentum space which gives -¢ ¢ is computed through a self-consistent calculation. The order parameters which characterize the phase are nearest-neighbor two particle pairing correlators s s v s s v s s v where v m are again the directions spanned by the three nearest neighbors. Following [19], we rewrite the order parameter as w u , where the u m are a basis for the 3  cyclical permutation group: A non-zero value of any w m coefficient (that is, of | | d ) signals a superfluid phase, while the particular coefficient renders information about the symmetry of that phase. As we will see below, the results that we obtain in the attractive regime are consistent with those of Poletti et al [19]. Along the vertical axis we find that the topological phase described by the winding number disappears for repulsive interactions around U t . As shown in figure 3(c), this transition is due to a spontaneous symmetry breaking of the z s expectation values, signaling the transition into a CDW where particles are polarized into one sublattice or the other. To confirm the phase separation or CDW order, we have studied a similar order parameter in the MPS simulation. In particular, we have computed the expectation value n n a b s s á ñ, which measures the coexistence of atoms in neighboring sites in the honeycomb lattice. For large U this value decreases continuously down to zero, confirming the separation of species in the lattice. Interestingly, from the point of view of the CDW order, this looks like a cross-over, but from the point of view of the winding number this looks like sharp phase transition into a disordered regime. Our MPS simulations cannot resolve if this is merely an artefact of poor resolution on the numerical side (for instance because the norm of v becomes so small that our computation of the winding number is inaccurate). We have also studied what happens for negative or attractive interactions. When U 0 < , the ground state configuration at zero chemical potential has a filling fraction larger than 1/2, that is more than one particle per unit cell. Despite this enlarged filling, the ordered phase still persists until U 1 = -, as evidenced by the winding number (figures 3(a) and (b)). At this critical interaction the lattice becomes perfectly filled ( figure 3(d)), forming a Mott insulator with two particles per site. This trivial configuration cannot be reproduced with the mean-field, because that wavefunction assumed half-filling. In order to study this region of strong attractions we have to use the BCS ansatz developed above (5). Figures 4(a)-(d) show the outcome of that ansatz. They reveal that the superconductor phase that was predicted for the honeycomb lattice [19] exists also in the full Haldane model. Most importantly, the symmetry of that phase is strongly dependent on the breaking of time-reversal symmetry in the non-interacting limit, i.e. of the topological phase of the U=0 Hamiltonian.

Discussion
Our study reveals that the topological insulator phase survives for a wide range of interaction. This phase is faithfully detected by the winding number. Given that earlier exact diagonalizations computed the Chern number in selected regimes of parameters [7,8], this leads us to conclude that the winding number can also reproduce the Chern number in interacting systems with moderate correlations. The ideas put forward in this work have very straightforward extensions to other models, including other topological insulators and composite systems [17], and other types of interactions, such as on-site Hubbard terms or long-range interactions. In all these systems it would be interesting to see whether the winding number may act as a precursor of other strongly correlated phases.
One may also wonder about the applicability of our results to the recent beautiful experiment demonstrating the Haldane model in an optical lattice [12]. This experiment uses laser assisted tunneling to implement t a b , , relying on the original lattice to supply t. Compared with our original proposal [14], it has the problem that the small overlap between neighboring sites would lead to a small value of U. In this case it would be more advantageous to implement other types of interactions, such as relying on two-level atoms and implementing on-site repulsion or attraction. Such models fall out of the scope of this work. Nevertheless, they could be studied using a generalization of our MPS and mean-field methods above, combined with the study of partial winding numbers [17] and see how they relate to the global topological structure of the bands.
Finally, we would like to emphasize the mean-field ansatz (8), which could be of interest to other contexts and models. Such a momentum-representation mean-field theory differs from other mean-field theories that have been developed in position space [29]. They connect to long-range interaction classical spin models that have been long analyzed in the literature. The connection between these models, the symmetries of the underlying topological insulator and how these are broken by the interaction terms could be a profitable avenue to understand the nature of the phase transitions that have been found in this study.