# Topological transitions in a model with

Particle-Hole symmetry, Pancharatnam-Berry Curvature and Dirac Points

###### Abstract

We study the topology and geometry of a fermionic model on the honeycomb lattice with spin-dependent hopping which breaks the time-reversal and charge-conjugation symmetries but preserves their composition. We show that in such a case the Zak phases are topological invariants that characterize the semi-metallic state at half-filling and determine the edge state structure. As the strength of the spin-dependent hopping varies, the model shows several Lifshitz transitions corresponding to creation and merging of multiple Dirac points. We discuss the possible realization of this model in cold atom systems and propose experimental signals of detecting the Dirac points and Pancharathnam-Berry curvature of the bands.

## I Introduction

There is much current interest in the topologically non-trivial phases of fermionic bands, since such phases have the potential of realizing quasi-particles with fractional quantum numbers and statistics. The topological phases of two dimensional insulators can be classified in terms of the sum of the Chern numbers of the occupied bands. The Chern numbers are the Pancharathnam-Berry (PB) curvature field integrated over the Brillioun zone. They are topological invariants as they do not change under smooth changes of the hamiltonian parameters that do not close the gap. The presence or absence of certain discrete symmetries leads to a more detailed classification of the possible phases schnyder2008; kitaev2009; lu2012.

The PB curvature characterizes how the phases of the single-particle wave-functions twist over the Brillouin zone. The semi-classical Sundaram-Niu equations sundaram1999; haldane2004 provide a nice physical interpretation of PB curvature as magnetic field in momentum space. In the presence of an external force, it induces the so called anomalous velocity perpendicular to the external force. This leads to a Hall conductance even in the absence of an external magnetic field, a phenomenon dubbed the anomalous Hall effect. The value of the Hall conductance is equal to the PB curvature field integrated over the occupied states. When the highest occupied energy band is partially filled the system is called a topological Fermi liquid haldane2004. When the Fermi level lies in the gap, the Hall conductance is the sum of the Chern numbers of the occupied bands and is hence quantized. Such systems are called Chern insulators.

The PB curvature becomes singular at the points where two bands touch due to the multivaluedness of the wave-function like at the Dirac points. Since the PB curvature is singular at the Dirac points, we refer to them as Dirac punctures (DP) in the Brillouin zone. The Chern-number of each band becomes ill-defined but the sum of the Chern numbers of the two bands is still well defined. If two bands touch and detach at DPs as some hamiltonian parameter is varied, the Chern number can redistribute between the bands leading to a transition between two different topological phases. Several models with Dirac points hassan2013b; montambaux2009; hasegawa2012; zhu2007 and non-zero PB curvatures jo2012; tarruell2012; aidelsburger2011; jaksch2003; sorensen2005; gerbier2010; cooper2011; cooper2011b have been studied.

In this paper we investigate the non-interacting limit of a model, the Kitaev-Hubbard model (KHUB) duan2003 presented by the hamiltonian,

(1) |

describes a system of fermions on a honeycomb lattice. labels the sites of a honeycomb lattice, the spin and the nearest neighbour link. There is a spin-independent hopping term with strength , a time-reversal breaking, spin-dependent hopping with strength and an onsite repulsion term of strength .

The KHUB model was originally proposed duan2003 as a way to realize Kitaev’s honeycomb model kitaev2006 in the large limit, and half filling. The model has been studied earlier in the regime hassan2013; hassan2013b. At half filling there is a chiral semi-metal phase at small and a Mott transition into an Algebraic Spin Liquid hassan2013 at large . At quarter filling there is a topological Fermi liquid phase and transition into Chern insulator phase hassan2013b. The transition occurs at at . As is increased, the critical value of decreases to at .

We investigate the non-interacting limit, in this work. There are topological transitions in this regime corresponding to creation and merging of DPs which occur at . It has been shown that these features persist even at non-zero values of jeanpaul2013. We study, in detail, the topology of the phases and the transition between them. We also discuss the proposed experimental realization of the model in cold atom systems duan2003 and point out possible experimental signals of the topological features.

