Interaction-dependent photon-assisted tunneling in optical lattices: a quantum simulator of strongly-correlated electrons and dynamical Gauge fields

(i) One-dimensional (1D) scheme: To introduce the main ideas in the simpler setting, let us start by considering a 1D Fermi–Hubbard model obtained from equation (2) for $\{{V}_{0,y},{V}_{0,z}\}\gg {V}_{0,x},$ such that only tunneling along the x-axis is relevant

$$\begin{eqnarray}&&{H}_{{\rm{FH}}}={H}_{{\rm{loc}}}+{H}_{{\rm{kin}}}+{V}_{{\rm{int}}}=\displaystyle \sum _{i,\sigma }{\epsilon }_{i,\sigma }{f}_{i,\sigma }^{\dagger }{f}_{i,\sigma }^{\phantom{\dagger }}-\displaystyle \sum _{i,\sigma }\left({t}_{x}{f}_{i,\sigma }^{\dagger }{f}_{i+1,\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)+\displaystyle \frac{1}{2}\displaystyle \sum _{i,\sigma }{U}_{\sigma \bar{\sigma }}{f}_{i,\sigma }^{\dagger }{f}_{i,\bar{\sigma }}^{\dagger }{f}_{i,\bar{\sigma }}^{\phantom{\dagger }}{f}_{i,\sigma }^{\phantom{\dagger }}.\end{eqnarray} \tag{ 4 }$$

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, figure 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 [58]. 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 ${\tilde{V}}_{0}\ll {V}_{0,x}$ (i.e. ${\tilde{t}}_{x}={\tilde{V}}_{0}{\rm{exp}}\{-\frac{{\pi }^{2}}{4}{({V}_{0,\alpha }/{E}_{{\rm{R}}})}^{1/2}\}\ll {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

$$\begin{eqnarray}&&{H}_{{\rm{mod}}}(t)=\displaystyle \sum _{i,\sigma }\displaystyle \frac{{\tilde{V}}_{0}}{2}\mathrm{cos}\left({\rm{\Delta }}{{kX}}_{i}-{\rm{\Delta }}\omega t+\varphi \right){n}_{i,\sigma }^{\phantom{\dagger }},{n}_{i,\sigma }^{\phantom{\dagger }}={f}_{i,\sigma }^{\dagger }{f}_{i,\sigma }^{\phantom{\dagger }},\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Delta }}k=({{\bf{k}}}_{1}-{{\bf{k}}}_{2})\cdot {{\bf{e}}}_{x}$ is the wavevector difference, ${X}_{i}=\frac{\lambda }{2}i$ stands for the minima of the original optical-lattice potential, ${\rm{\Delta }}\omega ={\omega }_{1}-{\omega }_{2}$ is the detuning of the laser beams, and φ is the relative phase with respect to the static optical lattice. By setting ${\rm{\Delta }}\omega \approx {U}_{\uparrow \downarrow }/r$ for a positive integer $r\in {\mathbb{Z}},$ the above Hubbard blockade for ${U}_{\uparrow \downarrow }\gg {t}_{\alpha }$ can be overcome through the absorption of r photons from the periodic driving (see figure 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 thisinteraction-dependent PAT, we move to the interaction picture with respect to ${U}_{0}(t)={\mathcal{T}}\left(\mathrm{exp}\{{\rm{i}}{\displaystyle \int }_{0}^{t}{\rm{d}}\tau ({V}_{{\rm{int}}}+{H}_{{\rm{mod}}}(\tau ))\}\right),$ such that the fermionic annihilation operators become

$$\begin{eqnarray}&&{U}_{0}(t){f}_{i,\sigma }{U}_{0}^{\dagger }(t)={{\rm{e}}}^{-{\rm{i}}{{tU}}_{\uparrow \downarrow }{n}_{i,\bar{\sigma }}}{{\rm{e}}}^{{\rm{i}}\displaystyle \frac{{\tilde{V}}_{0}}{2{\rm{\Delta }}\omega }\mathrm{sin}({\rm{\Delta }}{{kX}}_{i}-{\rm{\Delta }}\omega t+\varphi )}{f}_{i,\sigma }\\ &&={{\rm{e}}}^{-{\rm{i}}{{tU}}_{\uparrow \downarrow }{n}_{i,\bar{\sigma }}}\displaystyle \sum _{n\in {\mathbb{Z}}}{J}_{n}\left(\displaystyle \frac{\eta }{2}\right){{\rm{e}}}^{{\rm{i}}n({\rm{\Delta }}{{kX}}_{i}-{\rm{\Delta }}\omega t+\varphi )}{f}_{i,\sigma },\end{eqnarray} \tag{ 6 }$$

where we have gauged away an irrelevant phase by transforming the fermion operators⁴ . The second part of the equality is obtained after introducing the important parameter

$$\begin{eqnarray}&&\eta ={\tilde{V}}_{0}/{\rm{\Delta }}\omega ,\end{eqnarray} \tag{ 7 }$$

and using the Jacobi–Anger expansion for first-order Bessel functions ${J}_{n}(z),$ namely ${{\rm{e}}}^{{\rm{i}}z\mathrm{sin}\theta }={\displaystyle \sum }_{n\in {\mathbb{Z}}}{J}_{n}(z){{\rm{e}}}^{{\rm{i}}n\theta }$ [60]. For simplicity, we set ${\rm{\Delta }}{{kX}}_{i}=(\frac{1}{2}{\rm{\Delta }}k\lambda )i=\pi 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 figure 1(a)). However, this could be generalized to ${\rm{\Delta }}{{kX}}_{i}=(\frac{1}{2}{\rm{\Delta }}k\lambda )i=\pi 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 $\alpha =\pi /6,$ see figure 1(b)).

By substituting the expression (6) in the kinetic Hamiltonian ${H}_{{\rm{kin}}}(t)={U}_{0}(t){H}_{{\rm{kin}}}{U}_{0}^{\dagger }(t),$ one finds

$$\begin{eqnarray}&&{H}_{{\rm{kin}}}(t)=-\displaystyle \sum _{i,\sigma }\left({t}_{x,\bar{\sigma }}(t){f}_{i,\sigma }^{\dagger }{f}_{i+1,\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right),{t}_{x,\bar{\sigma }}(t)={t}_{x}{{\rm{e}}}^{-{\rm{i}}{{tU}}_{\uparrow \downarrow }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}f(t),\end{eqnarray} \tag{ 8 }$$

where we have introduced the population difference operator

$$\begin{eqnarray}&&{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}={n}_{i+1,\bar{\sigma }}-{n}_{i,\bar{\sigma }},\end{eqnarray} \tag{ 9 }$$

and a dynamical dressing function

$$\begin{eqnarray}&&f(t)=\displaystyle \sum _{n,m}{J}_{n}\left(\displaystyle \frac{\eta }{2}\right){J}_{m}\left(\displaystyle \frac{\eta }{2}\right){{\rm{e}}}^{-{\rm{i}}(n\pi i-m\pi (i+1))}{{\rm{e}}}^{-{\rm{i}}(n-m)\varphi }{{\rm{e}}}^{{\rm{i}}(n-m){\rm{\Delta }}\omega t}.\end{eqnarray} \tag{ 10 }$$

As announced earlier, the tunneling that connects single-occupied sites to doubly-occupied ones, yielding $\langle {\rm{\Delta }}{n}_{i+1,\bar{\sigma }}\rangle =\pm 1,$ is negligible in the absence of the driving ${\tilde{V}}_{0}=0.$ In this limit, the dressing function is $f(t)=1,$ such that the dressed tunneling can be neglected $\langle {t}_{x,\bar{\sigma }}(t)\rangle ={t}_{x}{{\rm{e}}}^{\mp {\rm{i}}{{tU}}_{\uparrow \downarrow }}\approx 0$ in a rotating-wave approximation for ${t}_{x}\ll {U}_{\uparrow \downarrow }$ (see figure 1(c)). By switching on the periodic driving ${\tilde{V}}_{0}\ne 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 $\pm {U}_{\uparrow \downarrow }=(n-m){\rm{\Delta }}\omega$ (figure 1(d)). In particular, by assuming that

$$\begin{eqnarray}&&{t}_{x},\delta {U}_{\uparrow \downarrow }=({U}_{\uparrow \downarrow }-r{\rm{\Delta }}\omega )\ll {U}_{\uparrow \downarrow }\approx r{\rm{\Delta }}\omega ,\end{eqnarray} \tag{ 11 }$$

we can neglect the majority of tunneling events using a similar rotating-wave argument, except for those that satisfy $n=r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}+m.$ Accordingly, the dressing function becomes simplified

$$\begin{eqnarray}&&f(t)=\displaystyle \sum _{m}{J}_{m}\left(\displaystyle \frac{\eta }{2}\right){J}_{m+r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left(\displaystyle \frac{\eta }{2}\right){{\rm{e}}}^{{\rm{i}}m\pi }{{\rm{e}}}^{-{\rm{i}}\pi r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}i}{{\rm{e}}}^{-{\rm{i}}r\varphi {\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}.\end{eqnarray} \tag{ 12 }$$

We further assume that the laser detuning is chosen in such a way that r is an even integer, and set ${{\rm{e}}}^{-{\rm{i}}\pi r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}i}=1$ for any population difference. Making use of the Neumann–Graf addition formula for Bessel functions [60], namely ${\displaystyle \sum }_{n\in {\mathbb{Z}}}{J}_{n}(z){J}_{n+\nu }(z){{\rm{e}}}^{{\rm{i}}n\theta }={J}_{\nu }\left(2| z\mathrm{sin}(\theta /2)| \right){{\rm{e}}}^{{\rm{i}}(\pi -\theta )\nu /2},$ we can express the PAT in terms of a single Bessel function

