Topological order in the Haldane model with spin-spin on-site interactions

Ultracold atom experiments allow the study of topological insulators, such as the noninteracting Haldane model. In this work we study a generalization of the Haldane model with spin-spin on-site interactions that can be implemented on such experiments. We focus on measuring the winding number, a topological invariant, of the ground state, which we compute using a mean-field calculation that effectively captures long range correlations and a matrix product state computation in a lattice with 64 sites. Our main result is that we show how the topological phases present in the noninteracting model survive until the interactions are comparable to the kinetic energy. We also demonstrate the accuracy of our mean-field approach in efficiently capturing long-range correlations. Based on state-of-the-art ultracold atom experiments, we propose an implementation of our model that can give information about the topological phases.


Introduction
Ultracold atoms trapped in optical lattices allow the simulation of rich physical models that cannot be easily recreated in solid state materials [1,2].Great experimental breakthroughs include the simulation of the Bose-Hubbard [3] and Fermi-Hubbard models [4,5], the implementation of fermionic and bosonic lattice models with artificial gauge fields, such as the Hofstadter model [6][7][8] or topological Floquet band models [9], and, of particular relevance for this work, the realization of the Haldane model [10].The Haldane model [11] is an example of a 2D topological insulator (TI), a linear model with a nontrivial topological invariant which, while being insulating in their bulk, supports protected currents along its edges [12].An open problem in the study of topological matter, and the focus of this work, is the interplay between interactions and the topological characterization of the system.An ambitious goal would be to understand why and how ultrastrong interactions can convert a topological insulator into a topologically ordered quantum state of matter, which features degenerate, locally indistinguishable ground states and topologically nontrivial excitations.However, before achieving this goal, it would be very interesting to focus on the first steps of this transition, studying the persistence of topologically nontrivial phases in the presence of interactions.Some work has already been done in the field of interacting topological insulators, such as the Haldane-Hubbard model in the Mott insulator limit [13][14][15][16][17][18][19], the Hubbard model with nearest neighbor interactions [20], the synthetic Creutz-Hubbard model [21] or a general TI in the context of the Fractional Quantum Hall Effect [22].Our novel contribution in this work is to study the topological phases of the Haldane model with spin 1/2 fermions for a broad range of on-site interactions.This model, which represents a straightforward extension of recent experiments with noninteracting fermionic atoms [10], is studied using two methods: (i) a mean-field approach in momentum space that effectively takes into account long range correlations and (ii) a matrix product state (MPS) ansatz in a 2D lattice [23,24] with 64 sites.We have found that the topological phases which appear in the noninteracting model extend to values of the interaction which are comparable to the kinetic energy of the atoms in the lattice.For stronger interactions, the ground state becomes a Mott insulator, where, in accordance to earlier work, a nontrivial spin order may appear [13][14][15][16][17][18][19].On account of our numerical results, we propose an ultracold atom experiment that can directly simulate our model, based on previous experiments with TIs [6,7,10] and the measurement of topological properties [25].
The article is structured as follows.In section 2 we describe a generalization of the Haldane model with spin-spin on-site interactions, how topology arises in the noninteracting limit and how it can be detected by measuring the winding number.In section 3 we begin our study of interactions, introducing a mean-field variational wavefunction that exactly reproduces the noninteracting ground state.This function is a product state in momentum state and may capture the long-range correlations needed to describe a topological phase with interactions.Using a global optimization procedure, we estimate the ground state energy and wavefunction within this ansatz.We find evidence of a topological phase for nonzero interactions, as well as a cross over into a Mott insulator and double occupancy regions for strongly repulsive and strongly attractive interactions, respectively.In section 4 we study the same problem, using the MPS ansatz spanned over a 2D honeycomb lattice.Heavy simulations with up to 64 lattice sites confirm the predictions of the mean-field ansatz, including the topological phase transition and cross-overs, and show the accuracy of this simple wavefunction when estimating the ground state energy.Section 5 discusses the modifications needed to implement our model using state-of-the-art experiments [10], including state preparation and detection of topological and trivial phases.We close this work with a brief summary and discussion in section 6.

A B
Figure 1.Scheme of the honeycomb lattice with sublattices A and B in colors blue and white, respectively.Parameters t, t a and t b stand for the hopping amplitudes and p for the external phase.represents the energy imbalance between sublattices and U is the on-site interaction energy between particles with opposite spin.

