Measuring topology in a laser-coupled honeycomb lattice: From Chern insulators to topological semi-metals

Ultracold fermions trapped in a honeycomb optical lattice constitute a versatile setup to experimentally realize the Haldane model [Phys. Rev. Lett. 61, 2015 (1988)]. In this system, a non-uniform synthetic magnetic flux can be engineered through laser-induced methods, explicitly breaking time-reversal symmetry. This potentially opens a bulk gap in the energy spectrum, which is associated with a non-trivial topological order, i.e., a non-zero Chern number. In this work, we consider the possibility of producing and identifying such a robust Chern insulator in the laser-coupled honeycomb lattice. We explore a large parameter space spanned by experimentally controllable parameters and obtain a variety of phase diagrams, clearly identifying the accessible topologically non-trivial regimes. We discuss the signatures of Chern insulators in cold-atom systems, considering available detection methods. We also highlight the existence of topological semi-metals in this system, which are gapless phases characterized by non-zero winding numbers, not present in Haldane's original model.

Abstract. Ultracold fermions trapped in a honeycomb optical lattice constitute a versatile setup to experimentally realize the Haldane model [Phys. Rev. Lett. 61, 2015Lett. 61, (1988]. In this system, a non-uniform synthetic magnetic flux can be engineered through laser-induced methods, explicitly breaking time-reversal symmetry. This potentially opens a bulk gap in the energy spectrum, which is associated with a non-trivial topological order, i.e., a non-zero Chern number. In this work, we consider the possibility of producing and identifying such a robust Chern insulator in the laser-coupled honeycomb lattice. We explore a large parameter space spanned by experimentally controllable parameters and obtain a variety of phase diagrams, clearly identifying the accessible topologically non-trivial regimes. We discuss the signatures of Chern insulators in cold-atom systems, considering available detection methods. We also highlight the existence of topological semi-metals in this system, which are gapless phases characterized by non-zero winding numbers, not present in Haldane's original model.

Introduction
Topological phases of matter have been a topic of great interest in condensed matter physics since the discovery of the integer quantum Hall effect [1]. They are characterized by transport properties -such as a quantized Hall conductivity -that depend on the topological structure of the eigenstates [2], and not on the details of the microscopic Hamiltonian. As a result, such properties are remarkably robust against external perturbations. Integer quantum Hall phases, the first topological insulating phases to be discovered [1], are realized by applying a large uniform magnetic field to a quasi-ideal two-dimensional electron gas, as formed in layered semiconductors structures.
The presence of a uniform magnetic field is not, however, a necessary condition to produce quantum Hall states, as first realized by Haldane [3]. He proposed a remarkably simple model on a honeycomb lattice, with real nearest-neighbor (NN) hopping and complex next-nearest neighbor (NNN) hopping mimicking the Peierls phases experienced by charged particles in a magnetic field. Although the magnetic flux through an elementary cell of the honeycomb lattice is zero, a staggered magnetic field present within this cell locally breaks time-reversal symmetry. Haldane showed that this model supports phases that are equivalent to integer quantum Hall phases: they correspond to insulators with quantized Hall conductivities, σ H = ν e 2 /h where e is the electron charge. In this manner, it is possible to generate a quantum Hall effect without a uniform external magnetic field. The integer ν = ±1 (depending on the particular values of the microscopic parameters) is a topological invariantthe Chern number -characteristic of the phase and robust with respect to small perturbations [4,2]. More recently, a more broad concept of topological insulators has emerged, classifying all possible topological phases for non-interacting fermions in terms of their symmetries [5,6]. In this modern terminology, the Haldane model belongs to the class A of Chern insulators, which are topologically equivalent to the standard quantum Hall states.
The Haldane model has not been directly realized in solid-state systems, due to the somewhat artificial structure of the staggered magnetic field. Interestingly, ultracold atomic gases [7,8] appear better suited to achieve this goal [9,10]. In recent years, many proposals have been put forward to realize artificial magnetic fields for ultracold atoms (see [11] for a review). Staggered fields are relatively easier to implement than uniform ones [12,13,14], and have already been realized in a square optical lattice [15]. Building on these ideas, Alba and coworkers [16] proposed a model very similar to Haldane's that could be realized with ultracold atoms. Their variant is based upon a state-dependent honeycomb optical lattice [8], where cold atoms in two different internal "pseudospin" states are localized at two inequivalent sites of the elementary cell. Additionally, laser induced transitions [12,17] between the nearest-neighbor sites lead to pseudospin-dependent hopping matrix elements containing phase factors, schematically depicted in Fig. 1. Furthermore, Alba and coworkers [16] suggested a measurement based on spin-resolved time-of-flight (ToF) experiments to identify topological invariants.
The present work provides a systematic analysis of the model proposed in Ref. [16] and identifies parameter regimes where Chern insulators emerge. The goals are: firstly, to serve as a detailed guide to possible experiments aiming at realizing such topological phases; and secondly, to discuss the subtle issue of identifying them through ToF methods. Ref. [16] focused on the very special case when one of the NNN hopping amplitudes was zero, and we find that such a system is semi-metallic, not a quantized Chern insulator. In this limit, where the bulk energy gap is closed, the Hall conductivity is no longer simply given by a Chern number σ H = ν e 2 /h, and therefore, this transport coefficient generally looses its topological stability. Yet, the ToF method seems to give a non-trivial signature in this regime, which is robust with respect to small variations of model parameters. The subtlety is that the ToF method of Ref. [16] actually measures a winding number [18], which only coincides with a topologically protected Hall conductivity when the energy gap is open. If this condition is met, the ToF method of Ref. [16] then produces a reasonable experimental measure of the topologically invariant Chern number, and we indeed verify its robustness when varying the system parameters. Interestingly, if the bulk gap is closed, we find that the winding number measured from a ToF absorption image might still depict a stable plateau when varying the microscopic parameters, under the condition that the Fermi energy is exactly tuned at the gap closing point. In this work, such gapless phases associated with a non-trivial winding number will be referred to as topological semimetals. Absent in the original Haldane model, they constitute intriguing topological phases, which can be created and detected in the laser-coupled honeycomb lattice.
In this work, several types of band structures and topological orders will therefore be present: (1) Chern insulating phases, i.e. gapped phases with non-trivial Chern numbers ν = ±1, (2) Topological semi-metals, i.e. gapless phases associated with a non-trivial winding number, (3) Standard semi-metals, i.e. gapless phases with the two bands touching at the Dirac points, as in graphene [19,20], and which are found at the transition between two topological phases. This paper is structured as follows. In Sect. 2, we introduce the model and discuss how the energy band topology can be characterized in terms of Chern numbers. We also discuss the magnetic flux configuration as a function of the model parameters, highlighting the time-reversal-symmetry breaking regimes. Sect. 3 presents the main results, where the phase diagrams are investigated as a function of the microscopic parameters. In Sect. 4, we examine the signatures of the ToF method [16], and compare its results when applied to a Chern insulator or to a semi-metallic phase, i.e., when the topological bulk gap is absent. We summarize the results in Sect. 5, and discuss an extension which implements the Kane-Mele model leading to Z 2 topological insulators [21].

The Hamiltonian
In the model introduced in Ref. [16], cold fermionic atoms are trapped in a honeycomb structure formed by two intertwined triangular optical lattices, whose sites are labeled by A and B respectively [cf. Fig. 1 (a)-(c)]. In the tight-binding regime -applicable for sufficiently deep optical potentials V A,B (x) -atoms are only allowed to hop between neighboring sites of the two triangular sublattices, which correspond to next-nearest neighbors of the honeycomb lattice (denoted n τ , m τ , with τ = A, B). The secondquantized Hamiltonian takes the form The hopping factors between NN and NNN sites of the honeycomb lattice are indicated by t A , t B and te iφ , with φ ≡ φ(m A , m B ) given by Eq. (3). The lattice spacing is a √ 3, and we set a = 1 in the main text, which defines our unit of length. (b) Three-beam laser configuration giving rise to the desired spin-dependent hexagonal lattice, which we describe for 40 K. In this vision, the lasers are detuned between the D1 and D2 lines of the 4S-4P transition, whereby the state-independent (scalar) light shift is zero. The remaining spin-dependent potential -an effective Zeeman magnetic field -is depicted in (c). The strength of the "same-spin" hopping, i.e. ta and t b , is governed by the choice of internal states: the pair |f = 9/2, m F = 7/2 and |f = 7/2, m F = 7/2 produce ta ≈ t b as they have opposite magnetic moments. In contrast, the choice |f = 9/2, m F = 9/2 and |f = 7/2, m F = 7/2 produces ta = t b . The effective Zeeman shift is plotted with a color scale where blue indicates the potential minima for pseudo-spin up atoms, forming the A sublattice; and red indicates the minima for pseudo-spin down atoms, forming the B sublattice. Not shown are an additional pair of Raman lasers, also in the ex−ey plane, that couple between the different sublattices (red and blue in (c)).
whereâ m A (b m B ) is the field operator for annihilation of an atom at the lattice site r m A (r m B ) associated with the A (B) sublattice, and where t A,B are the tunneling amplitudes. Furthermore, the two sublattices are coupled through laser-assisted tunneling, where hopping is induced between neighboring sites of the honeycomb lattice by a laser coupling the two internal states associated with each sublattice. This corresponds to tunneling processes linking nearest neighbors sites of the honeycomb lattice, denoted as m A , m B , which are described by the Hamiltonian Here, the phases φ(m A , m B ) generated by the laser fields are the analogs of the Peierls phases familiar from condensed matter physics [22,23], with r m A and r m B specifying the nearest neighboring sites of the hexagonal lattice. Following the approach of Jaksch and Zoller [12], these phases can be expressed in terms of the momentum p transferred by the laser-assisted tunneling as so that the phases have opposite signs for r A → r B and r B → r A hoppings (cf. Fig.  2(a) and Refs. [11,13]). Finally, the model also features an on-site staggered potential, described byĤ which explicitly breaks the inversion symmetry of the honeycomb lattice [3]. The total Hamiltonian, given bŷ is characterized by the hopping amplitudes (t, t A and t B ), the momentum transfer p = (p x , p y ), as well as the mismatch energy ε.
To eliminate the explicit spatial dependence of our Hamiltonian (5), we perform the unitary transformation

giving a transformed Hamiltonian
with new Peierls phases given bỹ The transformed Hamiltonian (7), featuring complex hopping terms along the links connecting NNN sites, is similar to the Haldane model [3], but with important differences highlighted in Sect. 2.4. Since r n A,B − r m A,B = δ µ − δ ν = ±a λ are the primitive lattice vectors of the honeycomb lattice (see Fig. 1), where µ, ν, λ = 1, 2, 3, the phasesφ in Eq. (8) no longer depend on the spatial coordinates. Therefore, the Hamiltonian (7) is invariant under discrete translations, [Ĥ tot , T 1,2 ] = 0 where T 1,2 ψ(r) = ψ(r + a 1,2 ), allowing us to invoke Bloch's theorem and reduce the analysis to a unit cell formed by two inequivalent sites A and B. In momentum space, the Hamiltonian takes the form of a 2 × 2 matrix, where and k = (k x , k y ) belongs to the first Brillouin zone (FBZ) of the system. We rewrite this Hamiltonian in the standard form whereσ is the vector of Pauli matrices, and d(k) has real-valued Cartesian components defined by The eigen-energies of the Hamiltonian (11) are E ± (k) = (k) ± d(k), where we introduced the "coupling strength" d(k) = |d(k)|.
Our Hamiltonian (11)-(13) differs from the expression derived in Ref. [16], where different Peierls phases were used ‡. Both models are exactly equivalent when t B = 0, where the system describes a semi-metal (cf. Section 3 and 4).
We now briefly describe how the energy spectrum changes with the parameters (t A , t B , ε) of the Hamiltonian (11)- (13), in the absence of momentum transfer p = 0. When t A,B = ε = 0, the band structure is that of graphene [19,20], namely, the spectrum is given by The two bands touch at zero energy for particular points K ± (the so-called Dirac points), where g(K ± ) = 0, and around which the spectrum is quasi-linear with momentum, E ± (k) ≈ ±v F |k|. We will still use the term "Dirac" points to denote K ± , even if the gap is open (in the vicinity of these points the excitations describe massive Dirac fermions). For ε = 0, a bulk gap ∆ ∝ ε opens at the Dirac points, where the gap width is defined as ∆ = min(E + ) − max(E − ) §. This is not a necessary condition to open a gap, as the NNN couplings t A , t B are also able to do so. For t A,B = 0 and ε = 0, the spectrum is now with D ± (k) = −(t A ± t B ) f (k). Next we note that |g(k)| 2 = 3 + 2f (k), showing that a gap ∆ ∝ |t A − t B | opens at the Dirac points due to NNN couplings. For finite momentum transfer p = 0, the energy spectrum leads to more complex spectral structures and phases, to be explored in Section 3. ‡ In Ref. [16], Peierls phases were considered to be of the form φ(m A , m B ) = p · (rm A − rm B ) instead of Eq. (3), cf. the Supplemental Material in Ref. [16]. We note that the correct form (3), used in the present work, corresponds to the synthetic Peierls phases that can be realized with cold atoms in optical lattices, following the method of Ref. [12]. In the following, we study the phases of non-interacting fermions in an opticallattice setup described by Hamiltonian (11)-(13). Such a system forms a metal (or a semi-metal) when the gap is closed ∆ = 0, and an insulator when ∆ > 0. In the latter case, we set the Fermi energy E F in the middle of the bulk gap. This classification in terms of the band structure is not exhaustive, and it must be completed by a description of the topological properties of this band structure. This is examined in the following Section 2.2. In addition, the properties of some peculiar semi-metals are also explored in this work (cf. Section 4).

The Chern number
When the two-band spectrum E(k) exhibits an energy gap ∆, one can define a topologically invariant Chern number [24], which encodes the topological order of the system. As shown in Ref. [4], the Chern number ν is equal to the transverse Hall conductivity, σ H = ν in units of the conductivity quantum, provided the Fermi energy is located in the bulk gap. The Chern number is given by the standard TKNN expression [4,25] where |u (−) (k) denotes the single-particle eigenstate associated with the lowest bulk band E − (k). The Berry's connection -or vector potential -A(k) is defined by This quantity, which defines the parallel transport of the eigenstates over the FBZ [24], also determines the topological order of the system [2]. The integration in Eq. (18) is taken over the FBZ, a two-torus denoted as T 2 , where the contribution due to any singularities of A(k) -to be discussed later on -should be excluded. It is convenient to parametrize the "coupling" vector d(k) in terms of the spherical angles θ ≡ θ(k) and φ ≡ φ(k), defined as where φ = π − arg g(k) for t > 0. In what follows we shall assume that t > 0, without loss of generality. In this representation, the Hamiltonian (11) takes the form The lowest eigenstate of (21) is given by and from Eqs. (19)- (22), we obtain an explicit expression for the Berry's connection, A crucial point to note is that the Berry's connection (23) has singularities at the points in k-space where d x (k) = d y (k) = 0 and d z (k) < 0. The condition d x (k) = d y (k) = 0, which coincides with the zeros of the function g(k), is always fulfilled at the special points K ± = ( 2π 3 , ± 2π 3 √ 3 ), where we have set a = 1. The second condition d z (K ± ) < 0 is only satisfied for certain values of the model parameters (t A,B , ε, p). In terms of the coupling vector d, the singularity takes place at the "South pole" where θ = π and φ is arbitrary, so that the state |u (−) is multivalued there. Note that this singularity can be removed locally by a gauge transformation, but not globally [26]. Moreover, we find that the phase φ = π−arg g(k) yields opposite vorticities at the two inequivalent Dirac points, where γ ± denotes closed loops around the two Dirac points K ± . If these singularities were absent, the integrand in Eq. (18) would constitute an exact differential form over the entire FBZ. In this trivial case, Stokes theorem would then ensure that the integral in Eq. (18) is zero, since this exact two-form is integrated over a closed manifold . To account for these singularities, Stokes theorem can be applied to a contour avoiding them [2,27]. In particular, the Chern number (18) can be written as a sum of integrals performed over the excluded singularities, i.e. by contributions from small circles of infinitesimal radius γ ± around the excluded Dirac points k = K ± at which A(k) is singular, Using Eq. (24) and taking into account the fact that cos θ(k) remains well-defined close to K ± , we find the simple expression for the Chern number which only involves the sign of the "mass" term d z (k) (13) at the two inequivalent Dirac points K ± . A detailed demonstration of Eq. (26), which further highlights the role played by the singularities, is presented in Appendix A. The important result in Eq. (26) shows that the Chern number ν can now be directly evaluated, without performing the integration over the FBZ in Eq. (17). From Eqs. (13)-(26), one can already deduce that non-trivial Chern numbers ν = 0 can only be obtained when d z (k) has opposite signs at the two inequivalent Dirac points K ± , which can only be achieved for p = 0. In the following Section 2.3, we give a physical interpretation in terms of effective magnetic fluxes and time-reversal-symmetry breaking. We also comment on pathological time-reversal symmetric configurations, which necessarily lead to a trivial topological order ν = 0.
To conclude this Section, we note that a Chern insulator is also characterized by current-carrying edge states that propagate along the edge of the system. This edge transport is guaranteed by the opening of a non-trivial bulk gap (∆, ν = 0, cf. Fig. 6 (b)), and it leads to the quantization of the Hall conductivity via the bulkedge correspondence [27]. The latter is observed through transport measurements The FBZ is a two-dimensional torus T 2 , which is a closed manifold. See also Appendix A. in solid-state experiments. In the cold-atom framework, such measurements are not convenient, as they would require atomic reservoirs coupled to the optical lattice. However, alternative methods, based on Bragg spectroscopy [28,29], have been proposed to extract and image these topological edge states [30]. We will use the appearance of chiral edge states later in this paper to strengthen the identification of Chern insulators (Section 3.2). They are obtained from the spectrum of Hamiltonian (5) in a finite geometry [27], as explained in the Appendix B.

Flux configurations and physical description of the model
In this Section, we examine the effects of the Raman-induced phases in Eq. (3) from a less formal point of view, by associating effective "fluxes" to these Peierls phases. First, one can evaluate the number of magnetic flux quanta penetrating each hexagonal plaquette , which yields (cf. Fig. 2 Therefore, in the absence of NNN hopping (i.e. t A,B = 0), the system has a trivial flux configuration Φ = 0 and remains invariant under time reversal. Importantly, when NNN hopping terms are introduced (i.e. t A,B = 0), triangular sub-plaquettes are penetrated by non-zero magnetic fluxes, explicitly breaking timereversal symmetry and potentially leading to QH phases [3]. Considering the subplaquettes formed by the A − B and A − A hoppings, illustrated in Fig. 2 (a), one finds that where we expressed the recoil momentum p = p 1 b 1 + p 2 b 2 in terms of the basic reciprocal lattice vectors b 1,2 , for which b j · a l = 2πδ jl . The sub-plaquettes formed by the A − B and B − B hoppings have a similar flux structure. Thus, the spacedependent Peierls phases (3) produce a flux configuration characterized by three local fluxes Φ 1,2,3 , and which is translationally invariant over the whole lattice (cf. Fig. 2 (a)). We also note that α Φ α = 0, which indicates that the total flux penetrating each hexagonal plaquette remains zero, as found above [3].
The system remains invariant under time reversal when H( where {Φ 1,2,3 } represents the flux configuration stemming from a given p. Besides the obvious case p=0, we find from Eq. (27), that this occurs: • if p 1 and p 2 are both integers, i.e. if p is a vector of the reciprocal lattice.
• if one of the components p 1 , p 2 or p 2 − p 1 is an even integer, and in particular, if p is collinear with one of the basis vectors b 1,2 . For example, when p 1 = 0 (resp. In these pathological "staggered flux" cases, the system remains invariant under time reversal and therefore topologically trivial (note that the number of magnetic flux quanta Φ α is only defined modulo 1). We can verify that these singular time-reversal configurations equally correspond to the condition As established in Eq. (26), the condition (28) naturally leads to a trivial Chern insulator ν = 0 when ∆ > 0, as expected for a time-reversal-invariant system exhibiting a gap. One can check that the condition (28) can be simply rewritten in terms of the vector p = whose solutions exactly reproduce the pathological cases listed above. When p satisfies the condition (29), the system is necessarily a trivial insulator or a standard semi-metal depending on the other parameters. In Section 3.1, we explore other values of the momentum recoil p, where trivial or non-trivial phases can be found depending on the specific values of the parameters t A , t B , ε.

Comparison with the Haldane model
We conclude this Section by comparing the laser-coupled honeycomb lattice (5), with the original Haldane model (cf. [3] and Fig. 2 (b)). In the latter, the hopping factor t 1 between NN sites of the honeycomb lattice is real, while NNN hoppings t 2 are multiplied by a constant phase factor e ±i2πφH (the sign being determined by the orientation of the path). Thus, in the Haldane model, the three small triangular subplaquettes illustrated in Fig. 2 (b) are all penetrated by the same flux Φ 1,2,3 = φ H , whereas the large central triangular plaquette is penetrated by a flux −3φ H . This leads to a staggered magnetic field configuration, with a vanishing total flux penetrating the hexagonal unit cell Φ( ) = 0. We stress that time-reversal symmetry is necessarily broken in the Haldane model, for any finite value of the phase φ H = 0. This important difference between the two models highlights the richness of the lasercoupled honeycomb lattice (11)- (13), where the flux configuration and the nature of the spectral gaps strongly depend on the orientation of the vector p entering the Peierls phases.

Phase diagrams for topological insulating phases
In this section, we perform a systematic characterization of the phase diagram. We set the nearest-neighbor tunneling amplitude to t = 1, thus effectively measuring all energies in units of t. The Cartesian components of the recoil momentum p x and p y are conveniently measured in units of K x and K y , which are the coordinates of the Dirac point K + , with K x = 2π/3 and K y = 2π/3 √ 3. Following the discussions in the preceding Section, we can expect three different phases : • a semi-metal (energy gap ∆ = 0), • an insulator (energy gap ∆ = 0) with trivial topology (ν = 0), • a Chern insulator (∆ = 0, ν = 0).
At this point, let us remind ourselves that the Chern number ν defined in Eq. (17) characterizes the topological order of insulating phases [2]. However, the expression in Eq. (26) could also be formally computed for a semi-metal configuration (∆ = 0), but in this case, the index ν cannot be associated with a robust and topologically protected Hall conductivity. This fact, which is crucial from the experimental detection point of view, is further elaborated in the next Section 4. In this Section, where the focus is set on Chern insulators, we are therefore looking for wide regions in parameter space where both ∆ and ν are non-zero. In Section 3.1, we consider how the system evolves as the recoil momentum p is varied without staggered potential (ε = 0). We examine further the role of anisotropy in the tunneling energies (t A = t B ) in Section 3.2, and finally the role of a staggered potential (ε = 0) in Section 3.3.

Recoil momentum
We first investigate the effects of the Raman recoil momentum p. Here, the staggered potential is set to ε = 0. The phase diagrams shown in Fig. 3 illustrate the appearance of topological phases as a function of the Cartesian components p x and p y , for several values of the tunneling rates t A,B . The areas corresponding to nontrivial topological phases, characterized by the Chern numbers ν = ±1, are indicated by blue and red colors, respectively. Green areas correspond to the trivial insulating phase ν = 0, and white areas signify the "undesired" metallic regime (∆ ≈ 0). The size of the bulk gap ∆ is simultaneously shown through the color intensity. Panel (a) shows the isotropic case with equal next-nearest neighbor tunneling amplitudes set to t B = t A = 0.3t. Here, non-trivial topological phases (ν = ±1) are generally separated by semi-metallic or metallic phases, and these topological regions depict triangular patterns. Panels (b) and (c) correspond to anisotropic cases where the hopping amplitude t B is reduced to t B = 0.2t in panel (b) and set to zero in panel (c). As the anisotropy increases,  Figure 3. Phase diagrams as a function of the recoil momentum components px and py. In all the figures, we set t = 1, ε = 0 and t A = 0.3t. In panel (a) we set t B = t A = 0.3t, and in panel (b) t B = 0.2t. The extreme case t B = 0 is shown in (c). The white regions correspond to metallic phases (i.e. vanishing of the gap ∆ ≈ 0), the blue and red regions correspond to topological phases with ν = ±1. The green regions correspond to trivial insulating phases ν = 0. The "resized" FBZ is indicated by a hexagon, which also serves to highlight the angle dependence with respect to the inverse lattice vectors. The size of the gaps is indicated by the intensity: the lightest shades denote areas where the gaps are ∆ < 0.1t and the brightest areas correspond to 1.5t < ∆ < 2t.
the metallic regions and trivial insulating phases progressively modify the non-trivial islands.
In the special case t B = 0, we find that all the regions that were non-trivial for t B > 0 reduce to semi-metals: when t B = 0, no topological insulating phase is found (contrary to what the Skyrmion behavior of the vector d would suggest [16], see Section 4). We stress that the semi-metallic behavior of the special case t B = 0 is found for the entire parameter space (i.e. for all t A , ε and p), and equally happens for the case t A = 0 and t B = 0. This subtle effect is highlighted in Fig. 6 (c), presented in Section 3.2, where the band structure E = E(k y ) clearly shows the indirect gap closing for the case t B = 0. This energy spectrum suggests that a small perturbation could open the bulk gap and lead to a Chern insulator. However, in Section 3.3, we show that the staggered potential does not open such a non-trivial gap in the case t B = 0. We therefore conclude that the condition t A,B = 0 should be satisfied to generate a robust Chern insulator.
The Hamiltonian is a periodic function of p, and the resulting periodicity of phases is conspicuous in the phase diagrams illustrated in Fig. 3. The central elementary lattice (the "resized" FBZ) cell is marked by a black hexagon ¶ in all panels of Fig. 3. We find that the most convenient non-trivial topological insulating phases (i.e. phases protected by the largest bulk gaps ∆ ∼ 2t) are found for p ∝ (sin N π/3, cos N π/3), where N is an integer. Therefore, setting p x = 0 potentially leads to topological phases with large bulk gaps, which is the most interesting situation for an experimental realization (cf. Section 4).
We now explore how the topological phases evolve as the laser recoil momentum p and the tunneling amplitudes t A,B are modified. As motivated above, we set p x = 0, and then compute the phase diagrams in the p y − t A plane. First of all, we investigate ¶ To be more precise, the arguments of the cosines and the complex exponentials in the Hamiltonian feature p/2, thus the "resized" FBZ is twice larger than the actual FBZ. The panels of Fig. 3 show rectangular regions containing exactly four Brillouin zones. the isotropic case t A = t B (the effects of anisotropy will be discussed in Section 3.2). The phase diagram presented in Fig. 4 indicates that in the realistic situation where t A ≈ t B , the sizes of the topological gaps ∆ are maximum for t A ≈ t B ≈ 0.3t, where ∆ ≈ 2t for p y ≈ 4K y (see also Fig. 3 (a)). Furthermore, this figure indicates that one should generally observe phase transitions between metallic and non-trivial topological phases as p y is varied. Importantly, we note that the system remains metallic (∆ = 0) when the "natural" hoppings t A,B are larger than the Raman-induced hopping t, in particular when t A ≈ t B ≈ t. In the following, we show that an anisotropy t A = t B , or the inclusion of a staggered potential ε = 0, can turn this metallic phase into a topological one.

Anisotropy
In Fig. 5, we show the phase diagram in the plane p y − t B for a large and fixed value of the tunneling rate t A = t. This important result shows that when t A ≈ t, the anisotropy |t A − t B | = 0 is necessary to open non-trivial topological gaps. This effect occurs for a relatively large range of the anisotropy, namely for t B ∈]0, t A ], and for specific values of the momentum p y . For larger anisotropy |t A −t B | > t, the topological phases are destroyed and only metallic and trivial insulators survive. Specific phase transitions between semi-metallic and Chern insulating phases, indicated in Fig. 5 by three successive dots, are further illustrated through the edge-state analysis, in Fig. 6. In panel 6 (b), one indeed observes the presence of topological edge states within the bulk gap, which is the hallmark of a Chern insulator, i.e. through the bulk-edge correspondence [27]. Finally, we note the robustness of the topological edge states within the semi-metal regime ∆ = 0, in Figs. 6 (a),(c), a fact which is further analyzed in Section 4.

Staggered potential
In this Section, we explore the effect of the staggered potential. In Fig. 7, we show the phase diagram as a function of the staggered potential strength ε and of the recoil momentum component p y , for several configurations of the tunneling amplitudes t A,B (we set p x = 0). First, we show the case t A = t B = t in Fig. 7 (a). In this situation, large metallic regions and small non-trivial islands are found in the phase diagram, which can already be anticipated from Fig. 4 for ε = 0. Interestingly, in the totally symmetric case, where t A = t B = t, the topological phases vanish for ε = 0, and they are thus separated along the ε axis + . These results indicate that the staggered + Note that the vanishing of the topological insulating phases for t A = t B = t and ε = 0 can be visualized in Fig. 4. potential is necessary to induce topological phases in this situation where t A = t B = t. However, for large values of the staggered potential, a trivial phase with ν = 0 is always privileged, in agreement with the general belief that such a staggered potential generically leads to trivial phases [3].
For t A = t B < t, non-trivial Chern insulating phases can be formed both with and without the staggered potential. In Fig. 7 (b), we illustrate the effects of the staggered potential for the optimized values of the tunneling rates t A = t B = 0.3t. Here, one observes two topological phases with ν = ±1, which are separated by a small metallic region. This result highlights that one should generally observe phase transitions between semi-metallic and non-trivial topological phases as p y is varied. On the other hand, varying the staggered potential to large values always privileges the transition to a trivial phase with ν = 0.
In the extreme case where t B = 0, shown in Fig. 7 (c), one finds the contour of the phase diagram presented in Ref. [16]. However, we stress that the two central regions featured in this diagram do not correspond to Chern insulating phases, as their corresponding bulk gap is closed. This indirect gap closing is further illustrated in Fig. 6 (c). In the next Section, we analyze this important point in more details.

The winding number and the ToF measurement: Chern insulators, topological semi-metals and Skyrmions
In the previous Section, we have identified the topological insulating phases that could be realized in our cold-atom system, when the parameters (t A,B , ε, p) are tuned in the gapped regimes ∆ > 0. In these situations, the Chern number (17) associated with the low-energy eigenstate |u − can be defined, and its experimental measure would witness a clear manifestation of non-trivial topological order. However, contrary to solid-state experiments where the Chern number is directly evaluated through a Hall conductivity measurement [1], it can only be observed indirectly in the cold-atom framework [16,31,32,33,34,10,35,30]. In this Section, we analyze in detail the topological orders which could be detected through a ToF experiment [16], and further discuss the role played by the bulk gap ∆ in this context. First of all, let us note that the Hamiltonian (11) can be associated with a topological (Pontryagin) winding number [18,36,37,38], which measures the number of times the unit vector n(k) = d(k)/d(k) covers the Bloch sphere S 2 as k evolves on the entire FBZ [18]. When w = 0, this leads to a Skyrmion configuration for the vector field n(k). As will be discussed later in this Section and depicted in Fig. 10, the Skyrmion configuration corresponds to a situation where the unit vector n(k) entirely covers the Bloch sphere once, which for the present model implies that the vector n(k) points in opposite directions (i.e. North and South poles) at the two inequivalent Dirac points, The winding number w characterizes the map n(k) : T 2 → S 2 defined in Eq. (13), and therefore, it is not necessarily related to the spectrum or eigenstates of the Hamiltonian (11) -contrary to the Chern number (17), which is a mathematical index associated with the state |u − [24].
In this work, a topological semi metal denotes a gapless phase ∆ = 0, characterized by a non-trivial winding number w = 0. The fate of the winding number w and its corresponding Skyrmion pattern will be discussed in Subsection 4.2, where these structures are shown to remain stable when ∆ = 0, as long as the gap does not close at the Dirac points. In fact, when the gap is open ∆ > 0, the Chern number (17) is exactly equal to the winding number (30) as can be demonstrated using Eqs. (18), (19) and (22) (cf. also Refs. [18,37,38]). As a corollary, the result in Eq. (31) can be easily deduced from Eq. (26). From the equivalence (32), we observe that the Chern insulating phases discussed in the previous Sections are characterized by a non-trivial winding number w = 0, and therefore, they are also associated with a Skyrmion pattern. In summary, measuring the winding number w in an experiment would allow to equally identify Chern insulators (∆ > 0) and topological semi-metals (∆ = 0).
As first observed in Ref. [16], the vector field n(k) could be detected through a ToF absorption image. From such data, one could then evaluate the winding number w, using a discretized version of Eq. (30). This detection method is based on the fact that n(k) can be expressed in terms of the momentum densities * The Chern number (17) and its corresponding fibre bundle structure [24] could also be formally defined when the gap is indirectly closed, such as in Fig. 6 (c). Thus, the Chern number ν and the winding number w are formally equivalent under the more general gap-opening condition [39]: E − (k) < E + (k) for all k ∈ FBZ. In the present model, this condition reads dz(K ± ) = 0. ρ A,B (k) associated with the two spin species A, B (cf. Fig. 1). Defining the regions for k ∈ K (−) and k / ∈ K (+) , Unfortunately, one cannot generally determine the regions K (±) in an experiment, unless the Fermi energy is exactly located in a bulk gap (in which case K (−) = FBZ and K (+) = ∅). Therefore, the vector field n(k) can only be approximately reconstructed from the data when both bands E ± (k) are partially filled. In fact, if we apply the relation ρ B (k) − ρ A (k) = n z (k) to every pixel of a ToF image , and discretize the expression (30) to evaluate the winding number from this data, we would experimentally measure the following quantity where µ, ν, λ = x, y, z, and where e x,y are the two unit vectors defined on the discretized FBZ. When the Fermi energy is set within a bulk gap, only the first sum contributes K (−) = FBZ , and the quantity w ToF converges towards the winding number w as the resolution of the grid is increased (cf. Appendix C). When the gap is closed, and if the Fermi energy is tuned such that the bulk bands E ± (k) are only partially filled, the quantity w ToF will generally deviate from the quantized value w. Consequently, the assumption of a perfectly filled lowest band (K (−) = FBZ and K (+) = ∅), as considered in the calculations of Ref. [16], is crucial in the case ∆ = 0. However, we indicate that this condition would be difficult to fulfill in an experiment, due to experimental imperfections and finite temperatures. We now illustrate this discussion in Subsections 4.1-4.2, where the signatures of Chern insulators and topological semi-metals are compared and commented. Let us finally remark that the quantity w ToF defined in Eq. (33) is strictly equivalent to the discretized expression for the Hall conductivity σ H , which is not necessarily quantized in the general case where the Fermi energy is not located in a bulk gap (cf. Appendix D).

The Chern insulators
When a spectral gap is opened, ∆ > 0, the winding number (30) is exactly equal to the Chern number (18). This potentially gives rise to a Chern insulator, as illustrated in Figs. 8 (a)-(b), where we compare how the energy gap ∆, the winding number w and the ToF measurement w ToF vary as a function of the recoil momentum p y . As expected from the topological property of the Chern number, we find that the ranges where ∆ > 0 and ν = w = ±1 lead to the clear plateaus depicted by the observable w ToF (p y ) ≈ ±1 (cf. Appendix C for a discussion on finite size effects). We also The other components nx,y(k) of the vector field could be obtained through similar measurements, combined with a rotation of the atomic states [16].  [16]. In all the figures, a vertical dashed line shows the value py = 4Ky used in Fig. 9.
demonstrate the robustness of these plateaus in Fig. 9 (a), where the quantity w ToF is computed as a function of the Fermi energy E F , and where a large plateau ∼ ∆ is observed for fixed values of the other parameters. Therefore, the ToF winding number w ToF shows a robust behavior, and exhibits a clear plateau when the Fermi energy lies in the band gap. In other words, the Chern insulating phase is characterized by a "quantized" winding number w ToF , which is protected against finite changes of the parameters through the existence of a topological bulk gap † †.
Let us stress the important fact that the Chern number ν in Eq. (18) no longer reflects the quantized Hall conductivity when the bulk gap is closed, in which case the system effectively describes a metallic phase. However, the topological order and Skyrmion patterns associated with the winding number w survive even when the bulk gap is closed, as we now explore in the next Subsection. † † The plateaus depicted by w ToF are strictly equivalent to Hall conductivity plateaus (cf. Appendix D). However, since the Hall conductivity is not measured in cold-atom experiments, we choose to represent the observable quantity w ToF in our plots, rather than σ H .

The topological semi-metals
First of all, we find that the ToF winding number w ToF , given by Eq. (33), can be robust even when the gap is closed. This effect is illustrated in Figs. 8 (c)-(d), where we compare how the energy gap ∆, the winding number w and the ToF winding number w ToF vary as a function of the recoil momentum p y . From Figs. 8 (c)-(d), we find that the winding number w displays non-trivial plateaus w = ±1, in regions where the bulk energy gap is closed ∆ = 0. In Fig. 8 (d), we precisely set the Fermi energy at the gap closing point E F = −0.25t, which can be determined from the spectra in Figs. 10(b)-(e). In this specific configuration, the observable winding number w ToF (p y ) depicts plateaus, and it converges towards the quantized value w ToF → w as the resolution of the grid is increased (cf. Appendix C). However, as the Fermi energy is tuned away from this ideal value, we find that the plateaus w ToF (p y ) ∼ ±1 progressively loose their robustness. This dramatic effect is illustrated in Fig. 9(b), where w ToF is computed as a function of the Fermi energy. Here, in contrast with the Chern insulator case shown in Fig. 9(a), the winding number w ToF is strongly parameter-dependent: it only reaches w ToF ≈ w = +1 at the specific Fermi energy E F = −0.25t (cf. Fig. 9 (b)). Consequently, the robust behavior of these topological semi-metals will only be observed if the Fermi energy is precisely tuned at the gap closing point (such as in Fig. 8 (d)). This important fact makes topological semi metals more challenging to detect than Chern insulating phases. We point out that the phase diagrams and Skyrmion configurations presented in Ref. [16], and reproduced in Figs. 8(c)-(d), only feature trivial insulators and "topological semi-metals", since all the computations were performed for the peculiar configuration t B = 0, whose corresponding spectrum remains gapless in the "non-trivial" regions (cf. also Fig. 10).
The topological semi-metal is an intriguing phase, and in this new context, an important question arises: If the topological winding number w remains stable for ∆ = 0 (cf. Fig. 8 (d)), under which condition is this quantity changing its value? We address this question by analyzing the energy spectrum E(k y ) together with the skyrmion pattern depicted by the vector n(k) in Fig. 10. Here, the parameters are the same as in Fig. 8 (d) and p y is varied between K y and 8K y , where transitions between The parameters are the same as in Fig. 8 (c)-(d), namely t A = 0.5t, t B = 0 and ε = −0.5t A . Note that the bulk gap is closed ∆ = 0 in all the figures except in (a). The x and y components of the normalized vector field n(k) = (dx(k), dy(k), dz(k))/d(k) are represented by red arrows for dz(k) > 0 and blue arrows for dz(k) < 0. The corresponding winding numbers w = 0, ±1 are also indicated. The location of the two Dirac points K ± are indicated by two vectical lines (top) and circles (below). A non-trivial winding number w = ±1 is clearly seen when the vector n(k) has covered the whole Bloch sphere once. Namely, when the North (n = +1z) and South (n = −1z) poles have been reached at the two inequivalent Dirac points K ± (cf. Eqs.(26)-(32)). Note that topological phase transitions w → w occur through a gap closing point, which is accompanied with the vanishing of dz(K D ), at the Dirac point K D . w = 0 ↔ w = +1, but also between w = +1 ↔ w = −1, are expected. From Fig. 10, we find that the winding number w, and its corresponding Skyrmion pattern, remain extremely stable as long as the energy bands do not touch at the Dirac points. When a direct gap closing occurs at the Dirac point K D , where K D denotes K + and/or K − , we observe a topological phase transition signaled by a change in the winding number w (cf. Fig. 10 (b) and (d)). In particular, we find that this direct gap closing is accompanied with the cancellation d z (K D ) = 0 at the band touching point K D .
In fact, the gap-closing condition d z (K D ) = 0 should necessarily be satisfied at the transition between different values of the winding number w, as can be deduced from the equivalence (32) and from the simple expression (26). The topological phase transitions are clearly visible on the Skyrmion patterns of Fig. 10, where a non-trivial winding number w = ±1 emerges when the vector n(k) has covered the whole Bloch sphere once. We remind the reader that this full covering of the Bloch sphere is achieved when the vector field n(k) reaches the North (n = +1 z ) and South (n = −1 z ) poles at the two inequivalent Dirac points K ± . In Fig. 10, we observe a radical change in the behavior of n(K ± ) as p y is varied. For example, for p y = K y ( Fig. 10 (a)), the vector field n(k) visits the North pole twice (i.e. at K + and K − ) but never the South pole (w = 0), while for p y = 4K y (Fig. 10 (c)) the vector visits the entire Bloch sphere once (w = 1). Between these two topologically different configurations, a direct gap closing occurs at k = K − for p y = 2K y (cf. Fig. 10 (b)), a singular situation where the gapless phase is equivalent to a standard semi-metal [19,20]. We note that transitions w = 0 ↔ w = ±1 require a single gap closing point (cf. Fig. 10 (b)), while transitions w = +1 ↔ w = −1 involve two gap-closing points (cf. Fig. 10 (d)). Therefore, in agreement with the equivalence (32), we observe that the topological phase transitions between different topological semi metals are of the same nature as the transitions between different Chern insulators, in the sense that both phenomena occur through direct gap closing (driven here by the control parameter p y ).
In summary, we conclude that the laser-coupled honeycomb lattice and the ToF method of Ref. [16] offers the possibility to explore the topological order of topological semi-metals, which survive in the absence of a band gap. However, we remind the reader that this detection scheme relies on the evaluation of the winding number w ToF , through a ToF measurement of the vector field n(k), which only converges towards the quantized value w for a complete filling of the lowest energy band E − (k). Thus, the experimental detection of topological semi-metals would constitute a subtle task, in the sense that the Fermi energy should be finely tuned in order to maximize the filling of the lowest band ( Fig. 8 (d)). Let us stress that the winding number can only take three possible values, w = 0, ±1. Therefore, an experimental plateau w ToF (p y ) ∼ ±1, stemming from a slightly incomplete filling of the band E − (k) and from finite size effects, would already provide an acceptable witness of non-trivial topological order.
We end this Section by observing that the transitions between topologically different semi-metals are driven by the laser recoil momentum p, and therefore, this interesting effect cannot be captured by the original Haldane model.

Conclusion
In this work, we explored the rich properties of the laser-coupled honeycomb lattice, which is described by the Hamiltonian (11)- (13). We demonstrated the existence of robust Chern insulators in this system, which can be reached in experimentally accessible regions of the large parameter space. In particular, we showed that the possibility of producing such non-trivial phases highly depends on the laser-coupling, through the orientation of the momentum transfer p and the effective (laser-induced) tunneling amplitude t. We showed that it is important to finely tune the ratios t A /t B and t A,B /t in order to open large and robust topological bulk gaps of the order ∆ ∼ 2t. We also discussed the role of the staggered potential ε, which is shown to be crucial in the fully symmetric regime t = t A,B , and which could also be used to drive transitions between topological phases of different nature. Importantly, we addressed the question of detectability in the context of the quest for robust Chern insulators, and we stressed the importance of identifying regimes corresponding to large bulk gaps. We showed that an experimental measure of the topological winding number (30)-(33), e.g. through ToF measurement [16], yields a strong signature for two types of topological phases: the Chern insulating phase and the topological semi-metal (a semi-metal characterized by a non-trivial winding number). Importantly, we showed that the detection of the topological semi-metal would require a delicate tuning of the Fermi energy, which privileges the search for Chern insulating phases from an experimental point of view.
The Chern insulator could alternatively be detected through the identification of chiral edge states, which are protected by the topological gap. A clear signature could be obtained, for example, using the shelving method described in Ref. [30]. From the spectra presented in Figs. 6(a),(c) and Fig. 10, we find that these topological edge states remain robust in the topological semi-metallic phase: the edge states can only disappear from the bulk gap through direct band-touching processes at the Dirac points. However, the experimental identification of these robust edge states for the semi-metallic regime, e.g. using the shelving method, remains an open question to be explored.
Let us finally end this work by mentioning the fact that this system could be directly extended to reproduce the spinful Kane-Mele model for Z 2 topological insulators [21] (see also its generalizations [40,41,42,43,44,45]). In this case, each triangular sublattice should trap atoms in two internal atomic states (Zeeman sublevels), yielding a "spin"-1/2 structure. These atoms should then be coupled independently by lasers in such a way that the tunneling operators, which are 2 × 2 matrices acting between NN sites n A and m B , have the form where σ Z acts on the "spins" and U (n A , m B ) = U † (m B , n A ). In this spinful honeycomb lattice configuration, non-trivial Z 2 topological phases featuring helical edge states [21,46,47], should be reached in the non-trivial regions identified in Section 3. Thus, the versatile laser-coupled honeycomb lattice is well suited for the exploration of two-dimensional topological phases with cold atoms [48,10,49,50,30,51,52].

Appendix A. Analytical calculation of the Chern number
In this Appendix, we provide a more detailed calculation of the Chern number (18), which further highlights the role played by the singularities at the Dirac points K ± . First, we express the Chern number as where the Berry's curvature F = F xy dk x ∧ dk y is a two-form associated with the Berry's connection A = A µ dk µ = i u (−) |∇ µ |u (−) dk µ through the exterior derivative F = dA, with Note that contributions due to any (gauge-dependent) singularities of A(k) should be excluded in the first-line of Eq. (A.1). Considering the gauge in which the lowest eigenstate of the Hamiltonian (21) is given by we find At this point, let us note that the Berry's curvature F is a gauge invariant quantity, which remains well defined over the entire FBZ. In contrast, the Berry's connection A depends on the gauge and can potentially possess singularities within the FBZ. In the present gauge, the singularities correspond to cos θ(k) = −1, which can only happen at a Dirac point K D , under the condition that d z (K D ) < 0 (cf. also main text).
We stress that such singularities, if present, could either take place at one or at two inequivalent Dirac points K D = K ± , depending on the model parameters (t A,B , ε, p) that determine the specific values of d z (K ± ). In the following, we will show that it is the number of singularity points that determines the non-triviality of the Chern number in Eq. (A.1). To do so, let us consider the following situations: Absence of singularities When the Berry's connection is regular over the entire FBZ, the Berry's curvature F = dA is an exact differential form. In this case, the Chern number (A.1) is given by the integral over a closed manifold (i.e. the two-torus T 2 ) of an exact differential form, which is trivial from Stokes theorem. In particular, when d z (K ± ) > 0 at the two inequivalent Dirac points, the Berry's connection (A.4) remains regular over the entire FBZ and the Chern number necessarily vanishes.
One singularity at a unique Dirac point Now suppose that the Berry's connection A in Eq. (A.4) is singular at one Dirac point, say K D = K + . This situation occurs when d z (K + ) < 0 and d z (K − ) > 0, which we now consider to be the case in this paragraph.
In the presence of such a singularity, the Berry's curvature F is no longer an exact differential form, and Stokes theorem cannot be applied globally over the torus. In order to compute the integral (A.1), we partition the FBZ into two complementary regions, R I and R II , whose common boundary ∂R is chosen to be a loop γ + encircling K + . Here, we define the region R II as the one that contains the Dirac point K + at which the singularity takes place (cf. Fig. A1 (b)). Then, we define specific gauges within each region [26,2,27], In this patchwork configuration, the Berry's connection A = {A I , A II } is a locallydefined quantity, which is now given by inside the regions R I and R II , respectively. We note that the gauge structures of the two individual regions are connected at the frontier ∂R = γ + through the gauge transformation (A.10) Furthermore, we note that the Berry's connection A II (k) is now regular at the Dirac point K + , where d z (K + ) < 0. Therefore, the locally-defined Berry's connection A = {A I , A II } is regular over the entire FBZ, and the integral (A.1) can now be computed by applying Stokes theorem to the two different regions [2,27] Therefore, when the singularity only takes place at the Dirac point K + , the Chern number is non-trivial and its value is directly related to the vorticity v + associated with this Dirac point [27].
In the opposite situation, where the singularity only takes place at the other Dirac point K D = K − , namely when d z (K + ) > 0 and d z (K − ) < 0, a similar calculation (with ∂R = γ − ) yields simultaneously removes the singularities at both Dirac points. In this case, Stokes theorem can be applied globally over the entire FBZ, leading to a zero Chern number (cf. the case with no singularity).
Synthesis From the results presented above, we conclude that the Chern number ν characterizing the topological order of our system can only take non-trivial values ν = ±1 when the Berry's connection (A.4) features a unique singularity inside the FBZ. Therefore, this non-trivial regime is reached when the function d z (k) has opposite signs at the two inequivalent Dirac points K ± . Then, from Eqs. (A.11)-(A.12), we finally obtain the result ν = 1 2 sign d z (K + ) − sign d z (K − ) , (A.14) already announced in Eq. (26).

Appendix B. Calculation of topological edge states in finite geometries
It is a standard procedure to determine the edge state structure by considering a semi-infinite system, namely using a cylindrical geometry in which periodic boundary conditions have only been applied to one spatial direction. Here, we assume that the system is closed along the y direction only, and we write the single-particle eigenfunctions as ψ A,B (r) = exp(ik y y) u A,B (r), where u A,B (r + a 3 ) = u A,B (r). Setting Ψ n = (u A (r n A ), u B (r n A − δ 2 )), where the index n labels the sites along the open direction (chosen along x here), we obtain the Harper-like equation EΨ n = DΨ n +RΨ n+1 +R † Ψ n−1 , D = −ε − 2t A cos p · a 3 /2 − k y (a 3 ) y −t e −iky(δ2)y + e −iky(δ1)y −t e iky(δ2)y + e iky(δ1)y ε − 2t B cos p · a 3 /2 + k y (a 3 ) y , R =     −t A e −ip·a1/2 e iky(a1)y + e −ip·a2/2 e iky(a2)y −te −iky(δ3)y 0 −t B e ip·a1/2 e iky(a1)y + e ip·a2/2 e iky(a2)y     . (B.1) The energy spectrum E = E(k y ), describing the bulk but also the edge states, can be obtained by solving the corresponding 2L × 2L Hamiltonian matrix numerically, where n = 1, . . . , L. We stress that the chiral edge states, identified with this method, could lead to clear signatures in an optical-lattice setup, even in the presence of an external confining trap [10,30,51]. The Hall conductivity is given by the Kubo formula [4] σ H = e 2 i V Eα<EF ∂ kx u α (k)|∂ ky u α (k) − (k x ↔ k y ), where V is the volume and where the sum α takes into account the contribution of all the occupied states |u α = |u ± (k) with energy E α = E ± (k) < E F . In our two-band system, the sum in Eq.(D.1) can be decomposed into two parts x (k) is the Berry's curvature associated with the state |u (±) (k) , and where we used the fact that F in the limit ∆k x,y → 0. When the Fermi energy is not located in a bulk gap, the Hall conductivity must be computed using the more general expression σ xy = e 2 h 1 2π which takes into account the fact that both bands E ± (k) could be partially filled.