Topological phase transitions in the non-Abelian honeycomb lattice

Ultracold Fermi gases trapped in honeycomb optical lattices provide an intriguing scenario, where relativistic quantum electrodynamics can be tested. Here, we generalize this system to non-Abelian quantum electrodynamics, where massless Dirac fermions interact with effective non-Abelian gauge fields. We show how in this setup a variety of topological phase transitions occur, which arise due to massless fermion pair production events, as well as pair annihilation events of two kinds: spontaneous and strongly-interacting induced. Moreover, such phase transitions can be controlled and characterized in optical lattice experiments.


I. INTRODUCTION
The ability to control systems at the quantum level offers an alternative and exciting route to deepen our understanding of nature [1]. In particular, the experimental design of manybody Hamiltonians provides a controlled method to study phenomena which were only believed to occur in the realm of condensed-matter physics. Ultracold atoms in optical lattices constitute a rich playground where the behavior of condensedmatter systems can be mimicked [2,3], such as the superfluid-Mott insulator transition for bosons [4] and fermions [5,6]. Often, the experimental versatility surpasses the possibilities of real materials, envisaging scenarios that were previously thought of as mere fiction. For instance, neutral bosonic or fermionic matter can be subjected to effective non-Abelian gauge fields [7,8]. In this work, we shall describe the properties of the non-Abelian honeycomb lattice, namely, a 2D free Fermi gas in a honeycomb lattice subjected to additional non-Abelian gauge fields. The bare Fermi gas can be considered to be a quantum-optical analogue of graphene [9], an interesting material with relativistic massless excitations that has been recently synthesized in a laboratory [10]. As shown below, the addition of external non-Abelian gauge fields induces a variety of topological phase transitions beyond the usual Landau symmetry-breaking paradigm [11], which can be observed in cold-atom experiments.
The very unusual properties of graphene, a single layer of carbon atoms packed in a honeycomb lattice, rely on the fact that the low-energy excitations display a linear dispersion relation [12], and are thus described by massless Dirac fermions [13]. Accordingly, it is possible to observe exotic effects in low-temperature table-top experiments, which usually belong to high-energy physics (see [14] and references therein). In this context, phenomena such as Klein tunneling [15,16], or the relativistic extension of Landau levels [17][18][19] have already been observed. Interestingly, such effects might lead to the experimental realization of fullyrelativistic Schrödinger cat states, the so-called Dirac cat states [20]. Besides, the transport properties of graphene are determined by the underlying relativistic excitations, which are responsible for the anomalous half-integer quantum Hall effect (QHE) [21][22][23][24][25][26]. These effects were also discussed in the context of ultracold atoms in honeycomb [9], and T 3 (rhombic) lattices [27].
In this work, we describe the novel effects that take place in fermionic optical lattices when external gauge fields are switched on. Let us note that additional Abelian gauge potentials [8,[28][29][30][31][32][33][34][35][36][37][38] lead to analogues of well-known effects, such as the Hofstadter butterfly [28], the Escher staircase [29], the integer QHE [39,40], the anomalous QHE in the honeycomb [41] and Kondo lattices [42], or even the fractional QHE [43]. When the external field is non-Abelian [7,44,45], a whole plethora of new phenomena arises. In this regime, it is possible to realize a quantum-optical spin Hall effect [46], to observe a modified Metal-Insulator transition [47], or to induce quasi-relativistic effects [48], such as Zitterbewegung [49]. It is also possible to observe non-Abelian effects in atom optics [50], to obtain an optical-lattice version of a spin field-effect transistor [51], to ascertain the absence of localization in disordered relativistic one-dimensional systems [52], or to study the effects of non-Abelian magnetic monopoles [53] and Yang Mills theories [54]. Such proposals [7,44] can also be used to produce effective spin-orbit interactions [55], to realize topological order [56,57], or to generalize the integer QHE to a non-Abelian scenario [58]. Let us finally remark that non-Abelian features also provide an interesting setup where Dirac fermions emerge in a square lattice [40,59], a regime previously restricted to staggered fields [60,61].
It is the purpose of the present article to study the effects of non-Abelian gauge fields on the emerging relativistic fermions that arise in the optical-lattice analogue of graphene. As discussed below, these fields lead to spontaneous and stimulated processes, where pairs of relativistic fermionantifermion are created or annihilated. Intriguingly, these interacting scattering processes occur in a non-interacting Fermi gas, which clearly shows that the background non-Abelian field has dramatic consequences on the emerging low-energy theory. Interestingly, these creation-annihilation processes can be understood in the light of topological quantum phase transitions [11], since the number of conical singularities and thus the topology of the Fermi surface is altered. Additionally, the creation-annihilation events occur between Dirac fermions with opposite topological charges [62], which allow us to further characterize the different topological phases. An emphasis is placed upon the detection of this effect from density measurements.
This paper is organized as follows, in Sec. II we review the properties of a 2D Fermi gas in a honeycomb optical lattice. In Sec. III, we introduce an additional non-Abelian gauge field, and describe its effect on the emerging relativistic fermions. Here, new phases with different Fermi surface topologies arise, and the gauge-induced phase transitions between them are presented in Sec. IV. In Sec. V, we present a detailed study of the anomalous quantum Hall effect that arises in the different topological phases. In Sec. VI, we discuss different techniques to experimentally detect and characterize these phases. Finally, Sec. VII contains our conclusions.

II. ABELIAN HONEYCOMB LATTICE
In this section, we review the peculiar properties of a free Fermi gas loaded in a honeycomb lattice (see Fig. 1(a)). In condensed matter, graphene provides an ideal realization of this system, where conduction electrons are tightly bound to the carbon atoms distributed along a honeycomb crystal lattice [12,14]. The recent isolation of graphene [10] reveals an experimental path towards truly 2D-materials which present distinctive features from their 3D-counterparts. In particular, graphene has the band structure of a semi-metal, since the valence and conduction bands touch at isolated points of the Brillouin zone. Moreover, the dispersion relation around such points is linear, and thus low-energy excitations behave as relativistic massless Dirac fermions. In accordance, this noteworthy material connects two different fields: relativistic quantum electrodynamics and condensed-matter physics.
A quantum-optical analogue of this 2D material can be achieved with fermionic atoms, such as 40 K or 6 Li [5,6], loaded in a honeycomb optical lattice [9] (see Fig. 1(a)). The corresponding optical dipole potential can be generated by three coplanar lasers in standing- [9] or traveling-wave configurations [63]. In fact, the bosonic superfluid-Mott insulator transition in both honeycomb and triangular 2D lattices has been recently observed [64]. In order to achieve the half-filling and non-interacting regime, pairwise interactions should be controlled by means of Feschbach resonances [65]. Let us remark that this setup offers an optical version of freestanding graphene, where dislocations, impurities, phonons, or curved defects are completely absent. Furthermore, the controllability of the hopping strengths and on-site energies, leads to anisotropic and massive regimes that go beyond the real material possibilities [9,46].
A. Band structure: Tight-binding model and conical singularities Let us review the band structure of fermions hopping along the honeycomb lattice [12] (cf. Fig. 1(a)), which has a two-atom basis with unit vectors a 1 = a 2 (3, , and lattice spacing a. Fermions in the Asublattice hop to three nearest neighbors of the B-sublattice , according to the following Hamiltonian Here, the fermionic operators a † i , b † i (a i , b i ) create (annihilate) a fermion on i−th site of the A-and B-sublattice respectively, and t stands for the nearest-neighbor hopping energy. Note that the hopping strength is controlled through the optical potential depth [63]. The translationally invariant Hamiltonian in Eq. (1) can be diagonalized in momentum space by the introduction of sublattice operators a k = ∑ r j ∈A e ikr j a j , and b k = ∑ r j ∈B e ikr j b j . Accordingly, the Hamiltonian becomes a and Ψ k = (a k , b k ) t is a spinor operator. Note that the Hamiltonian in Eq. (2) is expressed in terms of the sublattice flip operators σ + = | ↑ ↓ |, σ − = | ↓ ↑ |, with ↑ (↓) representing the A(B)-sublattice. This Hamiltonian is readily diagonalized, yielding the following band structure which is represented in Fig. 2(a). In this figure, we observe that the two energy bands touch at different points inside the first Brillouin zone, and therefore the Fermi surface at half filling is reduced to a discrete set of isolated points (see also Fig. 1(b)). Besides, a closer inspection shows that the energy dispersion around these singular points is conical E ∼ |k| (see Figure 1: (a) The honeycomb lattice can be described as two triangular Bravais lattices A, and B, spanned by the unit vectors a 1 and a 2 . The nearest neighbors of the B-sublattice are A-sublattice sites given by u 1 , u 2 , u 3 . (b) The reciprocal lattice is spanned by b 1 , b 2 , and the first Brillouin zone corresponds to a rotated hexagon, whose corners correspond to Dirac points (Note that only K, K are independent). Fig. 3(a)), which is a clear signature of the relativistic nature of the low-energy excitations at half filling. Let us now describe the content of figs. 3(a) and 3(b) in more detail. In Fig. 3(a), we schematically represent the distribution of Dirac cones inside the hexagonal first Brillouin zone. We find six conical singularities localized on the corners of the Brillouin zone, where the valence and conduction bands touch (see also Fig. 2(a)). Note that only two of such Dirac points are independent, the remaining being obtained from this pair by mere translations in the reciprocal unit basis b 1 = 2π 3a (1, The two independent Dirac points K, K are represented by opaque cones, whereas the remaining cones are transparent. Accordingly, we can state that two conical singularities arise in the band structure of the honeycomb lattice. This schematic figure is further supported by Fig. 3(b), where the contour of the conduction band in Eq. (3) has been represented. In this contour plot, the color scale is the following: dark blue colors represent Figure 2: Band structure and conical singularities of the Fermi gas in a honeycomb lattice. a) Graphene regime, where the bands touch at the corners of the hexagonal Brillouin zone. Accordingly, the fermion gas can be described as a semiconductor with vanishing gap or a semi-metal. b) Π-flux regime α = β = π/2, where the number of singularities is clearly raised with respect to graphene. c) Non-Abelian regime α = β = 0.2, where the valence and conduction bands split in two, and the structure of the conical singularities is modified as a consequence of the coupling to the gauge field. energies close to zero E F = 0, while light red colors correspond to higher energies, and thus to empty fermionic states. From this figure, we observe that six zero-energy points arise in the corners of the Brillouin zone, an hexagon represented with dashed lines. It is precisely at these zero-energy points that the conduction and valence bands touch, giving rise to the conical singularities formerly described. As shown in the next section, massless Dirac fermions emerge as the fundamental low-energy excitations around these singular points.