The model
The Haldane model is a Hamiltonian that describes the motion of particles in a honeycomb lattice with nearest-and next-to-nearest-neighbor hopping amplitudes.The topological nature of the model arises from the existence of complex hopping amplitudes, which in our work we place in the nearest-neighbor hopping.The honeycomb lattice is composed of two triangular lattices, which we denote A and B. Each site will be able to host up to two fermionic particles in state s ∈ {↑, ↓}, which we denote with creation operators a † i,s and b † i,s for the two sites of the i-th unit cell.The full model describing this system reads where H 0,s is the kinetic energy of each fermionic species, and U is the on-site interaction between distinguishable fermions.
The kinetic part has the usual form This noninteracting model contains a first-neighbor hopping, t, which is affected by a phase that grows linearly along the lattice direction p, φ ij = p (x i − x j ).There are also two next-to-nearest neighbor hopping amplitudes t a and t b , on the respective sublattices (see figure 1); and an on-site energy imbalance between lattices .In this work we have set: t = 1, t a = −t b = 0.1 and = 0.The hopping Hamiltonian for each spin population ( 2) is a two band model that can be written as an effective magnetic field acting in momentum space.Introducing the Fourier transformed operators a k,s and b k,s , and the pseudospin operators σx = a † b+b † a, σz = a † a − b † b and σy = −i[σ z , σ x ], we find The field B(k) determines a preferred orientation of the expectation value for the singleparticle operators σ i .The ground state Ψ 0 of H 0s is a half-filled state, for which the pseudospin field, S s (k), which is proportional to the expected value of the pseudospin operators, σ k,s , satisfies The direction and amplitude of the field B (k) and of the pseudospin S s are given by where v i and w i are the vectors that connect a site with its nearest-and next-to-nearestneighbors, respectively.The half-filled ground state without interactions is therefore completely determined by the direction of the pseudospin fields S s (k) of both spin populations over the Brillouin zone.An interesting property of the hopping Hamiltonian ( 2) is that different configurations of the vector field S s can be regarded as topologically distinct.More precisely, S s (k) may be regarded as a mapping from the Brillouin zone onto the sphere, |S s | = 1.These mappings may be classified by the number of times the pseudospin S s wraps around the sphere, what is called the winding number [26].
Whenever the band wraps one or more times around the sphere, the model is said to be topologically non-trivial.This results in a nonzero Chern number and the existence of topologically protected edge currents when the model is embedded in a finite domain with boundaries or holes.
A very important open question is what happens to the topological insulator model (1) when there are active interactions between particles, U = 0.In such situations, the single-particle picture is no longer valid, but the system may still exhibit topological properties, such as a nontrivial pseudospin pattern S s (k), which gives rise to different phases, or the existence of edge states, or, in the limit of infinitely strong interactions, perhaps the appearance of true topological order.
In this work we will analyze (1) using different variational methods -mean-field theory and matrix-product states-, with the goal of classifying the approximate ground states of the model.Our main tool in doing this will be the use of the winding number, ν, because it is a topological property that can be experimentally and numerically determined, without access to the wavefunctions or parallel transport.

Mean-field
We have seen above that the ground state at half filling and without interactions is uniquely determined by the pseudospin field S s (k) over the Brillouin zone, with one particle per unit cell in momentum space.It makes sense, therefore, to approximate the ground state at nonzero interactions by a product state in momentum space This ground state is constructed by placing one particle per spin in each of the momentum states, c † k,s . These modes are characterized by the fact that 1 2 0|c k,s σk,s c † k,s |0 = S s (k), where S s (k) is a vector of norm 1/2 and our variational parameter.All the information about the ground state is thus given by two independent pseudospin fields S ↑ and S ↓ , with uniform norm.The variational energy of this wave function then reads as the functional where N is the total number of unit cells in our lattice.We remark that, in the absence of interactions U = 0 this mean-field approach is exact.The form of the variational energy suggests that, when the interactions are weak enough compared to the kinetic energy terms, |U | |B|, the topological phase will remain unaffected, and the ground state will adopt a spin order determined by the effective magnetic field S k,↓ = S k,↑ ∝ −B k .For strong repulsive interactions U |B|, the mean field model suggests that the vector field develops an antiferromagnetic order, where all spin components point along the Z direction, S z ↓ = −S z ↑ .If we inspect the expression for the pseudospin we realize that this implies a phase separation between spin components, whereby different spin orientations sit on different sublattices, i.e particles ↑ and ↓ are placed on sublattices A and B, respectively, or vice versa.Moreover, because the transverse components are zero S x,y 0, the winding number vanishes in this phase.Similarly, for strong interactions, U −|B|, we expect a ferromagnetic order to develop, which amounts to having bunching of particles on the same lattice sites.
We have optimized numerically the variational energy for lattices with up to 400 sites, using a linear search algorithm that starts with a random distribution of vector fields.For each of the variational ground states, we have computed both a discrete approximation to the winding number (8), as well as an order parameter that detects phase separation and bunching.In the figure 2a we plot the winding number of the mean-field ground state of a lattice with 20 × 20 unit cells.We see that the topological phases survive until the strength of the interactions, |U |, is around three times bigger than the kinetic energy, t.In 2b we plot the order parameter coming from the expression of the variational mean-field energy, In regions which correspond to a topological phase, Π is strictly zero, but it breaks into two distinct phases when the topological order disappears.It then becomes negative in the case U > 0, when the interactions are repulsive, which means that particles of opposite spin sit on different sublattices, therefore signaling a charge density wave.When interactions are attractive, U < 0, the expectation value Π becomes positive, meaning that all particles sit on the same sublattice.

