An Optical-Lattice-Based Quantum Simulator For Relativistic Field Theories and Topological Insulators

We present a proposal for a versatile cold-atom-based quantum simulator of relativistic fermionic theories and topological insulators in arbitrary dimensions. The setup consists of a spin-independent optical lattice that traps a collection of hyperfine states of the same alkaline atom, to which the different degrees of freedom of the field theory to be simulated are then mapped. We show that the combination of bi-chromatic optical lattices with Raman transitions can allow the engineering of a spin-dependent tunneling of the atoms between neighboring lattice sites. These assisted-hopping processes can be employed for the quantum simulation of various interesting models, ranging from non-interacting relativistic fermionic theories to topological insulators. We present a toolbox for the realization of different types of relativistic lattice fermions, which can then be exploited to synthesize the majority of phases in the periodic table of topological insulators.


Introduction
In a seminal paper published in 1982 [1], R. P. Feynman discussed in great detail the problems connected with the numerical simulation of quantum systems. He envisaged a possible solution, the so-called universal quantum simulator, a quantummechanical version of the usual simulators and computers currently exploited in many applications of the "classical" world. If realized, such a device would be able to tackle many-body problems with local interactions by using the quantum properties of nature itself [2]. Interestingly, even without the advent of a fully universal quantum computer, the construction of small dedicated devices, also known as purpose-based quantum simulators, would already be of significant importance for our understanding of quantum physics. The basic idea is to engineer the Hamiltonian of the quantum model of interest in a highly-controllable quantum system, and to retrieve all the desired information with a measurement of its properties. Many research fields would eventually benefit from such devices, such as two-and three-dimensional many-body physics, non-equilibrium dynamics, or lattice gauge theories [3].
Recently, the scientific community is considering ultra-cold atoms as one of the most promising candidates for the realization of a wide variety of dedicated quantum simulations [4,5]. Indeed, these gases are genuine quantum systems where the available experimental techniques offer an impressive degree of control together with high-fidelity measurements, thus combining two fundamental requirements for a quantum simulator. Among the most recent experimental achievements, we would like to mention the observation of Anderson localization in disordered Bose-Einstein condensates (BEC) [6,7], the research on itinerant ferromagnetism with cold fermions [8], or the reconstruction of the equation of state of fermionic matter in extreme conditions, such as in neutron stars [9].
An important drawback in the applicability of cold atoms as quantum simulators is the difficulty of coupling their spatial degrees of freedom to external magnetic fields. This prevents a direct simulation of the quantum Hall physics [10], the controlled observation of whose extraordinary phenomenology would shed new light on quantum many-body theory. One way to overcome this problem is to dress the system with ingenious laser schemes, which mimic the effect of an external magnetic field, and thus allow the neutral atoms to behave as effectively charged particles [11]. This approach led recently to the realization of neutral BECs coupled to external effective magnetic and electric fields [12,13], or even with an effective spin-orbit coupling [14]. More generally, the scientific community has now realized that even in presence of an optical lattice, dressing cold gases with suitable optical and microwave transitions could push the experiments beyond the standard superfluid -Mott insulator transition, and significantly widen the spectrum of the models that are currently being simulated [15]. The possible applications of such optical-lattice-based quantum simulations are numerous and diverse, ranging from the realization of Abelian and non-Abelian static gauge fields [11,16,17,18,19] to that of quantum Hall states [20,21,22,23,24]; from the study of the anomalous quantum Hall effect [25,26] to the quantum spin Hall effect [27,28]; from three-dimensional topological insulators [29], to flat-band physics with a non-trivial topological order [30], or non-Abelian anyons [31]. Recently, a big effort has also been put in designing schemes where the exotic effects associated to relativistic quasiparticles, such as the Klein tunneling and the Zitterbewegung, arise in a controlled table-top experiment [32,33,34,35,36,37,38,39].
In this article, we elaborate on the idea of using a spin-independent bi-chromatic optical lattice dressed with suitable Raman transitions to simulate interesting noninteracting field theories of lattice fermions. We present a concrete proposal to create a three-dimensional optical lattice that traps a multi-species atomic gas, and to tailor arbitrary spin-dependent hopping operators. We have already shown how this setup could break the SU(2) invariance of the hopping rates for spin-1 atoms in spin-independent lattices, and how the simulation of systems subjected to three-body repulsion could benefit from it [40]. Here, we extend this idea further, and show that the same setup allows for the realization of hopping operators which modify the atomic hyperfine state. Combining this trapping scheme with Fermi gases, we show that this platform would open a new route towards the simulation of high-energy physics and topological insulators. This paper is organized as follows: in section 2, we describe qualitatively the idea of using an optical superlattice to realize a general hopping operator for a multi-species cold gas of alkalis. Further analysis and technical details are given in section 3, where we also present some numerical results that support the possibility of controlling a spin-flipping tunneling in this platform. The reader not interested in these technical details may skip this content without prejudicing the comprehension of the following sections. Some final remarks on the proposal are presented in section 4. In section 5, we discuss the possible applications of the described scheme, focusing on relativistic theories and topological insulators, and trying to give a list of the most interesting phenomena which could be explored. Finally, we present our conclusions in section 6.

