Interaction-Dependent Photon-Assisted Tunneling in Optical Lattices: A Quantum Simulator of Strongly-Correlated Electrons and Dynamical Gauge Fields

We introduce a scheme that combines photon-assisted tunneling by a moving optical lattice with strong Hubbard interactions, and allows for the quantum simulation of paradigmatic quantum many-body models. We show that, in a certain regime, this quantum simulator yields an effective Hubbard Hamiltonian with tunable bond-charge interactions, a model studied in the context of strongly-correlated electrons. In a different regime, we show how to exploit a correlated destruction of tunneling to explore Nagaoka ferromagnetism at finite Hubbard repulsion. By changing the photon-assisted tunneling parameters, we can also obtain a $t$-$J$ model with independently controllable tunneling $t$, super-exchange interaction $J$, and even a Heisenberg-Ising anisotropy. Hence, the full phase diagram of this paradigmatic model becomes accessible to cold-atom experiments, departing from the region $t\gg J$ allowed by standard single-band Hubbard Hamiltonians in the strong-repulsion limit. We finally show that, by generalizing the photon-assisted tunneling scheme, the quantum simulator yields models of dynamical Gauge fields, where atoms of a given electronic state dress the tunneling of the atoms with a different internal state, leading to Peierls phases that mimic a dynamical magnetic field.

Quantum many-body physics studies systems of interacting particles governed by the laws of quantum mechanics. This task becomes particularly challenging in a variety of contexts in which the interactions induce strong inter-particle correlations. For instance, this strongly-correlated behavior appears in condensed-matter models whenever the system cannot be divided into weakly-interacting parts, such that the whole cannot be understood as a sum of its parts and perturbative methods become futile [1]. This inherent complexity underlies the abundance of interesting phases of matter that emerge at different scales, but also the difficulty in understanding them from an original microscopic model (e.g. high-T c superconductivity [2]). The same occurs at much higher temperatures and densities, where quarks and gluons interact strongly, and lead to a variety of phases that defy our current understanding (e.g. quark matter [3]). In the opposite regime, that of extremely low temperatures and densities, ultracold atomic gases trapped by electromagnetic fields are gradually becoming a paradigm of strongly-correlated behavior in quantum many-body physics [4]. In contrast to the above condensed-matter and high-energy scenarios, ultracold atoms have a unique property: their microscopic properties can be fully characterized and controlled in experiments. This experimental control has reached such a status that the dream of exploiting a quantum system to understand the properties of a complex quantum many-body model (i.e. a quantum simulator [5]) is already an experimental reality [6].
Ultracold gases of neutral atoms can be trapped in periodic optical potentials obtained from the interference of laser beams. The dynamics of the atoms in these optical lattices resembles that of tightly-bound electrons in metals, such that this system can be considered to be a synthetic solid whose dimensionality and lattice structure can be experimentally tailored, while the nuisance of impurities, disorder, or other uncontrolled microscopic degrees of freedom present in real solids, is totally absent. Starting from this synthetic solid, it is possible to design a variety of quantum many-body models whose microscopic parameters can be experimentally characterized and controlled. For instance, the scattering of atoms leads to a short-range interaction that can be tuned all the way from weak to strong repulsion, such that the superfluid-insulator quantum phase transition of the bosonic [7][8][9] and fermionic [10][11][12] Hubbard models becomes accessible to experiments. For sufficiently strong repulsion, the half-filled Hubbard model leads to a Heisenberg antiferromagnet [13], which yields a playground for quantum magnetism with two-component bosonic atoms [14,15], and the starting point to study high-T c superconductivity with fermionic ones [16,17] upon controlled doping (i.e. inserting atomic vacancies with respect to the half-filled system).
In this article, we will combine this strongly-correlated behavior with external periodic drivings to obtain a flexible quantum simulator of quantum many-body models. In the context of optical lattices, there is a large body of relevant results regarding periodic drivings by modulations of the trapping optical potential. For instance, it is possible to periodically modulate the phase of the laser beams forming the optical lattice, as already demonstrated in experiments of chaotic dynamics with cold atoms [18,19]. Another possibility is to modify the detuning of these laser beams linearly in time, usually referred to as lattice acceleration, which leads to a linear gradient (i.e. constant force) in the lattice reference frame, and gives rise to Bloch oscillations [20]. From this perspective, the previous phase modulation [18,19] may also be interpreted as a periodic forcing. The combination of these two forces permitted probing the Wannier-Stark ladder spectrum [21], and testing the phenomenon of coherent destruction of tunneling in the absence of the gradient [22,23], and photon-assisted tunneling against the gradient [22,24]. We shall be particularly interested in such photon-assisted tunneling effect, whereby the atoms can tunnel in the presence of an energy penalty (i.e. the linear gradient) by absorbing photons from the external driving (i.e. the periodic phase modulation).
Photon-assisted tunneling (PAT) by phase modulation has also turned out to be a useful tool for quantum simulations. The dependence of the dressed tunneling on the modulation parameters has been used to drive the system across the superfluidinsulator transition [25,26], and to control the tunneling anisotropy of Bose-Hubbard models in triangular lattices leading to magnetic frustration [27,28]. A subject of research that has received considerable attention recently is the quantum simulation of orbital magnetism, whereby the atoms mimic the behavior of electrons in solids subjected to additional magnetic fields [6]. Since the atoms are neutral, one must design specific schemes to simulate the effect of artificial/synthetic magnetic fields [30,31], and PAT by phase modulation has also been exploited in this respect (see [29] for a recent review that also covers schemes that do not exploit PAT). When the phase modulation leads to an inhomogeneous periodic forcing [32], it is possible to dress the tunneling with an effective complex phase. Unfortunately, this simple proposal does not allow for the quantum simulation of synthetic magnetic fields [33], and alternative schemes have been considered. For instance, two-tone phase modulations lead to synthetic fluxes in arbitrary lattices [34], while the simpler single-tone phase modulations yield staggered fluxes in certain types of lattices [35][36][37].
A different possibility would be to abandon the periodic phase modulation, and investigate other types of drivings that can lead to the aforementioned synthetic Gauge fields. Instead of modulating the phase of the optical lattice, one can introduce a simple periodic driving by using a pair of slightly detuned Raman beams (i.e. a moving optical lattice), which has been considered in the context of PAT for trapped ions [38], ultracold atoms [39], and generic lattice models [40] that can be applied to a variety of contexts. For ultracold atoms, this moving optical lattice yields, in a certain regime, an inhomogeneous periodic modulation of the on-site energies of the effective Hubbard model, which can be exploited as a flexible PAT toolbox for quantum simulations of synthetic Gauge fields [38][39][40][41]. Here, the atoms tunnel in the presence of an energy penalty (i.e. again, a linear gradient) by absorbing photons from the external periodic driving (i.e. this time, the moving optical lattice), and acquire a Peierls phase that plays the role of a synthetic Gauge field, and depends upon the wavevectors of the Raman beams.
In this work, we explore a modification of this scheme by considering that the energy penalty can also be caused by the on-site Hubbard interactions, yielding a Hubbard blockade that inhibits the tunneling of atoms involving double occupation of a lattice site. The combination of this Hubbard blockade with the periodic driving by a moving optical lattice will induce an interaction-dependent photon-assisted tunneling. Let us note that the interplay of Hubbard interactions, linear gradients, and phase modulation of the optical lattice, has been shown to be responsible for interaction-shifted resonances in the PAT of Bose-Hubbard dimers [42] and chains [43]. Similar effects have been observed experimentally by considering a periodic modulation of the intensity of the optical-lattice laser beams [44,45], rather than the aforementioned phase modulation. This interactiondependent PAT can lead to new schemes to control effective magnetic Hamiltonians [44,45], or to methods that enhance the effects of three-body interactions [46]. We should also mention other proposals that are relevant for the particular subject of our work. These concern the engineering of density-dependent tunnelings by either combining laser-assisted schemes with statedependent lattices [47] in the spirit of the original proposal [30], or a periodic modulation of the Hubbard interactions [48][49][50][51].
In this work, we will show that the interaction-dependent PAT by a moving optical lattice offers a very flexible quantum simulator for paradigmatic models of strongly-correlated electrons, and can even allow for the quantum simulation of synthetic Gauge fields that are dynamical, in contrast to the static ones mentioned above.
This article is organized as follows. In Sec. II, we introduce the scheme to implement the interaction-dependent PAT with ultracold atoms in optical lattices, and derive a set of effective Hamiltonians that depend on the specific driving, lattice dimensionality, and fermionic/bosonic quantum statistics. The scope of the many-body phenomena that can be studied through these effective Hamiltonians is discussed in Sec. III. Finally, we present our conclusions and outlook in Sec. IV.

