{'fraction_non_alphanumeric': 0.03025729656045374, 'fraction_numerical': 0.008300674277646601, 'mean_word_length': 4.56968546637744, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulating Dirac fermions with Abelian and non-Abelian gauge fields in optical lattices

In this work we present an optical lattice setup to realize a full Dirac Hamiltonian in 2+1 dimensions. We show how all possible external potentials coupled to the Dirac field can arise from perturbations of the existing couplings of the honeycomb lattice model, without the need of additional laser fields. This greatly simplifies the proposed implementations, requiring only spatial modulations of the intensity of the laser beams. We finally suggest several experiments to observe the properties of the Dirac field in the setup.


Introduction
Lattice models have long been a powerful tool to study quantum mechanics and quantum field theories. They allow for computational treatment of analytically intractable problems, from high-Tc superconductors [1] to colour confinement in QCD [2]. However, the computational cost of dealing with these problems can be exorbitant [3]. In the spirit of Feynman [4], great effort has been devoted to simulating these models in intrinsically quantum systems, which can be more efficient in reproducing key characteristics of these problems. Indeed quantum simulators have managed to produce results in the study of relativistic quantum mechanics [5,6] and quantum phase transitions [7,8,9,10].
Optical lattices are one of the most promising candidates for quantum simulations, due to their scalability, tunability and their versatility in terms of geometry and dimensionality. They have already been successful in simulating quantum phase transitions of Bose-Hubbard and spin Hamiltonians and synthetic gauge theories [7,10,11,13], while further theoretical work has shown the possibility of addressing field theory and topological aspects of condensed matter phyisics [12,13,16,14,15]. In this work we propose a way to study in an optical lattice the 2+1 Dirac fields coupled to both Abelian and non-Abelian gauge fields. Our work connects with previous proposals for simulating Dirac fermions in various types of lattices [17,18,19], including massive fermions [20] and QED simulations using BEC for the bosonic fields [21]. The novelty of this work is that instead of artificially induced phases or bosonic baths, we rely on intensity modulations of the trap and of additional laser beams to give rise to all sorts of phenomena: from the appearance of synthetic electromagnetic fields and a Dirac mass to exotic flavour-coupling perturbations. Our idea is intimately connected to the way mechanical deformations induce effective magnetic fields in graphene [22,23,24], but makes use of the greater tunability of optical lattices to create a larger family of potentials.
The basic ingredient in this work is a setup consisting on two state-dependent triangular lattices, connected by a Raman laser. The spatial modulations of the intensity of the laser give rise to a position-dependent hopping which, as we show, is equivalent to a deformation of the lattice. The combination of this with other lattice perturbations allows us to produce effective Abelian fields (magnetic and electric potentials), non-Abelian fields (also known as flavour coupling terms), scalar fields and an effective mass. This setup is particularly interesting because it also helps us in the detection of the gauge fields, using the tools from Ref. [15], but the ideas put forward in this work can also be generalized to other lattices with Dirac fermions, such as the honeycomb lattice or the kagome lattice, which have been recently demonstrated in the lab [25].
This work is organized in a self-contained manner, evolving from well know developments (Sections 2-3) all the way up to the final implementation and detection. In section 2 we review how the Dirac equation arises from a tightbinding treatment of the honeycomb lattice. We pay special attention to the tools for deriving the continuum limit and to the appearance of an extra degree of freedom, usually called flavour [23]. In section 3 we analyse the effect of perturbations on the lattice, interpreting the result as a coupling between the Dirac field and scalar, Abelian and non-Abelian external gauge fields. In section 4 we present in detail our proposal for the experimental implementation of the Dirac model and of the external potentials using ultracold atoms in an optical lattice. After discussing a trapping scheme based on state-dependent optical lattices, section 5 shows how to implement the effective potentials from section 3 using intensity modulations of the laser beams. Section 4.3 elaborates on the experimental feasibility of the trapping scheme, estimating the Hamiltonian parameters and interaction strengths from band structure calculations. Following the setup, section 6 discusses different experimental protocols to characterize and measure the previous quantum simulations. Finally, in section 7 we discuss our proposal, relating it to recent proposals in the fields of optical lattices and graphene.

From tight-binding models to Dirac fields
Here we develop a theoretical framework in which 2D fields arise as the continuum approximation to a quadratic tight-binding model on a honeycomb lat- Figure 1: Two triangular sub-lattices displaced one with respect to the other to generate a honeycomb lattice. The vectors v 1 and v 2 are the hexagonal lattice generator vectors, whose length is the triangular lattice period, d = |v 2,3 |. In addition we define v 0 = 0 and v 3 = v 1 +v 2 . The unit cells are shown with dashed rectangles, and are labelled by m. Γ A and Γ B are the hopping parameters related with jumps between positions inside sub-lattices A and B while J i are the hopping parameters related with hexagonal jumps between sublattices, each one related with the corresponding v i , i = 0, 1, 2.
tice. Later in section 4 we show how these Hamiltonians describe non-interacting atoms in a strong periodic confining potential, but in this section the focus is on the abstract model. More precisely, we are interested in the derivation of honeycomb lattice band structure, whose excitations behave as relativistic Dirac particles, and how this procedure is modified to obtain a coupling between the Dirac particles and emergent gauge fields.

Tight-binding model for two coupled sub-lattices
Consider a system of fermionic particles in a strong confining periodic potential. We model these particles with fermionic creation and annihilation operators which describe the presence or absence of a fermion in a particular lattice site. The lattice is bipartite, which means that it can be divided into two disjoint sets of sites A and B, where the nearest neighbours of any A-site are all elements of B, and viceversa [cf. figure 1a)]. We consider nearestneighbour hoppings between different sublattices of a bipartite lattice, and also next-to-nearest-neighbour (or intra-lattice) hoppings. For the hexagonal spatial geometry shown in figure 1, our model Hamiltonian reads: where the J i are the hopping parameters related with the nearest neighbour hoppings; Γ A and Γ B model the tunneling effect inside sublattices A and B respectively; a † m and a m (b † m and b m ) are the creation and annihilation operators Right: zoom around the Dirac points where the energy bands take conical shape and the Hamiltonian of the system can be approximated as a Dirac Hamiltonian, being E ∝ |q|, equation (6).
for an atom in the position m of sublattice A (B), and the v's refer to the v vectors represented in figure 1b. The first line of equation (1) represents the hexagonal hoppings, and the second line represents the two kinds of triangular hoppings. Note also that J 0 , J 1 , J 2 , Γ A and Γ B have units of energy in this notation.