The Setup and The Idea
We consider the following atomic three-dimensional optical potential where x = (x 1 , x 2 , x 3 ), q = 2π/λ L (λ L is the wavelength of the laser), and where V 0 , ξ > 0 represent the potential amplitudes. The low-energy structure of this potential Figure 2. Sketch of the atomic structure of 40 K: from the electronic structure (L is the electronic angular momentum) to the fine structure (J = L + S, where S is the electronic spin) to the hyperfine structure (F = J + I, where I is the nuclear spin). The latter is drawn in the specific case of an external magnetic field present. The last box shows the optical spin-independent potential which traps equally all the hyperfine levels.
is a cubic array of main minima separated by "secondary" minima located in the middle of each lattice link (see figure 1). We note that additional higher-order minima are also present, but will not play any role in the phenomena discussed in this article. Due to the specific form of the potential in equation (1), the Hamiltonian can be divided into three independent terms, each one depending on one of the three couples of conjugate operators, {x i , p i } i∈1,2,3 . Consequently, the Bloch functions of the nth band with energy E n (p), can be written as ψ n,p (x) = j ψ n,pj (x j ). In order to discuss the effects occurring on the scale of one lattice site, Wannier functions can be introduced for each band w n,R (x) = 1 V e −iR·p ψ n,p (x)dp = w n,R1 (x 1 ) w n,R2 (x 2 ) w n,R3 (x 3 ).
Like Bloch functions, Wannier functions belonging to different bands form an orthonormal basis, and one can thus expand the Hamiltonian in such a basis. Since the Wannier functions are not eigenstates of the Hamiltonian, this expansion leads to a Hubbard model describing the tunneling of atoms between neighboring sites, together with a local on-site interaction coming from the scattering of the cold gas [41]. This setup can be used for the simulation of a lattice field theory, where the field operators are identified with the atomic creation-annihilation operators in the Wannier basis of the lowest energy band (i.e. the states localized in the main minima of the lattice). Conversely, higher energy bands provide auxiliary levels that shall be used as a resource to tailor the tunneling processes. The main result of this article is the claim that a complicated though not unfeasible combination of current technologies leads us to the realization of the following Hamiltonian Here, we are considering a multi-species fermionic scenario with many hyperfine levels of the same atom: c † rτ (c rτ ) creates (annihilates) a fermion with hyperfine spin τ localized in the main minima of the superlattice at r = m 1 a 1 + m 2 a 2 + m 3 a 3 , where m j ∈ {1...L j }, L j stands for the number of lattice sites along the x j axis, and a j is the lattice spacing in the j-th direction. The parameter t ν stands for the strength of the laser-assisted tunneling in theν direction, with with ν ∈ {a 1 , a 2 , a 3 }, which shall be described below. The operators U ν describe the tunneling from r to r + ν, and are . Sketch of a laser-assisted tunneling induced in the presence of a superlattice. Two physical hyperfine states belonging to the F = 9/2 manifold are connected via Raman couplings with the intermediate level of an auxiliary state belonging to the F = 7/2 manifold. If the coupling is detuned enough, the F = 7/2 level can be adiabatically eliminated: no population is left there and an effective coupling is engineered between neighboring sites. Left: scheme for a spin-preserving (i.e. diagonal) hopping. Right: scheme for a spin-flipping hopping.
a common feature in lattice gauge theories. We have also included an on-site Raman term Λ, of strength Ω, that induces a certain transition between the hyperfine states. Note that we use Gaussian units and = 1. The claimed possibility of engineering a wide range of hopping operators U ν , together with the state-of the-art control of the atomic interaction, makes already our system a versatile quantum simulator of lattice field theories. In this manuscript, we focus on non-interacting theories, which can be realized either with dilute systems, or by employing Feshbach resonances to lower the interaction strength (see e.g. [6,7]). We stress that interesting phenomena can also be observed in non-interacting gases when additional ingredients are introduced in their dynamics, such as disorder or the assisted-hopping processes discussed in this article.
Let us note that the control of the homogeneous tunneling for a single-species atomic gas is straightforward, and would not even require the superlattice (ξ = 0) [42]. Moving to a many-species case, one runs into the problem that a general hopping operator also entails terms flipping the atomic hyperfine spin (simply referred as spin in the following), which are not easily engineered. Here, we propose to realize such couplings by combining Raman transfers and a bi-chromatic superlattice (ξ = 0 in equation (1)). The proposal can be applied to all the alkalis notwithstanding their bosonic or fermionic nature. In the following, however, we shall focus in the fermionic scenario, which is best explained with the following practical example.
Let us consider an ultra-cold cloud of non-interacting 40 K atoms in the presence of a magnetic field of intensity B. Such field lifts the spin degeneracy within the two atomic hyperfine manifolds of the ground state, F = 9/2 and F = 7/2, according to the following relations (see also figure 2): where m F is the projection of the hyperfine spin along the quantization axis defined by the magnetic field, µ B is the Bohr magneton, g F is the hyperfine Landé Factor, and ∆ HF stands for the hyperfine splitting. These hyperfine levels are all trapped into the same spin-independent optical potential (1). Depending on the lattice theory we want to simulate, we select a subset of these hyperfine levels described theoretically by creation-annihilation operators in the lattice sites. We then identify such fields with the components of the lattice field theory to be simulated. This leads us to divide the hyperfine levels into two subsets: the subset of "physically meaningful" states, which belong to the hyperfine manifold F = 9/2, and the usually larger subset of auxiliary levels that shall be used to assist the tunneling and create the desired hopping operator.
Regarding the hopping operator in equation (2), we address each of its matrix elements [U ν ] τ τ separately. Given a matrix element (i.e. once we have identified the initial and final hyperfine levels to be connected by the assisted tunneling), we choose an auxiliary level belonging to the hyperfine manifold F = 7/2 trapped in the middle of the link. These levels provide intermediate "bus" states that shall be used as a resource to assist the tunneling as follows. The couplings between the atoms in the main sites, R 1 , and the "bus" states, R 2 , are realized via optical twophoton Raman processes transferring a net momentum q t . They have a mathematical expression proportional to the overlap integral of the initial and final Wannier functions: w * n2,R2 (x)e iqt·x w n1,R1 (x)dx. This integral is not zero because of the term e iqt·x , which is of course relevant only if 2π/|q t | is of the order of the lattice spacing. Since this regime cannot be achieved with microwave transitions, one is motivated to employ two-photon Raman transitions. Interestingly enough, it is possible to eliminate adiabatically the intermediate level and obtain an effective four-photon coupling between neighboring sites (see figure 3). We stress that different matrix elements can be engineered at the same time thanks to the magnetic-field splitting of the hyperfine levels (3): the involved atomic transitions become non-degenerate and can be individually addressed with different lasers. Furthermore, the use of coherent laser light for the Raman transitions entails the additional advantage of being able to deal with complex phases, and thus to realize complex gauge structures at will. The realization of the non-diagonal matrix elements requires the lattice to be slightly staggered, a technique discussed also in reference [17]. Summarizing, this proposal tries to exploit a hierarchy of energies characterizing atomic gases in optical lattices in order to assist the tunneling between neighboring sites with controlled adiabatic eliminations (see table 1).
The on-site spin-flipping Λ in equation (2) can be performed with standard technology based on microwave transitions, or Raman transitions carrying negligible momentum. Furthermore, these terms can also be exploited to correct spurious onsite couplings which may be induced by the laser scheme. Unfortunately, we note that there is no selection rule relying on the polarization properties of the light and the hyperfine moment of the atoms, which can be used to realize the different on-site and nearest-neighbors spin-flipping processes. The superimposed magnetic field cannot be aligned at the same time with the propagation vector of all the three lasers, aligned along the three Cartesian axes, which would be the case in which circularly polarized light could be exploited to induce controlled transitions. Conversely, we shall rely on the different Zeeman-shifted energies to selectively address the different couplings between the internal states.