B. Low-energy excitations: relativistic fermions
To study the low-energy excitations, we expand the Hamiltonian in Eq. (2) around the conical singularities located at . The long-wavelength excitations are obtained after expanding the momentum around each Dirac point, such as K, namely p =h(k − K), which leads to where α x = σ y , α y = −σ x , are the Dirac matrices, which are reduced to the usual Pauli matrices in a 2+1 Minkowski spacetime. This Hamiltonian is completely equivalent to the Dirac Hamiltonian of massless fermions, if we understand the Fermi velocity c = 3 2 tā h as the effective speed of light. The relativistic Hamiltonian in Eq. (4) readily accounts for the conical singularities E ∼ c|p| observed in Fig. 2(a). Additionally, it identifies the lattice geometry (i.e. the existence of two interpenetrating triangular sublattices) as the responsible of the spinor structure of the effective relativistic theory. Let us note that the low energy quantum field theory associated to the remaining Dirac point K also describes emerging massless fermions, but with a different representation of the Dirac matrices (i.e. satisfying the Clifford algebra {α j , α k } = 2δ jk ), namely, α x = σ y , α y = σ x . Therefore, two species of independent Dirac fermions N d = 2 emerge in the low-energy theory of a two-dimensional Fermi gas in the half-filled lattice.
The effective relativistic picture emerging at low energies in Eq. (4) should be compared with the exact energy bands in Eq. (3). In Fig. 4(a), we represent the exact band structure around the Dirac point K (transparent surfaces), and compare it with the conical energy dispersion predicted from the Dirac Hamiltonian (solid surfaces). Accordingly, we show that the description of low-energy excitations in terms of relativistic massless fermions is valid in a neighborhood of the Dirac points (i.e.δ k = k − K 1/a).

III. NON-ABELIAN HONEYCOMB LATTICE
In the previous section, we derived the effective low-energy quantum field theory for the half-filled fermionic honeycomb lattice, and showed that two species of massless Dirac fermions emerge at large wavelengths, each around a different Dirac point. The coupling of such massless fermions to Figure 3: Distribution of Dirac points with their associated topological charges along the first Brillouin zone for (a) Graphene P 1 , (c) R u 3 -reversed graphene P 2 , (e) R u 3 R u 1 -reversed graphene P 3 , (g) R u 1 -reversed graphene P 4 , (i) Π-flux P 5 , (k) Non-Abelian point α = β 1. Contour plot of the energy bands for (b) Graphene P 1 , (d) R u 3 -reversed graphene P 2 , (f) R u 3 R u 1 -reversed graphene P 3 , (h) R u 1 -reversed graphene P 4 , (i) Π-flux P 5 , (l) Non-Abelian point α = β 1.
external electromagnetic fields leads to interesting effects distinctive of relativistic field theories, such as the Klein paradox [16]. In this section, we shall focus on the fermion coupling to constant magnetic fields. In the Abelian scenario, this leads to the noteworthy anomalous integer quantum Hall effect, either in graphene [21][22][23][24][25][26], or honeycomb optical lattices [41]. As discussed below, this paradigm should be reconsidered when matter is coupled to non-Abelian gauge potentials. Let us remark that even if such fields are not available in condensed-matter experiments, it is possible to synthesize them in a controlled manner by means of laser-assisted tunneling processes in optical lattices [7], or dark-state methods [44]. Therefore, surpassing graphene-based materials, optical lattices present exotic avenues at the forefront of condensed matter, high-energy physics, and atomic physics.
In non-Abelian lattice gauge theory, a unitary matrix U i j is associated to the link connecting the lattice points r i → r j , which must be generated by the elements of the underlying Lie group. Consequently, the Abelian hopping Hamiltonian in Eq. (1) must be generalized to where a † iτ , b † iτ (a iτ , b iτ ) create (annihilate) a fermion at lattice site r i of the A-and B-sublattice, with a color-like component τ = 1, 2, ..., N corresponding to the fundamental representation of the Lie group. Such a situation has been depicted in Fig. 5 for the particular case of SU(2) fields (i.e. fermions have two colors τ = 1, 2 ). Here, we observe that the coupling of matter to non-Abelian fields induces unitary transformations along each hopping path: U 1 = e iατ x , U 2 = I, and U 3 = e iβ τ y , where τ x , τ y are Pauli matrices in the colorcomponents, and α, β represent the non-Abelian fluxes (see Fig. 5). In the momentum representation, the Hamiltonian in Eq. (5) and the spinor Ψ k = (a k1 , a k2 , b k1 , b k2 ) t includes pseudo-spin and color degrees of freedom, which correspond to the underlying sublattices and to the internal fermionic components. In the following sections, we study the pattern of emerging relativistic fermions associated to different gauge fluxes (α, β ). As shown below, the total number of conical singularities and their distribution within the Brillouin zone depends on the particular values of the gauge fluxes. We shall characterize different Abelian and non-Abelian phases in terms of these low-energy massless fermions.