In section II, we review the effects of the discrete symmetries on the PB curvature. We show that the PB phase vanishes for insulators with particle-hole symmetry at half filling. We then examine the topological features of bands with DPs and particle-hole symmetry in detail in section III. Here we explore the non-interacting hamiltonian of KHUB with periodic boundary condition (PBC) and show the topological transitions that characterize them. Further we study the model with open boundary conditions (OBC) and discuss the edge currents. In section IV we discuss the possible experimental realizations of the model and its topological features. Following the scheme proposed by Duan et. al.duan2003 we present a derivation of the spin-dependent periodic potential that leads to spin-dependent hopping terms in the hamiltonian. We show that the multiple DPs can be probed using Bloch-Zener oscillations and we present a semi-classical analysis to show how the PB curvature of the occupied bands can be detected. We summarise and discuss our results in section V.

## Ii Discrete symmetries and the PB curvature

As mentioned earlier, topological phases of insulators have been classified according to the presence or absence of certain discrete symmetries schnyder2008; kitaev2009; lu2012, namely time reversal symmetry (TRS), charge conjugation symmetry (CCS) and their composition which we call particle-hole symmetry (PHS). In this section, these symmetries are briefly reviewed and the constraints that are imposed on the PB curvatures are examined for number-conserving, non-interacting, 2-dimensional fermionic systems. The hamiltonian for such systems can be written as

(2) |

where goes over the Brillouin zone of a 2-dimensional Bravais lattice, is the single-particle hamiltonian and label the sublattice and spin indices. The single-particle hamiltonian of KHUB in Eq.(1), has corresponding to two sublattice and two spin orbitals in every unit cell. Setting , it can be written as

(3) |

where and are matrices in the sublattice and spin space respectively.

(4) | |||||

(5) |

where and . We denote the spectrum as, , ,

(6) |

In terms of these single particle eigen-functions, the PB vector potential, and curvature, are given by

(7) |

From the PB curvature, the Chern number

(8) |

can be computed. Now we discuss effect of the discrete symmetries on the energy bands and the PB curvature, one by one.

### ii.1 Time-reversal symmetry (TRS)

The time-reversal transformation replaces particles (holes) with momentum by particles (holes) with momentum . It is an anti-unitary transformation which we denote by ,

(9) |

All transition amplitudes are invariant under this transformation if there is a unitary matrix with such that,

(10) |

Under the time-reversal transformation, . Thus if it is a symmetry, then the Chern numbers, are all 0.

The KHUB satisfies the condition

(11) |

Thus for time-reversal symmetry to hold the condition in Eq.(10) needs to be satisfied for finite , implying that the matrix has to anti-commute with all the three Pauli matrices. Since such a matrix does not exist for any , the model in general is not TRS. But at two special points, with and with , where anti-commutes with and the model preserves TRS.

### ii.2 Charge conjugation symmetry (CCS)

The charge-conjugation transformation replaces particles with momentum by holes with momentum and vice-versa. It is unitary transformation in the many-body Hilbert space that we denote by ,

(12) |

All transition amplitudes are invariant under this transformation if there is a unitary matrix with such that,

(13) |

If the system has CCS, then all the single particle energies come in pairs with and . corresponds to the band index with negative of the energy of . The positive and negative energy bands have opposite Chern numbers. From Eq.(11) it follows that the KHUB has CCS only at (with ) and at (with ).

### ii.3 Particle-hole symmetry (PHS)

The particle-hole transformation which we refer to as the composition replaces particles with momentum by holes with momentum and vice versa. Note that the nomenclature is not uniform in the literature. For example Schnyder et. al. schnyder2008 refer to what we call CCS as the particle-hole symmetry and what we call PHS by “chiral” or “sublattice” symmetry.