Realization of Spin-Dependent Hopping Operators
In this technical section, we theoretically and numerically confirm the qualitative scheme presented above. We study two simple but important cases: the realization of diagonal and non-diagonal hopping operators for a two-species atomic gas. These can be considered as the main building blocks needed to realize any tunneling operator even in situations with more than two atomic species.

Coupling Between Different Hyperfine Manifolds
The most fundamental ingredient of this proposal is the possibility of using Raman processes to induce controlled atomic transitions between different hyperfine states of the electronic ground state L = 0 (L is the total electronic angular momentum). These transitions are realized with two lasers via adiabatic elimination of the electronically excited manifold L = 1. In the following, we address the atomic levels as |L, α, k , with α labeling the hyperfine degrees of freedom (see also figure 2 for some insights on the internal structure of 40 K), and k the quantum numbers of the center-of-mass wavefunction (in our case, the Wannier functions of the optical potential). As discussed in [40], the induced Raman coupling between the state |0, α, k and |0, α , k can be written as follows: This expression clearly factorizes the following contributions: • the time-dependence of the effective coupling and its effective frequency, which is the difference between the frequencies of the two lasers ω = ω 1 − ω 2 ; • the dependence on the center-of-mass degrees of freedom, S k k = k |e −i(p2−p1)·x |k , where p 1 and p 2 are the momenta of the two lasers; • the dependence on the initial and final internal states and on the polarization properties of light, Ω α α , which is a function of the dipole matrix elements between the initial (final) state and the excited levels.
Next, we specify (4) to the superlattice setup of section 2, i.e. we will consider Raman transitions in presence of lattices characterized by a Wannier function trapped in the middle of each link.