Matrix product state simulations
As a check of the validity of our mean-field approach we have used an MPS ansatz to compute the ground state of the interacting model and the winding number.MPS is a powerful method to compute the ground state of quantum lattice models that relies on the low scaling of the entanglement entropy in ground states in 1D and 2D.It describes a quantum state as a product of tensors on every lattice site where the dimension of the i j indices are the physical dimensions of the system at the j-th lattice site, and A i j [j] are the tensors corresponding to the same site, whose dimensions depend on the bipartite entanglement of the state around j.This ansatz is an exact representation of every quantum state for sufficiently large tensors.However, for practical applications one usually sets a maximum size, called bond dimension, χ, for every A i j [j], discarding configurations that contribute the least to total the ground state.
In our calculations we have constructed an MPS by mapping a lattice with N × N unit cells [23,24] to two consecutive 1D systems in a "snake" configuration, with a total of 2 × 2 × N 2 sites, one for each spin component [cf.figure 3].We have focused our efforts in diagonalizing a 4 × 4 lattice using a maximal bond dimension of χ = 120.Lattices smaller than 4 × 4 unit cells did not output correct results for the winding number, already in the noninteracting analytical solution.Bigger lattices demanded a too large bond dimension that exceeded our computational resources.
As hinted in our mean-field computation (figure 2), the nontrivial topological phases extend until the interactions are comparable to the kinetic energy, |U | ∼ 3t.For stronger interactions, the term U i∈A,B n i↑ n i↓ becomes dominant and nonlocal expectation values ( c † i c j , i = j) approach zero, where different ferro-or antiferromagnetic orders appear, which correspond to different distributions of the spin components in the two sublattices.We identified the double occupancy as an appropriate order parameter to measure these distribution effects In figure 4 we show the main results of our MPS simulations.Due to the big computational effort needed to get the MPS ground state at every point we have chosen to compute them only along some relevant lines in phase space, denoted by dashed lines in figure 4a.The winding numbers along these lines are shown in figures 4b and 4d.These winding numbers have been computed by extracting the expectation values c † i c j from the MPS and then Fourier transforming them into bigger lattices of 8×8, 12×12 or 16×16 unit cells in momentum space.This Fourier interpolation of the pseudospin field, allows us to make a more accurate determination of the winding number in momentum space.We only keep simulation results that have converged under this interpolation procedure.The selected values compare favorably with the mean-field predictions for large lattices 20 × 20, away from the phase transition.
Near the topological phase transitions, the winding number fails to converge in the MPS simulations [cf.gray shaded areas in figure 4b,d].However, this is not surprising.As a mean-field shows, it is difficult to determine accurately the winding number with just 4×4 pseudospins [cf.yellow line in figure 4b,d].Around these regions, the gap closes, and it is therefore very difficult for an algorithm to determine whether the pseudospin field maps into a closed or an open surface, respectively figures 2b and c.
In figure 4c we plot the double occupancy of the ground state for the MPS computation and the mean-field at the line φ/π = 0.5.The expectation value of this observable should approach 1/4 when the interactions have no effect on the ground state and the two topological insulators are uncoupled (2).However, D should become 1/2 in the strongly attractive phases and 0 in the strongly repulsive phases, where particles with opposite spin are either grouped or separated, respectively.We observe a good agreement between both simulations, specially at the plateau D = 1/4, where our results show a non-trivial topological phase.
Finally, figure 4e shows the ground state energy of the mean-field simulation, together with that of the MPS computation.The mean-field energy is in a remarkably good agreement with MPS for all parameters, and it is even slightly below for |U | ≤ |t|.While surprising, this can be justified by the fact that our momentum-space meanfield ansatz is nonlocal, and contains long range correlations that can be challenging to capture with the MPS, one-dimensional ansatz.