The particle-hole transformation defined above is a symmetry if

(14) |

The KHUB has PHS with at all values of . This symmetry is very common in condensed matter systems. It occurs in all bipartite lattices where the fermion hopping is only from one sublattice to the other. PHS implies that all the single particle levels come in pairs with and . The sum of the PB curvature over the positive and negative energy bands are equal to zero individually as the sum of the PB curvature over all the bands is identically zero. Hence if there is PHS, then

(15) |

Thus the total PB curvature vanishes for insulators with PHS at half-filling.

## Iii Topology of bands with DP and PHS

In this section, we consider the case when the highest negative energy band and the lowest positive energy band touch at DPs, where is an even integer. We denote the DPs by . We will show that for systems with PHS at half filling the PB curvature is given by,

(16) |

where is the PB flux passing through . Consequently, the Zak phaseszak1989; saket2013; delplace2011, defined as

(17) |

are topological invariants. These quantities are independent of the contour , provided it does not cross a DP. They are completely determined by the position and indices of the DPs. DPs lead to non-dispersive edge modes and we will show that the wave-vectors of these edge modes are determined by the Zak phases of the loops that wind around the Brillouin zone.

PHS implies that in the basis where is diagonal, the single-particle hamiltonian is of the form,

(18) |

In general the two blocks defined above can have different dimensions, say and , for example a bipartite lattice with different number of and lattice sites. However if , there will be zero eigenvalues at every , i.e. flat bands. While this may have interesting effects, we concentrate on the case so that we have an even number of bands, . It is then convenient to replace the index by a pair .

Using the fact that every matrix admits a singular value decomposition, we express as,

(19) |

where are unitary matrices and is a diagonal matrix, . The eigenvalues of the hamiltonian are then , the eigen-vectors being,

(20) |

We can write , where are matrices with unit determinant. The PB vector potential and curvature summed over all the negative energy bands can be computed to be,

(21) | |||||

(22) |

Thus can be non-zero only at points where has a vortex type singularity. From Eq.(19), we see that is the phase of . Since the matrix elements of are smooth functions of , can be multi-valued only at points where . These are precisely the DPs. Thus we have proved Eq.(16) showing that the PB curvature for systems with PHS at half filling is that of a set of vortices at the DPs. We also see that contains complete information of the topology of the system. The zeros of the determinant are the positions, , of the vortices. PB flux passing through is , where is the winding number of the phase of around it. We discuss the topological properties discussed above in the context of the KHUB in the following sections.

### iii.1 Topology of KHUB with PBC

We now apply the above developed formalism to the non-interacting hamiltonian KHUB which has PHS. It exhibits nontrivial properties hassan2013b such as topological Lifshitz transitions and non-zero Chern numbers which we discuss in detail here, with fixed at unity.

The first and the second band of our four band hamiltonian overlap in the range beyond which there is a non-zero gap between the bands at all . Our model also features multiple DPs whose number changes as a function of the spin-dependent hopping parameter . The transition points are seen at , and (FIG.(2)).

The location of the DPs can be determined from the energy spectra as the values at which the eigenvalues vanishes. At the DPs the wave-functions of the two sublattices decouple and we get

(23) |

We look for solutions in the direction, which imposes the condition on to be

(24) |

This condition is satisfied by for all . At , the graphene limit, doubly degenerate DPs are located at and summing up to a total of DPs. With increasing , this condition is satisfied by another value of . Thus for there are a total of DPs located at given by

(25) |

where the last two are related to through the underlying honeycomb lattice symmetry. At , six of these DPs merge in pairs at , and , leaving only those at . For , there are only 2DPs. At , six DPs emerge from and move away from each other in the Brillouin Zone with increasing . FIG.(3) shows the DPs for various . The merging and emerging of the DPs, previously discussed in other systems in degail2012; lim2012; fuchs2012; montambaux2009; wunsch2008, is a topological Lifshitz transition degail2012; lim2012.