Developing an Effective "6-Level Model"
Let us address the simulation a theory characterized by two-component fields. Following the discussion of section 2, we take two states of the F = 9/2 manifold of 40 K, for instance |9/2; m F = 7/2 and |9/2; m F = 9/2 , and map them into the theory to be simulated. Here and in the following subsections, we discuss the laserassisted hopping in the diagonal case (m F preserved while hopping) and non-diagonal case (m F flipped while hopping). For the diagonal case, we develop the "6-level model" depicted in figure 4. We consider one physically meaningful state, say |F = 9/2, m F = 9/2 , and one auxiliary state, say |F = 7/2, m F = 7/2 . Moreover, we consider different Wannier states for each of them, two localized in main sites (k = 1 and 3) and one in the intermediate link (k = 2). The model includes the effects of undesired couplings and additional levels, and its limitations, together with the approximations on which it relies, will be discussed at the end of the paragraph. We can identify the states with the short notation |F, k rather than with the longer previous one |0 α k . Below, we give an analytical estimate of the population transfer rate, whereas in the next subsections we present the numerical time-evolution for physically interesting cases.
The model is parametrized by six relevant couplings between the different Wannier functions S k k (see figure 4), whose properties are listed below. We exploit the existence of theorems which assure the possibility, in our case, of considering three real and exponentially localized Wannier functions w j (x), j ∈ {1, 2, 3} [43]. We write the parameters S k k factorizing out the space dependence of the coupling e iqt·xj , where x j is the position of the point around which the Wannier function w j (x) is localized, The parameters S 1,1 and S 2,2 describe two on-site couplings, whereas S 1,2 is the coupling between a main site and an intermediately trapped state (see figure 4). The last relation states that couplings between neighboring main sites are negligible. The relation between the other four overlap factors depends on the particular experimental situation. In this case, we are interested in the simplest scenario where a single Raman transition induces all these couplings, which leads us to In order to make this scheme simpler, we assume q t = 2q L , and thus e iqt·(x2−x1) = 1. As we will argue below, transferring a momentum which does not fulfill this requirement is not a problem since the resulting phase can be gauged away. The phase 2q t · x 1 can also be put to zero for the moment, since its role only becomes important when one needs to give a phase to different matrix elements. In the following, we will also consider situations where the coupling between the lattice sites 2 and 3 could be induced by lasers propagating in the opposite direction, where S 1,2 = S 2,1 = S 3,2 = S 2,3 . Taking these considerations into account, the Hamiltonian reads as follows (see figure 4 for the definitions of δ, ∆ and ω): Once we apply the unitary transformation the three levels |9/2, k become degenerate. In case the three inequalities |S i,j Ω|/(δ − d) 1 are fulfilled, it is possible to use second-order perturbation theory in order to develop an effective Hamiltonian describing the dynamics within the sub-manifold we are interested in, namely Remarkably enough, this Hamiltonian leads to the desired transfer rate of population from level |9/2, 1 to |9/2, 3 , and viceversa. The main contribution is the direct coupling A second contribution, which in our system will prove to be not-negligible, comes from a sort of "adiabatic elimination" of the level |9/2, 2 , namely Accordingly, we have derived the desired effective Hamiltonian where the Raman lasers assist the hopping of the physically meaningful F = 9/2 levels, after the auxiliary F = 7/2 bus states have been adiabatically eliminated. In the following sections, we shall address the range of validity of the approximations leading to this Hamiltonian, and compare it with the exact numerical investigation of the initial Hamiltonian (7). We want to stress here that even if the integrals in the definition (5) of the S kk can be complex numbers, this does not have any physical influence on this proposal. Indeed, even if the effective coupling between neighboring main sites −J was complex, its spatially uniform phase can be gauged away with a space-dependent unitary transformation (even in the case of periodic boundary conditions). Conversely, the non-uniform phase coming from the e iqt·x k factor, which arises when q t is not parallel to the direction of the hopping it assists, cannot be gauged away even in presence of open boundary conditions. Such a phase, which is not related to the fact that the integrals in (5) are complex, can be used to simulate an external uniform magnetic field [16,17]. Finally, we underline that in our setup, where the tunneling along each axis is induced by lasers propagating parallel to the axis itself, both complex phases can be gauged away. In order to simulate a magnetic field, therefore, one should move slightly away from this configuration and engineer a Raman coupling whose effective transmitted momentum does not run parallel to the links of the lattice. We will not consider this situation in this article because the models of interest in section 5 do not require such space-dependent phase.

Range of validity of the "6-Level Model"
The presented "6-level model" strongly relies on two approximations: (i) considering the bands of the lattice as being flat; (ii) neglecting delocalized higher-energy free states.
If these approximations are not justified for a given experimental configuration, spurious population transfers to next-neighboring sites would arise.
The approximation (i) is required to fulfill the core idea of the proposal, namely the adiabatic elimination of the intermediate level. This is demonstrated with a model which considers only a subset of the Hilbert space spanned by the real eigenstates of the Hamiltonian (Bloch functions), considering just three of their linear combinations (the Wannier functions w k=1 (x), w 2 (x) and w 3 (x)). This is equivalent to approximating the dispersion laws of the band as being flat, neglecting thus possible curvature effects, and is legitimated as long as the width of the band is much smaller than the detuning of the transition δ − d. In case the degeneracy of the Bloch functions cannot be assumed, all the Bloch functions should be considered in order to quantitatively estimate the spurious effects cited above. In general, this issue sets a trade-off for the relative depth ξ of the secondary lattice in (1): on one hand, a shallow lattice (ξ < 1) is desirable because the Wannier function of the intermediate minimum w k=2 (x) is not strongly localised and laser-induced transitions are favored (|S 1,2 | ∼ |S 1,1 |). On the other hand, the more the wavefunction is delocalized, the more the band bends, eventually  2) connects the |F = 9/2, m F = 9/2 (|F = 9/2, m F = 7/2 ) states to their auxiliary state. Detuning allows independent control of the hopping rates. Right: exact timeevolution of the "6-levels model" (7), showing the coherent population transfer between sites 1 and 3 of the spin state |F = 9/2, m F = 9/2 . The parameters used are listed in table 2. The inset shows the maximal populations of the six considered levels labelled |F, k as in (7) and shows that only a small fraction of the population is lost in auxiliary levels.
becoming parabolic at k = 0 with a bandwidth comparable to the detuning. In our numerical simulations we consider ξ = 1, which is a reasonable middle-way.
Regarding the issue (ii), higher-energy bands could become important in the presence of intense Raman transitions Ω and large detunings δ − d, which couple them to the lowest-band states. The presented analytical and numerical studies do not take into account these effects since they consider only three Wannier functions and effectively only two bands, even though including bands with localized Wannier functions would just imply a renormalization of the numerical coefficients S k,k . A different problem is the case of high-energy strongly parabolic bands, whose Wannier functions are not strongly localized. The effect of such states is not considered by our model, which is that of spreading population among many next-and furtherneighboring main sites. From an experimental point of view, we expect a trade-off to arise between a large detuning regime, allowing powerful lasers and strong effective couplings with noisy spurious population transfers, and a small detuning one, with clean but small couplings.