A. Wilson loop and non-Abelian regimes
The non-Abelian gauge fields that induce these conditional hoppings are defined through a generalized Peierls substitution U i j = V † i e iaA µ (r l ) V j , where the independent matrices V j reflect the local gauge invariance, and the gauge field A µ µ µ (r l ) is located at the link r l = 1 2 (r i + r j ) and directed towards µ µ µ = r j − r i . The particular choice U 1 = e iατ x , U 2 = I, and U 3 = e iβ τ y in Fig. 5, is locally equivalent to the following SU(2)-generated gauge field A = ∑ aµ A a µ (r l )τ a e µ In the continuum limit a → 0, and according to non-Abelian gauge theory, the associated strength tensor , describes a constant magnetic field along the z−axis B = 2αβ √ 3 τ z e z . On the lattice, these external fields are further characterized by the Wilson loop W = trU , where is the unitary transformation around an elementary plaquette. The Wilson loop represents a gaugeinvariant non-local observable related to the non-Abelian flux through an hexagonal plaquette, and is expressed as follows W (α, β ) = 2 cos 2 α + 2 cos 2β sin 2 α.
In Fig. 6, we depict the Wilson loop as a function of the external fluxes, and identify the Abelian (|W | = 2) and non-Abelian (|W | = 2) regimes. Figure 5: Scheme for the fermionic honeycomb lattice subjected to SU(2) gauge fields, where each hopping is dressed by U 1 = e iατ x , U 2 = I, and U 3 = e iβ τ y . We have used double-line links to indicate the underlying color degrees of freedom τ = 1, 2 associated to each fermion. The gauge fields on the links induce a color-dependent hopping between adjacent sites, which has been represented by squares.
• Π-flux point P 5 for (α, β ) = ( 1 2 π, 1 2 π). Note also that the lines connecting P 1 , P 2 , P 3 , P 4 correspond to Abelian phases (see Fig. 6). In the next section, we derive the band structure and low-energy massless fermions for these different regimes. We show that the four graphene-like points P 1 , P 2 , P 3 , P 4 have N d = 4 species of emerging Dirac fermions (notice the doubling due to color-degeneracy), and are thus topologically equivalent. Conversely, the Π-flux phase presents N d = 8 different gapless fermions, and thus a distinct Fermi surface topology. As shown in Sec. V, transport properties in this Π-flux phase, such as Hall conductances, clearly deviate from the graphene-like phases.

Non-Abelian regimes:
The remaining values of the gauge fluxes define non-Abelian phases with |W | < 2, whose lowenergy properties shall dramatically depend on the particular values of the fluxes. Below, we study the non-Abelian regime α = β 1 close to the Graphene point P 1 , where the low energy dispersion is linear and saddle-point in orthogonal directions. In this respect, the effective mass of quasiparticles at large wavelengths is highly anisotropic, and leads to novel features in the transport properties.

B. Effective low energy theories in Abelian regimes
In this section, we study the emerging quantum field theories of the non-Abelian honeycomb lattice at half-filling, centering our attention to the Abelian regimes introduced above.
(3) and shown in Fig. 2, with the peculiarity that each band is doubly-degenerate. As a consequence, the number of independent Dirac points is raised to N d = 4, fulfilling around the singularity K, and around K . As expected, the species of emerging massless fermions is doubled due to the color degrees of freedom.
R u 3 -reversed graphene: In the Abelian phase corresponding to α = 0, β = π (point P 2 in Fig. 6), the Hamiltonian can be obtained from graphene's analogue in Eq. (9), after reversing the hopping amplitude along u 3 (i.e. t → −t). The energy bands fulfill whose contour clearly shows that the pair of conical singularities are no longer at the corners of the Brillouin zone, but well inside it (see the contour energy in Fig.3(d)). The pair of conical singularities, now located at , correspond to the corner Dirac points of graphene, which have been translated in momentum space along 120 o direction. Such interesting transport of Dirac cones around the Brillouin zone is schematically represented in Fig. 3(c). Here, transparent cones represent the initial Dirac fermions for vanishing fluxes α = β = 0, and opaque cones represent the transported cones at α = 0, β = π. The direction of transport is shown by black arrows. The effective low energy theory at half-filling around these points is also captured by the relativistic Dirac equation, with color-degenerate massless fermions (i.e. a total of N d = 4 fermions), that behave according to the Hamiltonian around the singularity K, where a different representation of the Clifford algebra is chosen α x = −σ y , α y =σ x , wherẽ and lies in the complex unit circle. Around K , the low-energy Hamiltonian leads to α x =σ y , α y =σ x . R u 3 R u 1 -reversed graphene: Considering the gauge fluxes α = β = π corresponding to the point P 3 in Fig. 6, the Hamiltonian is analogous to graphene with the reversed hoppings t → −t along the directions u 1 , u 3 . The underlying energy bands are and lead to the conical singularities depicted in figs. 3(e) and 3(f). The pair of conical singularities lie at K = 2π , and correspond to the corner singularities translated in momentum space along a 180 o direction. At halffilling, the low energy excitations around such points are also relativistic massless fermions with α x = σ y , α y = σ x for K, and α x = σ y , α y = −σ x for K .
R u 1 -reversed graphene: Considering the gauge fluxes α = π, β = 0 corresponding to P 4 in Fig. 6, the Hamiltonian is analogous to graphene with the reversed hoppings t → −t along the direction u 1 . The underlying energy bands become represented in figs. 3(g) and 3(h). The pair of singularities lie 3 ), and correspond to the corner singularities translated in momentum space along a 60 o direction. At half-filling, the low energy excitations around such points are also relativistic massless fermions with α x = −σ y , and α y =σ x for K, whereσ j correspond to Eq. (15) for w = 1 2 (1 + i √ 3). Equivalently, the K emerging fermion leads to the Clifford algebra α x =σ y , and α y =σ x .
Let us finally remark the fact that all the Abelian phases described so far simply correspond to finite translations of the graphene points, and thus present the same number of emerging fermions N d = 4. Accordingly, the physical properties of the two-dimensional Fermi gas at half filling are expected to be analogous to graphene. The unique effect of a non-vanishing Abelian flux α, β is to transport the relativistic fermions around the Brillouin zone. Below, we study richer situations where the topology of the Fermi surface is modified by a variation of the number of gapless fermionic species. It is also important to note that these Abelian phases can also be achieved by simply tuning the anisotropy of the hopping parameters [9], without the need to insert external fields. Nonetheless, the phases to be discussed below cannot be achieved in this fashion, and constitute a non-trivial consequence of the applied gauge fields.
Π-flux point: Let us now focus on the Abelian regime α = β = π 2 , which corresponds to the Π-flux phase (i.e. lattice fermions pick up an overall minus sign by hopping around the elementary plaquette). This regime corresponds to the point P 5 in Fig. 6, where the fermionic Hamiltonian induces a color flip transformation when the fermion hops along directions u 1 and u 3 . The energy bands have been represented in Fig. 2(b), where one observes how the number of conical singularities is increased with respect to the previous Abelian phases. Let us note that the colordegeneracy is lifted due to the color-flip induced hopping. The distribution of conical singularities in the Brillouin zone can be better appreciated in figs. 3(i) and 3(j), where N d = 8 independent Dirac points appear in positions We conclude that the topology of the Fermi surface in the Π-flux phase is not equivalent to the remaining four Abelian phases described so far. Indeed, the number of Dirac fermions is distinct, which reveals the existence of processes, deeply rooted in the non-Abelian phase, where fermions scatter and pair-production (annihilation) is allowed. Such processes shall be discussed in detail in the forthcoming section, where a topological argument will clarify why such interacting events occur even when the Fermi gas is non-interacting.