Energy bands and Dirac field
The honeycomb Hamiltonian is recovered from equation (1) by making J 0 = J 1 = J 2 ≡ J and Γ A = Γ B = 0: To calculate the energy bands, we do a Fourier transform over the creation and annihilation operators as where k is the momentum, r m is the position of the site m, the integration limit Ω makes reference to the first Brillouin Zone (1BZ) and a similar expression is used for b. The diagonalization of the Hamiltonian in momentum space leads to the eigenvalues of the energy that form the energy band structure; namely which are represented in figure 2. There are six points in the 1BZ for which E = 0, but only two of them are non-equivalent and we name them K ± = (∓ 4π 3 , 0)/d, where d is the triangular lattice spacing [26]. Assuming only small energy perturbations around E = 0, the excitations will be confined to the so called Dirac cones. To find the effective theory for these low-energy and longwavelength excitations one makes an expansion around the Dirac points of E(k), introducing q = k − K ± and assuming that for the relevant states it remains "small" (i.e. |q| × d 1). Under these conditions, it is possible to work with the Dirac cones as if they extended over all values of q, with integrals in d 2 k over the 1BZ being replaced by an integral in d 2 q over the whole momentum space. By transforming the operators back to position space via inverse Fourier one gets an effective theory for the continuum fields {φ a (r), φ b (r)}, satisfying the Dirac-like Hamiltonian (cf. equation 15) where σ = (σ x , σ y ) are the usual Pauli matrices and the subscript ± stands for the choice of cone (K ± ).