Diagonal Hopping Operator
We now explicitly study the possibility of realising a diagonal tunneling operator. We numerically simulate the Hamiltonian (7) with a simple Runge-Kutta algorithm. We did not include in the simulation hyperfine states different from |F = 9/2, m F = 9/2 Table 2. Parameters used for the numerical simulation of the diagonal hopping in subsection 3.4. We list the numerical values of all the main parameters characterizing the atomic transitions and the Raman couplings. The first Raman coupling induces the hopping of the F = 9/2, m F = 9/2 whereas the second one addresses the F = 9/2, m F = 7/2 (such states were however not considered in the simulation).
Level: |F, m F , k Energy Parameters  and |F = 7/2, m F = 7/2 because they are strongly detuned from those we are considering. However, for completeness, we include the presence of a second Raman coupling which would be needed to induce the hopping of |F = 9/2; m F = 7/2 and check that it is unimportant. We show in figure 5 the numerical results. The realistic parameters used in this simulation are listed in table 2. The population is coherently transferred between two neighboring levels and only a negligible fraction is lost in auxiliary states. Regarding the validity of the "6-levels model", for the lattice considered here, the bandwidths of the two bands are respectively 0.2 kHz and 9.1 kHz, which should be compared with the considered detuning of 300 kHz. In these and the following simulations, the employed numerical values have only an illustrative purpose and other regimes could be considered.
Therefore, these results confirm the plausibility of our scheme to induce a laserassisted tunneling between the atoms sitting in the main minima of the optical lattice. To make the simulation toolbox reacher, we now address the possibility to control a spin-dependent hopping process.

Non-Diagonal Hopping Operator
In order to study the realisation of the non-diagonal hopping operator, we consider an enlarged 12-levels model, which is a generalization of the previous one taking into account more hyperfine states. We want now to transfer population between the manifolds |F = 9/2, m F = 9/2 and |F = 9/2, m F = 7/2 and consider as auxiliary states |F = 7/2, m F = 7/2 and |F = 7/2, m F = 5/2 .
A big issue which must be solved to engineer such a hopping is the arousal of undesired spin-flipping terms induced by the laser. In this paper we consider the possibility of staggering the lattice with an additional optical field, in order to lift the degeneracy between the different sites of the optical lattice, in the same fashion of [17]. Such a staggering can be done also in three dimensions since the cubic lattice is bipartite, and we consider staggering values of 10 − 15 kHz. Figure 6 sketches the experimental scheme we have in mind and shows the exact time evolution of the population transfer between the two levels of F = 9/2 in two neighboring sites. Interestingly enough, we show a flip of the Zeeman spin during the tunneling process, and thus obtain the promised spin-dependent hopping operator. The parameters of the simulation can be found in table 3.

From a Spin-Dependent Hopping Operator to a Quantum Simulator
In the previous section, we discussed how the superlattice geometry could be used to create non-trivial hopping operators on each link. Here we want to assemble these ingredients and discuss how to use them to engineer interesting quantum simulator in arbitrary dimensions, considering the advantages and disadvantages of the proposal.
First of all, we stress that the lasers needed to engineer the hopping along one direction must transfer momentum along that same direction. Therefore, just by controlling the beam propagation directions, we can tailor different tunneling operators along each axis. This is an important feature which will be largely exploited in the proposals listed in section 5. This kind of directionality selection rule is also responsible for avoiding the population of higher-order minima which do not lie on the edges of the unit lattice cell. Since they are not connected to the main minima by a line parallel to a cartesian axis, we do not consider any momentum transfer along such direction, and thus the formal orthogonality of the Wannier functions localized in those minima is never lifted.
Due to the very general formulation of the superlattice potential, the setup is well-suited also to work in two and one dimensions. Moreover in 1D it is possible to align the magnetic field splitting the hyperfine sublevels with the optical lattice. In this case, it is possible to use specific light polarization to selectively couple different atomic levels. As a short-term goal, it would be very interesting to understand what is the most interesting physics which could be simulated in a one-dimensional system, where the presence of more symmetries could lower the experimental intricacies. We partially address this question in section 5, where we argue that several onedimensional topological phases could be realized. In a three-(two-)dimensional case, with the magnetic field aligned with the diagonal of the cube (square), the polarization of the light cannot be exploited, only energy-based selection rules are reliable.
Unfortunately, energy-based selection rules are not enough to prevent, in the case of spin-flipping operators, spurious on-site spin-flipping couplings. We proposed to solve this issue by staggering the optical lattice, i.e. lifting the degeneracy of the lattice sites of some tens of kHz (see also the discussion in reference [17]). One should also mention that the diagonal hopping operators can still be engineered in presence of such staggering, with the only additional issue of using two Raman couplings (as in the non-diagonal case) to match the energy difference between sites.
Each of the matrix elements of the hopping operators is realized via an effective four-photon process. This means that the spin of the atom can be flipped of at most |∆m F | = 4: a careful analysis is needed in case one is interested in simulating a theory with more than 4 fields, because some hopping matrix elements might become avoided.
Finally, the description given in this proposal is essentially at the single-particle level where no many-body effects have been used. As a consequence, the proposal works both for bosons and fermions, which is a valuable result. The proposals of section 5 assume the possibility of switching off atomic interactions with Feshbach resonances. Nonetheless, it would be of the utmost interest to study how interactions could fit into this picture, and to explore the rich variety of models that could emerge. We leave this problem for future work.
We now conclude the part of the article devoted to the description of the experimental setup. We believe to have provided the relevant results supporting the initial claim that it was indeed possible to realize a system whose low-energy structure is described by Hamiltonian (2).