Experimental setup
Ultracold atoms trapped in optical lattices provide an excellent tool to study the topological properties of the Haldane model, in any of its flavors.This has been demonstrated by the recent experimental work in [10].A crucial step to realize topological order in the lattice, is to implement an effective gauge field on the phases of the hopping amplitudes when a particle moves through the lattice.There are two essential ways in which this is done.The first one is through laser assisted hopping [6,7], in which the lattice is divided into two sublattices which host atoms in different hyperfine states.A laser beam then transfers atoms between sublattices, while imparting a complex phase.A second way to make complex hoppings is to periodically shake the lattice [10,27].This creates an effective Hamiltonian such that particles acquire a complex phase as they hop between sites.This second method was the one used to recreate the Haldane model in reference [10].In this work the authors also measured the Berry phases acquired by the atom while moving through the Brillouin zone, an evidence of the topological order.
The model implemented by Jotzu et al. is topologically equivalent to the Hamiltonian that we have studied (1) in absence of interactions: both are related by a smooth, unitary deformation of the pseudospin field that does not change the topological invariants.Our goal would be then to add interactions to the experiment, so as to observe the survival of topological invariants and the phase transitions that have been evidenced in previous sections.In [10] the authors used 40 K atoms in two internal hyperfine states, |F, m F = |9/2, −9/2 and |9/2, −7/2 , but their interactions were suppressed using Feshbach resonances so that the effective model was a noninteracting version of (1).The introduction of interactions would "simply" require adiabatically moving away from the Feshbach resonance.As we have shown in previous sections, the ground state in the interacting regime |U | < |t|, is smoothly connected to U = 0 many-body state.We therefore expect that adiabatically increasing the interactions should not introduce significant heating.Regarding the experimentally achievable parameters, the authors report to first neighbor hoppings to be t F N /h = [−746(81), −527 (17), −527 (17)]Hz and the second neighbor hoppings t SN /h = [14, 14, 61]Hz, where h is Planck's constant.The ratio t SN /t F N ∼ 0.05 is comparable to our simulations, and the experimental parameters set a limit of interaction strength of about U ∼ 3t F N ∼ h × 1.5kHz.
Finally, to measure the winding of the system's ground state one has to measure first the pseudospin field, S (k).This can be done using TOF images [28,29], which reproduce the momentum density distributions in both sublattices, n a,b (k), and of every spin population.The winding number ( 8) is then straightforward to compute from the pseudospin field.

Summary and discussion
In this work we have studied the topological phase diagram of a topological insulator with spin-spin on-site interactions.To identify these topologically non-trivial phases we have estimated the winding number.Our main result is the observation that the topological phases from the noninteracting topological insulator model are extremely robust against the introduction of repulsive and attractive interactions, disappearing when interactions and hopping are comparable, |U |/t ∼ 3.This prediction is supported by mean-field and MPS calculations, both of which show a remarkably good agreement in the detection of the topological phases.We have further characterized the topological phase transitions using the double occupancy (14) and the product of the mean pseudospin fields in the z direction.Both order parameters showed some kind of discontinuity when the winding number changed, giving us a strong signal of a transition into a topologically trivial phase.We have also made a proposal to study these phase transitions using an existing, state-of-the-art implementation of the Haldane model [10], with the only requirement that interactions are reintroduced in the experiment.Finally, the remarkable agreement between our mean-field ansatz and the MPS simulations hint at the capacity of the former to capture long-range correlations and entanglement.Further work to improve the mean-field ansatz, and study time-dependent variations are in order.

Figure 2 .
Figure 2. Mean-field results of a lattice with 20×20 unit cells.(a) Topological phases in the mean-field simulation signaled by the winding number.(b) and (c) Map of the pseudospin field onto the Bloch sphere at points in phase space shown in (a).(b) represents a pseudospin configuration that wraps the whole sphere, therefore having |ν b | = 1; while the pseudospin configuration in (c) does not wrap around a whole sphere, thus ν c = 0. (d) Order parameter Π (12).

Figure 3 .
Figure3.Scheme of our "snake" MPS spanning over a 4 × 4 unit cell honeycomb lattice with spin.The order in which the physical indices are covered is indicated by the thick green line.The MPS goes first through the lattice with spin up and then through the lattice with spin down.

Figure 4 .
Figure 4. Results from the MPS computation with a 4 × 4 unit cell lattice at χ = 120.(a) Lines in phase space at which we have computed the winding number of the ground state.(b) Winding number at φ/π = 0.5.The lines represent mean-field results with different lattice sizes while points represent our MPS computation.The shaded areas are regions where the winding computation does not converge.(c) Double occupancy of the points at φ/π = 0.5.(d) Winding number at the line U/t = 1.The double occupancy at this line is constant, D = 0.25.The symbols are the same as in figure (b).(e) Energy of the mean field computation compared to the MPS energies.