In order to examine the Lifshitz transitions, we employ either the density of states or the thermodynamic consequences of the Fermi velocity, depending upon the transition point in question. The density of states does not change behaviour for the transition at . The Fermi velocity, which varies linearly with for and is thus expected to vanish at , remains non-zero and finite at that value. This should reflect in many of the thermodynamic properties of the system, and thus it can be used as a probe of this Lifshitz transitions.

The energy dispersion relation of the system for close to each of the DP is linear and is given by

(26) |

where and are small deviations away from the DP and and are constants dependent on . The density of states thus varies linearly with the energy for all except for the values at which the DPs merge and emerge. At , its behaviour changes sharply, with the dominant contribution varying as the square root of the energy. This is because, the dispersion relation takes the form

(27) |

at , and merging points for . On the other hand around the emerging point for , the dispersion relation takes the form

(28) |

thereby giving a constant and a linear contribution to the density of states. Very close to , however the constant term dominates. This sharp change in the density of states at and makes it possible to probe the Lifshitz transitions.

Information about the DPs as well as the PB curvature of the system, as shown earlier, can be obtained from the phase of . In FIG.(4), we plot this for various values of . For , there are eight distinct points around which the phase changes discontinuously by a value of corresponding to the DPs. On the other hand, there are only two such points at .

At the DPs, as shown earlier the PB curvature is singular. Applying a small staggered mass term to induce a gap at the DPs we compute the PB curvature of the second band at two values, shown in FIG.(5). The PB curvature of the second band shows a peak at the DPs.

Using the PB curvature obtain the Chern numbers, for and for . At half-filling, the total Chern number given by vanishes, implying that the Hall conductance also vanishes hassan2013; hassan2013b. Remarkably, even though the Chern number for the lowest and the highest bands are both equal to , the for these bands is not negative for all values of . This surprising result is true for the other bands as well with the signs flipped appropriately. FIG.(6) shows the PB phase as a function of filling for the lowest band, clearly depicting this behaviour. A consequence of this behaviour, can be seen in optical lattice experiments which shall be further discussed in section (IV.2.2). Thus the non-interacting KHUB with PBC shows intriguing topological character.

### iii.2 Topology of KHUB with OBC

The properties of the KHUB that were discussed till now are for PBC. We study the edge states in this model in a cylindrical geometry with zig-zag edges. There are zero energy edge states between the second and the third band and chiral edge states between the bottom two and the upper two bands. The number and the location of the zero-energy edge states in the quasi-momentum direction change as a function of which can be determined using the Zak phase zak1989; saket2013; delplace2011 around a closed contour. There have been proposals to probe these phases in optical lattices atala2012.

At the graphene limit, , there are continuous zero-energy edge states for each of the two spin species for , making a total of states. Here the Zak phase is for , and elsewhere. For values of between these two limits the edge states are not continuous, Fig:(7). The doubly-degenerate edge states in shift to , respectively, forming unique states and thus preserving the total number. On the other hand, for , there are unique continuous edge states for . Beyond again patches of zero-energy edge states occur.

The edge states carry current due to the breaking of time-reversal symmetry in the model. Using the Heisenberg equation of motion for the density operator and the density-current continuity equation, we can compute a general expression for the charge current between two sites on each of the , and links. From the current along the link

(29) |

the average charge current for cylindrical geometry can be computed

(30) | ||||

(31) |

Similarly the average current on the link can also be calculated.

Using these expressions, we find that the total average charge current explicitly involves the time-reversal breaking spin-dependent strength . Thus the existence of non-zero currents at the two edges of the system in FIG.(8) can be attributed to a non-zero . At there is no current at the edges, as expected. As a non-zero edge current appears and is initially negative at the left edge and positive at the right edge. However, by the signs of the currents on the two edges have flipped. The sign of the edge current also depends on which states are filled for . Since the chiral edge states have an opposite and larger contribution than those in the bulk, the edge current flips sign at quarter filling when the former begins to fill, as shown in FIG.(9).