Applications of the Quantum Simulator
In this section, we would like to discuss the use of the described setup as a quantum simulator. We will address a range of lattice field theories for relativistic fermions [44], and explore exotic phases of matter known as topological insulators [45]. In general, the task of a purpose-based quantum simulator is to realize a system described by an effective Hamiltonian H eff that reproduces faithfully the properties of the model to be simulated. In our case, this model corresponds to relativistic lattice fermions H rel , or topological insulators H top . The resource to accomplish such a task is the microscopic control over the superlattice setup, which we have argued previously to be described by the Hamiltonian (2), rewritten here for reading convenience, The main objective in the following subsections is to control and manipulate: (i) the optical lattice dimension D; (ii) the tunneling strengths t ν ; (iii) the spin-dependent hopping operators U ν ; (iv) the on-site Raman transitions Ω, Λ; such that the Hamiltonian of equation (2) simulates the desired physics, namely From a condensed-matter perspective, exotic phases are frequently associated to strongly-correlated regimes and many-body interactions. There are however distinguished exceptions to this paradigm, such as graphene [46] and topological insulators [45], where quadratic fermionic Hamiltonians contain a wealth of non-trivial phenomena. In the case of graphene, a two-dimensional layer of graphite, the lowenergy carriers can be described by emerging relativistic fermions without mass. On the other hand, topological insulators are exotic holographic phases with an insulating bulk, and a peculiar boundary that hosts robust conducting modes protected by topology arguments. In both cases, the transport properties of the material differ significantly from the standard solid-state theory. In the subsections below, we show how the Hamiltonian (2) serves as a versatile simulator of these two interesting phases of matter.

The Zoo of Relativistic Lattice Fermions
The properties of a relativistic spin-1/2 fermion with mass m are described by the famous Dirac Hamiltonian [47] H = dr Ψ(r) † H DI Ψ(r), where α ν , β are the so-called Dirac matrices fulfilling a Clifford algebra, {α ν , α µ } = 2δ νµ , {α ν , β} = 0, and c stands for the speed of light. Here, Ψ(r) is the N D -component fermionic field operator, where N D = 2 for one and two spatial dimensions, and N D = 4 for three spatial dimensions. Our objective now is to construct an effective Hamiltonian starting from equation (2) that closely resembles the relativistic field theory in equation (11). The underlying setup consists of a gas of ultracold 40 K atoms, which is a non-relativistic system; nonetheless we can design it as a quantum simulator of relativistic particles by exploiting the quantum statistics and a peculiar engineerable Fermi surface.
Naive Massless or Massive Dirac Fermions. The idea is to engineer translationally invariant hopping operators, U ν = e iφν Aν , according to the SU(N D ) group, where N D has been defined above. For the particular choices specified in table 4, one finds that the Hamiltonian in equation (2) in momentum space becomes where Ψ k is a multicomponent fermi operator that contains the different N D hyperfine levels involved in the simulation, and k is defined within the first Brillouin zone BZ. One readily observes that there are certain regimes, the so-called π-flux phases φ ν = π/2, where the energy spectrum develops N D = 2 D degeneracy points K d where the energy bands touch (K d ) = 0. Around these points K d = (d x π, d y π, d z π), where d ν ∈ {0, 1} is a binary variable, the low-energy excitations of the 40 K Fermi gas are described by the effective Hamiltonian where p d = k − K d represents the momentum around the degeneracy points, (α d ) ν = (−1) dν α ν are the Dirac matrices listed in table 4, and c = 2t x = 2t y = 2t z is the Fermi velocity that plays the role of an effective speed of light. Therefore, the Fermi surface of the half-filled gas consists of a set of isolated points, the so-called Dirac points, and the low-energy excitations around those points behave according to the Hamiltonian of massless Dirac fermions in equation (13). As occurs with graphene [46], Table 4. Quantum simulator of naive Dirac fermions. Each translationallyinvariant hopping operator Uν is expressed in terms of Pauli matrices {σx, σy, σz}, and a set of dimensionless fluxes {φ 1 , φ 2 , φ 3 }.
We also list a particular representation of the Dirac matrices αν , β for different spatial dimensions d, together with the important symmetries of the Hamiltonian that depend on Γ.
Let us note that we obtain an even number of relativistic-fermion species, each located around a different Dirac point (i.e. N 1 = 2 for one dimension, N 2 = 4 for two dimensions, and N 3 = 8 for three dimensions). This doubling of fermionic species is a well-known phenomena in lattice gauge theories [44], where the fermions in equation (13) would correspond to the so-called naive Dirac fermions [48]. As predicted by the Nielsen-Ninomiya theorem [49], this doubling cannot be avoided without breaking an underlying symmetry: • for D odd {H d DI , Γ 1 } = 0 is an involution known as chiral symmetry; is an antiunitary symmetry known as time-reversal (see table 4).
According to these results, we have a versatile quantum simulator of massless Dirac fermions in any spatial dimension D = 1, 2, 3. In D = 1, 2 they coincide with the famous Weyl fermions, whereas in D = 3 they contain a couple of Weyl fermions with opposed helicities. Note that this scheme can also be extended so as to simulate exotic Weyl fermions of any arbitrary spin s [50]. Besides, our quantum simulator also allows us to make these fermions massive, thus reaching the desired Hamiltonian in equation (11). The idea is to control the on-site Raman transitions such that Λ = β listed in table 4. In such case, the Rabi frequency plays the role of the mass mc 2 = 2Ω, and the effective Hamiltonian in equation (13) becomes Therefore, this quantum simulator can explore both the non-relativistic and the ultrarelativistic limits of the theory.
Wilson and Kaplan Fermions. From a lattice gauge theory perspective, the additional fermions around K d = 0 are spurious doublers that modify the physics at long wavelengths. A partial solution is to give the doublers a very large mass m K d c 2 , so that they effectively decouple from the low-energy physics of the Dirac fermion at K d = 0, namely m K d m K0 . We must find a way to engineer a momentumdependent mass that differs from the global on-site Raman mass discussed above. By combining the laser-assisted tunneling listed in table 4, with the additional terms U ν = ie iϕν β [29], the momentum-space Hamiltonian H in equation (12) becomes wheret ν are the additional laser-assisted tunneling strengths. Once more, for the π-flux phases ϕ ν = π/2, the effective Hamiltonian in equation (11) is modified into where m ν c 2 = 2t ν . Let us emphasize that our quantum simulator provides a complete control over the different masses, since m depends on the on-site Raman transition strengths, whereas m ν depends on the assisted-hopping strength, and thus on the laser power. In particular, when these parameters fulfill ν m ν = m (i.e. m x = m for D = 1, m x +m y = m for D = 2, and m x +m y +m z = m for D = 3), we obtain a single massless Dirac fermion at K d = 0, whereas the remaining doublers have been boosted to much higher energies. We thus achieve a quantum simulator of the so-called Wilson fermions of any spatial dimension D = 1, 2, 3 [51]. We note that this decoupling between a single massless Dirac fermion and its doublers is not in conflict with the Nielsen-Ninomiya theorem since the introduced mass terms explicitly break the aforementioned symmetries. This is particularly important in odd dimensions, where the theory does not preserve chiral symmetry, a fundamental concept in the standard model classifying right/left-handed particles Γ 1 Ψ = ±Ψ. To preserve such symmetries, the concept of Kaplan fermions arises [52], namely massless Dirac fermions bound to a lower-dimensional domain wall located at r * ⊥ where the Wilson mass gets inverted m K d = −|m| + 2|m|θ(r ⊥ − r * ⊥ ). Since we have a complete experimental access to the parameters of the Wilson mass, it is also possible to tune them such that ν m ν > m, and thus the mass m K d < 0 gets inverted, and one gets a lower-dimensional massless fermion bound to the region where this mass inversion takes place.
Let us close this subsection by underlying the versatility of our setup as a quantum simulator of a diverse set of relativistic lattice fermions. Not only can we implement massless Dirac fermions of any dimensionality, thus exploring their connection to Weyl fermions, but we can also control their mass. This leads us to the concept of massive Dirac fermions, and the notorious Wilson and Kaplan fermions dealing with the fermion doubling problem. Interestingly enough, the physics behind these high-energy particles is intimately related to the materials known as topological insulators [45], which are the subject of the following subsection. In fact, it is always possible to find a Kaplan-fermion representative within each class of topological insulators [53].