II. INTERACTION-DEPENDENT PHOTON-ASSISTED TUNNELING
In this section, we present a detailed proposal to combine PAT by periodic drivings with strong Hubbard interactions in experiments of ultracold alkali atoms in optical lattices. We show that by controlling (i) the atomic interactions by Feshbach resonances, and (ii) an additional moving optical lattice, one can exploit an interaction-dependent PAT to delve into interesting quantum many-body models that arise in the condensed-matter and high-energy scenarios.
The starting point is, as customary [4,52], a trapped atomic gas described in second quantization where we seth = 1 henceforth. Here, Ψ † σ (r), Ψ σ (r) create-annihilate atoms with mass m σ at the position r, and in the electronic state |σ corresponding to a particular energy level ε σ of the atomic groundstate manifold. To remain as general as possible, we consider that the components labelled by σ may correspond to the states of a bosonic gas, a fermionic one, or a mixture of both, which will determine the particular algebraic relations of the creation-annihilation operators.
We have introduced an optical trapping potential V ot (r r r) = ∑ α V 0,α sin 2 (kr α ) + 1 2 mω 2 t,α r 2 α that consists of: (i) A stateindependent periodic potential, where V 0,α are the ac-Stark shifts of independent pairs of retro-reflected laser beams of wavelength λ = 2π/k, which propagate along the axis α ∈ {x, y, z}, and are far detuned with respect to the excited atomic states. To obtain state-independent potentials, we assume that the detunings of the laser beams with respect to the excited states are much larger than the energy splittings in ε σ , and that the retro-reflected beams along each axis have parallel linear polarizations [53]. To obtain independent potentials along each axis, the pairs of interfering beams must have orthogonal polarizations, or detuned frequencies, with respect to other pairs of beams propagating along a different axis. Therefore, it is possible to tune the lattice depths V 0,α independently by controlling the beam intensities, which allows to tailor the effective dimensionality of the system. (ii) A harmonic trapping caused by a combination of the laser Gaussian profile and the retro-reflection scheme, where mω t,α λ 2 is the characteristic trapping energy assumed to be sufficiently weak ω t,α ω 0 σ ,α = 2E R,σ V 0,α /E R,σ , where E R,σ = k 2 /2m σ is the so-called recoil energy.
The final ingredient of the cold-atom Hamiltonian (1) is the s-wave scattering, which dominates at sufficiently low temperatures. This is described by a contact pseudo-potential V σ σ int (r − r ) = 4πa σ σ δ (r − r )/2µ σ σ characterized by the reduced masses µ σ σ = m σ m σ /(m σ + m σ ), and the scattering lengths a σ σ for the collisions of two atoms in the internal state |σ , σ . Such scattering lengths can be modified experimentally through an external magnetic field via the so-called Feshbach resonances [54]. In the following sections, we show how to exploit an interaction-dependent PAT as new tool to engineer quantum many-body Hamiltonians by tuning these scattering lengths appropriately in the presence of a weak moving optical lattice.
A. Scheme for a periodically-modulated ultracold Fermi gas Let us consider a single-species gas of fermionic atoms with two hyperfine states |↑ = |F, M , |↓ = |F , M , such that there is a unique mass m ↑ = m ↓ =: m and recoil energy E R,↑ = E R,↓ =: E R . We introduce the Wannier basis, where w(r − R i ) are the Wannier functions, and f i,σ are the fermionic operators that annihilate an atom of pseudospin σ = {↑, ↓} at the minima of an optical lattice potential R i labelled by the vector of integers i. We shall consider cubic optical lattices, although we note that the scheme detailed below can be directly applied to any other lattice geometry. In this basis, the general Hamiltonian (1) can be expressed in terms of the standard Fermi-Hubbard model [11], namely Figure 1: Scheme of the spin-independent PAT for fermions: (a) Laser scheme corresponding to the static optical lattice formed by the retro-reflected beams (red arrows) of frequency ω 2 , and the moving optical lattice formed by a slightly detuned Raman beam of frequency ω 1 (blue arrow) and the counter-propagating laser beam of frequency of frequency ω 2 (red arrow). All the laser beams have a linear polarization perpendicular to the plane of the figure, as illustrated by the filled circles. (b) Same as before, but considering that the detuned Raman beam (blue arrow) forms an angle with respect to the static optical lattice beams (red arrows). (c) Fermionic atoms in two hyperfine states |↑ (green circles), |↓ (orange circles) are trapped at the nodes of a static optical lattice potential (red lines). In the regime of strong s-wave scattering, atoms tunnel with strength t x between unoccupied sites centered at the energies ε ↑ , ε ↓ . Conversely, tunneling of one atom to an already-occupied site is inhibited by the large energy penalty t x U ↑↓ (i.e. Hubbard blockade). (d) Moving optical lattice potential (blue lines for snapshots of the wave traveling at speed v = ∆ω/∆k). The tunneling involving doubly-occupied sites can be reactivated when the atoms absorb r photons from the moving lattice, providing the required energy to overcome the interaction penalty. (e,f) The PAT contains diagonal (e) and off-diagonal (f) correlated events. The diagonal terms correspond to a dressed tunneling within the subspaces of single-or doubly-occupied sites, and is controlled by the Bessel function J 0 := J 0 (η). The off-diagonal tunneling connects the subspaces of singleand doubly-occupied sites, and is controlled by the Bessel function J r := J r (η) and dressed by the tunneling phase e ±irϕ . Also note that the residual dressed interaction δU ↑↓ is changed with respect to the original bare one U ↑↓ .
where we have introduced the unit vectors e α , and the notation σ = {↓, ↑} for σ = {↑, ↓}. Here, ε i,σ = ε σ + ∑ α 1 2 mω 2 t,α (R i,α ) 2 includes the hyperfine energies and the weak parabolic trapping potential, t α is the tunneling strength of atoms between neighboring potential wells along the α-axis, and U σ σ = U σ σ stands for the on-site interaction due to s-wave scattering, which only allows for interactions between fermions of a different state. As customary, we have neglected long-range tunnelings and interactions, which requires sufficiently deep optical lattices {V 0,x ,V 0,y ,V 0,z } E R . This can be justified considering that while longer range terms are exponentially suppressed with the distance by exp{−mω 0 σ ,α (R i,α − R j,α ) 2 /4}. Let us also note that, for ω t,α ω 0 σ ,α , the harmonic trapping does not modify the tunneling, but simply leads to a local term in the Wannier basis that has been incorporated in the local on-site energies ε i,σ of the Fermi-Hubbard model (2).
We consider the limit of very strong repulsion U ↑↓ t α , such that the bare tunneling events connecting single-occupied sites to doubly-occupied ones are energetically inhibited, as depicted in Fig. 1(c). We shall refer to this tunneling suppression as a Hubbard blockade by reminiscence of the Coulomb blockade that inhibits the sequential tunneling of electrons through quantum dots. The idea is to overcome this Hubbard blockade via the phenomenon of PAT (i.e. the fermions obtain the required energy for tunneling by absorbing photons from an external periodic driving). As shall be shown below, the tunneling of fermions between two lattice sites will depend on the density of fermions of the opposite pseudospin populating those sites, which shall be exploited to build a quantum simulator. We now discuss two possible periodic drivings that lead to such PAT, and organize our presentation by introducing the less demanding schemes first, adding more complexity gradually.
1. Two-component fermions in spin-independent moving optical lattices (i) One-dimensional scheme: To introduce the main ideas in the simpler setting, let us start by considering a one-dimensional (1D) Fermi-Hubbard model obtained from Eq. (2) for {V 0,y ,V 0,z } V 0,x , such that only tunneling along the x-axis is relevant As an external periodic driving, we consider a moving optical lattice stemming from a pair of non-copropagating laser beams along the x-axis. These beams are slightly detuned with respect to each other (i.e. traveling wave as opposed to the standing wave of the static optical lattice, Fig. 1(d)), but again far detuned with respect to the excited states (i.e. Raman beams). Moreover, they have the same linear polarization as the laser beams of the static optical lattice to ensure an spin-independent potential [53]. Since this moving lattice could induce a spurious tunneling due to the recoil kick imparted by the lasers, we assume that its intensity is much weakerṼ 0 V 0,x (i.e.t x =Ṽ 0 exp{− π 2 4 (V 0,α /E R ) 1/2 } t x in the Gaussian approximation). In this regime, the effect of the moving lattice is a periodic spin-independent modulation of the trapping frequencies of each potential well where ∆k = (k 1 − k 2 ) · e x is the wavevector difference, X i = λ 2 i stands for the minima of the original optical-lattice potential, ∆ω = ω 1 − ω 2 is the detuning of the laser beams, and ϕ is the relative phase with respect to the static optical lattice. By setting ∆ω ≈ U ↑↓ /r for a positive integer r ∈ Z, the above Hubbard blockade for U ↑↓ t α can be overcome through the absorption of r photons from the periodic driving (see Fig. 1(d)). To be more precise, as the driving comes from a two-photon ac-Stark shift, the process involves absorbing r photons from one laser beam and subsequently emitting them onto the other laser beam.
To provide explicit expressions for this interaction-dependent PAT, we move to the interaction picture with respect to U 0 (t) = T exp{i t 0 dτ(V int + H mod (τ))} , such that the fermionic annihilation operators become where we have gauged away an irrelevant phase by transforming the fermion operators [55]. The second part of the equality is obtained after introducing η =Ṽ 0 /∆ω, and using the Jacobi-Auger expansion for first-order Bessel functions J n (z), namely e iz sin θ = ∑ n∈Z J n (z)e inθ [56]. For simplicity, we set ∆kX i = ( 1 2 ∆kλ )i = πi, which can be achieved with laser beams of the standing and moving lattices of the same wavelength, and both propagating along the x-axis. In this configuration, it thus suffices to add a single laser beam detuned with respect to the optical-lattice laser beams (see Fig. 1(a)). However, this could be generalized to ∆kX i = ( 1 2 ∆kλ )i = πi/r, which may be relevant if the detuned Raman beam does not propagate along the x-axis, but makes some angle with respect to that axis (e.g. r=2 for an angle α = π/6, see Fig. 1 By substituting the expression (6) in the kinetic Hamiltonian H kin (t) = U 0 (t)H kin U † 0 (t), one finds where we have introduced the population difference operator ∆n i+1,σ = n i+1,σ − n i,σ , and a dynamical dressing function As announced earlier, the tunneling that connects single-occupied sites to doubly-occupied ones, yielding ∆n i+1,σ = ±1, is negligible in the absence of the drivingṼ 0 = 0. In this limit, the dressing function is f(t) = 1, such that the dressed tunneling can be neglected t x,σ (t) = t x e ∓itU ↑↓ ≈ 0 in a rotating-wave approximation for t x U ↑↓ (see Fig. 1(c)). By switching on the periodic drivingṼ 0 = 0, this tunneling becomes assisted by the harmonics of the dressing function that are close to resonance with the Hubbard interaction, namely for the integers fulfilling ±U ↑↓ = (n − m)∆ω ( Fig. 1(d)). In particular, by assuming that we can neglect the majority of tunneling events using a similar rotating-wave argument, except for those that satisfy n = r∆n i+1,σ + m. Accordingly, the dressing function becomes simplified We further assume that the laser detuning is chosen in such a way that r is an even integer, and set e −iπr∆n i+1,σ i = 1 for any population difference. Making use of the Neumann-Graf addition formula for Bessel functions [56], namely ∑ n∈Z J n (z)J n+ν (z)e inθ = J ν (2|z sin(θ /2)|) e i(π−θ )ν/2 , we can express the PAT in terms of a single Bessel function ; (e-f) initial state with a pair of atoms on the right well |0 1 , . We use the same criterion in all the figures: the dashed lines correspond to the tunneling for the un-drivenṼ 0 = 0 dimer, while the solid lines and symbols stand for the resonantly-driven dimer ∆ω = U ↑↓ /2 (i.e. two-photon assisted tunneling r = 2) withṼ 0 = 3∆ω corresponding to either the numerical solution of the exact Hamiltonian (4)-(5), or the numerical solution of the effective Hamiltonian (12)-(13), respectively. (a, b, c) Photon-assisted tunneling for the spin-up atom, and Pauli blockade of the spin-down atoms which cannot tunnel due to the exclusion principle. (d) Maximal arrival density of the spin-up atom max{ n 1,↑ (t) : 0 < t < 2π/t x }, which displays minima exactly at the zeros of the Bessel function J 0 (z 0,n ) = 0 (i.e. coherent destruction of tunneling). (e, f, g) Photon-assisted tunneling for both the spin-up and spin-down atoms. (h) Maximal arrival density of the spin-up atom max{ n 1,↑ (t) : 0 < t < 2π/t x }, which displays minima exactly at the zeros of the Bessel function J 2 (z 2,n ) = 0. In comparison to (d), we observe that the coherent destruction of tunneling depends on the background of spin-down atoms, leading to the correlated coherent destruction of tunneling exploited in Sec. III. which should be understood in terms of its Taylor series expansion.
The total time-evolution operator U(t) = U † 0 (t)e −it ∑ i δU ↑↓ n i↑ n i↓ e −iH eff t , can thus be expressed in terms of a time-independent Hubbard Hamiltonian of the form (2). However, the dressed tunneling strengths now depend on the density of fermions of the opposite pseudospin, and the residual Hubbard interaction depends on the resonance condition in Eq. (9), such that As announced in the introduction, the Hubbard-blockaded tunneling becomes activated through a PAT phenomenon, and leads to a density-dependent tunneling that can be written as follows where we have defined the hole number operators h i,σ = 1 − n i,σ . The first term in Eq. (13) describes the tunneling within the subspace of single-occupied sites H s , whereas the second one corresponds to tunneling within the subspace of doubly-occupied sites H d (see Fig. 1(e)). These subspaces can be described as two Hubbard sub-bands centered around ε s = 0 and ε d = δU ↑↓ . Finally, the third and fourth terms stand for tunneling events connecting the single-occupied to the doubly-occupied subspaces (see Fig. 1(f)). These four terms can be thus understood as diagonal and off-diagonal tunnelings. We note that a similar classification of the tunneling events of the original Hubbard model (4) can be performed by using H kin → ∑ i,σ (n i,σ + h i,σ )H σ ,i kin (n i,σ + h i,σ ). Let us emphasize, however, that the ratio of these diagonal/off-diagonal processes cannot be controlled, which contrasts the PAT Hamiltonians (13), where one can adjust the intensity of the moving optical latticeṼ 0 such that the ratio of the Bessel functions attains the desired value. This will be crucial to obtain a tunable t-J model with fully controllable parameters in Sec. III B. In Sec. III A, we will use this formulation to connect the effective model to the so-called bond-charge interactions, which leads to a quantum simulator of exotic Hubbard models. Moreover, the tunneling of one pseudospin acquires a complex phase that depends on the density of the other pseudospin, which will be crucial for the quantum simulation of dynamical Gauge fields in Sec. III C, when complemented with additional terms that allow us to control each pseudospin independently.
In order to test the validity of our derivations, we compare numerically the dynamics obtained from the effective Hamiltonian (12)- (13), and the periodically driven one (4)-(5) in the simplest setting: a Fermi-Hubbard dimer (see Figs. 2(a,e)). In , which does not display Hubbard blockade as the tunneling preserves the number of doubly-occupied sites. Nonetheless, the bare tunneling for the spin-up atoms (see dashed lines of Fig. 2(b)) is renormalized due to the periodic driving, as shown by the different population dynamics displayed by the solid lines (exact) and the symbols (effective). The excellent agreement between the solid lines and the symbols proves the validity of our derivations, and the accuracy of the interaction-dependent PAT Hamiltonian (12)- (13). Regarding the dynamics of the down-spin atoms, we note that these cannot tunnel due to the Pauli exclusion principle, as depicted in Fig. 2 , which suffers a Hubbard blockade as the tunneling must change the number of doubly-occupied sites. Hence, in the absence of the driving, the atomic tunneling is totally forbidden (see dashed lines of Fig. 2(f,g)). By switching on the driving, we observe that the tunneling of both spin-up and spin-down atoms is reactivated, as shown by the solid lines (exact) and the symbols (effective), which again show an excellent agreement supporting our analytical results.
Let us now address the phenomenon of correlated destruction of tunneling exploited in Sec. III for the quantum simulation of strongly-correlated models. According to Eq. (13), the tunneling is dressed by a different Bessel function depending on the particle-hole densities, and it can get coherently suppressed when the driving parameter η coincides with a zero of the corresponding Bessel function. In Fig.2(d), we observe this effect for the first pair of zeros, η = z 0,n with n = 1, 2, of the Bessel function J 0 (z 0,n ) = 0, which are displayed by the dashed dotted lines. We see how the maximal average population that reaches the left site of the Hubbard dimer vanishes when the driving ratio coincides with any of the zeros. In Fig. 2(h), we see that for a different particle-hole distribution, the coherent destruction of tunneling takes place at the zeros of a different Bessel function, namely η = z r,n for n = 1, 2 zeros of the Bessel function J r (z r,n ) = 0 for the chosen r = 2. Since the zeros of the two Bessel function do not coincide, we can independently suppress the tunneling correlated to a particular particle-hole distribution (i.e. correlated destruction of tunneling ), which will be relevant in in Sec. III.
We have so far presented numerical tests supporting the validity of the resonant PAT, ∆ω = U ↑↓ /r, such that the residual interactions of the dressed Fermi-Hubbard model (12) vanish δU ↑↓ = 0. However, our analytical results show that finite Hubbard interactions δU ↑↓ = (U ↑↓ − r∆ω) can be achieved by changing the velocity of the moving optical lattice ∆ω = U ↑↓ /r, which will be crucial for several quantum simulations in Sec. III. Let us test this result by numerically integrating an adiabatic evolution according to the effective (12)-(13) and periodically-driven (4)-(5) Hamiltonians for a half-filled dimer ( Fig. 3(a)). We study the evolution of the system (see Fig. 3(b)), for a slow ramp of the tunneling strength t x → t x (1 − (δt x )t) with a rate δt x . Initially, the dimer is prepared in the groundstate, which resembles a Fermi sea with the spin-up/down atoms delocalized along the dimer, corresponding to the groundstate of the Fermi-Hubbard dimer for t x δU ↑↓ . After the quench t f ≈ 1/δt x , the dimer should be in a spin singlet state corresponding to the groundstate of an antiferromagnetic Heisenberg model that arises for t x δU ↑↓ In Fig. 3(c), we represent the numerical results for the Heisenberg-singlet fidelity F s (t) = | HS|Ψ(t) | 2 as a function of the showing two different tunneling dressings that lead to the two-tone beating displayed in (b) for the initial state |Ψ 0 = |0 1 , ↑↓ 2 . As usual, the solid lines correspond to the numerical solution of the exact Hamiltonian (4)- (5) in the presence of an additional gradient, while the symbols stand for the numerical solution of the effective Hamiltonian (12) with the modified tunnelings (15) due to the presence of the gradient. (c) Scheme for PAT for states with single occupancies | ↓ 1 , showing two different tunneling dressings that yield the two-tone beating displayed in (d) for the initial state |Ψ 0 = | ↓ 1 , ↑ 2 . ramp time, and for different ramp rates. We observe that the fidelity approaches F s (1/δt x ) ≈ 1 for the very slow ramps, where the adiabatic evolution is expected to be more accurate. Once again, the good agreement between the effective (12)-(13) and periodically-driven (4)-(5) Hamiltonians, support our claim that one can study the effects of finite Hubbard interactions, and their interplay with the dressed PAT tunneling.
Before moving to the PAT in higher dimensions, let us mention that we could gain additional flexibility in the scheme by introducing an additional linear gradient, which may come from a lattice acceleration, an external electric field, or a magneticfield gradient. If we tune the gradient such that it coincides with the on-site interaction, we can generalize Eq. (12) by substituting According to this expression, the off-diagonal tunnelings connecting doubly-to single-occupied sites (i.e. fourth term in Eq. (15) represented in the upper panel of Fig. 4(a)), and single-to doubly-occupied sites (i.e. third term in Eq. (15) represented in the upper panel of Fig. 4(c)), depend on different Bessel functions. This leads to a two-tone beating in the tunneling dynamics, as shown in Figs. 4(b),(d), which also serve as tests of the validity of our analytical derivations.
(ii) Higher-dimensional scheme: The scheme presented above can be directly generalized beyond 1D. The static optical lattice should be modified such that it allows for tunneling along two As can be observed from Eqs. (6)- (8), to assist the tunneling along a given direction, it is crucial that the periodic modulation (5) has a phase that varies along that particular direction. Therefore, we would need to include additional moving optical lattices that propagate along the remaining axes, dressing the corresponding tunneling along two α = {x, y}, or three α = {x, y, z} directions. One may consider adding one independent detuned laser beam per axis, paralleling the construction of the one-dimensional case. Otherwise, one could simply tilt the laser beam of the one-dimensional case, such that it has a non-vanishing projection propagating along each axis. The former scheme would lead to independent moving lattices along each axis whose intensity and frequency can be tuned separately, whereas the latter would lead to a non-separable moving lattice that dresses all the different tunnellings with the same intensity and frequency, albeit one could play with the propagation angle.
For simplicity, we consider the first situation, such that the periodic driving is where we have introduced the labeling indexes i = (i x , i y ) for 2D, and i = (i x , i y , i z ) for 3D. As before, we have assumed that for 2D (Ṽ 0,x V 0,x , andṼ 0,y V 0,y ), and for 3D (Ṽ 0,x V 0,x ,Ṽ 0,y V 0,y , andṼ 0,z V 0,z ), such that the moving lattices do not Figure 5: Scheme of the spin-dependent PAT for fermions: (a) Laser scheme similar to Fig. 1, but with the linear polarization of the detuned Raman beam (blue arrow) rotated by an angle θ p with respect to the static lattice laser beams. This gives rise to a spin-dependent moving optical lattice. (b) Same as before, but considering that the detuned Raman beam (blue arrow) forms an angle with respect to the static optical lattice beams (red arrows). (c, d) The moving optical-lattice potential for each pseudospin (green lines for |↑ , and orange lines for |↓ ) reactivates the tunneling involving doubly-occupied sites. The PAT contains diagonal terms (e) corresponding to a spin-dependent dressed tunneling within the subspaces of single-or doubly-occupied sites controlled by the Bessel function J 0,σ := J 0 (η σ ). The off-diagonal terms (f) contain a spin-dependent tunneling connecting the subspaces of single-and doubly-occupied sites, which are controlled by the Bessel function J r,σ := J r (η σ ) and dressed by the spin-dependent tunneling phase e ±ir σ ϕ σ . modify the bare tunneling and only lead to a periodic modulation of the on-site energies. Each of these moving lattices assists the tunneling along a given direction, and does not interfere with the tunnelings along the remaining axes. Accordingly, the interaction-dependent PAT is a direct generalization of (12), which requires a parameter regime where ∆n i+e α ,σ = n i+e α ,σ − n i,σ , the bare tunnelings are approximated by Eq. (3), and the dressed tunnelings depend on where we have introduced η α =Ṽ 0,α /∆ω α . It is interesting to note that controlling the intensity difference of each moving optical lattice, we can tune the spatial anisotropy of the dressed tunnellings. The possibility of generalizing to 2D is especially interesting in the context of the t-J model, and its connection to high-T c cuprate superconductors, as outlined in Sec. III B.
We have thus seen that the interaction-dependent PAT with a moving optical lattice leads to effective Hubbard models of any dimensionality with dressed tunnelings that are density dependent. In the following section, we will show that by considering a state-dependent moving optical lattice, the PAT scheme becomes more flexible, which will allow us to target other quantum many-body models, in particular dynamical Gauge fields.

Two-component fermions in spin-dependent moving optical lattices
(i) One-dimensional scheme: Let us once again start with the less demanding case of 1D. We note that the far-detuned moving optical lattice can become state-dependent if the laser-beam polarizations are not collinear (see Fig. 5 (a,b)). This occurs even for detunings that are larger than the Zeeman and hyperfine splittings, as far as they do not exceed the fine-structure splitting [53]. For fermionic alkali atoms, this could turn to be incompatible with reaching ultracold temperatures, as the fine-structure splitting is rather small, and the residual photon scattering may become appreciable [31]. However, we stress that the moving lattice is by construction much weaker than the static spin-independent one. In fact, we can reduce the residual photon scattering by orders of magnitude by lowering the intensity of the moving-lattice laser beams, as far as their detuning is simultaneously lowered, such that the ratio η controlling the PAT (13) remains constant. We will thus assume that a conservative state-dependent moving optical lattice can be realized without increasing the photon scattering and heating the ultracold atomic gas.
Scheme for the spin-dependent PAT for the state |0 1 , ↑↓ 2 , showing that the spin-down atoms can be coherently frozen, as shown in (b) for V 0↑ = 3∆ω, and V 0↓ = 5.13∆ω such that J 2 (5.13) = 0. As usual, the solid lines correspond to the numerical solution of the exact Hamiltonian (4) with the spin-dependent periodic modulation (21), while the symbols stand for the numerical solution of the effective Hamiltonian (24) with the modified tunnelings (25). (c) Scheme for the spin-dependent PAT for the state |0 1 , ↑↓ 2 , showing that the spin-down atoms can be coherently frozen, as shown in (d) for V 0↑ = 5.13∆ω such that J 2 (5.13) = 0, and V 0↓ = 3∆ω.
In this case, we can generalize the driving (5) by including a state-dependent periodic modulation of the on-site energies where the spin-dependent amplitude must again fulfillṼ 0,σ V 0,x , and ϕ σ stands for a phase difference with respect to the static lattice that is generally state dependent [57,58] (see Figs. 5 (c,d)). Going back to the spin-independent scheme (5), the new modulation (20) can be achieved by rotating the polarization of the laser beam that is slightly detuned with respect to the staticlattice lasers [57]. In this case, the spin-dependent driving amplitudesṼ 0,σ can be tuned by controlling such an angle, or instead the direction of propagation of the laser beam with respect to the quantization axis [58]. Another possibility would be to resolve the hyperfine structure, such that one can exploit selection rules in the ac-Stark shifts. In fact, for pseudospins corresponding to the maximally-polarized Zeeman sublevels, it is possible to obtain optical lattices that selectively address a single pseudospin (i.e.Ṽ 0,↑ = 0,Ṽ 0,↓ = 0), or viceversa, as realized in ion-trap experiments [59]. This leads to a spin-dependent driving where the wavevector ∆k σ , detuning ∆ω σ , intensityṼ 0,σ , and relative phase ϕ σ can all be controlled independently for each pseudospin Paralleling the previous section, we will consider equal wavelengths of the static and moving optical lattices, such that ∆k σ X i = ( 1 2 ∆k σ λ )i = πi, although we remark again that the scheme also works for other propagation angles. Once the new periodic drivings (20)- (21) have been discussed, we can address the interaction-dependent PAT they give rise to. We shall use Eq. (21), as the results also encompass those related to the driving (20). By reproducing the steps that lead to the effective Hamiltonian (12) for the spin-independent driving, we find a parameter regime analogous to Eq. (9), namely and a new dressing function of the tunneling that becomes spin-dependent, namely where η σ =Ṽ σ  Fermi-Hubbard tetramer with parameters t x = t y = 0.1, U ↑↓ = 20, and subjected to a spin-dependent moving lattice with ϕ y,σ = 0, and ϕ x,σ = π/2, and ∆ω α,σ = U ↑↓ /2, ∀α, σ . (a) Scheme for the PAT for the state | ↑↓ 1 , 0 2 , 0 3 , 0 4 , where the frozen spin-down atom induces a π-flux in the tunneling of the spin-up atom, leading to the Aharonov-Bohm destructive interference in (b). (c,d) Populations for V 0↑ = ∆ω, and V 0↓ = 5.13∆ω such that J 2 (5.13) = 0, and the spin-down atom is frozen. As usual, the solid lines correspond to the numerical solution of the exact Hamiltonian (4) with the spin-dependent periodic modulation (21), while the symbols stand for the numerical solution of the effective Hamiltonian (28). Due to the interference, the spin-up atom cannot reach the opposite corner n 3,↑ (t) = 0. (e-h) Same as above, but for the state | ↑↓ 1 , ↓ 2 , ↓ 3 , ↓ 4 , where the spin-down atoms cannot tunnel due to the Pauli exclusion principle. For this density background, the tunneling phase of the spin-up atoms vanishes since there is no spin-down density gradient. Hence, this initial state does not lead to Aharonov-Bohm destructive interference, and the spin-up atom does indeed reach the opposite corner.
Remarkably, we find that the amplitude of the density-dependent tunneling can be controlled independently for each pseudospin and that the tunneling phase of one pseudospin depends on the density of the other pseudospin, which will be crucial for the quantum simulation of dynamical Gauge fields in Sec. III C. In order to benchmark these predictions, we study numerically a spin-dependent coherent destruction of tunneling in a Fermi-Hubbard dimer subjected to the spin-dependent moving optical lattice. According to Eq. (25), the dressed tunneling of a doubly-occupied half-filled dimer (see Fig. 6(a)) depends on the spin of the atom, such that the spin-up atoms tunneling is controlled by the Bessel function J r ↑ (η ↑ ), whereas the spin-down atoms tunneling depends on the Bessel function J r ↓ (η ↓ ). Therefore, by controlling the intensities of the the spin-dependent moving lattice, the dressed tunneling of the spin-down atoms can be coherently destructed J r ↓ (η ↓ ) = 0, while the spin-up atoms hop freely in the lattice J r ↑ (η ↑ ) = 0 (see Figs. 6(a)-(b)). Conversely, we can coherently freeze the spin-up atoms J r ↑ (η ↑ ) = 0, while the spin-down atoms hop freely in the lattice J r ↓ (η ↓ ) = 0 (see Figs. 6(c)-(d)). Let us emphasize the excellent agreement between our analytical description (symbols), and the exact dynamics of the periodic Hamiltonian (solid lines).
(ii) Two-dimensional scheme: The scheme presented above can be generalized beyond 1D. For the shake of concreteness, and for its particular interest in connection to the t-J model in Sec. III B, and the dynamical Gauge fields in Sec. III C, we will restrict to 2D (V 0,z {V 0,x ,V 0,y } E R ). The idea is to consider spin-dependent moving optical lattices along the x and y axes such that the relative phases of the moving lattices fulfill ϕ σ ,x ≥ 0, but ϕ σ ,y = 0. In this case, and after following the same steps as above in an analogous parameter regime t x ,t y , δU ↑↓ = (U ↑↓ − r σ ,α ∆ω σ ,α ) U ↑↓ ≈ r σ ,α ∆ω σ ,α = r σ ,α ∆ω σ ,α , ∀α, α , σ , σ one can derive the following effective Hamiltonian where η σ ,α =Ṽ σ 0,α /∆ω σ ,α . We thus see that when atoms tunnel along the x-axis, they acquire a dynamical phase that depends on the density of the other pseudospin, whereas they experience a vanishing phase when tunneling along the y-axis. This will be equivalent to the so-called Landau Gauge in Sec. III C, which is accompanied by a non-vanishing dynamical Wilson loop.
So far, all of our numerical tests have been independent of the phase of the moving optical lattices, as we have only addressed a Fermi-Hubbard dimer. By considering the simplest 2D case, a Fermi-Hubbard tetramer forming a square plaquette, we can already test numerically the predicted effect of the moving lattice phase, which according to Eq. (28), induces a density-dependent Peierls phase in the tunneling. To observe the effects of such a Peierls phase, we shall first exploit the above spin-dependent destruction of tunneling, or the Pauli exclusion principle, to freeze the dynamics of the spin-down atoms (see Figs. 6(a-b)). Then, the immobile spin-down atoms yield a density background that modifies the tunneling phase of the spin-up atoms. We explore this possibility in Fig. 7 for two different density distributions of the spin-down atoms, which lead to the presence/absence of an Aharonov-Bohm destructive interference in the tunneling dynamics of the spin-up atom. Besides confirming the validity of our analytical description (28) in a transparent scenario, let us note that by lifting the coherent destruction of tunneling, and allowing the spin-down atoms to tunnel, the Peierls phase will acquire its own dynamics, which will be crucial for the quantum simulation of dynamical Gauge fields in Sec. III C. At this point, the reader may skip the following sections, and move directly to the quantum simulator applications dealing with fermionic models in Sec. III.
where the Hamiltonian parameters coincide with the fermionic ones (3) under similar approximations. However, due to the different statistics, s-wave scattering among pairs of atoms with the same electronic state are now allowed (i.e. U ↑↑ ,U ↓↓ = 0), which gives more freedom for the interaction-dependent PAT. We shall discuss different regimes of interest for the quantum simulation of bosonic quantum many-body models, which can be achieved by tuning the Feshbach resonances.

Two-component hardcore bosons in spin-independent moving optical lattices
Let us start from the 1D Bose-Hubbard model obtained from Eq. (47) for {V 0,y ,V 0,z } V 0,x , such that only tunneling along the x-axis is relevant. In the hard-core limit, double occupancies of bosons with the same pseudospin are energetically forbidden (i.e. U ↑↑ ,U ↓↓ U ↑↓ t x ). In this limit, we can project out all states with sites occupied by more than one boson of the same pseudospin provided that the filling is n i,σ ≤ 1. By using the corresponding projector P s , we can map the bosonic creation-annihilation operators onto an su(2) spin algebra In such a hardcore limit, the Bose-Hubbard model (47) only contains the on-site Hubbard interactions for two bosons of opposite pseudospin, namelỹ Analogously, we must also project the spin-independent periodic modulation due to the moving optical lattice which, in the same regime as discussed for fermions (5), yields We can now proceed in analogy to the fermionic gas, as the difference between the fermionic operators and the hardcore-boson ones does not change any of the steps of the derivation. We move to the interaction picture with respect to the projected driving and Hubbard interactions U 0 (t) = T exp{i t 0 dτ(Ṽ int +H mod (τ))} , and follow the same steps as in the fermionic case to express the time-evolution operator where the dressed tunnelling strengths and phases are again density dependent. Given the one-to-one correspondence with the fermionic scheme, the numerical results to test the validity of our derivation by comparing the full and effective dynamics would not add anything different from Fig. 2, and we shall not include them here.
Since the hardcore-boson and the fermionic number operators have the same algebraic properties, we can rewrite this densitydependent tunnelling strength in complete analogy to the fermionic case where the hardcore hole operator ish i,σ = 1 −ñ i,σ . Therefore, the dressed tunnelings for hardcore bosons can be described pictorially by Figs. 1(e,f), which distinguish events that preserve/modify the double occupancy of the tunneling sites. Let us also note that the fermionic schemes for spin-dependent drivings in Sec. II A, and the generalization to higher dimensions, equally applies to the bosonic gas in this hardcore limit. Although we shall focus on the fermionic applications of the quantum simulator in Sec. III, we emphasize that all the quantum many-body models discussed there have a hardcore-boson counterpart, including the dynamical Gauge fields.

Two-component softcore bosons in spin-independent moving optical lattices
The objective of this section is to relax the hardcore constraint U ↑↑ ,U ↓↓ U ↑↓ t x , which forbids double occupancies of bosons with the same pseudospin. Let us, however, start by understanding the PAT of a single-component bosonic gas described by Eq. (47), but restricted to a single pseudospin (e.g. σ =↑). For notational convenience, we drop the pseudospin index, such that the Hamiltonian corresponds to the standard Bose-Hubbard model According to our discussion of the Hubbard blockade t x U, the tunnelling of bosons that changes the parity of the occupation number is energetically forbidden (see Fig. 8(a)). As customary, we activate this tunnelling by means of the periodic modulation  given by a moving optical lattice acting on the bosonic atoms (see Fig. 8(b)). Due to the different Hubbard interaction, which now involves a single pseudospin, one cannot use the expression in Eq. (6). Instead, we consider the interaction picture of a bond operator

again up to an irrelevant Gauge transformation, can be shown to be
where U 0 (t) = T exp{i t 0 dτ(V int + H mod (τ))} is the interaction-picture unitary. After defining the bosonic population difference operator ∆n i+1 = n i+1 − n i , and using the Jacobi-Auger expansion for each of the time-dependent exponentials, the expression of the kinetic energy, H kin (t) = U 0 (t)H kin U † 0 (t), becomes with exactly the same modulation function as in Eq. (8). We can then proceed by following the same steps as for the fermionic case, to find a parameter regime and finally obtain the effective interaction-dependent PAT Hamiltonian Here, we observe that the dressed tunnelling depends on the density of bosons of the same pseudospin at the sites involved in the tunnelling event (see Figs. 8(c,d)). Since the boson number per lattice site is not restricted anymore by the hardcore constraint, we cannot express it as the simple polynomial (34) quadratic in the density operators, but rather as a highly non-linear term.
Using the orthogonal projectors onto subspaces with a well-defined difference number of bosons P ∆n i+1 =± , this non-linearity becomes apparent where N b is the total number of bosons loaded in the optical lattice, and we have again assumed that r is an even integer. In analogy to the studies for the phase-modulation driving [42,43], we observe that there will be interaction-shifted resonances Initial state |0, 1 , which displays minima in max{ n 1 (t) : 0 < t < π/t x } exactly at the zeros of the Bessel function J 0 (z 0,n ) = 0. (a-b) Initial state |0, 2 , which displays minima in max{ n 1 (t) : 0 < t < π/t x } exactly at the zeros of the Bessel function J 2 (z 2,n ) = 0.
that correspond to the zeros of the Bessel functions for different density backgrounds. In fact, using our formalism, one could derive similar analytic expressions for a phase-modulation driving [42,43]. For the intensity-modulated lattices of the recent experiments [44,45], the situation is simpler as there can only be one particular occupation that is resonant with the periodic modulation of the tunneling matrix element, and no Bessel functions arise. The possibility of engineering Bose-Hubbard models with a density-dependent tunneling strength and phase by our PAT scheme can be considered as an alternative to the proposals based on periodic modulations of the s-wave scattering length [48,49]. Let us finally note that it is possible to generalize this scheme to higher dimensions, paralleling the fermionic case Sec. II A. In order to test the validity of our derivations, we have numerically compared the time-evolution predicted by either the effective Hamiltonian (40)-(41), or the periodically driven one (35)-(36) for a Bose-Hubbard dimer in the regime of interactiondependent PAT. In Fig. 9, we explore the real-time dynamics for different configurations of atoms in the initial state: (i) Nonblockaded regime: Figs. 9(a-b) represent the dynamics of the initial atomic configurations |0, 1 = b † 2 |vac , which does not display Hubbard blockade as it consists of a single atom. Yet, the bare tunneling (see dashed lines of Fig. 9(b)) is renormalized due to the periodic driving, as shown by the different population dynamics displayed by the solid lines (exact) and the symbols (effective). The same renormalization occurs for the configuration |1, 2 in Figs. 9(c-d), which is not blockaded as the overall occupation parity is conserved in the tunneling process. On top of the dressing of the tunneling, we observe a bosonic enhancement which leads to a doubled tunneling rate with respect to Fig. 9(b). (ii) Hubbard-blockaded regime: Figs. 9(e-h) represent the dynamics of the initial atomic configuration |0, 2 , which suffers a Hubbard blockade as the tunneling must change the overall occupation parity. Hence, in the absence of the driving, the atomic tunneling is totally forbidden (see dashed lines of Fig. 9(f)). By switching on the driving, we observe that the tunneling is reactivated, as shown by the solid lines (exact) and the symbols (effective). So far, all these simulations correspond to the resonant PAT, where the parameter regime (39) is achieved for δU = 0. In Figs. 9(g-h), we explore the off-resonant case δU > 0, and the possibility of describing the driving detuning as a residual Hubbard interaction. The agreement between the solid lines (exact) and the symbols (effective) in Fig. 9(h), shows that this is indeed the case. We observe that, as δU is increased, the periodic exchange of particles is gradually inhibited, as one would expect since the single-and doubly-occupied Hubbard bands become more and more separated in energy.
As a further numerical proof of the consistency of our effective description, let us explore the phenomenon of coherent destruction of tunneling. According to Eq. (41), the effective tunneling is dressed by a different Bessel function depending on the densities of the bosonic atoms. For instance, the PAT tunneling of Fig. 10(a) is controlled by J 0 (η), whereas the tunneling of Fig. 10(c) is controlled by J r (η). Therefore, whenever the driving parameter η coincides with a zero of the corresponding Bessel function, the tunneling should get coherently suppressed. In Fig.10(b), we observe this effect at η = z 0,n for n = 1, 2, 3 zeros of the Bessel function J 0 (z 0,n ) = 0, which are displayed by the dashed dotted lines. We see how the maximal average population that reaches the left site of the Hubbard dimer vanishes when the driving ratio coincides with any of the zeros. In Fig.10(d), we see that for a different atomic density distribution, the coherent destruction of tunneling takes place at the zeros of a different Bessel function, namely η = z r,n for n = 1, 2, 3 fulfilling J r (z r,n ) = 0 for the chosen r = 2.
Once the interaction-dependent PAT for the single-pseudospin bosons has been understood, and its validity has been checked numerically, we can turn our attention to the situation of two-pseudospin bosons without the hardcore constraint We shall be interested in the regime of U ↑↓ U ↑↑ ,U ↓↓ ,t x , where double occupancy of bosons of the same (different) pseudospin is allowed (forbidden). In this regime, as the intra-spin interactions do not blockade the tunneling, we only include the inter-spin interactions in the interaction picture U 0 (t) = T exp{i t 0 dτ(V int,↑↓ + H mod (τ))} , where a spin-independent moving optical lattice is applied to both pseudospins One can then see that, in the following parameter regime the dressed tunnelings will only depend on the density of atoms of the opposite pseudospin, as occurs for the hardcore bosons or the fermions. Therefore, the effective Hamiltonian becomes where we have introduced the residual interactionsŨ ↑↑ = U ↑↑ ,Ũ ↓↓ = U ↓↓ , andŨ ↑↓ =Ũ ↓↑ = δU ↑↓ . For ϕ = 0, we get an exotic Bose-Hubbard model with the analogue of the fermionic bond-charge interactions, whereby the dressed tunneling of one pseudospin depends on all the possible density backgrounds of the other pseudospin through the corresponding Bessel functions. The difference with respect to the fermionic bond-charge interactions (13) is the highly non-linear function of the bosonic densities, namely It is also worth commenting that, had we set U ↑↓ = U ↑↑ = U ↓↓ t x in the blockade conditions, the dressed tunneling would depend on the density of both pseudospins J r∆n i+1,σ (η) e −irϕ∆n i+1,σ → J r(∆n i+1,σ +∆n i+1,σ +1) (η) e −irϕ(∆n i+1,σ +∆n i+1,σ +1) , which might also be interesting regarding exotic Bose-Hubbard models. Let us note that, once again, this effective description matches perfectly the dynamics of the driven two-component Bose-Hubbard dimer (see Fig. 11). Once again, we could generalize to higher dimensions, or to spin-dependent drivings, and study the strongly-correlated models that arise. However, the properties particular for the restricted number of particles of fermions and hardcore bosons, which allow for instance to build a quantum simulator of dynamical Gauge fields (see Sec. III), cannot be generalized to the soft-core boson case. Let us finally note that, by considering the hardcore-boson limit on Eq. (45), one recovers the previous result (33).

C. Scheme for a periodically-modulated ultracold Fermi-Bose mixture
Let us now turn our attention to a mixture of bosonic and fermionic atoms, and consider that the two hyperfine states of the original formulation (1) Here, the Hamiltonian parameters now depend on the fermionic/bosonic nature of the atoms, such as the tunnelings which differ due to the different masses m f , m b through the corresponding recoil energies E R,f , E R,b . The Hubbard interactions  not only differ by the mass ratio, but also by the different scattering lengths for a boson-fermion a bf , and a boson-boson a bb collision. Once again, the fermion-fermion collisions are forbidden by the Pauli exclusion principle. We shall discuss different regimes of interest for the quantum simulation of Bose-Fermi quantum many-body models, which can be achieved by tuning the Feshbach resonances, and controlling these scattering lengths appropriately.
Let us focus on the 1D case. The moving optical lattices can be made dependent on the bosonic/fermionic atomic species, as the corresponding atoms usually have a very different atomic level structure. Therefore, in analogy to Eq. (21), we consider the following periodic driving where the amplitudesṼ 0,f/b wave vectors ∆k f/b , frequencies ∆ω f/b , and relative phases ϕ f/b , depend on the particular atomic species, and can be controlled experimentally. These periodic drivings will assist the tunneling against an energy penalty given by the Bose-Fermi interaction, which is very large where we have introduced the residual Bose-Fermi interactions δU bf , and two integers r b , r f ∈ Z that determine how many photons are absorbed from the periodic driving to overcome the energy penalty, and assist the tunneling.
In analogy with the two-component softcore bosons, only the interspecies Hubbard interactions inhibits the tunneling, such that the interaction picture must be U 0 (t) = T exp{i t 0 dτ(V int,bf + H mod (τ))} , leading to the time-evolution operator U(t) = U † 0 (t)e −it ∑ i δU bf n i,b n i,f e −iH eff t , with the following effective Bose-Fermi Hubbard Hamiltonian We thus observe that the bosonic/fermionic dressed tunnelings, both the amplitude and phase, depend on the fermionic/bosonic densities. Moreover, they can be independently tuned by controlling the parameters of the bosonic/fermionic moving lattice, which will be interesting for the quantum simulator of dynamical Gauge fields (see Sec. III). We also note that the generalization to 2D yields a Bose-Fermi analogue of the effective Hamiltonian (28).

III. APPLICATIONS OF THE PHOTON-ASSISTED TUNNELING QUANTUM SIMULATOR
We have introduced a scheme to implement an interaction-dependent PAT with trapped ultracold atoms. In this way, we have obtained a toolbox with different effective Hamiltonians that depend on the quantum statistics, dimensionality, and spindependent/spin-independent nature of the periodic driving used to assist the tunnelling. In this section, we will discuss how such a toolbox can be exploited for the quantum simulation of interesting problems in condensed matter and high-energy physics. Instead of focusing on a particular application, we have decided to stress the wide scope of the proposed quantum simulator by describing a variety of interesting quantum many-body models. All these problems share a common feature, they are described by complex quantum many-body models, either on the lattice or in the continuum, which still present open questions that defy the capabilities of existing numerical methods. We will describe the context of the particular models that can be targeted with the quantum simulator, and try to discuss the essence of the phenomena that they try to capture. Moreover, we will highlight what we believe are open questions of the models that have been studied in detail over the years, and point to future work required to understand the models that have not been explored in such a detail.

A. Bond-charge interactions and correlated electron-hole tunnelings
In the original derivation of the Hubbard model for electrons in transition metals [10], it was shown that a variety of longerrange terms also arise in the Hamiltonian. In addition to the tunnelling and the on-site density-density coupling, Coulomb interactions also yield nearest-neighbor terms that can be described as the repulsive interaction between: (i) charges localized at neighboring ions V, (ii) charges localized at an ion and a neighboring bond/link X, and (iii) charges localized at two neighboring bonds/links W. These bond-charge interactions can be responsible for a host of interesting effects in the context of charge density waves [61], ferromagnetism in itinerant electron systems [62], or alternative mechanisms of superconductivity based on electron holes [63]. However, it has been argued that the required values of X,W with respect to U,V to observe such effects are not likely to be realized in standard solid-state materials [64]. On the other hand, the possibility of controlling these terms in the synthetic solids offered by ultracold atoms in optical lattices has recently raised some interest in the community [65].
From the perspective of ultracold trapped atoms, one can evaluate these bond-charge terms for the contact interaction in Eq. (1). In the Wannier basis, one can introduce the bond operators, B σ i,j = f † i,σ f j,σ + f † j,σ f i,σ , to account for the state-dependent density located at the bond (i, j). Then, the bond-charge terms modify the standard Hubbard Hamiltonian (2) by introducing where we can express the different interaction strengths in the Gaussian approximation as follows From these expressions, we observe that the ratios V α /U,W α /U, and X α /U are vanishingly small for deep optical lattices [66]. An interesting possibility to reach regimes {V α ,W α , X α } ∼ U, where the bond-charge interactions can lead to new phases of matter, is to consider ultracold dipolar gases [65] or, as we show in this section, to exploit the PAT toolbox. At this point, we introduce an additional modification of the standard Hubbard model (2), whereby the tunnelling is not only modified by the density at separate sites (i.e. bond-charge interaction X α in Eq. (53)), but also by the density-density correlations Here,X α is the strength of this tunnelling, which cannot be obtained from any two-body interaction in the Wannier basis. In fact this term would require rather exotic three-body interactions, which cannot be realized even with ultracold dipolar gases [65]. These terms are interesting in a condensed-matter context, where they appear after reducing models with hybridized bands to single-band Hubbard models, as occurs in intermediate-valence solids [68], and the high-T c cuprates [69]. As shown below, the PAT toolbox can control all these terms in the effective cold-atom Hamiltonian.
1. Correlated electron-hole tunnelings: bond-ordered waves and triplet pairing In this section, we focus on a quantum simulator of the Fermi-Hubbard model with tunable ratios of t α /U, X α /U,X α /U, which may lead to very interesting many-body effects. Such a Fermi-Hubbard model H FH = H loc + H corr kin +V int can be rewritten in terms of asymmetric tunnelings that are correlated to the electron/hole occupation where we have introduced tunnelings in a hole-hole background t α hh = −t α , in an electron-hole background t α eh = −t α + X α , and in an electron-electron background t α ee = −t α +X α + 2X α . Let us focus for simplicity on the 1D limit of Eq. (56). The effective cold fermion Hamiltonian (12)- (13), which is obtained through the PAT by a spin-independent moving lattice, already contains these correlated particle-hole tunnelings. For the parameter regime fulfilling (9) with r = 2, and setting the moving lattice relative phase to ϕ = 0, one finds Accordingly, the correlated tunneling asymmetry can be tuned by modifying the intensity of the moving optical latticeṼ 0 , or its frequency ∆ω, such that η =Ṽ 0 /∆ω is varied. This effect could also be achieved through a periodic modulation of the s-wave scattering length as considered in other recent schemes [50,51]. The 1D Fermi-Hubbard model with the asymmetric correlated tunnelling strengths (57) hosts a variety of quantum many-body phases, contrasting the situation of the standard 1D Fermi-Hubbard model, where only insulating and Luttinger-liquid phases occur. We now comment on the phases that could be explored with cold atoms given the constraints imposed by the specific tunnelling rates (57). This model was initially studied for t x ee = t x hh < t x eh , where a groundstate with spin-density-wave order can be found for sufficiently strong repulsion at half-filling [74]. Later on, it was realized that for t x ee = t x hh > t x eh , a new density wave where charge alternates on the bonds (i.e. bond-ordered wave) can be stabilized for not too large interactions [75]. Interestingly, it has been recently shown that, by modifying the filling factor, the phase diagram of the model is much richer even for vanishing Hubbard interactions [51]. For instance, a superconducting phase with unconventional triplet pairing was identified. In addition to the interest of exploring these phases with cold atoms using the above PAT scheme (12)- (13), this could serve to benchmark the proposed quantum simulator, since these results rest on very accurate and efficient analytical and numerical methods that exist in 1D. Moreover, the quantum simulator could be used to study the fate of the predicted phases for different chemical potentials [51], and the appearance of new ones, as the Hubbard repulsion is switched on.
Once the quantum simulator has been verified, it would be very interesting to consider the 2D Fermi-Hubbard model with correlated tunnelings. Although the phase diagram is to the best of our knowledge mostly unknown, the results on the 1D model suggests that it can host a variety of new phases with respect to the standard 2D Fermi-Hubbard model, the understanding of which defies analytical and numerical methods. We believe that the possibility of finding interesting phases of matter, even above the stringent temperatures to observe magnetic ordering in the 2D Fermi-Hubbard model, is certainly worth exploring.
To introduce the topic of the following subsection, we emphasize that it is not possible to setX α = 0 without making X α = 0 simultaneously, as would be required to study solely the effects of bond-charge interactions (53). This also occurs for the schemes based on a periodic modulation of the s-wave scattering length [50,51]. In the following subsection, we shall show that this becomes possible by introducing an additional linear gradient in our scheme.

Bond-charge interactions: hole superconductivity and η-pairing
In this section, we focus on a quantum simulator of the Fermi-Hubbard model with bond-charge interactions, whose importance can be controlled by tuning the ratios of t α /U, X α /U. For simplicity, we focus on the 1D case, and consider the PAT by a spin-independent moving lattice (12) in the presence of an additional linear gradient, which leads to the dressed tunneling in Eq. (15). For the sake of concreteness, we consider a parameter regime fulfilling Eq. (9) for r = 2, and set ϕ = 0. By adjusting the intensity of the moving lattice toṼ 0 = 1.56∆ω (i.e. η = 1.56), such that terms like (55) in the effective Hamiltonian (12) vanishX = 0, we obtain where t = t x J 2 (η ), and the bond-charge interaction X = t x (J 0 (η )−J 2 (η )) alternates between neighboring sites. Other possible values of the moving-lattice intensity fulfillingX α = 0 are η = {4.89, 8.29, 11.53, . . .}, and correspond to solutions of the equation J 0 (η ) + J 4 (η ) − 2J 2 (η ) = 0. Moreover, one can also change the moving-lattice detuning such that r = {4, 6, 8, . . .}, and look for the solutions of J 0 (η ) + J 2r (η ) − 2J r (η ) = 0. This will yield different values of η , and allow for the tunability of the ratio of X/t, and the signs of X,t. This quantum simulator can explore the phenomenon of hole superconductivity [63], where the bond-charge interaction 0 < X t can be responsible for a superconducting phase even in the presence of a repulsive interaction δU ↑↓ > 0, as has been predicted using a mean-field approximation for 2D [70]. Since our effective Hamiltonian (58) can be generalized to higher dimensions following our results in Sec. II, the quantum simulator could test the correctness of such mean-field predictions [63]. From a broader perspective, the quantum simulator can explore the phase diagram of the model in different regimes, such as X > t, and study the effects of the bond-charge alternation in Eq. (58) Since the residual Hubbard interactions in (58) depend on the detuning of the photon-assisted scheme, one can also study the attractive case δU ↑↓ < 0. It has been shown [71] that a bond-charge interaction X = t stabilizes an η-pairing groundstate [72] for finite attraction, which for X = 0 only occurs for strictly infinite interactions δU ↑↓ → ∞ [73]. Although our quantum simulator cannot reach the exact condition of X = t (e.g. for η = 1.56, X/t ≈ 0.94), it has been argued that relaxing these conditions may still host the η-pairing groundstate, even if the exact methods of [71] cannot be applied any longer. It would be very interesting to explore this possibility with the quantum simulator, addressing the role of the finite residual gradient δU σσ in Eq. (58), the bond-charge alternation, and the possibility of achieving a regime X > t that is not feasible for the standard Hubbard model.
Before closing this section, we note that this can also be addressed with hardcore bosons, following our results in Sec. II.

B. High-T c superconductivity and itinerant ferromagnetism
To introduce the topic of this section, let us consider again the Fermi-Hubbard model (2) with additional bond-charge terms (53). The interaction between charges localized in neighboring bonds W α leads to a couple of terms: (i) a pair tunnelling for fermions of opposite pseudospin between neighboring pseudospin excitations, where we have introduced the spin-1/2 operators Regarding the direct exchange, the possibility of observing magnetic ordering of localized pseudospins in the cold-atom scenario would require a direct exchange J de α = 4W α that is much larger than the attainable values (54). The same occurs for transition metals, where the direct exchange cannot explain the appearance of magnetic ordering. As first pointed out by P. W. Anderson [13], magnetic ordering can also arise as a consequence of the strong Hubbard interactions that prevent conduction. Moreover, in the presence of doping, the interplay of this tendency towards magnetic ordering with the correlated dynamics of holes is one of the proposed mechanisms that could explain high-T c superconductivity [16].
In this section, we show that the effective Hamiltonian (18) for the PAT scheme with ϕ = 0 may open a new route to study models of itinerant quantum magnetism and high-T c superconductivity, provided that the limit of large repulsive interactions δU ↑↓ is considered. In this regime, the subspaces of single-occupied H s and doubly-occupied H d lattice sites become wellseparated Hubbard sub-bands. The kinetic energy in Eq. (18) contains terms that act within each of these sub-bands and terms that connect the two sub-bands, such as the operator K s→d : H s → H d expressed as or the operator K d→s : H d → H s , which can be expressed as We shall use this formulation to propose a quantum simulator that explores Nagaoka ferromagnetism even in the limit of finite repulsive interactions, and the full phase diagram of the t-J model in 1D and, more interestingly, in 2D.

Itinerant ferromagnetism: correlated destruction of tunnelling and Nagaoka ferromagnetism
For a strong Coulomb repulsion, electrons in undoped transition metals lower their energy by displaying antiferromagnetic ordering. In the standard Fermi-Hubbard model, the origin of this effect can be traced back to the so-called super-exchange interaction, whereby an antiferromagnetic spin pattern allows for virtual electron tunnelling between neighboring sites that lowers the kinetic energy [13]. The situation is utterly different for ferromagnetism, where the reliability of initial mean-field predictions of large ferromagnetic regions in the phase diagram is highly questionable [76]. One of the few rigorous results on the existence of ferromagnetism in the Fermi-Hubbard model is due to Y. Nagaoka [77], who showed that a single hole in a large class of half-filled Hubbard models can lead to a fully-polarized ferromagnetic groundstate for infinite repulsion.
The stability of this Nagaoka ferromagnet for different regimes has been a topic of recurrent interest in the literature. For two holes [78], the ferromagnet is no longer the groundstate. Nonetheless, for finite hole densities, the fully-polarized Nagaoka ferromagnet can be stable up to a critical hole doping [79]. Despite some initial discrepancy regarding its full polarization [80], more recent numerical results based on quantum Monte Carlo [81] and density-matrix renormalization group [82] agree with the above scenario [79]. Since the task of doping a transition metal with exactly one hole seems quite daunting, these results are crucial for an experimental realization of the Nagaoka ferromagnet. Another obstacle for the realization of a Nagaoka ferromagnet is the requirement of infinite repulsion. In a cold-atom context, alternatives exploiting a long-range double-exchange interaction in a two-band Hubbard model [83], or a large spin-imbalance in an optical lattice with a ladder structure [84] have been considered. We show below that our PAT scheme allows to access the physics of the infinite-repulsion Hubbard model, and thus Nagaoka ferromagnetism, even for finite Hubbard interactions.
Let us consider the effective kinetic energy obtained for the 2D or 3D scheme (18). By looking at the tunnelings in Eqs. (60)-(62), one notices that by modifying the intensity of the moving optical lattice such that J r (η ) = 0, while J 0 (η ) = 0, we can exactly cancel the terms that do not preserve the parity in the occupation number, namely Eqs. (61)- (62). For concreteness, we consider a parameter regime fulfilling (9) for r = 2, such that the moving lattices have the same intensity and detuning, and η =Ṽ 0,α /∆ω α = 5.135. Accordingly, only the parity-conserving tunnelling (60) is preserved. This effect is similar to the so-called coherent destruction of tunnelling [85], where the tunnelling of electrons subjected to a periodically-modulated force is totally suppressed for certain parameters of the force. In our case, the suppressed tunnelling is fully correlated to a particular electron-hole background (see Fig. 2), such that the effect can be understood as a correlated destruction of tunnelling.
For large but finite Hubbard repulsion δU ↑↓ , the Hubbard sub-bands are separated in energies and no term in the effective Hubbard Hamiltonian can connect them. By controlling the atomic filling factor such that n i,↑ + n i,↓ < 1, all the dynamics takes place within the single-occupation subspace, and is controlled by wheret α = t α J 0 (η ), and we have introduced the Gutzwiller projector onto the single-occupied sub-band P s = ∏ i (1 − n i,↑ n i,↓ ). Interestingly enough, the Hamiltonian (63) corresponds to the infinitely-repulsive Hubbard model supporting Nagaoka ferromagnetism. However, in our scheme, only a finite repulsion is required to allow for the adiabatic loading of the lower Hubbard sub-band. By varying the filling factor, the stability of the ferromagnetic phase and the full phase diagram can be explored experimentally. In particular, it can serve as a benchmark of the phase diagram presented [82], which is based on extrapolating numerical results for ladders with increasing number of legs, and has predicted an intermediate phase-separation region between the fully-polarized Nagaoka ferromagnet and the paramagnetic phase.

High-T c superconductivity: tunable t-J and t-XXZ models
In the previous section, we have considered an alternative route to access the physics of the limit of infinite repulsion in the hole-doped Fermi-Hubbard model (i.e. n i,↑ + n i,↓ < 1). However, a large but finite repulsion can also lead to very interesting physics. In this regime, the competition of the projected kinetic energy (63) with the antiferromagnetic super-exchange [13] leads to the so-called t-J model wheret is the tunnelling within the single-occupied subspace, andJ > 0 is the strength of the antiferromagnetic super-exchange interaction. In the case of square lattices, this Hamiltonian (64) has been considered as the canonical effective model in the theory of the high-T c cuprates by part of the scientific community [16,86,87]. In this context, the t-J model arises after mapping a more microscopic three-band Hubbard model of the cuprates [88] onto a single-band one [69,89]. By using the microscopic parameters of the three-band model [90], one finds thatJ/t ∼ 0.3 is the typical regime that can be realized in these materials. For the Fermi-Hubbard model with ultracold atoms (2), one can obtain an effective t-J model in the limit of very strong repulsioñ t U ↑↓ . Then, one findsJ = 4t 2 /U ↑↓ , which cannot attain values larger thanJ < 0.4t since the tunnelling must be at least t < 0.1U ↑↓ to allow for the perturbative process underlying the super-exchange. Therefore, the standard Fermi-Hubbard model is almost at the verge of the regime of importance for the hole-doped cuprates cupratesJ/t ∼ 0.3 [91]. Moreover, from a broader perspective, the rest of the rich phase diagram of the t-J model cannot be explored with these experiments.
We show below that our PAT scheme leads to a t-J model where the ratioJ/t can attain any desired value by controlling the intensity of the moving optical lattice. In this way, the full phase diagram of the t-J model may become accessible to cold-atom quantum simulators.
Let us consider the effective Hamiltonian obtained by the PAT scheme (18) for any dimensionality. We assume equal tunnelings t α =: t, driving parameters r α =: r, and driving ratios η α =: η, in all directions ∀α, and set ϕ α = 0. The limit of strong Hubbard repulsion in this case corresponds to tJ r (η) δU ↑↓ , where the parity-violating tunnelings described by the terms K s→d and K d→s in Eqs. (61)- (62) can only take place virtually. As in the standard Hubbard model [92,93], such virtual tunnelings can be calculated by a Schrieffer-Wolff-type unitary transformationH eff = e iS H eff e −iS , where S = −i(K s→d − K d→s )/δU ↑↓ is responsible for eliminating the energetically forbidden tunnelings to second order in the small expansion parameter ξ = tJ r (η) /δU ↑↓ . Considering the commutation properties of the different operators defined so far, one finds For hole-doping about half-filling, one can project onto the single-occupancy sub-band, which yields the aforementioned t-J model (64) with an additional densitydependent next-nearest-neighbor tunnelling where we have introduced the effective cold-atom parameters and the next-nearest-neighbor vectors u α n (e.g. for 1D u x 1 = 2e x , for 2D {u x 1 = 2e x , u x 2 = e x − e y , u y 1 = 2e y , u y 2 = e x + e y }). Let us note that this additional tunnelling requires that the target site is populated with a hole, as otherwise P s f † i,σ |Ψ = 0. Hence, close to half-filling n i ≈ 1, this term is reduced by a factor (1 − n i )/4 with respect to the Heisenberg super-exchange, and is typically neglected in the literature [86].
As announced above, we have obtained an effective t-J model where the ratio of the coupling constants (66), namelyJ/t = 4ξ J r (η)/J 0 (η), can be tuned by modifying the intensity of the moving optical lattice. Even if ξ ≤ 0.1 in the regime of validity of the t-J model, we can make J r (η) J 0 (η), such that exchange coupling is not required to be much smaller than the tunnelling as in the standard Fermi-Hubbard model. In this way, we can explore the full phase diagram of the t-J model.
In the 1D case, theoretical predictions about the phase diagram are supported by very accurate numerical methods [94]. Such results could serve to benchmark the accuracy of the proposed quantum simulator, which can be prepared in a regime corresponding to a metallic phase (i.e. a repulsive Luttinger liquid [95]), a gapless superconductor (i.e. an attractive Luttinger liquid [95]), or the so-called spin-gap phase [96], which consists of a quantum fluid of bound singlets with gapless density excitations, a gapped spin sector, and enhanced superconducting correlations (i.e. a Luther-Emery liquid [97]). Finally, for sufficiently strong super-exchange interactions, the antiferromagnetic order expels the doped holes leading to separated holerich and hole-poor regions (i.e. phase separation [98]). The richness of this phase diagram highlights the potential of our PAT scheme, and contrasts with the standard 1D Fermi-Hubbard model where only the repulsive Luttinger liquid can be achieved. For instance, the Luther-Emery liquid, which has eluded experimental confirmation so far, requiresJ ≈ 2.5t and thus lies out of the range of parameters that can be obtained from the repulsive Hubbard model. Moreover, our quantum simulator would allow to test the numerical results [99] predicting the disappearance of the phase separation in favor of enlarged superconducting and spin-gap regions, once the next-nearest-neighbor tunnelling terms in Eq. (65) are considered.
In the 2D case, a detailed understanding of certain regions of the phase diagram is still an open problem, and the subject of considerable debate. As emphasized in [100], theoretical predictions are difficult to verify due to (i) the absence of controlled analytical methods, and (ii) the limitation of numerical methods to small system sizes where finite-size effects can affect the predicting power. From this perspective, the proposed quantum simulator may eventually address some of the following open questions regarding the properties of the t-J model. Variational methods based on a resonating-valence-bond trial state [16] (i.e. a linear superposition of all possible configurations of singlet pairs with a weight that depends on the pairing symmetry), have predicted a rich phase diagram [86] with regions of (i) ferromagnetism, (ii) s-wave pairing, (iii) d-wave pairing, (iv) coexistent antiferromagnetism and superconductivity, and (v) phase separation. However, this variational approach introduces a certain bias through the choice of the particular set of ansatzs, and this compromises its reliability leading to considerable controversy in the community [2]. In particular, there are contradictory predictions for low dopings and not too large ratios of J/t, which turns out to be the regime of interest for the high-T c cuprates. For instance, the results of [86,102] contradict those of [98,101], regarding the onset of the phase separation in this low-doping region. There has also been some disagreement regarding the regions of superconductivity predicted by the variational approach [86], exact diagonalization [106], and quantum Monte Carlo [107]. From the perspective of the high-T c cuprates, addressing the conflicting predictions in [103,104] and [86,105] about the existence of stripe phases (i.e. inhomogeneous charge and spin distributions) in the t-J model is even more compelling, as these have been measured experimentally in the cuprates. If the t-J model is to function as a canonical model of the cuprates, as advocated in [16,86,87], it is important to settle this dispute and determine if it admits stripe phases. We believe that the proposed cold-atom experiment could be helpful in this respect.
Before closing this section, let us comment on two additional possibilities for the t-J model quantum simulator. The first, and most obvious one, is the possibility of controlling the spatial anisotropy of the parameters of the effective Hamiltonian (64) by simply exploiting the dependence of the dressed tunnelings on the lattice axes. Accordingly, the effective t-J model becomes anisotropic, such that the anisotropy of the tunnelingst α = t α J 0 (η α ), and the super-exchange couplingsJ α = 4t 2 α J r (η α ) 2 /δU ↑↓ , can be controlled through the intensities of the static and moving optical lattices along the different axes. Such anisotropy becomes especially interesting in the context of certain cuprate ladder compounds [108], which can be modeled by a number ∈ {1, . . . , n } of one-dimensional t-J chains that are coupled to each other by the transverse tunneling t and super-exchange coupling J . This defines the so-called rungs of the t-J ladder Hamiltonian H ladder tJ = P sH ladder tJ P s , wherẽ Such ladder Hamiltonians can be implemented in our quantum simulator if we supplement the above scheme (65) with additional static lattices along the y-axis with commensurate wavelengths with respect to the original lattice. For instance, combining two lattices with doubled wavelengthsλ y = 2λ y , leads to an array of decoupled ladders with n = 2 legs [109]. By adding further harmonics, one may create ladders with other numbers of legs, at least in principle (e.g. the first n harmonics of the Fourier series of the square-wave function yield an approximation to an array of decoupled n -legged ladders). Hence, the interactiondependent PAT scheme leads to Eq. (67) with tunable number of legs and Hamiltonian parameters together with the corresponding density-dependent next-to-nearest-neighbor tunneling (65), typically neglected for small dopings. These t-J ladders provide a very interesting interpolation between the well-understood 1D t-J model, and the more intriguing 2D case full of open questions of relevance to high-T c superconductivity. In fact, already for n = 2 legs, the doped holes tend to pair [110] developing a superconducting d-wave-like order [111]. In addition to the phases that also occur for the 1D t-J model [94], this d-wave superconductivity takes place in a wide region of the phase diagram [112], which includes the parameters relevant for the cuprates. Also, a very interesting even-odd effect reminiscent of the spin-gap presence/absence in the undoped system has been identified [113], whereby the hole d-wave pairing disappears for ladders with an odd numbers of legs. We note that some of these results [110] depend on ratios t /J < 1 that cannot be reached from standard ladder Hubbard models where t J . Likewise, the regime J t,t ,J, and its connection to short-range resonating valence bond states [114,115], cannot be reached from standard Hubbard models. It would be very interesting to test these predictions with our quantum simulator, which allows exploring all these parameter regimes.
Let us move to the last possibility of the t-J model quantum simulator: inducing a Heisenberg-Ising anisotropy in the superexchange interactions. This leads to a very interesting t-XXZ model described by the Hamiltonian H tXXZ = P sHtXXZ P s where the IsingJ z and flip-flopJ ⊥ interaction strengths are generally different, such that the t-J model is recovered whenJ ⊥ =J z . To achieve such an effective model with our quantum simulator, we must employ the PAT by a state-dependent moving optical lattice, which yields the effective Hamiltonian (24) in 1D, and (28) in 2D. We assume equal tunnelings t α =: t, spin-dependent driving parameters r σ ,α =: r σ , and ratios η σ ,α =: η σ , in all directions ∀α, and set ϕ σ ,α = 0. One may observe in Fig. 5(f) that second-order super-exchange interactions will depend on the spin configuration. In fact, we find such that the Heisenberg-Ising anisotropy ζ =J ⊥ /J z can be tuned all the way from small spin quantum fluctuations ζ → 0 (i.e. t-J z model), to the isotropic regime where spin quantum fluctuations play an important role ζ → 1 (i.e. t-J model). These Hamiltonians have been studied in the context of the propagation of a single hole in an antiferromagnetic matrix. For instance, in the t-J z model, the tunneling of the hole leaves behind a string of flipped Ising spins that costs an energy proportional to the string length, such that the holes are almost [116] localized to the site where they were doped [117]. The situation is considerably more complex as the spin fluctuations are switched on ζ > 0, and the t-J model is approached ζ → 1 [118], and some controversy regarding the limitations of the different analytical or numerical methods has been identified [119]. The possibility of controlling the amount of spin fluctuations in our quantum simulator, together with the possibility of using the high-resolution optics of quantum gas microscopes [120,121] to create localized holes and watch them propagate in real time, opens a very nice perspective in accessing this quantum many-body effect with ultracold atoms in optical lattices.
In the half-filled 1D case, it is well-known that the particle-hole excitations of the free tight-binding Hamiltonian can be mapped onto a pair of bosonic branches [95,130], which will play the role of the gauge degrees of freedom where we have introduce the quasi momentum q = 2πn/L for n ∈ Z + , the effective speed of light c ↓ = 2t ↓ eff a = t ↓ eff λ , and the bosonic operators for the particle-hole excitations b q,R (b q,L ) around the right (left) Fermi point. Moreover, we have assumed that the modifications with respect to half-filling coming from the weak parabolic trapping, encoded in ε i,↓ = ε ↓ + 1 2 mω 2 t,x X 2 i , can be accounted for using a local chemical potential δ µ i . Note that, in the continuum limit, this bosonic Hamiltonian becomes a scalar field theory with the energy zero set at ε ↓ , namely where the scalar field φ (x) = φ R (x) − φ L (x) is expressed in terms of the inverse Fourier transform of the bosonic operators b q,R , b q,L and their Hermitian adjoints, and π(x) is its canonically-conjugate momentum. Although there is no Gauge invariance in the massless scalar field theory (74), the dynamical bosonic field will interact with the fermionic matter, and can still lead to interesting phenomena. The particular form of the interaction brings us to the final identification (iii) the gauge field that dresses the tunnelling of the quantum matter in Eq. (71) becomes aA x (x i,i+1 ) = θ (n i+1,↓ − n i,↓ ), which amounts to the density difference of the 'Gauge' species. In order to express the gauge unitaries in terms of the bosonic particle-hole excitations, we use where we have again applied the continuum limit. To be consistent with such a limit, we note that the bare tunnelling of the fermionic quantum matter corresponds to the 1D version [131] of the so-called Kogut-Susskind fermions [124]. At half-filling, one obtains the massless Dirac field theory in the continuum limit, which minimally couples to a derivative of the scalar field where we have introduced the fermionic field operatorsΨ R (x) (Ψ L (x)) for the right (left) moving fermions with pseudospin σ =↑, the energy difference due to the different Zeeman shifts δ = ε ↑ − ε ↓ , and the effective speed of light c ↑ = t ↑ eff λ . We have obtained a peculiar relativistic theory of interacting quantum fields. Instead of the standard Yukawa coupling between scalar and fermionic fields ψφ ψ, we get a minimal coupling with the derivative of the scalar field ψγ 1 ∂ 2 x φ ψ. Moreover, the effective speeds of light of the scalar c ↓ and fermionic c ↑ particles can be tuned independently. Anyhow, scattering between the fermions will occur due to the exchange of scalar particles, such that the cold-atom experiment could be exploited to calculate scattering amplitudes in the spirit of [132]. However, it is not clear how one would create the initial incoming particles and measure the outgoing scattering probabilities in our scheme. A simpler goal would be to study collective properties of the model (76). For instance, for very large δ , fermion-fermion interactions will be mediated by the virtual exchange of scalar particles, such that the fermionic properties of the groundstate will be modified (e.g. correlation functions). Increasing the flux θ could lead to new phases departing from the Luttinger-liquid phase of the free Kogut-Susskind fermions. All these questions could be studied with the proposed quantum simulator.

Correlated topological insulators: Hofstadter-type models with dynamical Gauge fields
Topological insulators represent a family of holographic phases of matter that are insulating in the bulk and conducting at the boundaries [133]. These states are topologically different from trivial band insulators, as the bulk is characterized by a finite topological invariant that cannot be changed unless the bulk energy gap is closed (i.e. phase transition). Another difference occurs at the boundaries, where gapless edge states are responsible for the conductance. In the absence of a boundary energy gap, it is the chirality, or some additional symmetry of the problem, which underlies the robustness of the conductance. This is clearly exemplified by the so-called Hofstadter model [134], which describes fermions in a square lattice subjected to a perpendicular magnetic field, and displays the aforementioned bulk [135] and edge [136] properties. Similar phenomenology also arises in the Haldane model [137], which is a topological insulator in the same symmetry class (i.e. time-reversal symmetry breaking of the integer quantum Hall effect). Remarkably, other instances of topological insulators belonging to different symmetry classes have also been found, such as the Kane-Mele model [138] built from two time-reversed copies of the Haldane model, or the time-reversal Hofstadter model [139] built from two time-reversed copies of the Hofstadter model. and Bose-Fermi Hubbard models with generalized tunnelings whose strength depends on the atomic density through a Bessel function that is controlled by the intensity of the moving lattice, and whose Peierls phase depends on the atomic density through the phase of the moving lattice. We have argued that this effect can be exploited as a flexible tool to implement a variety of quantum simulations of quantum many-body models in the context of strongly-correlated electrons and high-energy physics. In particular, our scheme may allow to explore paradigmatic models, such as the t-J model, in regimes that were previously inaccessible to cold-atom experiments. Moreover, this proposal introduces a new perspective in the realization of dynamical Gauge fields departing from the lattice Gauge theory approach, which can also lead to interesting quantum many-body models.
Although we have focused on cold atoms, the scheme can be applied to other setups, provided that the relevant dynamics can be described by a lattice model, and one can control the periodic driving and the strong interactions. This might be the case of trapped-ion crystals, where the local vibrations and electronic states lead to a lattice model with bosonic and pseudospin degrees of freedom, and the interactions and periodic drivings are provided by their interaction with laser beams.A similar situation arises for superconducting circuits by considering photons in arrays of microwave cavities and superconducting qubits leading to the aforementioned lattice models.