Derivation of the Dirac Hamiltonian in position space
While the previous discussion reveals the main ingredients of the relativistic fields which arise from our discrete Hamiltonian, two concerns appear. One is the added difficulty of treating spatially dependent perturbations, caused by the two consecutive Fourier transforms. The other concern is the little attention paid to the fact that there are two intrinsically different cones, which in line with previous literature we call the two "flavours" of Dirac particles in the lattice. As we will show later, this flavour plays a role in the simulation of non-Abelian fields.
In order to extend the derivation of the Dirac dispersion relation to setups in which translational invariance is weakly broken, we rely on a continuum-field approach that bypasses the use of momentum space [27]. Starting with our tightbinding Hamiltonian (2), we approximate the Fock operators at each lattice site as the value of a continuous field defined over all space, but which varies so smoothly that it is approximately constant over each unit cell. Moreover, since we are interested in excitations around the K ± quasimomenta, these fields are the envelope of quasi-plane wavepackets around those pointŝ with an equivalent expression forb m . Here d 2 r is the unit cell area and the exponential containing the Dirac points K ± implements the desired wavepacket ansatz. Substituting these operators in the tight-binding model produces the following limit Hamiltonian where the Greek indices stand for the two possible flavour choices (± cones). Note that the index i = 0, 1, 2 runs over the three vectors v 0 , −v 1 , v 2 depicted in figure 1 and which connect different unit cells. This is so, because we are free to define the smoothly varying fields Ψ a and Ψ b on the middle point of the unit cell.
Since we have enforced the fields to be slowly varying, the presence of the rapidly oscillating exponential e ir(Kτ −K ζ ) imposes the condition τ = ζ and we end up with two flavour-decoupled integrals. A more quantitative argument would compare the τ = ζ term with one that oscillates with K + − K − . Since both are proportional to J, we only need to compare their ratio. Integrating by parts the oscillating term and using periodic or open boundary conditions we may write the rapidly oscillating term correction as which becomes small if the fields oscillate slowly, that is d|∇Ψ| |Ψ|.