A Toolbox for Topological Insulators
Topological insulators correspond to fermionic gapped phases of matter that are insulating in the bulk but allow a robust transport along the boundaries [45]. This Table 5. Periodic table of topological insulators. The underlying quadratic Hamiltonians H = αβ Ψ † α H αβ Ψ β + H.c., where α, β represent the lattice sites and the internal states of the fermion, can be classified according to the fundamental symmetries of time-reversal T , charge conjugation C, and the combination of both S = T C. The values T = 0, C = 0, S = 0, are assigned to Hamiltonians that break the symmetry, whereas T = ±1, C = ±1, S = 1 correspond to symmetry-preserving Hamiltonians, where T 2 = ±1, C 2 = ±1, S 2 = +1. There are six possible combinations for non-interacting fermionic Hamiltonians (and another four for pairing fermionic Hamiltonians), which lead to the classes listed in the first column. For each dimension d, there are three possible topological insulators among all these classes, and they are classified according to the integer Z or binary Z 2 nature of a certain topological invariant. In the column labelled QS, we list the particular instances that can be simulated with our super-lattice based quantum simulator.

Class
Name robustness is due to the existence of gapless edge excitations which are protected against disorder by topological arguments (i.e. they avoid Anderson localization [54], and thus transport charge even in the presence of strong disorder). The paradigmatic example of a topological insulator is the integer quantum Hall effect (IQHE) [45], a two-dimensional electron gas subjected to a strong magnetic field that displays a robust quantization of the transverse conductivity σ xy = n · e 2 /h, where n ∈ Z is related to the topological invariant known as the Chern number [56]. In this case, chiral electrons bound to the one-dimensional edges of the sample avoid back-scattering processes, and are thus immune to disorder [57]. Remarkably enough, the IQHE is only one instance of a large list of topological insulating phases. Each class can be characterized by a set of discrete fundamental symmetries [58], and a certain topological invariant, see table 5. Note that we have excluded the topological superfluids from this table, since their quantum simulation would require a pairing mechanism, and thus goes beyond the scope of this work [59,60]. We now discuss how our quantum simulator can reproduce the properties of many of these fascinating phases of matter following two possible strategies.
Bottom-Up Approach. In this case, one designs the ultracold-atom Hamiltonian so that it simulates a particular model belonging to the desired class of topological insulators. Therefore, a different experiment would be required for each class-oriented simulator. Two representative, yet reasonably simple examples are the Su-Schrieffer-Hegger model of polyacetilene [61], which is related to the D = 1 BDI topological insulator, or the π−flux phase of the fermionic Creutz ladder [62], which is related to the D = 1 AIII topological insulator. The former can be simulated by using a one-component Fermi gas in a one-dimensional dimerized optical superlattice, thus obtaining where δ quantifies the different tunneling strength between superlattice sites [35]. On the other hand, the Creutz ladder is described by where K, M are tunneling strengths, and θ is a magnetic flux piercing the ladder. This requires two Zeeman sublevels to be assigned to the fermion species a n , b n , and a one-dimensional laser-assisted tunneling U a1 = diag{e −iθ , e iθ },Ũ a1 = iσ x , together with the Raman on-site operator of strength M [63]. It is possible to continue this approach, proceeding thus to higher dimensions and different topological classes. Prominent examples would be the honeycomb timereversal breaking Haldane model [64,65] for the D = 2 topological insulator in class A, or the time-reversal Kane-Mele model in the honeycomb lattice in class AII [66], or other optical-lattice geometries [27,67,68]. Rather than following this route, we shall explore a different approach that is better suited to the superlattice-based simulator introduced above. Indeed, we shall argue that this quantum simulator allows the reproduction of most the topological phases in table 5.
Dimensional-Reduction Approach. In this case, the starting point is the quantum simulator of D-dimensional Kaplan fermions in equation (16). Depending on the particular choice of Dirac matrices, the inverted-mass regime shall correspond to a different class of topological insulators. Besides, in some situations, a dimensional reduction [53] that amounts to the increase of the optical-lattice depth in one direction, connect us to a different lower-dimensional class. We rewrite the full Hamiltonian where the Wilson mass is m K d = m − ν (−1) dν m ν , and where the Dirac matrices α d ν , β shall be selected so that the T , C, S symmetries are explicitly broken or preserved. This translationally-invariant Hamiltonian preserves these symmetries when the following conditions are met C : S :  (19) that directly lead to several classes of topological insulators. Each class, characterized by the discrete symmetries T , C, S, where T 2 = Θ T Θ * T = ±1, and C 2 = Θ C Θ * C = ±1. Besides, each class has a Wilson-fermion representative with a particular choice of the Clifford algebra αν , β that depends on the dimensionality D and the corresponding symmetries. We also highlight the topological insulators that can be obtained by dimensional reduction from a parent Hamiltonian, such as AII, D = 3 → AII, D = 2, or A, D = 2 → AIII, D = 1. We also list the unitary matrices Θ T , Θ C involved in the definition of the discrete symmetries.
where Θ T , Θ C are some unitary matrices. In table 6, we list the symmetry properties of different Kaplan-fermion Hamiltonians. It is important to note that these symmetries might correspond to the exact symmetries in nature (e.g. when considering the hyperfine levels {|F, m F , |F, −m F }, the time-reversal symmetry given by θ T = iσ y exactly correspond to time-reversal symmetry in nature). Conversely, these symmetries might otherwise correspond to the algebraic properties of the effective Hamiltonian. Let us emphasize, however, that as far as the disorder respects such symmetries, the robustness of the edge excitations is guaranteed. It would be of the greatest interest to design disorder breaking or preserving such symmetries, generalizing the studies on Anderson localization [6,7]. In table 6, we have listed the different topological insulators that can be simulated with our scheme. As shown in [29] for the particular case of three-dimensional AII insulators, the laser parameters can be controlled so that an odd number of Wilson masses are inverted. This mass-inversion occurs through a gap-closing point, and thus a quantum phase transition between a normal band insulator, and a topological one occurs. This new phase is characterized by an odd number of massless fermionic excitations (i.e. massless Dirac fermions) bound to the boundaries of the system, and protected by a topological invariant. In the three-dimensional case, this corresponds to an axion term that modifies the response of the system according to the so-called axion electrodynamics [69]. Remarkably, table 6 contains all the relevant information to explore the exotic properties of different topological insulators in a superlattice based experiment with ultracold atoms.