C. Effective low energy theories in a non-Abelian regime
In order to identify the effects of the non-Abelian fields, we focus in the regime α = β 1, where the topology of the Fermi surface can be worked out analytically. We note that a variety of different and interesting non-Abelian phases also occur for different fluxes (α, β ), leaving their numerical treatment for the following sections. The energy band structure can be expressed in terms of the following functions These energy bands have been represented in Fig.2(c) for the particular case of α = 0.2, where one readily observes that the non-Abelian features of the external fields dramatically modify the corresponding conical singularities. The valence and conduction bands each split into two components, and the conical singularities are displaced from the corners of the Brillouin zone (see also Fig. 3(k) and 3(l)). However, we note that the two initial (two-fold degenerate) Dirac points give birth to four (non-degenerate) independent Dirac cones, so that N d = 4 remains constant in this non-Abelian regime. The exact position of the displaced Dirac points depends on the non-Abelian flux α as dictated by the following expressions which reduce to the Abelian points for α = 0 as expected. Once the exact location of the non-Abelian Dirac points has been identified, it is possible to linearize the Hamiltonian in Eq. (6) and obtain the continuum theory describing lowenergy excitations. In the particular case of p =hδ k = h(k − K 1 ), we get the free massless Dirac Hamiltonian with a perturbation ∆H describing spin-and color-flip processes that arise due to the non-Abelian fluxes Here, we have also defined with the following parameters η ± = (ξ ± 1 − ξ 2 )e iπ/4 , and ξ = (−1 + √ α 2 + 2α 4 )/2(1 + α 2 ). From the expressions above, one can check that in the vanishing flux limit α → 0, the perturbation also vanishes ∆H → 0, and we recover the Abelian graphene theory. Let us remark here that the non-Abelian perturbation is essential to account for the displacement of the Dirac points. As follows from Fig. 4(b), the dispersion relation around the singular points is not simply conical. Indeed, at low energies, the spectrum becomes linear in one direction, and displays a saddle point in the transverse direction. In this respect, we can expect a dramatic modification of the massless fermions properties due to their coupling to these exotic non-Abelian fields.
We conclude that the effect of the synthetic gauge fluxes is two-fold: they modify the Fermi surface by transporting the Dirac points according to Eq. (27), and they are also responsible of a different low-energy theory in Eq. (28). Indeed, the specific form of the perturbation in Eq. (29) suggests the existence of certain emerging gauge fields that couple to the bare massless fermions and give rise to this two-fold effect. A similar effect was explicitly shown in graphene-like structures with a Kekule distortion [66,67].

IV. TOPOLOGICAL QUANTUM PHASE TRANSITIONS
In the previous sections, we have thoroughly described the pattern of Dirac points and emerging relativistic fermions for different Abelian and non-Abelian phases. These phases differ by the total number of Dirac fermions and their distribution within the Brillouin zone. Below, we characterize the topological properties associated to the pattern of massless fermions by means of a topological charge (i.e. topological invariant) µ = ±1 defined for each conical singularity [62]. Moreover, the topological properties of the Fermi surface for non-interacting Fermi systems fully determine the underlying quantum order, the paradigm of order in quantum states beyond the Landau symmetry-breaking description [11,68]. In this respect, the pattern of topological charges classifies the phases introduced above within different quantum-ordered universality classes having the same symmetry.
We also study how a modification the gauge fluxes (α, β ) connects different quantum orders (i.e. different patterns of topological charges). In this respect, non-Abelian gauge fields induce a quantum phase transition (QPT) (i.e. a phase transition at T = 0 driven by quantum fluctuations) between different vacua with the same symmetry, and cannot be classified by the Landau symmetry-breaking paradigm. This class of QPTs have been previously studied in gapped systems such as the fractional QHE [69], or gapless strongly correlated systems, such as certain types of spin liquids [70,71]. In the sections below, we show how such QPT occurs in the realm of Fermi gases subjected to non-Abelian fields, which can be probed in a cold atom experiment.