These rich topological properties of the model motivates the study of realizing the model and the properties in optical lattice experiments.

## Iv Experimental Realization

We study the realization of the model and the above mentioned topological features in cold atom experiments in the following sections. The model can be realized in optical lattice systems as shown by Duan et al.duan2003. We systematically derive the hamiltonian in the next section. We then discuss methods of probing the DPs using Bloch-Zener oscillations and we propose a method to detect the PB curvature in optical lattice systems.

### iv.1 Spin-dependent hopping in the honeycomb lattice

The KHUB is a model on the honeycomb lattice with spin-dependent hopping. There are multiple methods of obtaining spin-independent hopping on the honeycomb lattice loon2010. For example, the required honeycomb lattice can be generated by three intersecting laser beams at an angle of between them duan2003; loon2010. In this section, we systematically derive the spin-dependent hopping on the honeycomb lattice using the method suggested by Duan et. al. duan2003 which is different from that discussed earlier mazza2012.

Most fermionic optical lattice experiments are performed using atoms leblanc2006. In the absence of external magnetic field the and the levels of potassium each split into two hyperfine levels. Two of the hyperfine energy levels of are much lower in energy compared to levels of . FIG. (10) is a schematic of the three level system formed by the low levels of and a level of .

The lower two energy levels are separated from the third by a gap which is orders of magnitude largerleblanc2006 than the hopping parameter seen in typical optical lattice experiments mazza2012. Two blue de-tuned laser beams and excite virtual transitions between the first and the third () and the second and the third () levels respectively. Since these virtual transitions are fast compared to the hopping of the atoms, a local microscopic hamiltonian can be written as

(32) |

Here represents the energies of and . is the atomic creation operator for the -th energy level at . The last two terms in the above expression arise due to the interaction of the atom with the laser beams. The wavelengths of the laser beams is such that it only causes transitions from the energy levels and to the energy level with and representing the strength of these transitions respectively and represents the electromagnetic field. The effective two-level system can be obtained by integrating out the third energy level. The action in the path integral formalism for the three level system can be written as

(33) |

Integrating the third energy level we obtain the effective action of the two-level system as

(34) |

where the matrix is given as

(35) |

Since the beams are monochromatic, the electromagnetic fields can be written as , where is the frequency of transition from the -th, , energy level to the third energy level. This is clearly seen in FIG.(10). Thus the effective new matrix is given by

(36) |

The two low-lying energy levels labelled and can be represented by pseudo-spin indices and respectively. The effective potential seen by the pseudo-spins is given as

(37) |

When laser beams of varying intensities and directed along the wave-vectors are incident on this effective two-level atom, the field can be written as

(38) |

with the condition . This implies that the fields arising due to different laser beams are independent. The effective potential which depends on the spin becomes

(39) |

Thus the spin dependent potential can be varied by tuning the laser beams and .

The spin-dependent potential on the honeycomb lattice can be generated by tuning three lasers with different strengths oriented along the three directions , and , at an angle of relative to one another duan2003; loon2010, FIG.(1). The laser beam is sufficient to generate a spin dependent coupling on the link, . The lasers and directed along and with a relative phase difference are required to generate the couplings along these directions. Thus we have and respectively for the and directions. We now write the potential as a sum of spin-independent and spin-dependent parts,

(40) |

where is the effective space dependent magnetic field generated by the laser beams. Its individual components can be written as

(41) | ||||

(42) | ||||

(43) |

where is the wave-vector along and respectively. The spin independent potential can be written as

(44) |

This shows that the spin-independent part cannot be tuned individually as it is coupled to . So we add an additional spin-independent potential which can be tuned without affecting the spin-dependent part. Now the wave-function of the atom in the spin-dependent honeycomb lattice potential follows the time-independent Schrödinger equation of the form

(45) |