Conclusions
In this article, we have presented a concrete proposal for the realization of laserassisted tunneling in a spin-independent optical lattice trapping a multi-spin atomic gas. Remarkably enough, it is possible to tailor a wide range of spin-flipping hopping operators, which opens an interesting route to push the experiments beyond the standard superfluid -Mott insulator transition. The scheme we have presented combines bi-chromatic lattices and Raman transfers, to adiabatically eliminate the states trapped in the middle of each lattice link. These states act as simple spectators that allow us to assist the tunneling of atoms between the main minima of the optical lattice. This mechanism is clearly supported by our numerical simulations of the timeevolution of the atomic population between the different optical-lattice sites. Even if we focus on fermionic 40 K, we stress that the ingredients of this proposal do not rely on the atomic statistics, and could be thus used for all the alkalis.
We believe that such a device could have important applications in the quantum simulation of non-interacting lattice field theories, which are characterized, in their discrete version, by on-site and nearest-neighbor hopping Hamiltonians. Once the fields of the theory to be simulated are mapped into the atomic hyperfine states, the desired operators correspond to population transfers between such levels. The former can be realized by standard microwaves, whereas the latter might be tailored with the laser-assisted schemes described here.
Even though interactions are at the heart of a plethora of interesting effects, non-interacting fermionic theories already encompass a number of phenomena whose experimental realization would be of the greatest interest. In the second part of the article, we analyzed interesting physical models which could benefit from our proposal. In particular, we focused on relativistic field theories, and showed that there is a zoo of relativistic lattice fermions that can be addressed with this platform. Besides, we presented a toolbox to design particular assisted tunneling processes that lead us to the physics of topological insulators. Remarkably enough, this quantum simulator turns out to be extremely versatile, since most of the phases of the periodic table of topological insulators can be addressed.
Finally, let us comment on the possible combination of this proposal with the control of interactions already achieved in cold-atom gases. This might eventually boost experiments into regimes where classical numerical simulations fail, which we leave here as an outlook for future work. In particular, the problem of robustness of topological orders (classified for non-interacting theories) with respect to interactions is one of the most important challenges of the modern condensed matter [70]. We believe that a direct combination of our setup with Feshbach resonances will provide important insights into this unsolved question.