A. Topological characterization: winding numbers
The topology of the Fermi surface determines the type of effective theory and emergent symmetries at low energies. In order to classify the different universality classes in Fermi systems with gapless excitations, a topological winding number can be defined for the general Hamiltonian H = ∑ k Ψ † k H k Ψ k , where H k was defined in Eq. (6). This Hamiltonian fulfills a particle-hole symmetry {Γ, H k } = 0, where Γ = diag(I 2 , −I 2 ), and allows us to define the following winding number where C is a curve defined in momentum space where the Hamiltonian H k presents no zeros [62]. Such a winding num-ber, or the equivalent version where h k = −t ∑ j e −iku j U † j , has three important properties: Dirac fermions: The zeros of H k , and thus conical singularities, are detected inside the region delimited by the curve C, by means of a non-vanishing winding number µ = 0.
Homotopic invariant: The winding number is a homotopic invariant, and is thus independent with respect to distortions of the loop C as far as H −1 k exists. Topological invariant The winding number is a topological invariant, and is thus independent with respect to local perturbations of the Hamiltonian At half-filling, the Fermi surface is reduced to a set of points and we can characterize its topology by the pattern of conical singularities and their corresponding topological charges (K j , µ j ). As discussed in the next section, varying the gauge fluxes (α, β ) along a non-Abelian path γ : (α 0 , β 0 ) → (α f , β f ), such that |W | γ = 2, we observe that different quantum orders are connected by quantum phase transitions. We have first determined the pattern of topological charges in each phase, and the corresponding universality classes: Graphene universality class: The quantum order in this equivalence class is determined by two positive topological charges µ = +1, and two negative topological charges µ = −1 within the first Brillouin zone (see Fig. 3(a)). There are thus N d = 4 Dirac points, which can be distributed in different positions in momentum space (see figs. 3(a)-3(h)), maintaining the same quantum numbers and thus belonging to the same universality class. We can thus conclude that all the graphene points P 1 , P 2 , P 3 , P 4 describe the same quantum order.
Π-flux universality class: The quantum order in this equivalence class is determined by four positive topological charges µ = +1, and four negative topological charges µ = −1 within the first Brillouin zone (see Fig. 3(i)). Therefore, a total of N d = 8 species of gapless fermions arise.
As discussed so far, the different quantum orders are all determined by a distinct even number of gapless fermions. The adiabatic variation of the external gauge fluxes will give rise to processes, where the number of massless fermions is modified. Note however, that the total number of massless excitations shall always remain even, a distinctive feature of all the phases discussed above. A different prescription occurs for topological insulators with time-reversal invariance, where the number of massless excitations can be generally odd [72][73][74][75][76][77][78][79][80][81]. We believe that it would be interesting to study this type of topological phase transitions in such systems.

B. Abelian path: absence of a quantum phase transition
In this section, we study the Fermi surface topology as the gauge fluxes are modified along an Abelian path γ : (α 0 , β 0 ) → (α f , β f ) such that |W | γ = 2. According to the previous section, this path only joins regimes belonging to the same uni-versality class, that of graphene. Therefore, no phase transition occurs, and modifying the gauge flux only transports the N d = 4 topological charges in momentum space. We focus on the following Abelian path γ 1 (t) = (tπ, 0), with t ∈ [0, 1], connecting the 1 st graphene point P 1 to the 4 th graphene point P 4 , as represented in Fig. 6. In figs. 7(a)-7(h), the distribution of conical singularities along the Abelian path γ 1 is represented.
In these figs. 7(a)-7(h), the Brillouin zone has been depicted by a dashed hexagon, and the conical singularities correspond to the dark-blue points where the valence and conduction bands touch. This set of contour figures represent the conduction band energy obtained by diagonalizing the Hamiltonian in Eq. (5) for different fluxes along the path γ 1 . It clearly follows that, starting from graphene P 1 in Fig. 7(a), the initial doubly-degenerate topological charges are split into a pair of Dirac points that travel in opposite directions (see Fig. 7(b)). Note however that the total number of Dirac points N d = 4 and topological charges µ = {−1, −1, +1, +1} are always preserved from Fig.7(a) to Fig.7(h). Having the same pattern of topological charges, the different regimes along the Abelian path all belong to the same universality class, and have the same quantum order. The effect of a varying gauge flux along this path is simply the transport of topological charges in momentum space. Therefore, no QPT connecting two different quantum orders exists, a fact that also holds for any other Abelian path (e.g. that joining P 2 to P 3 in Fig. 6).

C. Non-Abelian path: topological quantum phase transition
From the previous section we conclude that the Abelian regime must be abandoned if different quantum orders are to be connected. Below, we show that a non-Abelian flux quench provides a rich scenario where scattering processes occur, involving the creation or the annihilation of massless Dirac fermions. Moreover, the pattern of topological charges is modified accordingly. Let us emphasize the exotic nature of these interaction processes, which occur in a non-interacting gas. As shown below, only scattering of fermions with opposite topological charges is allowed. Moreover, we can draw an interesting analogy for two types of topological QPTs, which can be driven by a stimulated or spontaneous creationannihilation of opposite topological charges.