$$\begin{eqnarray}&&{t}_{x,\bar{\sigma }}(t)={t}_{x}{{\rm{e}}}^{-{\rm{i}}t\delta {U}_{\uparrow \downarrow }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}{J}_{r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left(\eta \right){{\rm{e}}}^{-{\rm{i}}r\varphi {\rm{\Delta }}{n}_{i+1,\bar{\sigma }}},\end{eqnarray} \tag{ 13 }$$

which should be understood in terms of its Taylor series expansion.

The total time-evolution operator $U(t)={U}_{0}^{\dagger }(t){{\rm{e}}}^{-{\rm{i}}t{\displaystyle \sum }_{i}\delta {U}_{\uparrow \downarrow }{n}_{i\uparrow }{n}_{i\downarrow }}{{\rm{e}}}^{-{\rm{i}}{H}_{{\rm{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 equation (11), such that

$$\begin{eqnarray}&&{H}_{{\rm{eff}}}=\displaystyle \sum _{i,\sigma }{\epsilon }_{i,\sigma }{f}_{i,\sigma }^{\dagger }{f}_{i,\sigma }^{\phantom{\dagger }}-\displaystyle \sum _{i,\sigma }\left({t}_{x}{J}_{r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left(\eta \right){{\rm{e}}}^{-{\rm{i}}r\varphi {\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}{f}_{i,\sigma }^{\dagger }{f}_{i+1,\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\&&+\displaystyle \frac{1}{2}\displaystyle \sum _{i,\sigma }\delta {U}_{\sigma \bar{\sigma }}{f}_{i,\sigma }^{\dagger }{f}_{i,\bar{\sigma }}^{\dagger }{f}_{i,\bar{\sigma }}^{\phantom{\dagger }}{f}_{i,\sigma }^{\phantom{\dagger }}.\end{eqnarray} \tag{ 14 }$$

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

$$\begin{eqnarray}&&{J}_{r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left(\eta \right)={J}_{0}\left(\eta \right){h}_{i,\bar{\sigma }}{h}_{i+1,\bar{\sigma }}+{J}_{0}\left(\eta \right){n}_{i,\bar{\sigma }}{n}_{i+1,\bar{\sigma }}+{J}_{r}\left(\eta \right){n}_{i,\bar{\sigma }}{h}_{i+1,\bar{\sigma }}+{J}_{r}\left(\eta \right){h}_{i,\bar{\sigma }}{n}_{i+1,\bar{\sigma }}.\end{eqnarray} \tag{ 15 }$$

where we have defined the hole number operators ${h}_{i,\bar{\sigma }}=1-{n}_{i,\bar{\sigma }}.$ The first term in equation (15) describes the tunneling within the subspace of single-occupied sites ${{\mathcal{H}}}_{{\rm{s}}},$ whereas the second one corresponds to tunneling within the subspace of doubly-occupied sites ${{\mathcal{H}}}_{{\rm{d}}}$ (see figure 1(e)). These subspaces can be described as two Hubbard sub-bands centered around ${\epsilon }_{{\rm{s}}}=0$ and ${\epsilon }_{{\rm{d}}}=\delta {U}_{\uparrow \downarrow }.$ Finally, the third and fourth terms stand for tunneling events connecting the single-occupied to the doubly-occupied subspaces (see figure 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}_{{\rm{kin}}}\;\to \;$ ${\displaystyle \sum }_{i,\sigma }\left({n}_{i,\bar{\sigma }}+{h}_{i,\bar{\sigma }}\right){H}_{{\rm{kin}}}^{\sigma ,i}({n}_{i,\bar{\sigma }}+{h}_{i,\bar{\sigma }}).$ Let us emphasize, however, that the ratio of these diagonal/off-diagonal processes cannot be controlled, which contrasts the PAT Hamiltonians (15), where one can adjust the intensity of the moving optical lattice ${\tilde{V}}_{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 section 3.2. In section 3.1, 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 section 3.3, 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 (14), (15), and the periodically driven one (4), (5) in the simplest setting: a Fermi–Hubbard dimer (see figures 2(a), (e)). In figure 2, we explore the real-time dynamics for different configurations of atoms in the initial state. (i) Pauli-blockaded regime: figures 2(b), (c) represent the dynamics of the initial atomic configuration $| {\downarrow }_{1},\uparrow {\downarrow }_{2}\rangle ={f}_{1\downarrow }^{\dagger }{f}_{2\uparrow }^{\dagger }{f}_{2\downarrow }^{\dagger }| 0\rangle ,$ 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 figure 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 (14), (15). In particular, it shows that provided the constrains in equation (11) are carefully fulfilled by the system paramaters, terms beyond the rotating-wave approximation leading to equation (12), within the single-band approximation, do not lead to additional errors departing from the desired target Hamiltonian evolution. Regarding the dynamics of the down-spin atoms, we note that these cannot tunnel due to the Pauli exclusion principle, as depicted in figure 2(c). (ii) Hubbard-blockaded regime: figures 2(f), (g) represent the dynamics of the initial atomic configurations $| {0}_{1},\uparrow {\downarrow }_{2}\rangle ={f}_{2\uparrow }^{\dagger }{f}_{2\downarrow }^{\dagger }| 0\rangle ,$ 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 figure 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.

Zoom In Zoom Out Reset image size

Let us now address the phenomenon of correlated destruction of tunneling exploited in section 3 for the quantum simulation of strongly-correlated models. According to equation (15), 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 figure 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 figure 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 section 3.

We have so far presented numerical tests supporting the validity of the resonant PAT, ${\rm{\Delta }}\omega ={U}_{\uparrow \downarrow }/r,$ such that the residual interactions of the dressed Fermi–Hubbard model (14) vanish $\delta {U}_{\uparrow \downarrow }=0.$ However, our analytical results show that finite Hubbard interactions $\delta {U}_{\uparrow \downarrow }=({U}_{\uparrow \downarrow }-r{\rm{\Delta }}\omega )$ can be achieved by changing the velocity of the moving optical lattice ${\rm{\Delta }}\omega \ne {U}_{\uparrow \downarrow }/r,$ which will be crucial for several quantum simulations in section 3. Let us test this result by numerically integrating an adiabatic evolution according to the effective (14), (15) and periodically-driven (4), (5) Hamiltonians for a half-filled dimer (figure 3(a)). We study the evolution of the system (see figure 3(b)), for a slow ramp of the tunneling strength ${t}_{x}\to {t}_{x}(1-(\delta {t}_{x})t)$ with a rate $\delta {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}\gg \delta {U}_{\uparrow \downarrow }.$ After the quench ${t}_{{\rm{f}}}\approx 1/\delta {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}\ll \delta {U}_{\uparrow \downarrow }$

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{0}\rangle \approx | {\rm{FS}}\rangle =\displaystyle \frac{1}{2}\left({f}_{1\uparrow }^{\dagger }+{f}_{2\uparrow }^{\dagger }\right)\left({f}_{1\downarrow }^{\dagger }+{f}_{2\downarrow }^{\dagger }\right)| 0\rangle \longrightarrow | {{\rm{\Psi }}}_{{\rm{f}}}\rangle \approx | {\rm{HS}}\rangle =\displaystyle \frac{1}{\sqrt{2}}\left({f}_{1\uparrow }^{\dagger }{f}_{2\downarrow }^{\dagger }-{f}_{1\downarrow }^{\dagger }{f}_{2\uparrow }^{\dagger }\right)| 0\rangle ,\end{eqnarray} \tag{ 16 }$$

In figure 3(c), we represent the numerical results for the Heisenberg-singlet fidelity ${{\mathcal{F}}}_{{\rm{s}}}(t)={| \langle {\rm{HS}}| {\rm{\Psi }}(t)\rangle | }^{2}$ as a function of the ramp time, and for different ramp rates. We observe that the fidelity approaches ${{\mathcal{F}}}_{{\rm{s}}}(1/\delta {t}_{x})\approx 1$ for the very slow ramps, where the adiabatic evolution is expected to be more accurate. Once again, the good agreement between the effective (14), (15) 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.

Zoom In Zoom Out Reset image size

At this point, it is worth commenting on the effect of higher excited bands that would be present in the optical-lattice setup, but are not contained in the single-band approximation implicit to equation (2), and the rest of our treatment. The periodic modulation may also assist inter-band transition by multi-photon resonances where $n{\rm{\Delta }}\omega ={\rm{\Delta }}E,$ where $n\in {\mathbb{Z}}$ and Δ E is the energy gap between the lowest and some higher band. To avoid such processes, one must ensure that these resonances are avoided for the lowest-lying bands where the number of absorbed photons n can be the lowest. Eventually, this parameter choice may lead to a resonance with a much higher band for $n\ll 1,$ but provided that ${\rm{\Delta }}E\ll {\tilde{V}}_{0},$ the population transferred will be exponentially slower than the inter-band tunneling t_inter-band$\sim {t}_{x}{({\rm{\Delta }}E/{\tilde{V}}_{0})}^{n}\ll {t}_{x}.$ One must thus make sure that these inter-band population transfer is much slower than the time required for the experiment, which becomes essential for the cases that deal with the slower super exchange (16). In addition to these possible errors in the simulation, one also has to consider the effect of working in a different periodically-modulated picture, which can lead to micro-motion contributions at the driving frequency that can alter the experimental measurements [45].

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 magnetic-field gradient. If we tune the gradient such that it coincides with the on-site interaction, we can generalize equation (14) by substituting

$$\begin{eqnarray}&&{J}_{r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left(\eta \right)\to {J}_{r(1+{\rm{\Delta }}{n}_{i+1,\bar{\sigma }})}\left(\eta \right)={J}_{r}\left(\eta \right){h}_{i,\bar{\sigma }}{h}_{i+1,\bar{\sigma }}+{J}_{r}\left(\eta \right){n}_{i,\bar{\sigma }}{n}_{i+1,\bar{\sigma }}+{J}_{0}\left(\eta \right){n}_{i,\bar{\sigma }}{h}_{i+1,\bar{\sigma }} \\&&+{J}_{2r}\left(\eta \right){h}_{i,\bar{\sigma }}{n}_{i+1,\bar{\sigma }}.\end{eqnarray} \tag{ 17 }$$

According to this expression, the off-diagonal tunnelings connecting doubly- to single-occupied sites (i.e. fourth term in equation (17) represented in the upper panel of figure 4(a)), and single- to doubly-occupied sites (i.e. third term in equation (17) represented in the upper panel of figure 4(c)), depend on different Bessel functions. This leads to a two-tone beating in the tunneling dynamics, as shown in figures 4(b), (d), which also serve as tests of the validity of our analytical derivations.

Zoom In Zoom Out Reset image size

(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 (${V}_{0,z}\gg \{{V}_{0,x},{V}_{0,y}\}\gg {E}_{{\rm{R}}}$) or three ($\{{V}_{0,x},{V}_{0,y},{V}_{0,z}\}\gg {E}_{{\rm{R}}}$) directions. As can be observed from equations (6)–(10), 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 $\alpha =\{x,y\},$ or three $\alpha =\{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 1D 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

$$\begin{eqnarray}&&{H}_{{\rm{mod}}}(t)=\displaystyle \sum _{{\bf{i}},\sigma }\displaystyle \sum _{\alpha }\displaystyle \frac{{\tilde{V}}_{0,\alpha }}{2}\mathrm{cos}\left({\rm{\Delta }}{k}_{\alpha }{R}_{{\bf{i}},\alpha }-{\rm{\Delta }}{\omega }_{\alpha }t+{\varphi }_{\alpha }\right){n}_{{\bf{i}},\sigma }^{\phantom{\dagger }},\end{eqnarray} \tag{ 18 }$$

where we have introduced the labeling indexes ${\bf{i}}=({i}_{x},{i}_{y})$ for 2D, and ${\bf{i}}=({i}_{x},{i}_{y},{i}_{z})$ for 3D. As before, we have assumed that for 2D (${\tilde{V}}_{0,x}\ll {V}_{0,x},$ and ${\tilde{V}}_{0,y}\ll {V}_{0,y}$), and for 3D (${\tilde{V}}_{0,x}\ll {V}_{0,x},$ ${\tilde{V}}_{0,y}\ll {V}_{0,y},$ and ${\tilde{V}}_{0,z}\ll {V}_{0,z}$), such that the moving lattices do not 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 (14), which requires a parameter regime

$$\begin{eqnarray}&&{t}_{x},{t}_{y},{t}_{z},\delta {U}_{\uparrow \downarrow }=({U}_{\uparrow \downarrow }-{r}_{\alpha }{\rm{\Delta }}{\omega }_{\alpha })\ll {U}_{\uparrow \downarrow }\approx {r}_{x}{\rm{\Delta }}{\omega }_{x}={r}_{y}{\rm{\Delta }}{\omega }_{y}={r}_{z}{\rm{\Delta }}{\omega }_{z},\end{eqnarray} \tag{ 19 }$$

and yields the following effective Hamiltonian

$$\begin{eqnarray}{H}_{{\rm{eff}}} & = & \displaystyle \sum _{{\bf{i}},\sigma }{\epsilon }_{{\bf{i}},\sigma }{f}_{{\bf{i}},\sigma }^{\dagger }{f}_{{\bf{i}},\sigma }^{\phantom{\dagger }}-\displaystyle \sum _{{\bf{i}},\alpha }\displaystyle \sum _{\sigma }\left({t}_{\alpha }{J}_{{r}_{\alpha }{\rm{\Delta }}{n}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\bar{\sigma }}}\left({\eta }_{\alpha }\right){{\rm{e}}}^{-{\rm{i}}{r}_{\alpha }{\varphi }_{\alpha }{\rm{\Delta }}{n}_{{\bf{i}}+{{\bf{e}}}_{\alpha },\bar{\sigma }}}{f}_{{\bf{i}},\sigma }^{\dagger }{f}_{{\bf{i}},+{{\bf{e}}}_{\alpha },\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{{\bf{i}},\sigma }\delta {U}_{\sigma \bar{\sigma }}{f}_{{\bf{i}},\sigma }^{\dagger }{f}_{{\bf{i}},\bar{\sigma }}^{\dagger }{f}_{{\bf{i}},\bar{\sigma }}^{\phantom{\dagger }}{f}_{{\bf{i}},\sigma }^{\phantom{\dagger }},\end{eqnarray} \tag{ 20 }$$

where ${\rm{\Delta }}{n}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\sigma }={n}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\sigma }-{n}_{{\bf{i}},\sigma },$ the bare tunnelings are approximated by equation (3), and the dressed tunnelings depend on

$$\begin{eqnarray}{J}_{{r}_{\alpha }{\rm{\Delta }}{n}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\bar{\sigma }}}\left({\eta }_{\alpha }\right) & = & {J}_{0}\left({\eta }_{\alpha }\right){h}_{{\bf{i}},\bar{\sigma }}{h}_{{\bf{i}}+{{\bf{e}}}_{\alpha },\bar{\sigma }}+{J}_{0}\left({\eta }_{\alpha }\right){n}_{{\bf{i}},\bar{\sigma }}{n}_{{\bf{i}}+{{\bf{e}}}_{\alpha },\bar{\sigma }}\\ & & +{J}_{{r}_{\alpha }}\left({\eta }_{\alpha }\right){n}_{{\bf{i}},\bar{\sigma }}{h}_{{\bf{i}}+{{\bf{e}}}_{\alpha },\bar{\sigma }}+{J}_{{r}_{\alpha }}\left({\eta }_{\alpha }\right){h}_{{\bf{i}},\bar{\sigma }}{n}_{{\bf{i}}+{{\bf{e}}}_{\alpha },\bar{\sigma }},\end{eqnarray} \tag{ 21 }$$

where we have introduced ${\eta }_{\alpha }={\tilde{V}}_{0,\alpha }/{\rm{\Delta }}{\omega }_{\alpha }.$ 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}_{{\rm{c}}}$ cuprate superconductors, as outlined in section 3.2.

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.

(i) 1D 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 figure 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 [58]. 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 [33]. 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 (15) 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.

Zoom In Zoom Out Reset image size

In this case, we can generalize the driving (5) by including a state-dependent periodic modulation of the on-site energies

$$\begin{eqnarray}&&{H}_{{\rm{mod}}}(t)=\displaystyle \sum _{i,\sigma }\displaystyle \frac{{\tilde{V}}_{0,\sigma }}{2}\mathrm{cos}\left({\rm{\Delta }}{{kX}}_{i}-{\rm{\Delta }}\omega t+{\varphi }_{\sigma }\right){n}_{i,\sigma }^{\phantom{\dagger }},\end{eqnarray} \tag{ 22 }$$

where the spin-dependent amplitude must again fulfill ${\tilde{V}}_{0,\sigma }\ll {V}_{0,x},$ and φ_σ stands for a phase difference with respect to the static lattice that is generally state dependent [61, 62] (see figures 5 (c), (d)). Going back to the spin-independent scheme (5), the new modulation (22) can be achieved by rotating the polarization of the laser beam that is slightly detuned with respect to the static-lattice lasers [61]. In this case, the spin-dependent driving amplitudes ${\tilde{V}}_{0,\sigma }$ can be tuned by controlling such an angle, or instead the direction of propagation of the laser beam with respect to the quantization axis [62]. 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. ${\tilde{V}}_{0,\uparrow }=0,$ ${\tilde{V}}_{0,\downarrow }\ne 0$), or vice versa, as realized in ion-trap experiments [63]. This leads to a spin-dependent driving where the wavevector Δ k_σ, detuning Δω_σ, intensity ${\tilde{V}}_{0,\sigma },$ and relative phase φ_σ can all be controlled independently for each pseudospin

$$\begin{eqnarray}&&{H}_{{\rm{mod}}}(t)=\displaystyle \sum _{i,\sigma }\displaystyle \frac{{\tilde{V}}_{0,\sigma }}{2}\mathrm{cos}\left({\rm{\Delta }}{k}_{\sigma }{X}_{i}-{\rm{\Delta }}{\omega }_{\sigma }t+{\varphi }_{\sigma }\right){n}_{i,\sigma }^{\phantom{\dagger }}.\end{eqnarray} \tag{ 23 }$$

Paralleling the previous section, we will consider equal wavelengths of the static and moving optical lattices, such that ${\rm{\Delta }}{k}_{\sigma }{X}_{i}=(\frac{1}{2}{\rm{\Delta }}{k}_{\sigma }\lambda )i=\pi i,$ although we remark again that the scheme also works for other propagation angles. Once the new periodic drivings (22), (23) have been discussed, we can address the interaction-dependent PAT they give rise to. We shall use equation (23), as the results also encompass those related to the driving (22). By reproducing the steps that lead to the effective Hamiltonian (14) for the spin-independent driving, we find a parameter regime analogous to equation (11), namely

$$\begin{eqnarray}&&{t}_{x},\delta {U}_{\uparrow \downarrow }=({U}_{\uparrow \downarrow }-{r}_{\sigma }{\rm{\Delta }}{\omega }_{\sigma })\ll {U}_{\uparrow \downarrow }\approx {r}_{\uparrow }{\rm{\Delta }}{\omega }_{\uparrow }={r}_{\downarrow }{\rm{\Delta }}{\omega }_{\downarrow },\end{eqnarray} \tag{ 24 }$$

and a new dressing function of the tunneling that becomes spin-dependent, namely

$$\begin{eqnarray}&&{f}_{\sigma }(t)=\displaystyle \sum _{m}{J}_{m}\left(\displaystyle \frac{{\eta }_{\sigma }}{2}\right){J}_{m+{r}_{\sigma }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left(\displaystyle \frac{{\eta }_{\sigma }}{2}\right){{\rm{e}}}^{{\rm{i}}m\pi }{{\rm{e}}}^{-{\rm{i}}\pi {r}_{\sigma }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}i}{{\rm{e}}}^{-{\rm{i}}{r}_{\sigma }{\varphi }_{\sigma }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}.\end{eqnarray} \tag{ 25 }$$

where ${\eta }_{\sigma }={\tilde{V}}_{0}^{\sigma }/{\rm{\Delta }}{\omega }_{\sigma }.$ In this case, the detunings are chosen such that r_σ is an even integer for both pseudospins, such that we can thus set ${{\rm{e}}}^{-{\rm{i}}\pi {r}_{\sigma }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}i}=1.$ Using the Neuman-Graff addition formula once again, we find that

$$\begin{eqnarray}&&{H}_{{\rm{eff}}}=\displaystyle \sum _{i,\sigma }{\epsilon }_{i,\sigma }{f}_{i,\sigma }^{\dagger }{f}_{i,\sigma }^{\phantom{\dagger }}-\displaystyle \sum _{i,\sigma }\left({t}_{x}{J}_{{r}_{\sigma }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left({\eta }_{\sigma }\right){{\rm{e}}}^{-{\rm{i}}{r}_{\sigma }{\varphi }_{\sigma }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}{f}_{i,\sigma }^{\dagger }{f}_{i+1,\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\&&+\displaystyle \frac{1}{2}\displaystyle \sum _{i,\sigma }\delta {U}_{\sigma \bar{\sigma }}{f}_{i,\sigma }^{\dagger }{f}_{i,\bar{\sigma }}^{\dagger }{f}_{i,\bar{\sigma }}^{\phantom{\dagger }}{f}_{i,\sigma }^{\phantom{\dagger }},\end{eqnarray} \tag{ 26 }$$

Remarkably, we find that the amplitude of the density-dependent tunneling can be controlled independently for each pseudospin

$$\begin{eqnarray}&&{J}_{{r}_{\sigma }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left({\eta }_{\sigma }\right)={J}_{0}\left({\eta }_{\sigma }\right){h}_{i,\bar{\sigma }}{h}_{i+1,\bar{\sigma }}+{J}_{0}\left({\eta }_{\sigma }\right){n}_{i,\bar{\sigma }}{n}_{i+1,\bar{\sigma }}+{J}_{{r}_{\sigma }}\left({\eta }_{\sigma }\right){n}_{i,\bar{\sigma }}{h}_{i+1,\bar{\sigma }}+{J}_{{r}_{\sigma }}\left({\eta }_{\sigma }\right){h}_{i,\bar{\sigma }}{n}_{i+1,\bar{\sigma }},\end{eqnarray} \tag{ 27 }$$

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 section 3.3. 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 equation (27), the dressed tunneling of a doubly-occupied half-filled dimer (see figure 6(a)) depends on the spin of the atom, such that the spin-up atoms tunneling is controlled by the Bessel function ${J}_{{r}_{\uparrow }}({\eta }_{\uparrow }),$ whereas the spin-down atoms tunneling depends on the Bessel function ${J}_{{r}_{\downarrow }}({\eta }_{\downarrow }).$ 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}_{\downarrow }}({\eta }_{\downarrow }^{\star })=0,$ while the spin-up atoms hop freely in the lattice ${J}_{{r}_{\uparrow }}({\eta }_{\uparrow }^{\star })\ne 0$ (see figures 6(a), (b)). Conversely, we can coherently freeze the spin-up atoms ${J}_{{r}_{\uparrow }}({\eta }_{\uparrow }^{\star })=0,$ while the spin-down atoms hop freely in the lattice ${J}_{{r}_{\downarrow }}({\eta }_{\downarrow }^{\star })\ne 0$ (see figures 6(c), (d)). Let us emphasize the excellent agreement between our analytical description (symbols), and the exact dynamics of the periodic Hamiltonian (solid lines).

Zoom In Zoom Out Reset image size

(ii) Two-dimensional scheme: The scheme presented above can be generalized beyond 1D. For the sake of concreteness, and for its particular interest in connection to the t-J model in section 3.2, and the dynamical Gauge fields in section 3.3, we will restrict to 2D (${V}_{0,z}\gg \{{V}_{0,x},{V}_{0,y}\}\gg {E}_{{\rm{R}}}$). The idea is to consider spin-dependent moving optical lattices along the x- and y-axes

$$\begin{eqnarray}&&{H}_{{\rm{mod}}}(t)=\displaystyle \sum _{{\bf{i}},\sigma }\displaystyle \sum _{\alpha }\displaystyle \frac{{\tilde{V}}_{0,\alpha }^{\sigma }}{2}\mathrm{cos}\left({\rm{\Delta }}{k}_{\sigma ,\alpha }{r}_{{\bf{i}},\alpha }^{0}-{\rm{\Delta }}{\omega }_{\sigma ,\alpha }t+{\varphi }_{\sigma ,\alpha }\right){n}_{{\bf{i}},\sigma }^{\phantom{\dagger }},\end{eqnarray} \tag{ 28 }$$

such that the relative phases of the moving lattices fulfill ${\varphi }_{\sigma ,x}\geqslant 0,$ but ${\varphi }_{\sigma ,y}=0.$ In this case, and after following the same steps as above in an analogous parameter regime

$$\begin{eqnarray}&&{t}_{x},{t}_{y},\delta {U}_{\uparrow \downarrow }=({U}_{\uparrow \downarrow }-{r}_{\sigma ,\alpha }{\rm{\Delta }}{\omega }_{\sigma ,\alpha })\ll {U}_{\uparrow \downarrow }\approx {r}_{\sigma ,\alpha }{\rm{\Delta }}{\omega }_{\sigma ,\alpha }={r}_{{\sigma }^{\prime },{\alpha }^{\prime }}{\rm{\Delta }}{\omega }_{{\sigma }^{\prime },{\alpha }^{\prime }},\forall \alpha ,{\alpha }^{\prime },\sigma ,{\sigma }^{\prime }\end{eqnarray} \tag{ 29 }$$

one can derive the following effective Hamiltonian

$$\begin{eqnarray}{H}_{{\rm{eff}}} & = & -\displaystyle \sum _{{\bf{i}}},\sigma \left({t}_{x}{J}_{{r}_{\sigma ,x}{\rm{\Delta }}{n}_{{\bf{i}}+{{\bf{e}}}_{x},\bar{\sigma }}}\left({\eta }_{\sigma ,x}\right){{\rm{e}}}^{-{\rm{i}}{r}_{\sigma ,x}{\varphi }_{\sigma ,x}{\rm{\Delta }}{n}_{{\bf{i}}+{{\bf{e}}}_{x},\bar{\sigma }}}{f}_{{\bf{i}},\sigma }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{x},\sigma }^{\phantom{\dagger }}\right.\\ & & +\left.{t}_{y}{J}_{{r}_{\sigma ,y}{\rm{\Delta }}{n}_{{\bf{i}}+{{\bf{e}}}_{y},\bar{\sigma }}}\left({\eta }_{\sigma ,y}\right){f}_{{\bf{i}},\sigma }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{y},\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)+\displaystyle \sum _{{\bf{i}}}\delta {U}_{\uparrow \downarrow }{n}_{{\bf{i}},\uparrow }{n}_{{\bf{i}},\downarrow },\end{eqnarray} \tag{ 30 }$$

where ${\eta }_{\sigma ,\alpha }={\tilde{V}}_{0,\alpha }^{\sigma }/{\rm{\Delta }}{\omega }_{\sigma ,\alpha }.$ 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 section 3.3, 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 equation (30), 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 figures 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 figure 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 (30) 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 section 3.3. At this point, the reader may skip the following sections, and move directly to the quantum simulator applications dealing with fermionic models in section 3.

Zoom In Zoom Out Reset image size

Let us start from the 1D Bose–Hubbard model obtained from equation (49) for $\{{V}_{0,y},{V}_{0,z}\}\gg {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}_{\uparrow \uparrow },{U}_{\downarrow \downarrow }\gg {U}_{\uparrow \downarrow }\gg {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 $\langle {n}_{i,\sigma }\rangle \leqslant 1.$ By using the corresponding projector ${{\mathcal{P}}}_{{\rm{s}}},$ we can map the bosonic creation-annihilation operators onto an ${\mathfrak{su}}(2)$ spin algebra

$$\begin{eqnarray}&&{{\mathcal{P}}}_{{\rm{s}}}{b}_{i,\sigma }{{\mathcal{P}}}_{{\rm{s}}}=| {0}_{i,\sigma }\rangle \langle {1}_{i,\sigma }| =:\;{\tilde{b}}_{i,\sigma }^{\phantom{\dagger }},{{\mathcal{P}}}_{{\rm{s}}}{b}_{i,\sigma }^{\dagger }{{\mathcal{P}}}_{{\rm{s}}}=| {1}_{i,\sigma }\rangle \langle {0}_{i,\sigma }| \\&&=:\;{\tilde{b}}_{i,\sigma }^{\dagger },{{\mathcal{P}}}_{{\rm{s}}}{b}_{i,\sigma }^{\dagger }{b}_{i,\sigma }^{\phantom{\dagger }}{{\mathcal{P}}}_{{\rm{s}}}=| {1}_{i,\sigma }\rangle \langle {1}_{i,\sigma }| =:\;{\tilde{n}}_{i\sigma }.\end{eqnarray} \tag{ 32 }$$

In such a hardcore limit, the Bose–Hubbard model (49) only contains the on-site Hubbard interactions for two bosons of opposite pseudospin, namely

$$\begin{eqnarray}{\tilde{H}}_{{\rm{hBH}}} & = & {\tilde{H}}_{{\rm{loc}}}+{\tilde{H}}_{{\rm{kin}}}+{\tilde{V}}_{{\rm{int}}}:= {{\mathcal{P}}}_{{\rm{s}}}{H}_{{\rm{HB}}}{{\mathcal{P}}}_{{\rm{s}}}=\displaystyle \sum _{i,\sigma }{\epsilon }_{i,\sigma }{\tilde{b}}_{i,\sigma }^{\dagger }{\tilde{b}}_{i,\sigma }^{\phantom{\dagger }}\\ & & -\displaystyle \sum _{i,\sigma }\left({t}_{x}{\tilde{b}}_{i,\sigma }^{\dagger }{\tilde{b}}_{i+1,\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)+\displaystyle \frac{1}{2}\displaystyle \sum _{i,\sigma }{U}_{\sigma \bar{\sigma }}{\tilde{b}}_{i,\sigma }^{\dagger }{\tilde{b}}_{i,\bar{\sigma }}^{\dagger }{\tilde{b}}_{i,\bar{\sigma }}^{\phantom{\dagger }}{\tilde{b}}_{i,\sigma }^{\phantom{\dagger }}.\end{eqnarray} \tag{ 33 }$$

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

$$\begin{eqnarray}&&{\tilde{H}}_{{\rm{mod}}}(t):= {{\mathcal{P}}}_{{\rm{s}}}{H}_{{\rm{mod}}}(t){{\mathcal{P}}}_{{\rm{s}}}=\displaystyle \sum _{i,\sigma }\displaystyle \frac{{\tilde{V}}_{0}}{2}\mathrm{cos}\left({\rm{\Delta }}{{kX}}_{i}-{\rm{\Delta }}\omega t+\varphi \right){\tilde{n}}_{i,\sigma }^{\phantom{\dagger }}\end{eqnarray} \tag{ 34 }$$

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)={\mathcal{T}}\left(\mathrm{exp}\{{\rm{i}}{\displaystyle \int }_{0}^{t}{\rm{d}}\tau ({\tilde{V}}_{{\rm{int}}}+{\tilde{H}}_{{\rm{mod}}}(\tau ))\}\right),$ and follow the same steps as in the fermionic case to express the time-evolution operator $U(t)={U}_{0}^{\dagger }(t){{\rm{e}}}^{-{\rm{i}}t\displaystyle \sum _{i}\delta {U}_{\uparrow \downarrow }{\tilde{n}}_{i\uparrow }{\tilde{n}}_{i\downarrow }}{{\rm{e}}}^{-{\rm{i}}{\tilde{H}}_{{\rm{eff}}}t},$ in terms of

$$\begin{eqnarray}&&{\tilde{H}}_{{\rm{eff}}}=\displaystyle \sum _{i,\sigma }{\epsilon }_{i,\sigma }{\tilde{b}}_{i,\sigma }^{\dagger }{\tilde{b}}_{i,\sigma }^{\phantom{\dagger }}-\displaystyle \sum _{i,\sigma }\left({t}_{x}{J}_{r{\rm{\Delta }}{\tilde{n}}_{i+1,\bar{\sigma }}}\left(\eta \right){{\rm{e}}}^{-{\rm{i}}r\varphi {\rm{\Delta }}{\tilde{n}}_{i+1,\bar{\sigma }}}{\tilde{b}}_{i,\sigma }^{\dagger }{\tilde{b}}_{i+1,\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\&&+\displaystyle \frac{1}{2}\displaystyle \sum _{i,\sigma }\delta {U}_{\sigma \bar{\sigma }}{\tilde{b}}_{i,\sigma }^{\dagger }{\tilde{b}}_{i,\bar{\sigma }}^{\dagger }{\tilde{b}}_{i,\bar{\sigma }}^{\phantom{\dagger }}{\tilde{b}}_{i,\sigma }^{\phantom{\dagger }},\end{eqnarray} \tag{ 35 }$$

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 figure 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 density-dependent tunnelling strength in complete analogy to the fermionic case

$$\begin{eqnarray}&&{J}_{r{\rm{\Delta }}{\tilde{n}}_{i+1,\bar{\sigma }}}\left(\eta \right)={J}_{0}\left(\eta \right){\tilde{h}}_{i,\bar{\sigma }}{\tilde{h}}_{i+1,\bar{\sigma }}+{J}_{0}\left(\eta \right){\tilde{n}}_{i,\bar{\sigma }}{\tilde{n}}_{i+1,\bar{\sigma }}+{J}_{r}\left(\eta \right){\tilde{n}}_{i,\bar{\sigma }}{\tilde{h}}_{i+1,\bar{\sigma }}+{J}_{r}\left(\eta \right){\tilde{h}}_{i,\bar{\sigma }}{\tilde{n}}_{i+1,\bar{\sigma }},\end{eqnarray} \tag{ 36 }$$

where the hardcore hole operator is ${\tilde{h}}_{i,\bar{\sigma }}=1-{\tilde{n}}_{i,\bar{\sigma }}.$ Therefore, the dressed tunnelings for hardcore bosons can be described pictorially by figures 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 section 2.1, 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 section 3, we emphasize that all the quantum many-body models discussed there have a hardcore-boson counterpart, including the dynamical Gauge fields.

The objective of this section is to relax the hardcore constraint ${U}_{\uparrow \uparrow },{U}_{\downarrow \downarrow }\gg {U}_{\uparrow \downarrow }\gg {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 equation (49), but restricted to a single pseudospin (e.g. $\sigma =\uparrow$). For notational convenience, we drop the pseudospin index, such that the Hamiltonian corresponds to the standard Bose–Hubbard model

$$\begin{eqnarray}&&{H}_{{\rm{BH}}}={H}_{{\rm{loc}}}+{H}_{{\rm{kin}}}+{V}_{{\rm{int}}}=\displaystyle \sum _{i}{\epsilon }_{i}{b}_{i}^{\dagger }{b}_{i}^{\phantom{\dagger }}-\displaystyle \sum _{i}\left({t}_{x}{b}_{i}^{\dagger }{b}_{i+1}^{\phantom{\dagger }}+{\rm{H.c.}}\right)+\displaystyle \sum _{i}\displaystyle \frac{U}{2}{n}_{i}({n}_{i}^{\phantom{\dagger }}-1).\end{eqnarray} \tag{ 37 }$$

According to our discussion of the Hubbard blockade ${t}_{x}\ll U,$ the tunnelling of bosons that changes the parity of the occupation number is energetically forbidden (see figure 8(a)). As customary, we activate this tunnelling by means of the periodic modulation

$$\begin{eqnarray}&&{H}_{{\rm{mod}}}(t)=\displaystyle \sum _{i}\displaystyle \frac{{\tilde{V}}_{0}}{2}\mathrm{cos}\left({\rm{\Delta }}{{kX}}_{i}-{\rm{\Delta }}\omega t+\varphi \right){n}_{i}^{\phantom{\dagger }},\end{eqnarray} \tag{ 38 }$$

given by a moving optical lattice acting on the bosonic atoms (see figure 8(b)). Due to the different Hubbard interaction, which now involves a single pseudospin, one cannot use the expression in equation (6). Instead, we consider the interaction picture of a bond operator ${B}_{i,i+1}={b}_{i}^{\dagger }{b}_{i+1}^{\phantom{\dagger }}+{b}_{i+1}^{\dagger }{b}_{i}^{\phantom{\dagger }}$ which, again up to an irrelevant Gauge transformation, can be shown to be

$$\begin{eqnarray}&&{U}_{0}(t){B}_{i,i+1}{U}_{0}^{\dagger }(t)=\\&&{{\rm{e}}}^{-{\rm{i}}{tU}({n}_{i+1}-{n}_{i}+1)}{{\rm{e}}}^{-{\rm{i}}\displaystyle \frac{{\tilde{V}}_{0}}{2{\rm{\Delta }}\omega }\mathrm{sin}({\rm{\Delta }}{{kX}}_{i}-{\rm{\Delta }}\omega t+\varphi )}{{\rm{e}}}^{{\rm{i}}\displaystyle \frac{{\tilde{V}}_{0}}{2{\rm{\Delta }}\omega }\mathrm{sin}({\rm{\Delta }}{{kX}}_{i+1}-{\rm{\Delta }}\omega t+\varphi )}{b}_{i}^{\dagger }{b}_{i+1}^{\phantom{\dagger }}+{\rm{H.c.}},\end{eqnarray} \tag{ 39 }$$

where ${U}_{0}(t)={\mathcal{T}}\left(\mathrm{exp}\{{\rm{i}}{\displaystyle \int }_{0}^{t}{\rm{d}}\tau ({V}_{{\rm{int}}}+{H}_{{\rm{mod}}}(\tau ))\}\right)$ is the interaction-picture unitary. After defining the bosonic population difference operator ${\rm{\Delta }}{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}_{{\rm{kin}}}(t)={U}_{0}(t){H}_{{\rm{kin}}}{U}_{0}^{\dagger }(t),$ becomes

$$\begin{eqnarray}&&{H}_{{\rm{kin}}}(t)=-\displaystyle \sum _{i,\sigma }\left({t}_{x}(t){b}_{i}^{\dagger }{b}_{i+1}^{\phantom{\dagger }}+{\rm{H.c.}}\right),{t}_{x}(t)={t}_{x}{{\rm{e}}}^{-{\rm{i}}{tU}({\rm{\Delta }}{n}_{i+1}+1)}f(t),\end{eqnarray} \tag{ 40 }$$

with exactly the same modulation function as in equation (10).

Zoom In Zoom Out Reset image size

We can then proceed by following the same steps as for the fermionic case, to find a parameter regime

$$\begin{eqnarray}&&{t}_{x},\delta U=(U-r{\rm{\Delta }}\omega )\ll U\approx r{\rm{\Delta }}\omega ,\end{eqnarray} \tag{ 41 }$$

and finally obtain the effective interaction-dependent PAT Hamiltonian

$$\begin{eqnarray}&&{H}_{{\rm{eff}}}=\displaystyle \sum _{i}{\epsilon }_{i}{b}_{i}^{\dagger }{b}_{i}^{\phantom{\dagger }}-\displaystyle \sum _{i}\left({t}_{x}{J}_{r({\rm{\Delta }}{n}_{i+1}+1)}\left(\eta \right){{\rm{e}}}^{-{\rm{i}}r\varphi ({\rm{\Delta }}{n}_{i+1}+1)}{b}_{i}^{\dagger }{b}_{i+1}^{\phantom{\dagger }}+{\rm{H.c.}}\right)+\displaystyle \sum _{i}\displaystyle \frac{\delta U}{2}{n}_{i}({n}_{i}-1).\end{eqnarray} \tag{ 42 }$$

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 figures 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 (36) quadratic in the density operators, but rather as a highly nonlinear term. Using the orthogonal projectors onto subspaces with a well-defined difference number of bosons ${{\mathcal{P}}}_{{\rm{\Delta }}{n}_{i+1}=\pm {\ell }},$ this nonlinearity becomes apparent

$$\begin{eqnarray}&&{J}_{r({\rm{\Delta }}{n}_{i+1}+1)}\left(\eta \right)=\displaystyle \sum _{{\ell }^{\prime} =0}^{{N}_{{\rm{b}}}+1}{J}_{r{\ell }^{\prime} }(\eta ){{\mathcal{P}}}_{{\rm{\Delta }}{n}_{i+1}+1=\pm {\ell }^{\prime} }=\displaystyle \sum _{{\ell }^{\prime} =0}^{{N}_{{\rm{b}}}+1}{J}_{r{\ell }^{\prime} }(\eta )\displaystyle \prod _{{\ell }\ne {\ell }^{\prime} }\displaystyle \frac{{{\ell }}^{2}-{({\rm{\Delta }}{n}_{i+1}+1)}^{2}}{{{\ell }}^{2}-{\ell }{^{\prime} }^{2}},\end{eqnarray} \tag{ 43 }$$

where ${N}_{{\rm{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 [46, 47], we observe that there will be interaction-shifted resonances 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 [46, 47]. For the intensity-modulated lattices of the recent experiments [48, 49], 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 [53, 54]. Let us finally note that it is possible to generalize this scheme to higher dimensions, paralleling the fermionic case section 2.1.

In order to test the validity of our derivations, we have numerically compared the time-evolution predicted by either the effective Hamiltonian (42), (43), or the periodically driven one (37), (38) for a Bose–Hubbard dimer in the regime of interaction-dependent PAT. In figure 9, we explore the real-time dynamics for different configurations of atoms in the initial state: (i) Non-blockaded regime: figures 9(a), (b) represent the dynamics of the initial atomic configurations $| 0,1\rangle ={b}_{2}^{\dagger }| {\rm{vac}}\rangle ,$ which does not display Hubbard blockade as it consists of a single atom. Yet, the bare tunneling (see dashed lines of figure 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\rangle$ in figures 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 figure 9(b). (ii) Hubbard-blockaded regime: figures 9(e)–(h) represent the dynamics of the initial atomic configuration $| 0,2\rangle ,$ 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 figure 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 (41) is achieved forδ U = 0. In figures 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 figure 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.

Zoom In Zoom Out Reset image size

As a further numerical proof of the consistency of our effective description, let us explore the phenomenon of coherent destruction of tunneling. According to equation (43), the effective tunneling is dressed by a different Bessel function depending on the densities of the bosonic atoms. For instance, the PAT tunneling of figure 10(a) is controlled by ${J}_{0}(\eta ),$ whereas the tunneling of figure 10(c) is controlled by ${J}_{r}(\eta ).$ Therefore, whenever the driving parameter η coincides with a zero of the corresponding Bessel function, the tunneling should get coherently suppressed. In figure 10(b), we observe this effect at $\eta ={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 figure 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{H}_{{\rm{BH}}}={H}_{{\rm{loc}}}+{H}_{{\rm{kin}}}+{V}_{{\rm{int}}}=\displaystyle \sum _{i,\sigma }{\epsilon }_{i,\sigma }{b}_{i,\sigma }^{\dagger }{b}_{i,\sigma }^{\phantom{\dagger }}-\displaystyle \sum _{i,\sigma }\left({t}_{x}{b}_{i,\sigma }^{\dagger }{b}_{i+1,\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\&&+\displaystyle \frac{1}{2}\displaystyle \sum _{i}\displaystyle \sum _{\sigma ,{\sigma }^{\prime }}{U}_{\sigma {\sigma }^{\prime }}{b}_{i,\sigma }^{\dagger }{b}_{i,{\sigma }^{\prime }}^{\dagger }{b}_{i,{\sigma }^{\prime }}^{\phantom{\dagger }}{b}_{i,\sigma }^{\phantom{\dagger }}.\end{eqnarray} \tag{ 44 }$$

We shall be interested in the regime of ${U}_{\uparrow \downarrow }\gg {U}_{\uparrow \uparrow },{U}_{\downarrow \downarrow },{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)={\mathcal{T}}\left(\mathrm{exp}\{{\rm{i}}{\displaystyle \int }_{0}^{t}{\rm{d}}\tau ({V}_{\mathrm{int},\uparrow \downarrow }+{H}_{{\rm{mod}}}(\tau ))\}\right),$ where a spin-independent moving optical lattice is applied to both pseudospins

$$\begin{eqnarray}&&{H}_{{\rm{mod}}}(t)=\displaystyle \sum _{i,\sigma }\displaystyle \frac{{\tilde{V}}_{0}}{2}\mathrm{cos}\left({\rm{\Delta }}{{kX}}_{i}-{\rm{\Delta }}\omega t+\varphi \right){n}_{i,\sigma }^{\phantom{\dagger }}.\end{eqnarray} \tag{ 45 }$$

One can then see that, in the following parameter regime

$$\begin{eqnarray}&&{t}_{x},{U}_{\uparrow \uparrow },{U}_{\downarrow \downarrow },\delta {U}_{\uparrow \downarrow }=({U}_{\uparrow \downarrow }-r{\rm{\Delta }}\omega )\ll {U}_{\uparrow \downarrow }\approx r{\rm{\Delta }}\omega ,\end{eqnarray} \tag{ 46 }$$

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

$$\begin{eqnarray}&&{H}_{{\rm{eff}}}=\displaystyle \sum _{i,\sigma }{\epsilon }_{i,\sigma }{b}_{i,\sigma }^{\dagger }{b}_{i,\sigma }^{\phantom{\dagger }}-\displaystyle \sum _{i,\sigma }\left({t}_{x}{{\rm{J}}}_{r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left(\eta \right){{\rm{e}}}^{-{\rm{i}}r\varphi {\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}{b}_{i,\sigma }^{\dagger }{b}_{i+1,\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\&&+\displaystyle \frac{1}{2}\displaystyle \sum _{i}\displaystyle \sum _{\sigma ,{\sigma }^{\prime }}{\tilde{U}}_{\sigma \sigma ^{\prime} }{b}_{i,\sigma }^{\dagger }{b}_{i,{\sigma }^{\prime }}^{\dagger }{b}_{i,{\sigma }^{\prime }}^{\phantom{\dagger }}{b}_{i,\sigma }^{\phantom{\dagger }},\end{eqnarray} \tag{ 47 }$$

where we have introduced the residual interactions ${\tilde{U}}_{\uparrow \uparrow }={U}_{\uparrow \uparrow },$ ${\tilde{U}}_{\downarrow \downarrow }={U}_{\downarrow \downarrow },$ and ${\tilde{U}}_{\uparrow \downarrow }={\tilde{U}}_{\downarrow \uparrow }=\delta {U}_{\uparrow \downarrow }.$ 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 (15) is the highly nonlinear function of the bosonic densities, namely

$$\begin{eqnarray}&&{J}_{r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left(\eta \right)=\displaystyle \sum _{{\ell }^{\prime} =0}^{{N}_{{\rm{b}}}}{J}_{r{\ell }^{\prime} }(\eta ){{\mathcal{P}}}_{{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}=\pm {\ell }^{\prime} }=\displaystyle \sum _{{\ell }^{\prime} =0}^{{N}_{{\rm{b}}}}{J}_{r{\ell }^{\prime} }(\eta )\displaystyle \prod _{{\ell }\ne {\ell }^{\prime} }\displaystyle \frac{{{\ell }}^{2}-{({\rm{\Delta }}{n}_{i+1,\bar{\sigma }})}^{2}}{{{\ell }}^{2}-{\ell }{^{\prime} }^{2}}.\end{eqnarray} \tag{ 48 }$$

It is also worth commenting that, had we set ${U}_{\uparrow \downarrow }={U}_{\uparrow \uparrow }={U}_{\downarrow \downarrow }\gg {t}_{x}$ in the blockade conditions, the dressed tunneling would depend on the density of both pseudospins ${J}_{r{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left(\eta \right){{\rm{e}}}^{-{\rm{i}}r\varphi {\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\to {J}_{r({\rm{\Delta }}{n}_{i+1,\sigma }+{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}+1)}\left(\eta \right){{\rm{e}}}^{-{\rm{i}}r\varphi ({\rm{\Delta }}{n}_{i+1,\sigma }+{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}+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 figure 11).

Zoom In Zoom Out Reset image size

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 section 3), cannot be generalized to the soft-core boson case. Let us finally note that, by considering the hardcore-boson limit on equation (47), one recovers the previous result (35).

In this section, we focus on a quantum simulator of the Fermi–Hubbard model with tunable ratios of ${t}_{\alpha }/U,{X}_{\alpha }/U,{\tilde{X}}_{\alpha }/U,$ which may lead to very interesting many-body effects. Such a Fermi–Hubbard model ${H}_{{\rm{FH}}}={H}_{{\rm{loc}}}+{H}_{{\rm{kin}}}^{{\rm{corr}}}+{V}_{{\rm{int}}}$ can be rewritten in terms of asymmetric tunnelings that are correlated to the electron/hole occupation

$$\begin{eqnarray}&&{H}_{{\rm{kin}}}^{{\rm{corr}}}=\displaystyle \sum _{{\bf{i}},{\boldsymbol{\alpha }}}\displaystyle \sum _{\sigma }\left({t}_{{\rm{hh}}}^{\alpha }{h}_{{\bf{i}},\bar{\sigma }}{h}_{{\bf{i}}+{{\bf{e}}}_{\alpha },\bar{\sigma }}+{t}_{{\rm{eh}}}^{\alpha }({n}_{{\bf{i}},\bar{\sigma }}{h}_{{\bf{i}}+{{\bf{e}}}_{\alpha },\bar{\sigma }}+{h}_{{\bf{i}},\bar{\sigma }}{n}_{{\bf{i}}+{{\bf{e}}}_{\alpha },\bar{\sigma }})+{t}_{{\rm{ee}}}^{\alpha }{n}_{{\bf{i}},\bar{\sigma }}{n}_{{{\bf{ie}}}_{\alpha },\bar{\sigma }}\right){f}_{{\bf{i}},\sigma }^{\dagger }{f}_{{\bf{i}},+{{\bf{e}}}_{\alpha },\sigma }^{\phantom{\dagger }}+{\rm{H.c.}},\end{eqnarray} \tag{ 58 }$$

where we have introduced tunnelings in a hole–hole background ${t}_{{\rm{hh}}}^{\alpha }=-{t}_{\alpha },$ in an electron–hole background ${t}_{{\rm{eh}}}^{\alpha }=-{t}_{\alpha }+{X}_{\alpha },$ and in an electron–electron background ${t}_{{\rm{ee}}}^{\alpha }=-{t}_{\alpha }+{\tilde{X}}_{\alpha }+2{X}_{\alpha }.$

Let us focus for simplicity on the 1D limit of equation (58). The effective cold fermion Hamiltonian (14), (15), which is obtained through the PAT by a spin-independent moving lattice, already contains these correlated particle–hole tunnelings. For the parameter regime fulfilling (11) with r = 2, and setting the moving lattice relative phase to φ = 0, one finds

$$\begin{eqnarray}&&{t}_{{\rm{hh}}}^{x}=-{t}_{x}{J}_{0}(\eta ),{t}_{{\rm{eh}}}^{x}=-{t}_{x}{J}_{2}(\eta ),{t}_{{\rm{ee}}}^{x}=-{t}_{x}{J}_{0}(\eta ).\end{eqnarray} \tag{ 59 }$$

Accordingly, the correlated tunneling asymmetry can be tuned by modifying the intensity of the moving optical lattice ${\tilde{V}}_{0},$ or its frequency Δω, such that $\eta ={\tilde{V}}_{0}/{\rm{\Delta }}\omega$ is varied. This effect could also be achieved through a periodic modulation of the s-wave scattering length as considered in other recent schemes [55, 56].

The 1D Fermi–Hubbard model with the asymmetric correlated tunnelling strengths (59) 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 (59). This model was initially studied for ${t}_{{\rm{ee}}}^{x}={t}_{{\rm{hh}}}^{x}\lt {t}_{{\rm{eh}}}^{x},$ where a groundstate with spin-density-wave order can be found for sufficiently strong repulsion at half-filling [78]. Later on, it was realized that for ${t}_{{\rm{ee}}}^{x}={t}_{{\rm{hh}}}^{x}\gt {t}_{{\rm{eh}}}^{x},$ a new density wave where charge alternates on the bonds (i.e. bond-ordered wave) can be stabilized for not too large interactions [79]. 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 [56]. 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 (14), (15), 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 [56], 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 set ${\tilde{X}}_{\alpha }=0$ without making X_α = 0 simultaneously, as would be required to study solely the effects of bond–charge interactions (55). This also occurs for the schemes based on a periodic modulation of the s-wave scattering length [55, 56]. In the following subsection, we shall show that this becomes possible by introducing an additional linear gradient in our scheme.

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}_{\alpha }/U,{X}_{\alpha }/U.$ For simplicity, we focus on the 1D case, and consider the PAT by a spin-independent moving lattice (14) in the presence of an additional linear gradient, which leads to the dressed tunneling in equation (17). For the sake of concreteness, we consider a parameter regime fulfilling equation (11) for r = 2, and set φ = 0. By adjusting the intensity of the moving lattice to ${\tilde{V}}_{0}=1.56{\rm{\Delta }}\omega$ (i.e. ${\eta }_{\star }=1.56$), such that terms like (57) in the effective Hamiltonian (14) vanish $\tilde{X}=0,$ we obtain

$$\begin{eqnarray}{H}_{{\rm{eff}}} & = & \displaystyle \sum _{i,\sigma }({\epsilon }_{i,\sigma }+\delta {U}_{\sigma \bar{\sigma }}i){f}_{i,\sigma }^{\dagger }{f}_{i,\sigma }^{\phantom{\dagger }}-\displaystyle \sum _{i,\sigma }\left({{tf}}_{i,\sigma }^{\dagger }{f}_{i+1,\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\ & & +\displaystyle \sum _{i,\sigma }{{XB}}_{i,i+1}^{\sigma }\left({n}_{i+1,\bar{\sigma }}-{n}_{i,\bar{\sigma }}\right)+\displaystyle \frac{1}{2}\displaystyle \sum _{i,\sigma }\delta {U}_{\sigma \bar{\sigma }}{f}_{i,\sigma }^{\dagger }{f}_{i,\bar{\sigma }}^{\dagger }{f}_{i,\bar{\sigma }}^{\phantom{\dagger }}{f}_{i,\sigma }^{\phantom{\dagger }}.\end{eqnarray} \tag{ 60 }$$

where $t={t}_{x}{J}_{2}\left({\eta }_{\star }\right),$ and the bond–charge interaction $X={t}_{x}({J}_{0}\left({\eta }_{\star }\right)-{J}_{2}\left({\eta }_{\star }\right))$ alternates between neighboring sites. Other possible values of the moving-lattice intensity fulfilling ${\tilde{X}}_{\alpha }=0$ are ${\eta }_{\star }=\{4.89,8.29,11.53,\ldots \},$ and correspond to solutions of the equation ${J}_{0}\left({\eta }_{\star }\right)+{J}_{4}\left({\eta }_{\star }\right)-2{J}_{2}\left({\eta }_{\star }\right)=0.$ Moreover, one can also change the moving-lattice detuning such that $r=\{4,6,8,\ldots \},$ and look for the solutions of ${J}_{0}\left({\eta }_{\star }\right)\;+\;$ ${J}_{2r}\left({\eta }_{\star }\right)\;-\;$ $2{J}_{r}\left({\eta }_{\star }\right)=0.$ This will yield different values of ${\eta }_{\star },$ 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 [67], where the bond-charge interaction $0\lt X\ll t$ can be responsible for a superconducting phase even in the presence of a repulsive interaction $\delta {U}_{\uparrow \downarrow }\gt 0,$ as has been predicted using a mean-field approximation for 2D [74]. Since our effective Hamiltonian (60) can be generalized to higher dimensions following our results in section 2, the quantum simulator could test the correctness of such mean-field predictions [67]. 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 equation (60).

Since the residual Hubbard interactions in (60) depend on the detuning of the photon-assisted scheme, one can also study the attractive case $\delta {U}_{\uparrow \downarrow }\lt 0.$ It has been shown [75] that a bond-charge interaction X = t stabilizes an η-pairing groundstate [76] for finite attraction, which for X = 0 only occurs for strictly infinite interactions $\delta {U}_{\uparrow \downarrow }\to \infty$ [77]. Although our quantum simulator cannot reach the exact condition of X = t (e.g. for ${\eta }_{\star }=1.56,$ $X/t\approx 0.94$), it has been argued that relaxing these conditions may still host the η-pairing groundstate, even if the exact methods of [75] 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 $\delta {U}_{\sigma \bar{\sigma }}$ in equation (60), 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 section 2.

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 [80]. One of the few rigorous results on the existence of ferromagnetism in the Fermi–Hubbard model is due to Nagaoka [81], 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 [82], 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 [83]. Despite some initial discrepancy regarding its full polarization [84], more recent numerical results based on quantum Monte Carlo [85] and density-matrix renormalization group [86] agree with the above scenario [83]. 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 [87], or a large spin-imbalance in an optical lattice with a ladder structure [88] 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 (20). By looking at the tunnelings in equations (62)–(64), one notices that by modifying the intensity of the moving optical lattice such that ${J}_{r}({\eta }_{\star })=0,$ while ${J}_{0}({\eta }_{\star })\ne 0,$ we can exactly cancel the terms that do not preserve the parity in the occupation number, namely equations (63), (64). For concreteness, we consider a parameter regime fulfilling (11) for r = 2, such that the moving lattices have the same intensity and detuning, and ${\eta }_{\star }={\tilde{V}}_{0,\alpha }/{\rm{\Delta }}{\omega }_{\alpha }=5.135.$ Accordingly, only the parity-conserving tunnelling (62) is preserved. This effect is similar to the so-called coherent destruction of tunnelling [89], 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 figure 2), such that the effect can be understood as a correlated destruction of tunnelling.

For large but finite Hubbard repulsion $\delta {U}_{\uparrow \downarrow },$ 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 $\langle {n}_{{\bf{i}},\uparrow }+{n}_{{\bf{i}},\downarrow }\rangle \lt 1,$ all the dynamics takes place within the single-occupation subspace, and is controlled by

$$\begin{eqnarray}&&{H}_{{\rm{eff}}}={{\mathcal{P}}}_{{\rm{s}}}\displaystyle \sum _{{\bf{i}},\alpha }\displaystyle \sum _{\sigma }\left(-{\tilde{t}}_{\alpha }{f}_{{\bf{i}},\sigma }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right){{\mathcal{P}}}_{{\rm{s}}},\end{eqnarray} \tag{ 65 }$$

where ${\tilde{t}}_{\alpha }={t}_{\alpha }{J}_{0}({\eta }_{\star }),$ and we have introduced the Gutzwiller projector onto the single-occupied sub-band ${{\mathcal{P}}}_{{\rm{s}}}={\displaystyle \prod }_{{\bf{i}}}(1-{n}_{{\bf{i}},\uparrow }{n}_{{\bf{i}},\downarrow }).$ Interestingly enough, the Hamiltonian (65) 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 [86], 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.

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. $\langle {n}_{{\bf{i}},\uparrow }+{n}_{{\bf{i}},\downarrow }\rangle \lt 1$). However, a large but finite repulsion can also lead to very interesting physics. In this regime, the competition of the projected kinetic energy (65) with the antiferromagnetic super-exchange [13] leads to the so-called t-J model

$$\begin{eqnarray}&&{H}_{{tJ}}={{\mathcal{P}}}_{{\rm{s}}}{\tilde{H}}_{{tJ}}{{\mathcal{P}}}_{{\rm{s}}},{\tilde{H}}_{{tJ}}=\displaystyle \sum _{{\bf{i}},\alpha }\displaystyle \sum _{\sigma }\left(-\tilde{t}{f}_{{\bf{i}},\sigma }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)+\displaystyle \sum _{{\bf{i}},\alpha }\tilde{J}\left({{\boldsymbol{S}}}_{{\bf{i}}}\cdot {{\boldsymbol{S}}}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}}}-\displaystyle \frac{1}{4}{n}_{{\bf{i}}}{n}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}}}\right),\end{eqnarray} \tag{ 66 }$$

where $\tilde{t}$ is the tunnelling within the single-occupied subspace, and $\tilde{J}\gt 0$ is the strength of the antiferromagnetic super-exchange interaction. In the case of square lattices, this Hamiltonian (66) has been considered as the canonical effective model in the theory of the high-${T}_{{\rm{c}}}$ cuprates by part of the scientific community [16, 90, 91]. In this context, the t-J model arises after mapping a more microscopic three-band Hubbard model of the cuprates [92] onto a single-band one [73, 93]. By using the microscopic parameters of the three-band model [94], one finds that $\tilde{J}/\tilde{t}\sim 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 repulsion $\tilde{t}\ll {U}_{\uparrow \downarrow }.$ Then, one finds $\tilde{J}=4{\tilde{t}}^{2}/{U}_{\uparrow \downarrow },$ which cannot attain values larger than $\tilde{J}\lt 0.4\tilde{t}$ since the tunnelling must be at least $\tilde{t}\lt 0.1{U}_{\uparrow \downarrow }$ 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⁶ $\tilde{J}/\tilde{t}\sim 0.3$. At this point, it should be mentioned that reaching the low temperatures required to observe the effect of the antiferromagnetic super exchange at equilibrium is a great challenge that has only been achieved recently [95]. However, from a broader perspective, the rest of the rich phase diagram of the t-J model cannot be explored with these experiments. A possibility to attain tunability over these parameters, while at the same time controlling and homogeneous atomic doping from half-filling, would be to consider composite-fermion quasiparticles whose tunneling corresponds to a correlated tunneling of a boson–fermion pair in a Bose–Fermi mixtures [96]. We show below that our PAT scheme leads to a standard fermionic t-J model where the ratio $\tilde{J}/\tilde{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 (20) for any dimensionality. We assume equal tunnelings ${t}_{\alpha }=:\;t,$ driving parameters ${r}_{\alpha }=:\;r,$ and driving ratios ${\eta }_{\alpha }=:\;\eta ,$ in all directions $\forall \alpha ,$ and set φ_α = 0. The limit of strong Hubbard repulsion in this case corresponds to ${{tJ}}_{r}\left(\eta \right)\ll \delta {U}_{\uparrow \downarrow },$ where the parity-violating tunnelings described by the terms ${K}_{{\rm{s}}\to {\rm{d}}}$ and ${K}_{{\rm{d}}\to {\rm{s}}}$ in equations (63)–(64) can only take place virtually. As in the standard Hubbard model [97, 98], such virtual tunnelings can be calculated by a Schrieffer–Wolff-type unitary transformation ${\tilde{H}}_{{\rm{eff}}}={{\rm{e}}}^{{\rm{i}}S}{H}_{{\rm{eff}}}{{\rm{e}}}^{-{\rm{i}}S},$ where $S=-{\rm{i}}({K}_{{\rm{s}}\to {\rm{d}}}-{K}_{{\rm{d}}\to {\rm{s}}})/\delta {U}_{\uparrow \downarrow }$ is responsible for eliminating the energetically forbidden tunnelings to second order in the small expansion parameter $\xi ={{tJ}}_{r}\left(\eta \right)/\delta {U}_{\uparrow \downarrow }.$ Considering the commutation properties of the different operators defined so far, one finds that ${\tilde{H}}_{{\rm{eff}}}={H}_{{\rm{loc}}}+{K}_{0}+{V}_{{\rm{int}}}+([{K}_{{\rm{s}}\to {\rm{d}}},{K}_{{\rm{d}}\to {\rm{s}}}]+[{K}_{{\rm{s}}\to {\rm{d}}},{K}_{0}]+[{K}_{0},{K}_{{\rm{d}}\to {\rm{s}}}])/\delta {U}_{\uparrow \downarrow }+{\mathcal{O}}({\xi }^{2}).$ For hole-doping about half-filling, one can project onto the single-occupancy sub-band, which yields the aforementioned t-J model (66) with an additional density-dependent next–nearest–neighbor tunnelling

$$\begin{eqnarray}&&{\tilde{H}}_{{\rm{eff}}}={H}_{{tJ}}+{{\mathcal{P}}}_{{\rm{s}}}{\rm{\Delta }}H{{\mathcal{P}}}_{{\rm{s}}},{\rm{\Delta }}H \\&&=\displaystyle \sum _{{\bf{i}},\sigma }\displaystyle \sum _{\alpha ,n}\displaystyle \frac{J}{4}\left(-{f}_{{\bf{i}},\sigma }^{\dagger }{n}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\bar{\sigma }}^{\phantom{\dagger }}{f}_{{\bf{i}}+{{\bf{u}}}_{n}^{{\boldsymbol{\alpha }}},\sigma }^{\phantom{\dagger }}+{f}_{{\bf{i}},\bar{\sigma }}^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\sigma }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\bar{\sigma }}^{\phantom{\dagger }}{f}_{{\bf{i}}+{{\bf{u}}}_{n}^{{\boldsymbol{\alpha }}},\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right),\end{eqnarray} \tag{ 67 }$$

where we have introduced the effective cold-atom parameters

$$\begin{eqnarray}&&\tilde{t}={{tJ}}_{0}(\eta ),\tilde{J}=4{t}^{2}{J}_{r}{(\eta )}^{2}/\delta {U}_{\uparrow \downarrow },\end{eqnarray} \tag{ 68 }$$

and the next–nearest–neighbor vectors ${{\bf{u}}}_{n}^{\alpha }$ (e.g. for 1D ${{\bf{u}}}_{1}^{x}=2{{\bf{e}}}_{x},$ for 2D $\{{{\bf{u}}}_{1}^{x}=2{{\bf{e}}}_{x},{{\bf{u}}}_{2}^{x}={{\bf{e}}}_{x}-{{\bf{e}}}_{y},{{\bf{u}}}_{1}^{y}=2{{\bf{e}}}_{y},{{\bf{u}}}_{2}^{y}={{\bf{e}}}_{x}+{{\bf{e}}}_{y}\}$). Let us note that this additional tunnelling requires that the target site is populated with a hole, as otherwise ${{\mathcal{P}}}_{{\rm{s}}}{f}_{{\bf{i}},\sigma }^{\dagger }| {\rm{\Psi }}\rangle =0.$ Hence, close to half-filling $\langle {n}_{i}\rangle \approx 1,$ this term is reduced by a factor $(1-\langle {n}_{i}\rangle )/4$ with respect to the Heisenberg super-exchange, and is typically neglected in the literature [90].

As announced above, we have obtained an effective t-J model where the ratio of the coupling constants (68), namely $\tilde{J}/\tilde{t}=4\xi {J}_{r}(\eta )/{J}_{0}(\eta ),$ can be tuned by modifying the intensity of the moving optical lattice. Even if $\xi \leqslant 0.1$ in the regime of validity of the t-J model, we can make ${J}_{r}(\eta )\gg {J}_{0}(\eta ),$ 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 [99]. 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 [100]), a gapless superconductor (i.e. an attractive Luttinger liquid [100]), or the so-called spin-gap phase [101], 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 [102]). Finally, for sufficiently strong super-exchange interactions, the antiferromagnetic order expels the doped holes leading to separated hole-rich and hole-poor regions (i.e. phase separation [103]). 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, requires $\tilde{J}\approx 2.5\tilde{t}$ 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 [104] predicting the disappearance of the phase separation in favor of enlarged superconducting and spin-gap regions, once the next–nearest–neighbor tunnelling terms in equation (67) 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 [105], 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 [90] 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}_{{\rm{c}}}$ cuprates. For instance, the results of [90, 107] contradict those of [103, 106], 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 [90], exact diagonalization [111], and quantum Monte Carlo [112]. From the perspective of the high-${T}_{{\rm{c}}}$ cuprates, addressing the conflicting predictions in [108, 109] and [90, 110] 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, 90, 91], 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 (66) 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 tunnelings ${\tilde{t}}_{\alpha }={t}_{\alpha }{J}_{0}({\eta }_{\alpha }),$ and the super-exchange couplings ${\tilde{J}}_{\alpha }=4{t}_{\alpha }^{2}{J}_{r}{({\eta }_{\alpha })}^{2}/\delta {U}_{\uparrow \downarrow },$ 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 [113], which can be modeled by a number ${\ell }\in \{1,...,{n}_{{\ell }}\}$ of 1D 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}_{{tJ}}^{{\rm{ladder}}}={{\mathcal{P}}}_{{\rm{s}}}{\tilde{H}}_{{tJ}}^{{\rm{ladder}}}{{\mathcal{P}}}_{{\rm{s}}},$ where

$$\begin{eqnarray}{\tilde{H}}_{{tJ}}^{{\rm{ladder}}} & = & \displaystyle \sum _{i,\sigma }\displaystyle \sum _{{\ell }=1}^{{n}_{{\ell }}}\left(-\tilde{t}{f}_{i,{\ell },\sigma }^{\dagger }{f}_{i+1,{\ell }\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)+\displaystyle \sum _{i}\displaystyle \sum _{{\ell }}\tilde{J}\left({{\boldsymbol{S}}}_{i,{\ell }}\cdot {{\boldsymbol{S}}}_{i+1,{\ell }}-\displaystyle \frac{1}{4}{n}_{i,{\ell }}{n}_{i+1,{\ell }}\right)\\ & & +\displaystyle \sum _{i,\sigma }\displaystyle \sum _{{\ell }=1}^{{n}_{{\rm{l}}}}\left(-t^{\prime} {f}_{i,{\ell },\sigma }^{\dagger }{f}_{i,{\ell }+1\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)+\displaystyle \sum _{i}\displaystyle \sum _{{\ell }}J^{\prime} \left({{\boldsymbol{S}}}_{i,{\ell }}\cdot {{\boldsymbol{S}}}_{i,{\ell }+1}-\displaystyle \frac{1}{4}{n}_{i,{\ell }}{n}_{i,{\ell }+1}\right).\end{eqnarray} \tag{ 69 }$$

Such ladder Hamiltonians can be implemented in our quantum simulator if we supplement the above scheme (67) 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 ${\tilde{\lambda }}_{y}=2{\lambda }_{y},$ leads to an array of decoupled ladders with ${n}_{{\ell }}=2$ legs [114]. 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}_{{\ell }}$-legged ladders). Hence, the interaction-dependent PAT scheme leads to equation (69) with tunable number of legs and Hamiltonian parameters

$$\begin{eqnarray}&&\tilde{t}={t}_{x}{J}_{0}({\eta }_{x}),\tilde{J}=\displaystyle \frac{4{t}_{x}^{2}{J}_{r}{({\eta }_{x})}^{2}}{\delta {U}_{\uparrow \downarrow }},t^{\prime} ={t}_{y}{J}_{0}({\eta }_{y}),J^{\prime} =\displaystyle \frac{4{t}_{y}^{2}{J}_{r}{({\eta }_{y})}^{2}}{\delta {U}_{\uparrow \downarrow }},\end{eqnarray} \tag{ 70 }$$

together with the corresponding density-dependent next-to-nearest-neighbor tunneling (67), 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}_{{\rm{c}}}$ superconductivity. In fact, already for ${n}_{{\ell }}=2$ legs, the doped holes tend to pair [115] developing a superconducting d-wave-like order [116]. In addition to the phases that also occur for the 1D t-J model [99], this d-wave superconductivity takes place in a wide region of the phase diagram [117], 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 [118], whereby the hole d-wave pairing disappears for ladders with an odd numbers of legs. We note that some of these results [115] depend on ratios $t^{\prime} /J^{\prime} \lt 1$ that cannot be reached from standard ladder Hubbard models where $t^{\prime} \gg J^{\prime} .$ Likewise, the regime $J^{\prime} \gg \tilde{t},t^{\prime} ,\tilde{J},$ and its connection to short-range resonating valence bond states [119, 120], 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 super-exchange interactions. This leads to a very interesting t-XXZ model described by the Hamiltonian ${H}_{{tXXZ}}={{\mathcal{P}}}_{{\rm{s}}}{\tilde{H}}_{{tXXZ}}{{\mathcal{P}}}_{{\rm{s}}}$

$$\begin{eqnarray}&&{\tilde{H}}_{{tXXZ}}=\displaystyle \sum _{{\bf{i}},\alpha }\displaystyle \sum _{\sigma }\left(-{\tilde{t}}_{\sigma }{f}_{{\bf{i}},\sigma }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}},\sigma }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\&&+\displaystyle \sum _{{\bf{i}},\alpha }\left({\tilde{J}}_{\perp }\left({S}_{{\bf{i}}}^{x}{S}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}}}^{x}+{S}_{{\bf{i}}}^{y}{S}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}}}^{y}\right)+{\tilde{J}}_{z}\left({S}_{{\bf{i}}}^{z}{S}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}}}^{z}-\displaystyle \frac{1}{4}{n}_{{\bf{i}}}{n}_{{\bf{i}}+{{\bf{e}}}_{{\boldsymbol{\alpha }}}}\right)\right),\end{eqnarray} \tag{ 71 }$$

where the Ising ${\tilde{J}}_{z}$ and flip-flop ${\tilde{J}}_{\perp }$ interaction strengths are generally different, such that the t-J model is recovered when ${\tilde{J}}_{\perp }={\tilde{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 (26) in 1D, and (30) in 2D. We assume equal tunnelings ${t}_{\alpha }=:\;t,$ spin-dependent driving parameters ${r}_{\sigma ,\alpha }=:\;{r}_{\sigma },$ and ratios ${\eta }_{\sigma ,\alpha }=:\;{\eta }_{\sigma },$ in all directions $\forall \alpha ,$ and set ${\varphi }_{\sigma ,\alpha }=0.$ One may observe in figure 5(f) that second-order super-exchange interactions will depend on the spin configuration. In fact, we find

$$\begin{eqnarray}&&{\tilde{t}}_{\sigma }={{tJ}}_{0}({\eta }_{\sigma }),{\tilde{J}}_{\perp }=\displaystyle \frac{4{t}^{2}{J}_{{r}_{\uparrow }}({\eta }_{\uparrow }){J}_{{r}_{\downarrow }}({\eta }_{\downarrow })}{\delta {U}_{\uparrow \downarrow }},{\tilde{J}}_{z}=\displaystyle \frac{2{t}^{2}\left({J}_{{r}_{\uparrow }}{({\eta }_{\uparrow })}^{2}+{J}_{{r}_{\downarrow }}{({\eta }_{\downarrow })}^{2}\right)}{\delta {U}_{\uparrow \downarrow }},\end{eqnarray} \tag{ 72 }$$

such that the Heisenberg–Ising anisotropy $\zeta ={\tilde{J}}_{\perp }/{\tilde{J}}_{z}$ can be tuned all the way from small spin quantum fluctuations $\zeta \to 0$ (i.e. t-J_z model), to the isotropic regime where spin quantum fluctuations play an important role $\zeta \to 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 [121] localized to the site where they were doped [122]. The situation is considerably more complex as the spin fluctuations are switched on $\zeta \gt 0,$ and the t-J model is approached $\zeta \to 1$ [123], and some controversy regarding the limitations of the different analytical or numerical methods has been identified [124]. 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 [125, 126] 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.

Let us focus on the 1D case, and consider the PAT of a Fermi gas by a spin-dependent moving lattice (26) in a parameter regime fulfilling equation (24) for r_σ = 2, and $\delta {U}_{\uparrow \downarrow }=0.$ Let us note that the same can be obtained for a Bose–Fermi mixture (54), provided that the hardcore constraint is considered. In order to obtain a quantum simulator of quantum matter coupled to dynamical gauge fields (73), one needs to find a particular set of parameters such that: i the tunnelling amplitude (27) becomes a simple c-number, and (ii) the tunnelling phase in equation (26) is non-vanishing only for one of the pseudospin states. In order to fulfill (i), the ratio of the moving-lattice intensities with respect to the detuning must fulfill ${J}_{0}\left({\eta }_{\star }^{\sigma }\right)={J}_{2}\left({\eta }_{\star }^{\sigma }\right),$ which can be achieved with any of the following values ${\eta }_{\star }^{\sigma }\in \{4.89,8.29,11.53,\ldots \}.$ By direct substitution in equation (27), one finds that the dressed tunnelling strength becomes the desired c-number ${t}_{x}{J}_{{r}_{\sigma }{\rm{\Delta }}{n}_{i+1,\bar{\sigma }}}\left({\eta }_{\star }^{\sigma }\right)={t}_{x}{J}_{0}\left({\eta }_{\sigma ,\star }\right)=:\;{t}_{{\rm{eff,}}}^{\sigma }.$ In order to fulfill (ii), it suffices to set the relative phase of the moving optical lattice ${\varphi }_{\downarrow }=0$ for the pseudospins that will play the role of the Gauge field. Hence, the tunnelling phases in equation (26) fulfill ${\varphi }_{\downarrow }=0,$ and $2{\varphi }_{\uparrow }=:\;-\theta \ne 0,$ such that the atoms with pseudospin $| \downarrow \rangle$ play the role of the dynamical Gauge field for the atoms with pseudospin $| \uparrow \rangle$ (see figure 7), namely

$$\begin{eqnarray}&&{H}_{{\rm{eff}}}=\displaystyle \sum _{i}{\epsilon }_{i,\downarrow }{f}_{i,\uparrow }^{\dagger }{f}_{i,\downarrow }^{\phantom{\dagger }}-\displaystyle \sum _{i}{t}_{{\rm{eff}}}^{\downarrow }\left({f}_{i,\downarrow }^{\dagger }{f}_{i+1,\downarrow }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\&&-\displaystyle \sum _{i}{t}_{{\rm{eff}}}^{\uparrow }\left({f}_{i,\uparrow }^{\dagger }{{\rm{e}}}^{{\rm{i}}\theta \left({n}_{i+1,\downarrow }^{\phantom{\dagger }}-{n}_{i,\downarrow }^{\phantom{\dagger }}\right)}{f}_{i+1,\uparrow }^{\phantom{\dagger }}+{\rm{H.c.}}\right)+\displaystyle \sum _{i}{\epsilon }_{i,\uparrow }{f}_{i,\uparrow }^{\dagger }{f}_{i,\uparrow }^{\phantom{\dagger }}.\end{eqnarray} \tag{ 74 }$$

This is the general result of this section. The PAT scheme has allowed us to build an effective Hamiltonian where the atoms with one of the pseudospins hop freely in the lattice and play the role of a dynamical 'Gauge' field for the atoms with the remaining pseudospin. This must be contrasted with other interesting proposals [51, 54], where the tunnelling Peierls phase for atoms depends on their own density, such that the roles of matter and Gauge fields cannot be distinguished. We believe that these type of Hamiltonians (74), and their straightforward generalization to higher dimensions, to other lattices, to hardcore bosons, or to slightly modified PAT schemes⁷ , will lead to several interesting many-body phenomena that deserve to be studied in further detail. To illustrate this richness, we describe a particular example that leads to an interesting relativistic quantum field theory in the continuum.

The effective Hamiltonian (74) corresponds exactly to the structure of the dynamical gauge field theory in equation (73) if we make the following identifications: (i) The fermionic quantum matter is represented by the atoms with pseudospin $| \uparrow \rangle ,$ namely ${\psi }_{i}:= {f}_{i,\uparrow }.$ (ii) the Gauge degrees of freedom will be some collective low-energy excitations of the atoms with pseudospin $| \downarrow \rangle .$ 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 [100, 134], which will play the role of the gauge degrees of freedom

$$\begin{eqnarray}&&{K}_{\mathrm{eff},\downarrow }=\displaystyle \sum _{i}{\epsilon }_{i,\downarrow }{f}_{i,\downarrow }^{\dagger }{f}_{i,\downarrow }^{\phantom{\dagger }}-{t}_{{\rm{eff}}}^{\downarrow }\displaystyle \sum _{i}\left({f}_{i,\downarrow }^{\dagger }{f}_{i+1,\downarrow }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\to {H}_{{\mathcal{G}}} \\&&=\displaystyle \sum _{q}({\epsilon }_{\downarrow }+{c}_{\downarrow }q)\left({b}_{q,{\rm{R}}}^{\dagger }{b}_{q,{\rm{R}}}^{\phantom{\dagger }}+{b}_{q,{\rm{L}}}^{\dagger }{b}_{q,{\rm{L}}}^{\phantom{\dagger }}\right),\end{eqnarray} \tag{ 75 }$$

where we have introduce the quasi momentum $q=2\pi n/L$ for $n\in {{\mathbb{Z}}}^{+},$ the effective speed of light ${c}_{\downarrow }=2{t}_{{\rm{eff}}}^{\downarrow }a={t}_{{\rm{eff}}}^{\downarrow }\lambda ,$ and the bosonic operators for the particle–hole excitations ${b}_{q,{\rm{R}}}$ (${b}_{q,{\rm{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 ${\epsilon }_{i,\downarrow }={\epsilon }_{\downarrow }+\frac{1}{2}m{\omega }_{{\rm{t}},x}^{2}{X}_{i}^{2},$ can be accounted for using a local chemical potential $\delta {\mu }_{i}.$ Note that, in the continuum limit, this bosonic Hamiltonian becomes a scalar field theory with the energy zero set at ${\epsilon }_{\downarrow },$ namely

$$\begin{eqnarray}&&{H}_{{\mathcal{G}}}=\int \displaystyle \frac{{\rm{d}}x}{2\pi }{c}_{\downarrow }\left(\displaystyle \frac{1}{2}{\left(\pi (x)\right)}^{2}+\displaystyle \frac{1}{2}{\left({\partial }_{x}\phi (x)\right)}^{2}\right),\end{eqnarray} \tag{ 76 }$$

where the scalar field $\phi (x)={\phi }_{{\rm{R}}}(x)-{\phi }_{{\rm{L}}}(x)$ is expressed in terms of the inverse Fourier transform of the bosonic operators ${b}_{q,{\rm{R}}},$ ${b}_{q,{\rm{L}}}$ and their Hermitian adjoints, and $\pi (x)$ is its canonically-conjugate momentum. Although there is no Gauge invariance in the massless scalar field theory (76), 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 equation (73) becomes ${{aA}}_{x}({x}_{i,i+1})=\theta ({n}_{i+1,\downarrow }-{n}_{i,\downarrow }),$ 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

$$\begin{eqnarray}{{aA}}_{x}({x}_{i,i+1}) & = & a\theta \displaystyle \frac{{n}_{i+1,\downarrow }-{n}_{i,\downarrow }}{a}={a}^{2}\theta {\partial }_{x}\left({\psi }^{\dagger }(x)\psi (x)\right)\to {A}_{x}(x)\\ & = & \theta a{\partial }_{x}^{2}\phi (x),{U}_{i,i+1}=1+{\rm{i}}{{aA}}_{x}(x)+{\mathcal{O}}({a}^{2}{A}_{x}^{2}),\end{eqnarray} \tag{ 77 }$$

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 [135] of the so-called Kogut-Susskind fermions [129]. At half-filling, one obtains the massless Dirac field theory in the continuum limit, which minimally couples to a derivative of the scalar field

$$\begin{eqnarray}{H}_{{\rm{RQFT}}} & = & {H}_{{\mathcal{G}}}+\displaystyle \int \displaystyle \frac{{\rm{d}}x}{2\pi }\left({\tilde{{\rm{\Psi }}}}_{{\rm{L}}}^{\dagger }(x)(\delta +{c}_{\uparrow }(-{\rm{i}}{\partial }_{x}-{A}_{x}(x))){\tilde{{\rm{\Psi }}}}_{{\rm{L}}}^{\phantom{\dagger }}(x)\right.\\ & & +\left.{\tilde{{\rm{\Psi }}}}_{{\rm{R}}}^{\dagger }(x)(\delta +{c}_{\uparrow }(+{\rm{i}}{\partial }_{x}+{A}_{x}(x))){\tilde{{\rm{\Psi }}}}_{{\rm{R}}}^{\phantom{\dagger }}(x)\right),\end{eqnarray} \tag{ 78 }$$

where we have introduced the fermionic field operators ${\tilde{{\rm{\Psi }}}}_{{\rm{R}}}(x)$ (${\tilde{{\rm{\Psi }}}}_{{\rm{L}}}(x)$) for the right (left) moving fermions with pseudospin $\sigma =\uparrow ,$ the energy difference due to the different Zeeman shifts $\delta ={\epsilon }_{\uparrow }-{\epsilon }_{\downarrow },$ and the effective speed of light ${c}_{\uparrow }={t}_{{\rm{eff}}}^{\uparrow }\lambda .$ We have obtained a peculiar relativistic theory of interacting quantum fields. Instead of the standard Yukawa coupling between scalar and fermionic fields $\bar{\psi }\phi \psi ,$ we get a minimal coupling with the derivative of the scalar field $\bar{\psi }{\gamma }^{1}{\partial }_{x}^{2}\phi \psi .$ Moreover, the effective speeds of light of the scalar ${c}_{\downarrow }$ and fermionic ${c}_{\uparrow }$ 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 [136]. 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 (78). 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.

Topological insulators represent a family of holographic phases of matter that are insulating in the bulk and conducting at the boundaries [137]. 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 [138], which describes fermions in a square lattice subjected to a perpendicular magnetic field, and displays the aforementioned bulk [139] and edge [140] properties. Similar phenomenology also arises in the Haldane model [141], 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 [142] built from two time-reversed copies of the Haldane model, or the time-reversal Hofstadter model [143] built from two time-reversed copies of the Hofstadter model.

Paralleling the effect of interactions in the quantum Hall effect [144], one expects that even more exotic phases of matter will appear when considering the effect of correlations in the above models [145]. So far, the typical route to introduce such correlation effects has been to include the effect of on-site and nearest-neighbor Hubbard interactions in the Haldane [146], Kane–Mele [147], or time-reversal Hofstadter [148] models. These studies show that the topological features are robust to interactions, but no other exotic phases such as topological Mott insulators [149] (i.e. interaction-induced bulk gap and protected edge states) or topological fractional insulators [150] (i.e. fractional excitations and protected edge states) were found. In this section, we consider explicitly the 2D PAT by a spin-dependent moving lattice (28), and discuss how this scheme may be used to explore a new type of correlation effects introduced by substituting the fixed background gauge field by a dynamical one in the standard [138], and time-reversal [143], Hofstadter models. Given the recent realizations of both models with non-interacting cold atoms [44], we believe that future experiments will be able to explore the full phase diagram, and the possibility of finding more exotic phases brought by the interactions with the dynamical Gauge field.

(i) Hofstadter model in a dynamical Gauge field.— Let us consider the 2D scheme (28) leading to the effective Hamiltonian (30). We generalize our prescription for the 1D case, and set for ${r}_{\sigma }^{\alpha }=2,$ and $\delta {U}_{\uparrow \downarrow }=0.$ In order to make the dressed tunnelling amplitude a c-number, we tune again the ratio of the moving-lattice intensities with respect to the detunings to ${\eta }_{x}^{\sigma }={\eta }_{y}^{\sigma }={\eta }_{\star }^{\sigma }\in \{4.89,8.29,11.53,\ldots \},$ such that ${J}_{0}\left({\eta }_{\star }^{\sigma }\right)={J}_{2}\left({\eta }_{\star }^{\sigma }\right).$ We can thus define the effective tunnelling amplitudes along the x- and y-axes as ${t}_{{\rm{eff}},x}^{\sigma }:= {t}_{x}{J}_{0}\left({\eta }_{\sigma ,\star }\right),$ and ${t}_{{\rm{eff}},y}^{\sigma }:= {t}_{y}{J}_{0}\left({\eta }_{\sigma ,\star }\right).$ Additionally, we need to control the relative phases of the moving lattices, such that only one of the pseudospins develops a non-vanishing Peierls phase ${\varphi }_{\downarrow ,x}=0,$ but $2{\varphi }_{\uparrow ,x}\;:-\theta \ne 0.$ Accordingly, the effective Hamiltonian (30) becomes a Hofstadter model in the Landau gauge for the atoms with pseudospin $| \uparrow \rangle$

$$\begin{eqnarray}&&{H}_{{\rm{eff}}}={H}_{{\mathcal{G}}}-\displaystyle \sum _{{\bf{i}}}\left({t}_{{\rm{eff}},x}^{\uparrow }{{\rm{e}}}^{{\rm{i}}\theta \left({n}_{{\bf{i}}+{{\bf{e}}}_{x},\downarrow }-{n}_{{\bf{i}},\downarrow }\right)}{f}_{{\bf{i}},\uparrow }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{x},\uparrow }^{\phantom{\dagger }}+{t}_{{\rm{eff}},y}^{\uparrow }{f}_{{\bf{i}},\uparrow }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{y},\uparrow }^{\phantom{\dagger }}+{\rm{H.c.}}\right).\end{eqnarray} \tag{ 79 }$$

Here, the synthetic Gauge field according to equation (73) corresponds to $a{\bf{A}}({x}_{{\bf{i}},{\bf{i}}+{{\bf{e}}}_{x}})=\theta ({n}_{{\bf{i}}+{{\bf{e}}}_{x},\downarrow }-{n}_{{\bf{i}},\downarrow }){{\bf{e}}}_{x},$ and thus depends on the atoms with the remaining pseudospin $| \downarrow \rangle ,$ which evolve under the free tight-binding Hamiltonian

$$\begin{eqnarray}&&{H}_{{\mathcal{G}}}=-\displaystyle \sum _{{\bf{i}}}\left({t}_{{\rm{eff}},x}^{\downarrow }{f}_{{\bf{i}},\downarrow }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{x},\downarrow }^{\phantom{\dagger }}+{t}_{{\rm{eff}},y}^{\downarrow }{f}_{{\bf{i}},\downarrow }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{y},\downarrow }^{\phantom{\dagger }}+{\rm{H.c.}}\right).\end{eqnarray} \tag{ 80 }$$

The same can be obtained for a Bose–Fermi mixture (54), provided that the hardcore constraint is considered, such that the hardcore bosons play the role of the Gauge field. Let us note that the dynamical Peierls phase cannot be gauged away due to its inhomogeneity, and thus corresponds to a non-trivial 'Gauge' field. Alternatively, one can compute the Wilson loop operator around a square plaquette,

$$\begin{eqnarray}&&{W}_{\circlearrowleft}={U}_{{\bf{i}}+{{\bf{e}}}_{y},{\bf{i}}}{U}_{{\bf{i}}+{{\bf{e}}}_{x}+{{\bf{e}}}_{y},{\bf{i}}+{{\bf{e}}}_{y}}{U}_{{\bf{i}}+{{\bf{e}}}_{x},{\bf{i}}+{{\bf{e}}}_{x}+{{\bf{e}}}_{y}}{U}_{{\bf{i}},{\bf{i}}+{{\bf{e}}}_{x}}={{\rm{e}}}^{{\rm{i}}\theta ({n}_{{\bf{i}},\downarrow }-{n}_{{\bf{i}}+{{\bf{e}}}_{x},\downarrow }-{n}_{{\bf{i}}+{{\bf{e}}}_{y},\downarrow }+{n}_{{\bf{i}}+{{\bf{e}}}_{x}+{{\bf{e}}}_{y},\downarrow })}={{\rm{e}}}^{{\rm{i}}{\oint }_{\circlearrowleft}{\bf{A}}(x)\cdot {\rm{d}}{\bf{l}}}\end{eqnarray} \tag{ 81 }$$

and check that it is a non-trivial operator for $\theta \in (0,2\pi ).$ As the 'Gauge' fields commute at any position, one can regard the Hamiltonian (79) as a dynamical Abelian Hofstadter model. The quantum simulator will be able to explore this interesting model, and the fate of the Hofstadter quantum Hall phase, in the light of the aforementioned correlated topological insulators.

(ii) Time-reversal Hofstadter model in a dynamical gauge field.— We can now generalize the above construction to a time-reversal invariant situation, such as having two copies of the Hofstadter model subjected to anti-parallel magnetic fields. In our case, this is straightforward if we tune the relative phases as follows $2{\varphi }_{\downarrow ,x}=-2{\varphi }_{\uparrow ,x}=2\theta \ne 0.$ Moreover, we also set ${t}_{{\rm{eff}},\alpha }^{\uparrow }={t}_{{\rm{eff}},\alpha }^{\downarrow }=:\;{t}_{{\rm{eff}},\alpha }$ by controlling the moving optical lattices appropriately. Accordingly, the effective Hamiltonian (30) becomes a doubled Hofstadter model in the Landau gauge

$$\begin{eqnarray}{H}_{{\rm{eff}}} & = & -\displaystyle \sum _{{\bf{i}}}\left({t}_{{\rm{eff}},x}{{\rm{e}}}^{+{\rm{i}}\theta \left({n}_{{\bf{i}}+{{\bf{e}}}_{x},\downarrow }-{n}_{{\bf{i}},\downarrow }\right)}{f}_{{\bf{i}},\uparrow }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{x},\uparrow }^{\phantom{\dagger }}+{t}_{{\rm{eff}},y}{f}_{{\bf{i}},\uparrow }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{y},\uparrow }^{\phantom{\dagger }}+{\rm{H.c.}}\right)\\ & & -\displaystyle \sum _{{\bf{i}}}\left({t}_{{\rm{eff}},x}{{\rm{e}}}^{-{\rm{i}}\theta \left({n}_{{\bf{i}}+{{\bf{e}}}_{x},\uparrow }-{n}_{{\bf{i}},\uparrow }\right)}{f}_{{\bf{i}},\downarrow }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{x},\downarrow }^{\phantom{\dagger }}+{t}_{{\rm{eff}},y}{f}_{{\bf{i}},\downarrow }^{\dagger }{f}_{{\bf{i}}+{{\bf{e}}}_{y},\downarrow }^{\phantom{\dagger }}+{\rm{H.c.}}\right),\end{eqnarray} \tag{ 82 }$$

where atoms of any pseudospin play the role of the dynamical 'Gauge' field for the atoms with the remaining pseudospin, and the magnetic fluxes are opposite for each pseudospin. Accordingly, the time-reversal symmetry should invert the flux $\theta \to -\theta ,$ and flip the pseudospin ${f}_{{\bf{i}},\uparrow }\to {f}_{{\bf{i}},\downarrow }.$ Hence, the quantum simulator can explore the fate of the fate of the time-reversal Hofstadter quantum spin-Hall phase due to the presence of strong correlations, and dynamical effects of the Gauge field.