where is the total potential. The wave-function of the atoms can be expanded in terms of atomic orbitals which are localized in the sublattice of the -th triangular Bravais lattice, that is

(46) |

This allows us to compute the hopping parameter as

(47) |

Eq. (2), with the hamiltonian Eq.(3), is applicable when the nearest neighbours provide the dominant contributions to the hopping. Thus we need for the and links and for the link with . Once we obtain the spin-dependent hopping the onsite interactions between the atoms in the optical lattice systems can be created by Feshbach resonanceloon2010. At very large interactions both the Kitaev-Heisenberg model chaloupka2010; duan2003; hassan2013b; hassan2013 and the isotropic Kitaev model kitaev2006 can be obtained. Thus the anisotropic Kitaev model can be obtained by tuning the strengths of individual laser beams, and this model can be further mapped onto the toric-code hamiltonian kitaev2003.

### iv.2 Probing the topological properties of the model

#### iv.2.1 Locating the DPs: Bloch oscillations

Recently an experimental method to probe the DPs using Bloch-Zener oscillations breid2006; kolovsky2003 was suggested by Tarruell et al. tarruell2012, whereas a detailed method of numerically simulating such oscillations was discussed by Uehlinger et al. uehlinger2013. Here we probe the DPs in our model using Bloch-Zener oscillations.

At time , the tight binding Kitaev-Hubbard hamiltonian with a staggered potential and in the presence of a harmonic trap is given by

(48) |

Here is the strength of the staggered onsite potential, and are the strengths of the harmonic trap in the and directions, while and represent the spatial coordinates of the lattice site which are measured in terms of the lattice parameter .

We calculate the -particle many body ground state for this hamiltonian and evolve it using the total hamiltonian . Here the interaction term is that of an external force field of magnitude along on the lattice, and is given by

(49) |

where is the position vector of the lattice site. The Schrödinger evolution is

(50) |

where is measured in terms of the Bloch oscillation time period . We choose lattice sites in order to prevent the cloud from ever hitting the boundary. At every time step we measure the projection of the Fourier-transformed many-body density matrix on to the density matrix of the single-particle bands in the presence of the staggered potential,

(51) |

Here is the single-particle eigen-state of the -th band of the non-interacting hamiltonian with staggered mass. It is possible to project the density matrix because we have assumed that the trap potential varies slowly so that the single-particle bands do not change in the presence of the trap.

We also compute probability amplitude per particle

(52) |

where is the number of particles in the system. The quasi-momentum distribution of the particles clearly shows a sudden reduction in the density when a DP is encountered which gives us a method of probing them in experiments.

We study the quasi-momentum distribution for and . In FIG.(11) the quasi-momentum probability amplitude (in the orthogonal coordinates ) of particles in the second band for various instances during one Bloch oscillation is plotted. The parameters are Hz, , Hz, Hz, Hz and . The probability amplitude initially localized around the origin, moves along the direction and encounters a DP at location at time .

On further evolution the DPs at is probed finally returning to the center of the Brillouin zone after one oscillation. There is a transfer of particles to the higher bands close to the DPs as seen in FIG.(13) where we have plotted and . The two peaks in the figure corresponds to the DPs seen from the quasi-momentum distribution. The inset shows the probabilities for individual bands, with transitions occurring between each of the successive bands.

In contrast with the above situation, we show in FIG.(12) the Bloch-Zener oscillation of particles at keeping the rest of the parameters unchanged. At this and in this direction the system has four DPs, all of which are probed at different times by the cloud. There is a transfer of amplitude when the state passes through each of the four DPs.

In FIG.(13) we plot and for . The lower band only shows two peaks even though this system has four DPs. This is because we have been unable to resolve the transfer at each of the four DPs to high accuracy which is obtained only when there is a reasonable fraction of the cloud in the second band at time . This requires a large number of particles in the cloud, increasing its width in momentum space. This decreases the resolution of the DPs, and can only be circumvented by increasing the size of the system, which is limited by currently available computational power to us.