Stimulated creation-annihilation of Dirac fermions:
Here, we study a non-Abelian path γ 2 that joins the 1 st graphene point P 1 to the Π-flux regime P 5 , and then to the 3 rd graphene point P 3 (i.e. γ 2 (t) = (tπ,tπ), with t ∈ [0, 1], as represented in Fig. 6). This path is clearly non-Abelian, since |W | γ 2 = 2 except at t = 1 2 where the Π-flux regime is reached, and, as shown below, the pattern of topological charges determining the underlying universality class is abruptly modified for certain critical fluxes (α c , β c ).
In figs. 8(a)-8(h), we observe how the initial doublydegenerate Dirac points ( Fig. 8(a)) are split into a pair of fermions which move along the vertical axis in opposite directions ( Fig.8(b)). Remarkably, in Fig. 8(c) we can observe how additional Dirac fermions are being created along the borders of the Brillouin zone. A closer inspection of the energy band dispersion around k = 2π 3a (−1, −1 √ 3 ) (i.e. the region where these gapless fermions are created), reveals that there is a process of induced or stimulated pair production (see figs. 9(a)-9(d)), where a couple of massless fermions is created in the presence of a third gapless mode. Interestingly enough, the pair of created fermions carry opposite topological charges µ = ±1, so that the total winding number is conserved during the scattering process (see also Fig.9(e)). Note that the production of fermionic gapless modes occurs at T = 0, and thus a QPT between different quantum orders takes place. In the same manner symmetry-breaking phase transitions lead to gapless Goldstone bosons, QPTs driven by fluctuations of quantum order can give rise to gapless fermionic excitations. Also, analogous to the symmetry protection of Goldstone modes, these massless fermions are protected by quantum order (i.e. topological invariant charges µ = ±1).
The quantum phase transition taking place around α c = β c = π 4 , connects two different phases with N d = 4 and N d = 8 low energy Dirac points, so that the total topological charge is conserved. Let us also note that there is a critical point at α c = β c = 3π 4 where the inverse process takes place. In this case, two massless fermions with opposite charges are annihilated and connect the phases with N d = 8 and N d = 4 fermions. It is also interesting to remark that at criticality α c = β c = π 4 , one can show that the massless excitations are highly anisotropic, being relativistic in one direction, and nonrelativistic along the transverse axis.
In figs. 11(a)-11(d) we study the energy band landscape as the gauge fluxes are swept across the critical point α c = 2β c = 0.63π. Initially, the low energy excitations are clearly gapped and for low enough energies we have a local vacuum. Then, as the gauge fluxes increase, we observe how the local gap vanishes and two massless fermions appear. The topological charge distribution is schematically shown in Fig. 11(d), where the total topological charge is once more conserved. Let us note that for α c = 2β c = 0.7π (see Fig. 11(f)), the timereversed scattering process takes place. Here, we observe a pair fermion-antifermion being annihilated due to their opposite charges. These events of spontaneous pair creation or annihilation connect two distinct quantum orders N d = 4 → 8 by means of a non-symmetry breaking quantum phase tran-  Fig. 8(c), and the reverse annihilation process (time-reversed event) corresponding to Fig. 8(f). sition. A similar spontaneous creation-annihilation process is also illustrated in fig. 16 (see Sec. V). In the following sections, we address the issue of characterizing the above topological phases through accessible measurement techniques in the optical-lattice setup.
Let us close this section by comparing the topological QPT presented in this work, to the semi-metal to insulator transition in bare graphene-like systems [9,[82][83][84][85][86][87]. In graphene, a modulation of the hopping anisotropy induced by strain can lead to the merging of Dirac points, which turns graphene into an insulator. An insulating gap can also be dynamically induced by quenching the sublattice energy imbalance [86]. Conversely, in the non-Abelian honeycomb lattice, the merging of Dirac points does not lead to a global gap, and the semimetallic phase is always preserved. Accordingly, this opticallattice setup offers a broader spectrum of topological QPTs, where both spontaneous and stimulated creation-annihilation of Dirac fermions occur.

V. THE QUANTUM HALL MEASUREMENTS
The presence of several distinct Dirac points within the first Brillouin zone leads to dramatic consequences in the physical properties of the system. In the low-energy regime, these relativistic effects are adequately described by the underlying effective Dirac theories. An interesting effect to be considered in this context is the so-called anomalous quantum Hall effect (AQHE) [21,22], which has been discovered in graphene sheets subjected to external magnetic fields [23][24][25][26]. This effect is characterized by the peculiar sequence of plateaus in the transverse conductivity of the system as a function of the chemical potential: in the low-energy regime, where the spectrum exhibits the so-called relativistic Landau levels [17][18][19], the Hall conductivity is quantized according to σ xy = ± 4 h (ν + 1 2 ), where ν is an integer, h is Planck's constant and we set e = 1. The Hall measurements show steps of four that can be traced back to the presence of two independent two-fold degenerate Dirac points N d = 4, and a half-integer anomaly due to the zero-energy mode associated to the selfconjugate lowest Landau level.
In previous sections, it has been shown how several Dirac points are induced in our cold-atom system. In particular, we have shown that the number of independent Dirac cones N d is not necessarily conserved as the non-Abelian fluxes (α, β ) are varied. The aim of this section is to show that the number of emerging massless fermions N d can be inferred from Hall measurements. In particular, we show that the Hall plateaus at zero Fermi energy E F = 0 fulfill and thus provide a Dirac fermion witness. Moreover, the robustness of this witness is associated to the system zeromodes, and is guaranteed by the Hamiltonian particle-hole symmetry and topological index theorem [21]. The relation of the Hall conductivity to the number of Dirac fermions in Eq. (33) can be conjectured from the theoretical results in [21][22][23][24]. Below, we confirm this conjecture by numerically computing the Chern numbers associated to the energy band structure, and thus estimating the corresponding Hall conductivities. Furthermore, this conductivity can be experimentally accessed in cold atoms by the so called Streda formula [88].
In this section, we shall first review how such measurements can be performed in cold-atom experiments by analyzing the density profiles. Then, we show how the AQHE taking place in this system gives us crucial information on the topological QPTs that occur as the non-Abelian fluxes are varied. In particular, we discuss under what conditions the Hall measurements are indeed witnesses of such transitions.

A. The anomalous quantum Hall effect and the density profiles
In order to provoke the quantum Hall effect in our system, one has to open the central gap. This can be achieved by subjecting the system to an additional Abelian gauge field (i.e. an artificial magnetic field). Such synthetic magnetic fields can be produced for neutral atoms by rotating the lattice [33][34][35][36][37], by considering laser-induced methods [8,[28][29][30][31][32], or even by immersing the system in a rotating BEC [38]. The presence of an external Abelian gauge field modifies the hopping operators introduced in Sect. III, which now take the form where Φ is the effective magnetic flux quanta per unit cell. A clear manifestation of the external magnetic field can already be found by computing the energy spectrum: when the flux is given by the ratio of two integers, namely Φ = p/q, the initial band structure is split into several subbands. For such rational Abelian fluxes, we note that the Hamiltonian takes the form of a 4q × 4q matrix. Therefore the number of gaps highly depends on the integer q, but also on the specific values of the non-Abelian fluxes α and β . When plotted as a function of the flux Φ, these numerous subbands form a fractal structure which reminds the well-known Hofstadter butterfly [89], as depicted in Fig. 12 for the Abelian case α = β = 0. In the cold-atom framework, the conductivity tensor measures the response of the system to an external perturbation, such as a lattice acceleration, which generates a current. In the presence of an additional magnetic field, a Hall current is generated and this transverse transport is dictated by the Hall conductivity σ H . When the Fermi energy E F lies in a spectral gap, the Hall conductivity can be expressed as a sum of topological invariants where N ch (E n ) is the Chern number associated to the filled energy band E n , which is necessarily an integer [90]. Consequently, the sum in Eq. (35) characterizes the different topological phases associated to the spectral gaps. These quantum Hall phases can be reached by varying the Fermi energy E F , namely, by controlling the atomic filling factor. The computation of the Chern numbers can be achieved numerically through diverse techniques [40,58,91] and leads to the complete description of the QHE in our system. In Fig. 13 (a), we show the Hall plateaus in the Abelian case α = β = 0, as recently discovered in graphene. One clearly distinguishes the Dirac regime, situated in the range E F ∈ [−1, 1], where the anomalous double steps and a half-integer plateau are observed. It should be stressed that the range of this relativistic regime, where the dispersion relation is linear, is known to be bounded by singularities in the density of states (see Fig. 13 (b)).
Although optical-lattice experiments involving atomic currents are on their way, it is preferable to exploit another method to evaluate the Hall conductivity. Indeed, several works have recently proposed that the quantized Hall conductivity could be measured through an elegant analysis of the density profiles [88,[92][93][94]. Furthermore, this efficient method takes into account the presence of a harmonic potential V harm (r) that generally confines the atoms within the optical lattice. If this external potential varies slowly compared to the lattice spacing, one may use the local-density approximation (LDA) and define a local Fermi energy as The atomic density is then simply obtained through the expression n(r) = dε D(ε) Θ E F (r) − ε , where D(ε) is the density of states of the uniform system. Each time the Fermi energy E F (r) falls in a gap, the density n(r) depicts a plateau: there is a one-to-one correspondence between the spectral gaps represented in Fig. 12 and the plateaus of the density profiles n(r). We illustrate this phenomenon in Fig. 14, where two density profiles for Φ 1 = 1/157 and Φ 2 = 2/157 are represented, and where the correspondence with the gaps in Fig. 12 (b) has been highlighted with symbols. In this figure, we focus on the relativistic region around E F = 0 and set α = β = 0. Note that in the following, the density per plaquette is studied as a function of the Fermi energy n(E F ), being density profiles n(r) simply obtained by considering the LDA relation in Eq. (36). In the upper pannel of Fig. 14, a typical density profile n(r) obtained for Φ 1 = 1/157 is shown: the relativistic regime can already be identified by an unusual behavior (see inside the red rectangle). Having identified the unusual plateaus in the density profiles around E F , one may now evaluate the Hall conductivity through the Streda formula [88], which takes the form First of all, let us show how the Streda relation (37) allows one to detect the AQHE at the Graphene point P 1 (α = β = 0). Comparing the two density profiles in Fig. 14 Fig. 13 (a), the Hall conductivity undergoes the sequence σ H h ∈ {..., −10, −6, −2, 2, 6, 10, ...}, and therefore agrees with the relativistic prediction σ H = ± 4 h (ν + 1 2 ). Moreover, this example shows how a precise analysis of the density profile measurements obtained for two close values of the flux Φ allows one to observe the AQHE proper to the relativistic Figure 15: Two non-Abelian paths starting from the Graphene point (path I) and from the π-flux point (path II). The background contour corresponds to |W (α, β )| regime. Furthermore, these density profiles constitute a powerful tool to determine the number of Dirac cones situated inside the first BZ, since the latter is directly related to the jumps observed in the anomalous Hall sequences N d = 4. This important result will be exploited in the following section to study how the density profiles are modified as one varies the non-Abelian fluxes α and β . At this point, it is worth noticing that the density profiles, and thus the Hall plateaus, are not affected as one travels along an Abelian path (see Sect. IV B). The latter fact is in agreement with the fact that N d = 4 remains constant along such a path and that no topological QPT is thus observed.

B. The Hall plateaus along a non-Abelian path
In this section, we consider the two different paths within the α − β space in Fig. 15. The first path starts from the Π-flux point (α = β = π/2) and reaches α = 0.9(3π/4), β = π/2 in four steps (see the purple dots in Fig. 15). The second path starts from the graphene point (α = β = 0) and reaches α = β = π/4 in six steps (see the blue dots in Fig. 15). Both paths teach us in what interesting manner the density profiles and the Hall plateaus are modified as the number of Dirac cones N d changes.
Path I.-Let us first study the path starting at the Π-flux point (see path I in Fig. 15). In this path, the parameter β = π/2 is constant, while α is varied from π/2 to 0.9(3π/4). The corresponding energy spectra at vanishing flux Φ = 0 are represented in Fig. 16, where one observes how the eight initial isotropic Dirac points are progressively deformed as the flux is varied. As one travels between the Point 2 (α = 3π/5) and the Point 3 ( α = 0.87(3π/4)), one clearly observes that four Dirac points approach and finally merge at Point 4 α = 0.9(3π/4) (see Fig. 16 (b) and (c)). Indeed, at criticality α c = 3π/4, the massless fermions merge into a highly anisotropic excitation, the latter being linear in one direction and parabolic in the transverse. Below, we explain how this topological phase transition is detected through density measurements and the AQHE. Note that we focus on very low fluxes Φ, which allows one to discuss the Hall plateaus in terms of the energy spectra in Fig. 17.
• Point 1: α = π/2. We first compute σ H through the density profiles at the Π-flux point. The van Hove singularities (VHS) in the density of states (DOS) delimit the relativistic regime between E V HS = ±0.4, where we observe steps of eight in the Hall conductivity σ H = ± 8 h (ν + 1 2 ). This numerical computation emphasizes the presence of the N d = 8 massless Dirac fermions, as observed in Fig. 16 (a).
• Point 2 to Point 4: α = 3π/5 → 0.9(3π/4). As the non-Abelian flux α is varied, we observe a transition from Consequently, the density profile measurements performed along this path witness the topological phase transition characterized by the change N d = 8 → 4.
The above analysis confirms the aforementioned conjecture announced in Eq. (33), and thus proves that density measurements do indeed witness the topological phase transition.
Path II.-Let us then focus on the second path (see path II in Fig. 15). The corresponding energy spectra at vanishing flux Φ = 0 are represented in Fig. 17, where one observes that the two initial 2-degenerate Dirac points are progressively split into four non-degenerate singularities, and at some point become anisotropic (see Fig. 17 (d) and (f)). We study below • Point 2, α = β = π/16: as the degeneracy is lifted, four (non-degenerate) Dirac cones are present in the Brillouin zone. Our numerical results show that the Hall plateaus exhibit steps of four σ H = ± 4 h (ν + 1 2 ) in the range E ∈ [−0.11, 0.11], steps of two in the range E ∈ [−1, −0.11) (0.11, 1], and steps of one elsewhere |E| > 1. These three regimes are delimited by VHS in the DOS at E VHS 1 = ±0.11, and at E VHS 2 = ±1. This analysis suggests the existence of three different regimes, namely, a relativistic regime with four isotropic Dirac cones, an intermediate relativistic regime that is characterized by the reminiscent of the two initial cones, and finally a non-relativistic regime where the linear dispersion vanishes. Additionally, we have numerically confirmed that the central plateau at E F = 0 (i.e. σ H = −2 → 2) is robust, which is in accordance with symmetry and topological considerations on the corresponding zero-mode [21]. Therefore, the zero-energy Hall conductivity is an ideal witness of the number of emerging massless fermions.
• Point 3 to Point 6, α = π/8, β = π/16 → α = β = π/4: along this path, one observes the interesting fact that the density profiles, and thus the Hall plateaus, are affected by the anisotropic character of the Dirac cones. Hence, unusual behavior is observed in the relativistic regime where the conductivity evolves with steps of two in the anisotropic regime (cf. Figs. 17(d) and (f)), whereas it evolves with steps of four in the isotropic regime (cf. Figs. 17(c) and (e)). One may conclude that the anisotropic cones play a peculiar role in the anomalous quantum Hall effect: these anisotropic singularities do not contribute to the Hall conductivity for E F = 0.

VI. EXPERIMENTAL CHARACTERIZATION OF TOPOLOGICAL PHASES
In the previous section, we have described in detail a method of detection of the phenomena predicted in this paper using density profile measurements. This is the easiest and the most convenient method, and with current experimental setups can be realized with the accuracy of, say, 5 percent. It is reasonable to expect that by optimizing this method, accuracies of 1 percent are achievable. However, since the detection in quantum many-body systems is one of the most challenging tasks, we discuss in this section other possible measurements strategies to detect topological phases and Dirac fermions. We present the state-of-art of all methods and discuss their precision. All the detection methods discussed below have at least some of the following properties: • Spatial resolution. That means that the methods allow to resolve spatial (or equivalently Fourier transformed momentum) characteristics of the system.
• Color resolution. Note that the internal degrees of freedom used for decoding the color states will correspond to internal spin, or pseudo-spin states. Many methods will allow to resolve color states spectroscopically.
• Non-demolition property. That means that the methods will affect the measured quantum mechanical systems in a "minimal fashion", and can be even repeated to allow for time resolution. Only few methods have this property, such as contrast imaging of atomic density and spin polarization spectroscopy (see below).
Detection methods can be divided into two groups: those which address directly the properties of Dirac fermions, and those that have more general character and provide general properties of the system. a. Direct detection of Dirac fermions These methods can be further classified as: i) global detection mechanisms, that allow us to infer the total number of massless fermions (this is the reason for calling them global), and ii) local detection mechanism, that should allow us to map the whole Fermi surface, and thus locate each of the conical singularities.
Among the global detection mechanisms we can list: • Measurement of the atom density as a function of the chemical potential: Following [9], a measurement of the atom density close to zero chemical potential depends on the number of Dirac points in the Fermi surface. In Ref. [9] it was proposed to detect masslessto-massive Dirac particles transition; here we propose to use it to detect the change of number of the Dirac fermions. In principle, this method requires nothing but a precise measurement of the density, as in the case of the method discussed in the previous section. The seminal steps toward the precise density measurements were achieved quite recently in the first observations of the wedding cake structure of the Mott insulator -superfluid transition in the optical lattice loosely confined in an harmonic trap [95,96]. Perhaps the most recent stateof-art paper is Ref. [97] in which validation of the quantum simulator of MI-SF transition is realized. Here, the most accurate Quantum Monte Carlo simulations are directly compared with experimental data; in this way accuracies of 1 percent in density determination are reached.
The local detection mechanisms are obviously more demanding; among them we can point out: • Atomic angle-resolved photoemission spectroscopy ARPES method: a possible candidate for the adaptation of ARPES techniques to optical lattice setups, is based on using Raman transitions to a side level that is decoupled from the system [98]. This is a very powerful method that allows to determine directly the spatially resolved single particle density matrix, i.e. indirectly to map the Fermi surface. This new measurement technique has been realized for ultracold atom gases in Ref. [99]; like photoemission spectroscopy for electronic materials, it directly probes low energy excitations and thus can reveal excitation gaps and/or pseudogaps with accuracy less than 5 percent. Atomic ARPES is a destructive method, but in principle allows for color resolution (by splitting the relevant states in an external fields and using the spectral resolution), i.e. it allows to measure ρ st (x, x) = Ψ s †(x)Ψ t (x) , where s,t are the color indices.
• Direct band mapping is an alternative approach [100][101][102]. In this method one releases adiabatically the optical potentials to measure directly the quasi-momentum distribution. This can be combined with contrast imaging, so that can be combined with other measurements, such as density measurement in the same experimental run (see for example [6]). This method allows to map higher bands independently, which in principle would determine positions of the Dirac points directly. Also, it can incorporate the color resolution.
b. Indirect detection of Dirac fermions Apart from these direct methods, there is a whole variety of general detection methods that will provide additional information and characteristics of the system. Although such measurements will not give information about the number of Dirac points, they will characterize additional interesting aspects of the system: • Density imaging: Of course, the density measurements are the necessary ingredient of the direct methods based on Streda formula, or dependence on chemical potential, discussed above. Here, we consider them in general sense as a source of additional information. The density imaging allows for color resolution, and can be done destructively using absorptive imaging, or nondestructively using contrast imaging. It is used in the band mapping technique, but can be also used if we turn the optical potentials suddenly off to measure the momentum distribution, which can give a particularly interesting characterization of the considered states of the system.
• Noise interferometry: this very powerful (although destructive) method, proposed in Ref. [103], can be used to measure directly the two-point density-density correlation functions, which can be expressed in terms of Fourier transforms in the (quasi) momentum space of occupations of the corresponding single particle states.
Originally applied to bosons [104], it is very useful to characterize fermionic states [105]. It gives indirect information about density of states close to the Fermi energy.
• Hall conductivity measurements: this is discussed in detail in Sec. V, where we provide very detailed examples for predicted results of such measurements. This method is in the first place an indirect method. Still, since the anomalous quantum Hall effect is highly affected by the Fermi surface topology at E F = 0, the Hall conductivity, determined from density profiles measurements as in Refs. [88,[92][93][94], provides direct information about the change in the number of Dirac fermions.
• Spin polarization spectroscopy: This method (see [106] and references therein) may be applied to measure Fourier component of the spin density and their fluctuations. The methods is based on the quantum Faraday effect, i.e. coupling the polarized light beam to atoms, mapping the state of atoms onto light, and then "reading" off the information about atoms by measuring the fluctuations of the light beam using homodyne detection. It combines color and spatial resolution with non-demolition character, and is thus very promising and powerful. This method is very well developed experimentally (for a review see [107]), and the experiments with spatial, or, better to say, momentum resolution have been initiated in the group of E. Polzik. Similarly, as noise interferometry, it gives indirect information about density of states close to the Fermi energy.
• Momentum-resolved Bragg spectroscopy: this method was pioneered by W.D. Phillips [108] and W. Ketterle [109,110] groups. It allows to measure the structure factor, spin structure factor and other momentum and frequency resolved correlation functions. It also gives direct information of the dispersion relation of lowenergy (quasi) particle excitations. The current stateof-art of this method is presented in Ref. [111], where the accuracies of few percents are reached.
• Methods involving cavity QED: these ideas were mostly developed by H. Ritsch and co-workers ( [112,113], see also [114]), but have not yet been realized experimentally. They are somewhat similar to ARPES, with cavity-enhancing of the signal. In principle they are aimed on measuring the spatial density-density correlations by coupling the systems to the cavity mode, and measuring statistics of the light outgoing from the cavity.
Last, but not least, we should mention enormous progress in detecting atomic density in situ and in expanding clouds of a single atom resolution [115][116][117]. These methods will add on unprecedented precision to many of the above discussed detection schemes.

VII. CONCLUSIONS
In this article we have discussed the optical-lattice analogue of graphene subjected to additional non-Abelian gauge fields. We have shown how a long-wavelength perspective captures the relativistic nature of the excitations responsible for transport properties of the system. Besides, when the effective fields are switched on, the topology of the Fermi surface is modified, and novel phases with gapless fermionic excitations that differ from graphene arise. Such topological phases have been fully characterized by a topological charge assigned to each of the massless fermions.
We have also discussed how an adiabatic variation of the external fluxes along a non-Abelian path can be interpreted as a topological phase transition between different quantum orders. We have identified two processes responsible for such phase transitions, namely, spontaneous and stimulated scattering events where the number of gapless fermions is modified.
We have demonstrated that the creation and annihilation processes not only modify the Fermi surface's topology but also lead to interesting generalizations of the anomalous QHE. Moreover, we proved that a measure of this effect through atomic density profiles gives a direct signature of the number of Dirac fermions. In this sense, we emphasized that the diverse topological phases are characterized by topological charges (winding numbers), but also by a specific sequence of the underlying Chern numbers.
These systems provide a very challenging, but at the same time rich and interesting playground for applying various known detection techniques. Some of them provide direct information about Dirac fermions and topological phase transition, the other provide additional, but highly relevant characteristics. Despite the fact that quite a lot of methods exist, it is desirable to search for new ones especially suited for the systems considered here (attempting for instance to measure topological invariants directly).
We discussed here the physics of single particles. Obviously, it opens an avenue towards studies of even more challenging physics of interacting ultracold atoms in non-Abelian gauge fields, such as the fractional QHE or the Fermi colorsuperconductivity.
Summarizing, the effect of non-Abelian gauge fields on the emerging relativistic theory is far from trivial. Not only do such fields transport the topological charges attached to these Dirac points, but they also lead to scattering processes where pairs of fermions with opposite topological charges are created or annihilated, as a result of the transition between different quantum orders. The possibility to observe such exotic phenomena, and extensions thereof, in an optical lattice table-top experiment seems to be a challenging and rewarding task where interdisciplinary effects of condensed-matter, high-energy, and atomic physics can be observed.