At this point one may expandΨ(r +
Our resulting Hamiltonian is therefore diagonal both in position and in flavour space but couples the internal degrees of freedom a-b via an off-diagonal operator with only one nontrivial term: Since 1 + e −iKτ ·v2 + e iKτ ·v1 vanishes due to the definition of the Dirac cones, the coupling term simplifies to Introducing a momentum operator q = −i∇, our coupling term becomes C τ ∼ τ q x + iq y and the Hamiltonian spatial density is therefore: which is the expression of the massless Dirac Hamiltonian with Fermi velocity c = √ 3Jd/2 . Note that changing flavour is equivalent to changing the sign of q x in the Hamiltonian.

Extending the Dirac Hamiltonian by modifying the lattice parameters
We now generalize the free Dirac Hamiltonian to include a variety of external fields [28] H D = c α · p + βmc 2 + βV cov .
Here m is the mass of the particles and V cov is the most general covariant potential containing scalar, vector, matrix, pseudoscalar, pseudovector and pseudotensor fields. The Dirac matrices α i = γ 0 γ i and β = γ 0 , are defined in terms of the generators of the Clifford group, γ µ . In 2+1 dimensions this is a set of 2 × 2 matrices satisfying the anticommutation relations {γ µ , γ ν } = 2η µν with η µν = diag(1, −1, −1), where µ, ν = 0, 1, 2. The choice α = σ = (σ x , σ y ) and β = σ z corresponds directly to the Hamiltonian in equation (6). The Dirac Hamiltonian then has the form where A 0 (r) and A(r) are the usual scalar and vector potential that give rise to observable electric and magnetic fields, and V (r) is a scalar potential that mimics the effect of an imposed mass; e is an effective charge. Note that since the fields A and A 0 do not have any dynamics, we may assume e = 1, for convenience. Let us now show how to recover all terms in equation (15) by slightly perturbing the tight-binding model.

Generating a mass term
The simplest term that we can add to Hamiltonian equation (2) is an energy difference between atoms in sublattices A and B: which in the continuum limit becomes Comparing this with equation (15), we see that the position-independent term proportional to σ z can be identified with an effective mass.

Abelian potentials
We can extend our model by starting from equation (2) and smoothly modifying the hoppings as In the continuum limit this renders In the notation of the previous section, this is equivalent to a change in the "coupling term" between pseudospins: therefore allowing us to formally derive an Abelian external potential where

Scalar fields
Consider now perturbations which are diagonal in the internal space to the tight-binding Hamiltonian, such as those given by intralattice hoppings in equation (1): Following the continuum limit performed in section 3.1 the Hamiltonian in momentum space around the K τ point can then be written as or purposefully rewritten (reabsorbing the constant i e −iKτ vi flavour-dependent factor) as which has both a term proportional to I representing an electric potential, , and a term proportional to σ z which contributes to the effective mass.

Flavour-coupling perturbations
We have seen that spatial variations of the hopping elements emerge as an Abelian external field. These variations have to be small, but it should be stated that, in order to keep flavours decoupled, they also have to be slowly varying. Otherwise, if there is a modulation with wavevector comparable to the order (K + − K − ), the flavour-coupling terms τ = ζ in equation (8) do not vanish. What we suggest now is to introduce small perturbations whose wavelength is comparable to the lattice constant and thus bridge the difference in momentum between cones. For instance with a constant coupling strength χ x,i,m = χ x . These rapid oscillations only allow the survival of terms that couple different cones, cancelling all terms inside the same cone. The most general perturbation of this sort is where χ x is the spatial dependence of the slow envelope that surrounds our perturbation (26). Actually, this envelope can be "remodulated" in order to make flavour coupling also weakly spatial-dependent. Moreover, we can also introduce in equation (26) sine terms, (r) = χ y sin [r · (K + − K − )], which make the flavour-coupling term complex.

The complete Hamiltonian
Since there exist no pseudo-potentials (γ 5 = 1 in 2 + 1 dimensions), our idea of perturbing the tight-binding model parameters allows us to reconstruct all possible external potentials of the Dirac equation, plus an additional coupling between different types of particles. With all these elements we have the following effective single-particle Hamiltonian:

Trapping of atoms in an optical honeycomb lattice
In this section we introduce optical lattices and explain how to construct a honeycomb lattice with two state-dependent triangular optical lattices. We discuss how to implement a tight-binding Hamiltonian and propose an experimental setup to obtain the Dirac field and the external fields using solely Raman and detuned lasers.

From two triangular sub-lattices to a honeycomb lattice
Recent advances in the development of high-aperture objectives and their integration in optical traps open the door to the generation of almost arbitrary two-dimensional potential landscapes for ultracold atoms. The basic idea is that off-resonant light may be used to tightly confine atoms in the maxima or minima of intensity, recreating sophisticated lattice models [29]. While until now those minima and maxima were generated through the interference of multiple laser beams [7,29], a novel paradigm consists on shaping and organizing those intensity profiles by simply projecting sophisticated images on the two-dimensional focal plane of a lens. The first experiments along this line have reproduced the usual square lattice quantum simulations [30] and also demonstrated the first triangular lattices [31].
In this work we are particularly interested in the last of those setups, which combines two triangular lattices [31] in the same plane. In this experiment the trapping laser beams are first diffracted by a holographic mask with a triangular pattern, selecting the first diffraction orders, which are then collected by a powerful lens to create the imprinted intensity pattern at its focal plane. The relative phases of the three diffracted beams may be independently controlled in a way that allows the displacement of the resulting triangular lattice (cf. figure 3). When this procedure is applied to two laser beams that differ in frequency or polarization, it becomes possible to produce two triangular lattices, A and B, which coexist on the same plane and have a tunable relative separation. The result is the original setup introduced in figure 1 in an abstract way, where now A and B are physically implemented by an optical potential. In order to jump from two independent triangular lattices to the setup introduced in figure 1, we need means to introduce the couplings between different lattice sites, that is hoppings from one lattice to another, or within the same lattice. The most direct way to implement this in an optical system is to use the two lattices A and B to trap atoms of the same species but in different internal states, |a and |b respectively. We also need a way to rotate between |a and |b , which can be implemented via Raman transitions. In a situation with all these ingredients two kinds of hoppings are allowed for the atoms: on the one hand an atom in sub-lattice A (B) can tunnel between positions in its own sub-lattice, as contemplated in our model by the hopping parameter Γ A (Γ B ). This parameter can be controlled by increasing or decreasing the intensity I A (I B ) of the laser beam generating each sublattice. On the other hand a Raman-induced change in the internal state of an atom from |a to |b (from |b to |a ) will make the

Realization of two state-dependent triangular sub-lattices
How do we implement in practice the ideas from the previous subsection, and in particular the coupling between lattices? The engineering and control of statedependent lattices is a mature technology [29,32], which nevertheless requires some careful control of the atomic states and decoherence. We will briefly describe how this works for fermionic alkaline atoms and how this integrates with the projected lattices scheme.
As sketched in figure 4, it is possible to find a wavelengthλ l falling between the D1 and D2 lines of an alkaline atom for which polarised light only traps atoms in one of the ground states manifolds. More precisely, the fine-structure energy levels of alkaline atoms are denoted by |L, J, m J , where L is the electron angular momentum quantum number, J the combined orbital and spin momentum and m J is its projection along the quantization axis. When we illuminate with σ + polarised light at a frequency ω = 1 2 (ω D1 + ω D2 ), the ac-Stark shift that it induces on the 0, 1 2 , − 1 2 cancels due to the positive and negative contribution of off-resonant D1 and D2 transitions. The result is that circularly polarized σ ± light can only trap atoms in the 0, 1 2 , ± 1 2 states, respectively. In practice, however, the situation is more subtle because atoms also have some hyperfine structure, induced by the coupling between the electronic and nuclear angular momenta. Let us focus on the fermionic species 6 Li, in line with previous proposals [31]. Out of the hyperfine ground-states |F, m F HF =  If V + and a V − are the intensity distributions that result by illuminating with light in the σ + and σ − polarisations, the states |a and |b will feel the ac-Stark potentials V A = 2 3 V + + 1 3 V − and V B = V − . As shown in figure 5, the combined potentials can look like two displaced triangular lattices by choosing the adequate relative phases between the diffracted beams. Furthermore, the relative depths can be controlled by changing the ratio V + /V − , but this may require a new tuning of the relative phases.
The last ingredient is a coupling between both lattices. This can be done using a Raman laser that couples both internal states, |a and |b . When an atom in state |a is affected by the laser, it will switch state, but since the energy must be conserved, this implies tunneling to a neighboring site in the other lattice. This qualitative description assumes that atoms in the A and B sublattices are confined to the lowest energy band. A very important question is whether this tight-binding approximation is compatible with the Raman laser. This will be discussed in the following section.

Band structure calculations
Our experimental proposal relies on the possibility to combine two different lattices and couple them via Raman assisted tunneling. The feasibility of this procedure has been experimentally demonstrated in a superlattice experiment [33]. It has also been discussed at length in the various works that suggest implementing gauge fields via photon-assisted tunneling, beginning with the seminal paper by Jaksch and Zoller [12], and also in a related work on coupled 1D lattices [34]. Despite this, it is very illustrative to do a quantitative discussion of the lattice parameters involved in this setup, with the aim of clarifying what hopping, Γ A,B and coupling strengths, J i , can be achieved, and what is the limit of weak perturbations that we will rely on later in the manuscript.
Our basic tool in this discussion is the expansion of the field operator in terms of Wannier wavefunctions [29]. We will assume that the two triangular lattices are defined using a similar potential, V triang (x), which for the sake of concreteness we choose This potential gives rise to two single-particle Hamiltonians, with a relative displacement given by v 4 (see figure 1): whose eigenstates are the Bloch waves ψ k (x). The Wannier functions are sums over these Bloch states in a Brillouin zone, w(x) = 1 |B| 1/2 ψ k (x)d 2 k and we assume that they have the same shape for both lattices. Furthermore, we assume that there is a strong confining harmonic potential in the z-direction such that the atoms are restricted to its vibrational ground state with an approximate length scale ( /mω) 1/2 ∼ d.
In the low energy limit that is usual for these experiments, the fermionic field may be approximated by a linear combination of these localized stateŝ In this limit, the motion of particles is described by a tight-binding Hamiltonian, where the intralattice hopping amplitude is defined as the nearest-neighbor term of the Hamiltonian using this expansion, or and similarly for Γ B . If we also include a Raman coupling among lattices, characterized by a Rabi frequency Ω, we may still use the tight-binding approximation provided that this frequency is much smaller than the gap between the lowest energy band and the first excited band of the triangular lattice, |Ω| ∆E. If this is the case, we will have that the interlattice coupling can be estimated as Note that for equation (32) to remain valid, Ω has to be small when compared with V 0 , not with the actual hoppings. This means that the ratio between J and Γ A,B is not fixed, and that one hopping does not have to be small with respect to the other one.
We have computed numerically the Wannier functions for a triangular lattice setup using a discretization of the Brillouin zone with 100 modes and expanding the Bloch wave with up to 625 modes. Using this we have estimated the integrals corresponding to the hoppings, and also to the on-site and nearest-neighbor interactions, given by where C denotes the surface of the unit cell. The results are shown in Fig. 6, where we used the recoil energy E R = 4π 2 2 /2md 2 as unit of energy to ease the comparison, typical values for the lattice spacing (400−600 nm), and interaction strength values for alkaline fermions which range from g/d 3 0.01E R ( 40 K) to g/d 3 0.1E R ( 6 Li) [35]. Note how the ratio J/Ω is at least 100 times larger than Γ A,B /E R for reasonable values of the potential depth. This means that even if we have to impose Ω V 0 , we still can reach a regime in which J Γ A,B . Note also that while the interlattice term U AB decreases rapidly, there is still a window of values where it might be comparable to the influence of hopping, opening the door to experiments with gauge fields and interactions.

Perturbing the hopping parameters
In the previous discussion we presented a feasible setup to obtain a honeycomb tight-binding Hamiltonian with neutral atoms in an optical lattice. In this section we take the scheme one step further so that the hopping parameters can be locally changed to incorporate mass terms and pseudofields to the Dirac Hamiltonian. We propose to optically control the Dirac field physics much in the way mechanical strains have been suggested in graphene sheets [22,23,24]. Two proposals to implement variations are presented, the first of which is best in line with the setup shown in the previous section, the other being a more formal approach to general hopping calculations across potential wells.

Spatial dependence of Raman intensity
Interlattice hoppings are assisted by an external laser beam inducing Raman transitions between the internal states of the atoms, thus creating a transition amplitude which scales with the intensity of the beam. By changing the spatial dependence of this amplitude one can overprint a first-neighbours hopping perturbation i (r). This method is the easiest to implement in this setup, but is very limited and unable to act on intralattice (next-to-nearest-neighbour) transitions. Therefore it is the chosen method to add Abelian vector potentials and non-Abelian coupling terms.

Lattice distortions
A more general analysis can be made by considering spatial distortions of the lattice, that is, either through a relative displacement of one site with respect to its neighbor or by adjusting the width of the individual potential wells. The displacement may be achieved by controlling the relative phases of the conforming beams through the phase modulators P1 and P2 in figure 3. The width of the potential wells can be modified by a change in the height of the confining potential, using a mask which modulates the intensity of the laser beams. We now develop a simple model for calculating the influence of these two modifications in the hopping parameters.
We assume that the ideal (unperturbed) hopping parameters between wells behave as if lattice sites were spatially separated harmonic oscillators in an original reference frame r . The lattice distortion described locally by the transformation r = Ar renders the following single-well Hamiltonian: where the frequency of the site trap is also changed ω →ω(r), due to modulations of the intensity of the confining beam (the kinetic energy terms are not transformed since they are expressed in "real" space, while the potential term is meant to look like a perturbed harmonic oscillator). The localized Wannier function in the well is where B = √ A T A, σ 2 = /mω and N = πσ 2 / |B|. Notice that B is positivedefinite by construction and σ is r-dependent. If we assume no vibrational levels can be excited, the hopping parameter between first-neighbour wells separated by a vector a is simply proportional to the overlap between wavefunctions This will be valid as long as variations are smooth and small, corresponding to the wavefunctions being good approximations to the ground state of the potential well and the restriction to the lowest vibrational level. Therefore the dislocation must be small: A I + C → G = A T A I + (C + C T ), defining C as our generator. Let us call the unperturbed vector connecting both wells a 0 , so that the actual vector becomes a = A −1 a 0 = (I − C)a 0 . The bilinear matrix expands as B = √ G I + (C + C T ) I + 1 2 (C + C T ). Substituting these values yields: This result can be readily interpreted. There are two contributions: one (C + C T ) which modifies the hopping due to the change in distance between wells, and another (ω(r)) which refers to the intensity profile of the confining laser beam. We can therefore tune these parameters, e.g. displacing or rotating one sublattice on top of the other, to allow for variations of the tight-binding hoppings. As explained before, all these changes are used to simulate Abelian and non-Abelian external gauge fields.

Experimental detection
The most direct consequence of the appearance of Abelian and non-Abelian external fields is the distortion of the energy bands, ranging from the movement of the Dirac cones in momentum space to the appearance of a gap in the spectrum. Two methods are proposed here to observe these changes: the first one is the measurement of the momenta after removing the confining potential so as to explore the form of the energy bands; the second one focuses on measuring how these fields affect the dynamics of a single particle (or small group of particles) moving in the lattice.

Time of flight images
It is possible to probe the atomic population within the first Brillouin zone (BZ) of an optical lattice. The procedure consists of the following steps: (i) adiabatically switching off the lattice potentials so that the atoms quasimomenta are converted into real momenta, (ii) then letting the atoms expand freely during a certain time of flight and (iii) finally taking an absorption image of the expanded cloud [36].
In our setup the situation is a bit more complicated. Let us denote by |+, k and |−, k the eigenstates associated with upper and lower energy bands, E + (k) and E − (k), of the tight-binding Hamiltonian (2). In a ground state with a Fermi energy slightly above zero, we would have many atoms in the lower band sector occupying the whole of the first BZ and only few atoms in the upper band, concentrated around Dirac points (See figure 7a). Despite the difference, by means of absorption images one can not distinguish between |+ and |− states as both have components corresponding to the |a ,|b internal states of the atoms.
In order to picture the Dirac cones we need a method that discriminates between energy bands, say, transforming all the |+ states into |a 's and all the |− states into |b 's, while preserving the momentum, k. The adiabatic theorem provides us the way for doing this. Our starting point is the honeycomb lattice Hamiltonian in momentum spacẽ written with the pseudospin structure u k = (a k , b k ), the couplings f (k) = 1+e ik·v1 +e −ik·v2 , and the Pauli ladder operators σ ± . Note how the Hamiltonian   is a composition of commuting terms for each value of the momentum, k. We will adiabatically distort all terms, following the route in figure 8, which consists on first adding a mass term, ∼ mσ z , and then decreasing J down to zero. The protocol maps the two eigenstates of the initial Hamiltonian, |+ and |− , to |a and |b , accurately. This technique would allow for the experimental observation of the effects of an Abelian potential described in section 3, as represented in figure 7b, where the Dirac cone displacement manifests as a deformation of the |a -momentum distribution. In general all other effects, which are summarized in Table 1, could be observed by similar methods.

Spin textures in time-of-flight images
The bipartite nature of the honeycomb lattice allows us to separately probe the atomic population densities for each hyperfine state. We can compare these values at each point of the Brillouin zone to obtain a field S z (k) = a † k a k − b † k b k . Moreover, an adiabatic protocol such as the one described in the previous subsection or in-flight Raman-assisted internal state rotations provide a way of measuring the "rotated" fields S . This observable vector field on the Brillouin zone (S(k) = (S x , S y , S z )) is a powerful tool to analyze some distinct features of the ground state: the presence of a gap in the energy band structure (S z (K ± ) = 0), the characteristic state differences between cones or even the topological nature of the ground state [15]. In figure 9 we show a simulation of such a measurement in two distinct cases: figures 9a, 9b, 9c feature the value of the phase φ = tan −1 (S y /S x ) for different values of the Abelian potential A. The displacement of the cones provoked by A along the Brillouin zone is apparent, as it is the fact that the cones are vortices in the S x , S y vector field. Figure 9d illustrates how a gap or effective mass created by a non-zero value of Γ A can be observed by measuring S z . Once again, the Brillouin zone and the S x , S y vector field are depicted to illustrate the vortex effect of the Dirac cone. We compare the change in σx by an exact simulation of the particle dynamics with the Landau-Zener formula (purple). We denote by t = +(−)∞ the times where the particle is still (already) away from the cone.

Few particle dynamics
The method in the previous subsection gives us access, among other things, to the best known experimental observable in condensed matter physics, which is the density of states. It has been the easiest to measure and therefore has become the default choice in optical lattice simulations of solid state physics. However, the distinctive characteristics of the excitations in a half-filled hexagonal lattice may be worth an extra effort: trying to observe the behaviour of a group of particles with well-defined momentum obeying a complete Dirac equation will reveal some of the distinctive features of the simulated fields.
As an example we suggest using the Klein tunneling effect [37] to probe and measure the energy gaps between the two bands. Let us add a uniform electric field pointing along one direction, say x, described by a linearly growing potential, V x. As explained in Ref. [38], the effect of this potential is to accelerate particles, continuously increasing their momenta in time, until the particle reaches a boundary of the Brillouin zone. We have set an adimensional time scale t = (time) × J/ . At this point two things may happen. If the particle is far away from a Dirac singularity, or the minimum gap between energy bands (the effective mass m) is large compared to the acceleration, mc 2 V, the particle will simply reappear through the opposite side of the Brillouin zone, reversing its velocity and performing the so called Bloch oscillations. However, if the particle hits against the proximities of the K ± points and the gap is small, mc 2 V, the particle will experience a Landau-Zener process and jump to the opposite energy band, maintaining its group velocity.
As shown in figure 10, we have simulated numerically this process for different masses of the Dirac field. As a signature of the jump between bands we simply use the expectation value of k · σ, which is directly correlated to band excitation probability but easier to measure. Note how for the gapless phase, m = 0, there is a perfect jump (Figure 10a), and how the probability is well approximated by the Landau-Zener formula (Figure 10b), which allows us to reverse-engineer the experiment and fit the value of m.
How would this be implemented in an experiment? The idea would be to do the same optical lattice setup with a small number of fermions cooled down to the lowest value of the momenta, k = 0, in a honeycomb lattice that implements the desired Dirac Hamiltonian, with or without effective mass. One would then activate an electric field along the direction w, for a certain time t. After this time one would measure the state of the atoms, or more precisely the expectation value σ . By changing the duration of the field and its intensity, and monitoring the changes in k · σ, one would be able to reconstruct not only the Klein effect but also the whole spin texture of the bands, as discussed in the previous subsection.

Conclusions
In this work we have presented an experimental proposal to simulate Dirac fermions interacting with an effective gauge field. The fermionic component of the model is obtained by trapping atoms in two triangular optical lattices that are connected to form a honeycomb lattice [26]. As explained in the manuscript, the low-energy excitations of such a model may be described using an effective theory that consists of two flavors of non-interacting Dirac fermions. The gauge fields, on the other hand, arise from perturbations of the atom dynamics, such as lattice distortions, short-and long-wavelength modulations of the hopping amplitudes and state-dependent external potentials. Such perturbations are particularly easy to implement using our setup, which consists of two independent lattices.
Our work connects both with recent developments in the field of graphene [22,23,24] and with the field of quantum simulation of synthetic gauge fields [12,13,14,16], but with various advantages. On the implementation side, the optical lattice setup allows for a single-site resolution and a local customization of the potentials which is hard to think of in solid state implementations. Moreover, the use of two atomic species in a bipartite lattice introduces new measurement possibilities, such as the direct observation of the fermionic fields [See section 6.2], or the study of state-dependent Bloch oscillations [Section 6.3]. In comparison with other ultracold atom proposals, while we still rely on the use of assisted tunneling, ours is a static and straightforward setup, where tunneling is implemented by a simple optical field without hopping unitaries, and which nevertheless allows for realistic values of the couplings.
The proposal and the study in this work are also of theoretical interest. The work with optical lattices allows us to make an accurate and rigorous connection between the microscopic theory of trapped atoms and the simulated quantum field theory. While a similar work has been done for strained and curved graphene sheets [23], the optical lattice setup allows us to compute from first principles how the microscopic changes in the optical potential lead to "strain" and hopping distortions. In some cases, as in section 5.1 a simple modulation of the Raman lasers translates into a similar modulation of the hoppings, giving rise to the gauge fields. In other cases, we may concentrate on geometric deformations of the trapping potential and rigorously work out how they affect the hopping matrices, in a theoretically pleasant and flexible way. We foresee that the same tools developed in this manuscript will also help solve an open problem in the graphene world, which is the relation between microscopic deformations of the honeycomb lattice and the appearance of an effective metric and curvature. Instead of introducing an effective spin connection by hand, as it is currently done for carbon layers [23], we expect to be able to derive the same metric directly from deformations of the optical potential.
Finally, it is also remarkable the fact that our models allows the introduction of on-site and nearest-neighbor atomic interactions at no cost and with realistic values [Section 4.3]. At this point the model stops being a single-particle theory and becomes numerically intractable, entering the regime in which quantum simulation provides both new problems and interesting answers.