No DPs are encountered when we apply the force field in the direction. The quasi-momentum distribution for the second band corresponds to the second Brillouin zone in the momentum distribution obtained in optical lattice experiments. Thus from the above discussion, we see that the DPs can be probed in such experiments. direction or the

#### iv.2.2 Signal for PB curvature: Rotating condensates

In this section we propose a method of measuring the PB curvature in optical lattice experiments different from methods proposed previously price2012; atala2012. The hamiltonian in optical lattice experiment is generated as given in Eq. (48). Then using time of flight experiments which are the standard methods of imaging in optical lattice experiments, we study the density profiles of the atomic cloud. Here all the external potentials are switched off suddenly following which the atomic cloud is allowed to expand ballistically, the atoms now behaving as free particles. Images of the cloud at various intervals of time shows how the density distribution changes as the cloud expands. If the system initially had non-zero PB curvature, then the cloud along with expansion, rotates thereby probing the PB curvature.

To understand the density profiles we first review the Sundaram-Niusundaram1999 equations (SNE) and then the Thomas-Fermi approximation for a many-fermion system to address the problem of rotating condensates.

The Sundaram-Niu equations sundaram1999 govern the classical dynamics of a wave packet restricted to the band with momentum space width small compared to that of the Brillioun zone and a real space width small compared to the applied external field. These wave packets therefore have a width in real space that is large compared to the lattice spacing but small compared to the scale of the variation of the external fields. The equation of motion of the Bloch electron for two dimensional systems in the absence of magnetic field are

(53) | |||||

(54) |

where is the energy in the absence of external fields, is the PB field and the external potential. Thus induces a force-dependent anomalous velocity.

The Sundaram-Niu equations describe the wave-packet dynamics for a single particle. In optical lattice experiments, the many body dynamics of the atoms needs to be considered for which we use the Thomas-Fermi approximation. This approximation assumes that the ground state is described by a phase-space particle density that incorporates the Pauli exclusion principle. The number of fermions in a phase-space volume around the point , can be written as

(55) |

where is the Fermi energy level and is the single particle hamiltonian given by

(56) |

The initial particle density of the cloud can be calculated from Eq.(55). We now allow the cloud to expand freely and by using time of flight experiments the density at various times can be computed using Liouville’s theorem.

The Thomas-Fermi approximation can be extended to multiple bands if the bands are well separated and if the applied external potential varies slowly enough to prevent interband transitions. Thus, the total phase space density is the sum of the phase space densities for first band and second bands. Additionally the Thomas-Fermi approximation fails when there are DPs in the system.

The above formalism can be applied to our multi-band model since the bands are well separated. A slowly varying external potential is applied and the inversion symmetry is broken to open up a gap at the DPs. We use the same parameters as in uehlinger2013. The system is confined by a rotationally invariant external harmonic potential , where

(57) |

and , and are the dimensionless spatial coordinates of the lattice sites. Here nm is the wavelength of the laser beam, rad/s is the trapping frequency of the potential, is the mass of the atoms and Hz is the nearest neighbour hopping parameter. A staggered mass term Hz is added.

FIG.(14) shows the variation of the number density of the cloud with distance for both and . At it reaches a maximum since the trap potential is zero there. Away from , the total energy of each of the bands increases by , reducing the number of occupied states lying below and thus decreasing the density. At there is a plateau where the density has a constant value of . The plateau occurs due to the energy gap between the second and the first band. The size of the plateau decreases depending on the size of the gap between the two bands. Below , the first and the second bands overlap and hence we see a continuous change in the density as seen from the plot for .

In FIG.(15), we see that the PB phase of the occupied bands is zero at the center , since the Chern number of the contributing bands are equal and opposite. Away from the centre the trap potential reduces the number of occupied states below , decreasing the PB phase as a consequence. A minima is reached when the only contributing band is the lowest band after which the PB phase increases and reaches a value of zero when all the bands are empty. At , the PB phase also shows a plateau similar to that of the density whereas at the PB phase changes continuously. At , there are kinks where the PB phase is lesser than and greater than .

The kink appears because the PB phase of the band is not entirely positive when the Chern number is and vice-versa. To illustrate this we plot, in FIG.(6), the PB phase of first band as a function of the filling factor.

From the Thomas Fermi approximation, the velocity of the cloud can be computed and is as shown in FIG.(16). The velocity is zero where the PB phase is zero. The direction of the velocity at the boundary changes due to the change in the sign of the PB phase as shown in FIG.(6).

From the velocity, the total angular momentum density per particle can be computed and is found to be for and for . This value is large compared to the value obtained in bosonic optical lattice experiments chevy2000, and thus should be observable in such experiments.

The evolution of the cloud density after the traps have been switched off can be computed using Liouville’s theorem. At various times the cloud rotates as it expands as seen in FIG.(17). This rotation strongly indicates a non-zero PB curvature. The cloud also inherits the hexagonal structure of the underlying honeycomb lattice as it expands.

## V Conclusions and Discussion

To summarize, we have analyzed a fermionic model on the honeycomb lattice with spin-dependent hopping that breaks time-reversal symmetry but preserves the particle-hole (chiral/sub-lattice) symmetry, the Kitaev-Hubbard model. The model has DPs and is a semi-metal at half filling. We show that the particle-hole symmetry implies that the total PB curvature vanishes everywhere on the Brillioun zone at half filling. Consequently, the generalized Zak phases are topological invariants that are wholly determined by the positions and chiralities of the DPs. We show that all this information about the topology is contained in the determinant of a matrix defined in Eq.(5). We also numerically show that the structure of the non-dispersive edge states are determined by these topological invariants.

Multiple DPs exist in this model and as the strength of the spin-dependent hopping parameter, , is varied, topological Lifshitz transitions occur where the DPs are created and merge. At , the model is same as graphene and there are 4 DPs. As soon as changes from zero, there is a transition from to DPs. As increases the DPs migrate over the Brillioun zone and pairwise merge at , and , resulting in a transition DPs at . At , there is again a transition to DPs which now emerge from . The signal of these transitions can be seen in the density of states and in the edge state structure. The effect of broken time reversal symmetry of this model can be seen from the existence of the charge currents at the edges of the system. We observe that the edge currents change sign near the transitions.

A scheme to realise spin-dependent hopping in cold atom systems was proposed by Duan et al duan2003. We have analyzed their scheme and have provided a systematic derivation of the low energy effective spin-dependent periodic potential which leads to the spin-dependent hopping matrix elements in the tight-binding model.

Finally, we have examined experimental signals of the topological and geometric features of the model realized in cold atom systems. We have shown that Bloch-Zener oscillations in our system probes the location of the DPs and can hence be used to observe the creation, migration and the merging process. The confining trap breaks the particle-hole symmetry and makes the PB curvatures of the bands observable. We have shown that the effect of this is seen in the rotation of the expanding cloud when the trap is removed. For realistic atomic and trap parameters, we have shown that this provides a clear signal in time of flight experiments.

In conclusion, while the Kitaev-Hubbard model was initially proposed duan2003 to realize Kitaev’s honeycomb model at large and half filling, our work hassan2013; hassan2013b shows that it contains rich physics at all parameter ranges. At half filling, intermediate and large , there are phases with magnetic order and a transition to an algebraic spin liquid hassan2013. At small to intermediate it shows a very interesting phenomena of creation, migration and merging of DPs jeanpaul2013.

We have also shown that at quarter and three quarter filling, there is a Chern insulating phase which indicates that there could be fractional anomalous quantum Hall states when the bands are partially filled. We feel that all these theoretical results strongly motivate attempts to physically realize the model.

Acknowledgement: We thank Mukul S. Laad and David Sénéchal for useful discussions.