Artificial gauge fields and topology with ultracold atoms in optical lattices

The dynamics of charged particles in a periodic potential can be described using the tight-binding approximation. The description is based on the standard orbitals of the atoms, with the assumption that the overlap between wavefunctions of electrons on neighboring sites is small. This is a valid assumption for low temperatures, where all particles occupy the lowest energy level. The corresponding Hamiltonian takes the following form

$$\begin{eqnarray}&&{\hat{H}}_{0}=-J\displaystyle \sum _{m,n}\phantom{\rule{}{0.25ex}}({\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}+{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}}\phantom{\rule{}{0.25ex}}),\end{eqnarray} \tag{ 1 }$$

where ${\hat{a}}_{m,n}^{\dagger }$ and ${\hat{a}}_{m,n}$ are the creation and annihilation operators on site (m, n) respectively, m is the site index along x and n the one along y. Electrons in the lattice can tunnel between neighboring sites with strength J, which determines the kinetic energy of the system. These hopping matrix elements are modified according to the Peierls substitution [62], if an external magnetic field B = ∇ × A is applied, here A is the vector potential. As a result, tunneling is accompanied by phases ${\phi }_{m,n}^{i}=-{{qA}}_{m,n}^{i}/{\rm{\hslash }}$, i = {x, y}, commonly known as Peierls phases (figure 1) and the tight-binding Hamiltonian is modified accordingly as

$$\begin{eqnarray}\hat{H}= & -J\displaystyle \sum _{m,n}\phantom{\rule{}{0.25ex}}({{\rm{e}}}^{{\rm{i}}{\phi }_{m,n}^{x}}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}\\ & \qquad \quad +{{\rm{e}}}^{{\rm{i}}{\phi }_{m,n}^{y}}{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}}\phantom{\rule{}{0.25ex}}).\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

The Peierls phases are a manifestation of the Aharonov–Bohm phase experienced by a charged particle moving in a magnetic field

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{AB}}=-\displaystyle \frac{q}{{\rm{\hslash }}}{\oint }_{{ \mathcal C }}{\bf{A}}\cdot {\rm{d}}{\bf{r}}=-2\pi \,{{\rm{\Phi }}}_{B}/{{\rm{\Phi }}}_{0}=: -2\pi \alpha ,\end{eqnarray} \tag{ 3 }$$

where Φ_B is the magnetic flux through the area enclosed by the contour ${ \mathcal C }$ [63]. Equivalently one can define the magnetic flux per lattice unit cell in units of the magnetic flux quantum as

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{1}{2\pi }{\rm{\Phi }},\quad \mathrm{with}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\rm{\Phi }}=({\phi }_{m,n}^{x}+{\phi }_{m+1,n}^{y}-{\phi }_{m,n+1}^{x}-{\phi }_{m,n}^{y}).\end{eqnarray} \tag{ 5 }$$

For simplicity we will refer to Φ as the flux per unit cell of the lattice or the flux per plaquette. The reason for this will become apparent later.

In the presence of the vector potential A the discrete translation symmetry of the underlying lattice is no longer determined by the lattice translation operators. Even if the magnetic field B is homogeneous and invariant under discrete lattice translations, the vector potential giving rise to the complex tunnel matrix elements in (2) is generally not invariant under this translation. This prevents a direct application of the well-known Bloch theorem. Instead the new symmetries of the system are described by so-called MTOs [64–66], which are a combination of lattice translation and gauge transformation and can be expressed in the following form:

$$\begin{eqnarray}{\hat{T}}_{x}^{M} & = & \displaystyle \sum _{m,n}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}{{\rm{e}}}^{{\rm{i}}{\theta }_{m,n}^{x}},\\ {\hat{T}}_{y}^{M} & = & \displaystyle \sum _{m,n}{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}{{\rm{e}}}^{{\rm{i}}{\theta }_{m,n}^{y}}.\end{eqnarray} \tag{ 6 }$$

In order to find a complete set of commuting operators we require $[{\hat{T}}_{i}^{M},\hat{H}]=0$ and $[{\hat{T}}_{x}^{M},{\hat{T}}_{y}^{M}]=0$ [51, 67]. Let us consider a translation around a super-cell with dimensions k × l. For a homogeneous flux distribution one can show that the operators fulfill the following relation:

$$\begin{eqnarray}&&{{\rm{e}}}^{-{\rm{i}}{k}{l}{\rm{\Phi }}}{({\hat{T}}_{x}^{M})}^{k}{({\hat{T}}_{y}^{M})}^{l}={({\hat{T}}_{y}^{M})}^{l}{({\hat{T}}_{x}^{M})}^{k}.\end{eqnarray} \tag{ 7 }$$

Hence, for rational values of α = p/q ($p,q\in {\mathbb{Z}}$) the commutator $[{({\hat{T}}_{x}^{M})}^{k},{({\hat{T}}_{y}^{M})}^{l}\phantom{\rule{}{2.0ex}}]$ vanishes if

$$\begin{eqnarray}&&{kl}{\rm{\Phi }}=2\pi p\displaystyle \frac{{kl}}{q}\mathop{=}\limits^{!}2\pi \times \nu ,\quad \nu \in {\mathbb{Z}}.\end{eqnarray} \tag{ 8 }$$

This relation can be fulfilled by various super-cells, however, the smallest one is uniquely defined by kl = q. This cell is called magnetic unit cell and its area is determined by ${A}_{{\rm{MU}}}={{qa}}^{2}$, which is q times the area of the normal unit cell of the underlying lattice. It contains q non-equivalent sites, which gives rise to the fractal energy spectrum as a function of the applied field. The operators ${({\hat{T}}_{x}^{M})}^{k}\equiv {\hat{M}}_{x}^{k}$, ${({\hat{T}}_{y}^{M})}^{l}\equiv {\hat{M}}_{y}^{l}$ and $\hat{H}$ defined in equation (2) form a complete set of commuting operators and one can find simultaneous eigenstates Ψ_m,n by formulating a generalized Bloch theorem based on the magnetic translation symmetries:

$$\begin{eqnarray}{\hat{M}}_{x}^{k}\ {{\rm{\Psi }}}_{m,n} & = & {{\rm{e}}}^{{{\rm{i}}{k}}_{x}{ka}}\ {{\rm{\Psi }}}_{m,n},\\ {\hat{M}}_{y}^{l}\ {{\rm{\Psi }}}_{m,n} & = & {{\rm{e}}}^{{{\rm{i}}{k}}_{y}{la}}\ {{\rm{\Psi }}}_{m,n},\end{eqnarray} \tag{ 9 }$$

where k = (k_x, k_y) is defined within the first magnetic Brillouin zone (FBZ): −π/(ka) ≤ k_x < π/(ka), −π/(la) ≤ k_y < π/(la).

It is interesting to mention here, that the area of the magnetic unit cell is uniquely defined by α, its dimensions, however, are not, as illustrated in figure 2 for α = 1/4. A detailed computation of the MTOs for α = 1/4 and staggered flux distributions can be found in [51].

Zoom In Zoom Out Reset image size

The Peierls phases in Hamiltonian (2) depend on the particular gauge chosen for the description. However, physical observables do not depend on the gauge and one can always choose to work with a certain gauge that simplifies the calculations. In this chapter we will use the Landau gauge ϕ_m,n = (−Φ n, 0), where Hamiltonian (2) transforms to

$$\begin{eqnarray}\hat{H}= & -J\displaystyle \sum _{m,n}({{\rm{e}}}^{-{\rm{i}}{\rm{\Phi }}n}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}\\ & \qquad +{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}}),\end{eqnarray} \tag{ 10 }$$

commonly known as Harper–Hofstadter Hamiltonian [36–38]. In this gauge only tunneling along the x-direction is accompanied by Peierls phases, while hopping in the perpendicular direction remains real. The energy spectrum of this Hamiltonian exhibits a fractal self-similar structure as a function of the flux α, which is known as Hofstadter's butterfly [38] and is shown in figure 3(a).

Zoom In Zoom Out Reset image size

For rational values of the flux α = p/q the most convenient choice of the magnetic unit cell in the Landau gauge is the one oriented along the y-direction with dimensions (1 × q) · a², because the MTOs take the simple form

$$\begin{eqnarray}&&{\hat{M}}_{x}^{1}=\displaystyle \sum _{m,n}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n},\quad {\hat{M}}_{y}^{q}=\displaystyle \sum _{m,n}{\hat{a}}_{m,n+q}^{\dagger }{\hat{a}}_{m,n}.\end{eqnarray} \tag{ 11 }$$

They correspond to usual lattice translation operators, with the only difference that the translation along y is performed over an increased distance, namely q lattice sites. The magnetic unit cell contains a flux Φ_MU = p × 2π. In accordance with the generalized Bloch theorem (9), one can make the following ansatz to solve the Schrödinger equation (SE)

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{m,n}=\displaystyle \sum _{m,n}{{\rm{e}}}^{{{\rm{i}}{k}}_{x}{ma}}{{\rm{e}}}^{{{\rm{i}}{k}}_{y}{na}}{\psi }_{n},\quad \quad {\psi }_{n+q}={\psi }_{n}.\end{eqnarray} \tag{ 12 }$$

The wavefunction Ψ_m,n is expanded in single-particle on-site wavefunctions ψ_{i, j} ≡ ψ_{i, j} ∣i, j〉, and k_x, k_y are defined in the range $-\pi /a\leqslant {k}_{x}\lt \pi /a$ and $-\pi /({qa})\leqslant {k}_{y}\lt \pi /({qa})$. Indeed applying the MTOs defined in equation (11) one obtains

$$\begin{eqnarray}{\hat{M}}_{x}^{1}\ {{\rm{\Psi }}}_{m,n} & = & {{\rm{\Psi }}}_{m+1,n}={{\rm{e}}}^{{{\rm{i}}{k}}_{x}a}\ {{\rm{\Psi }}}_{m,n},\\ {\hat{M}}_{y}^{q}\ {{\rm{\Psi }}}_{m,n} & = & {{\rm{\Psi }}}_{m,n+q}={{\rm{e}}}^{{{\rm{i}}{k}}_{y}{qa}}\ {{\rm{\Psi }}}_{m,n},\end{eqnarray} \tag{ 13 }$$

justifying the wavefunction ansatz equation (12). Inserting the latter into the SE associated with Hamiltonian (10),

$$\begin{eqnarray}E{{\rm{\Psi }}}_{m,n} & = & -J({{\rm{e}}}^{-{\rm{i}}{\rm{\Phi }}n}\ {{\rm{\Psi }}}_{m+1,n}+{{\rm{e}}}^{{\rm{i}}{\rm{\Phi }}n}\ {{\rm{\Psi }}}_{m-1,n}\\ & & +{{\rm{\Psi }}}_{m,n+1}+{{\rm{\Psi }}}_{m,n-1}),\end{eqnarray} \tag{ 14 }$$

one obtains an infinite number of relations

$$\begin{eqnarray}E{\psi }_{n} & = & -J[2\cos ({k}_{x}a-{\rm{\Phi }}n){\psi }_{n}\\ & & +{{\rm{e}}}^{{{\rm{i}}{k}}_{y}a}{\psi }_{n+1}+{{\rm{e}}}^{-{{\rm{i}}{k}}_{y}a}{\psi }_{n-1}].\end{eqnarray} \tag{ 15 }$$

These equations can be recast into a q-dimensional eigenvalue problem

$$\begin{eqnarray}&&E({\bf{k}})\left(\begin{array}{c}{\psi }_{0}\\ {\psi }_{1}\\ \vdots \\ {\psi }_{q-1}\end{array}\right)=H({\bf{k}})\left(\begin{array}{c}{\psi }_{0}\\ {\psi }_{1}\\ \vdots \\ {\psi }_{q-1}\end{array}\right),\end{eqnarray} \tag{ 16 }$$

where H(k) is a q × q matrix of the following form:

$$\begin{eqnarray}H({\bf{k}})=-J\left(\begin{array}{ccccc}{h}_{0} & {{\rm{e}}}^{{{\rm{i}}{k}}_{y}a} & 0 & \cdots & {{\rm{e}}}^{-{{\rm{i}}{k}}_{y}a}\\ {{\rm{e}}}^{-{{\rm{i}}{k}}_{y}a} & {h}_{1} & {{\rm{e}}}^{{{\rm{i}}{k}}_{y}a} & \cdots & 0\\ 0 & {{\rm{e}}}^{-{{\rm{i}}{k}}_{y}a} & {h}_{2} & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ {{\rm{e}}}^{{{\rm{i}}{k}}_{y}a} & 0 & 0 & \cdots & {h}_{q-1}\end{array}\right),\end{eqnarray} \tag{ 17 }$$

with ${h}_{q}=2\cos ({k}_{x}a-q{\rm{\Phi }})$. For $\alpha \in {\mathbb{Z}}$ there is no external magnetic field and one recovers a single energy band, with the tight-binding dispersion of a simple 2D square lattice

$$\begin{eqnarray}&&E({\bf{k}})=-2J\cos ({k}_{x}a)-2J\cos ({k}_{y}a).\end{eqnarray} \tag{ 18 }$$

The corresponding bandwidth is E_bw = 2 × 4J. A rational flux per plaquette α = p/q leads to a splitting of this band into q subbands, E_μ(k), μ = {1, ..., q } resulting in a fractal spectrum (figure 3(a)). Figures 3(b), (c) show two exemplary bandstructures calculated for two different values of α. For irrational fluxes one obtains an infinite number of energy levels, which form a Cantor set [38].

The Hofstadter model is an important paradigmatic model for investigating the intriguing properties of topological phases of matter [3, 5]. It has motivated many theoretical and experimental studies also in artificial quantum systems [28]. It explicitly breaks TR symmetry and supports a series of QH and potentially fractional QH states [68–71]. QH insulators are characterized by a quantized Hall conductance [1, 2], which is directly related to the Chern number [7] of the occupied energy bands.

The Hall conductance σ_H is typically measured by sending a constant current through the material and observing the voltage difference that builds up across the sample perpendicular to the current. For low enough temperatures all energy levels below the Fermi energy E_F are filled. Assuming that the Fermi energy lies between two bands the conductance is given by the sum of the contributions from all occupied bands E_μ < E_F

$$\begin{eqnarray}&&{\sigma }_{{\rm{H}}}=\displaystyle \frac{{e}^{2}}{h}\displaystyle \sum _{{E}_{\mu }\lt {E}_{{\rm{F}}}}{\nu }_{\mu },\end{eqnarray} \tag{ 19 }$$

where ν_μ is the Chern number of the μth band.

This is very interesting because QH systems are insulating in the bulk. Indeed the conductivity defined above results from the conducting edge modes, which are gapless and chiral. Each edge mode contributes one quantum of conductance e²/h, such that the magnitude of the Hall conductivity is only determined by the number of edge modes, which in turn is directly related to the Chern number of the bulk. Hence, the existence of edge states can be seen as a manifestation of the bulk topological order [7]. This correspondence explains the robustness of the Hall conductivity against small perturbations, since they do not change the topological properties of the energy bands determined by ν_μ

$$\begin{eqnarray}{\nu}_{\mu} = & \frac{{\rm{i}}}{2\pi}{\int}_{{\rm{FBZ}}}\left(\left\langle \left.\frac{\partial {u}_{\mu}({\bf{k}})}{\partial {k}_{x}}\right|\frac{\partial {u}_{\mu}({\bf{k}})}{\partial {k}_{y}}\right\rangle\right.\\&\qquad\qquad-\,\left.\left\langle\left. \frac{\partial {u}_{\mu}({\bf{k}})}{\partial {k}_{y}}\right|\frac{\partial {u}_{\mu}({\bf{k}})}{\partial {k}_{x}}\right\rangle \right){{\rm{d}}}^{2}k\\=& \frac{1}{2\pi}{ \int}_{{\rm{FBZ}}}{{\rm{\Omega}}}_{\mu}({\bf{k}}){{\rm{d}}}^{2}k.\end{eqnarray} \tag{ 20 }$$

Here $| {u}_{\mu }({\bf{k}})\rangle$ are the periodic eigenfunctions, which are solutions of the eigenvalue equation (16), Ω_μ(k) is the Berry curvature and the integral is performed over the first Brillouin zone (FBZ) [72]. The topological properties of the bulk are usually defined for an infinite system and its relation to the topological properties of the edge modes is commonly known as bulk-edge correspondence [8–10]. In section 7 we will discuss the first observation of a Chern number with charge-neutral particles [73].

Several properties of the Hofstadter model are worth mentioning here. Since the sub-bandstructure for a certain value of α emerges from a topologically trivial tight-binding band, the Chern number of the total band necessarily vanishes, ${\sum }_{\mu }{\nu }_{\mu }=0$. One can further show that the particle-hole symmetry inherent to the model, results in a symmetric distribution of Chern numbers around zero energy (figures 3(b), (c)). Their values can be computed analytically based on a Diophantine equation [74, 75] or numerically for non-degenerate and degenerate bands following the work by Fukui et al [76].

Time-periodic Hamiltonians $\hat{H}(t+T)=\hat{H}(t)$ can be treated using Floquet's theorem [93–95], which states that the evolution of the system after one period T = 2π/ω can be described by an effective time-independent Hamiltonian ${\hat{H}}_{F}$. The difficulty consists in finding an analytic expression for it and moreover in designing a time-dependent system whose effective Hamiltonian exhibits the desired properties. Within the scope of this tutorial we restrict ourselves to the high-frequency limit, where the frequency associated with T is large compared to the other energy scales, ℏω ≫ J, such that the Floquet Hamiltonian can be computed in a perturbative manner.

We start with a reminder on Floquet theory. In operator notation the time-dependent SE is

$$\begin{eqnarray}&&{\rm{i}}{\rm{\hslash }}\displaystyle \frac{\partial }{\partial t}\hat{U}(t,{t}_{0})=\hat{H}(t)\hat{U}(t,{t}_{0}),\qquad \hat{U}({t}_{0},{t}_{0})={\mathbb{1}},\end{eqnarray} \tag{ 21 }$$

where $\hat{U}(t,{t}_{0})$ is the unitary time-evolution operator, that evolves a state $| \psi ({t}_{0})\rangle$ at t = t₀ to $| \psi (t)\rangle =\hat{U}(t,{t}_{0})| \psi ({t}_{0})\rangle$. In general it can be expressed as

$$\begin{eqnarray}&&\hat{U}(t,{t}_{0})={{ \mathcal T }}_{t}\exp \left[-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}{\int }_{{t}_{0}}^{t}\hat{H}({t}^{{\prime} }){\rm{d}}{t}^{{\prime} }\right],\end{eqnarray} \tag{ 22 }$$

where ${{ \mathcal T }}_{t}$ is the time-ordering operator, i.e. a short notation of an infinite series of commutator relations. For periodic systems the evolution operator fulfills [96]

$$\begin{eqnarray}&&\hat{U}({nT},0)={[\hat{U}(T,0)]}^{n}={[\hat{U}(T)]}^{n},\mathrm{with}\ n\in {\mathbb{N}},\end{eqnarray} \tag{ 23 }$$

where $\hat{U}(T,0)\equiv \hat{U}(T)$ is the evolution operator over one period T. The long-time behavior of the system can be described stroboscopically with $\hat{U}(t)$ at times t = nT [93, 96–99] using an effective time-independent Floquet Hamiltonian ${\hat{H}}_{F}$

$$\begin{eqnarray}&&\hat{U}(T)=\exp \left[-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}T{\hat{H}}_{F}\right],\end{eqnarray} \tag{ 24 }$$

where ${\hat{H}}_{F}$ is a Hermitian matrix.

Note that this definition of ${\hat{H}}_{F}$ that describes the stroboscopic dynamics is not unique. There are several Hamiltonians that share the same spectral and topological properties, which only differ by a gauge transformation. The description is useful for a first analysis of the problem, however, for a complete discussion of experimental results the full-time dynamics has to be taken into account, since it may lead to additional oscillations in the observables [100–102]. This additional dynamics is sometimes called micro-motion. Its importance is nicely illustrated based on the example of Paul traps, where the motion of the particle in the time-dependent potential can be partitioned into a slow and a fast part, where the latter corresponds to the micro-motion, which oscillates at the frequency of the driving potential [103–105]. To treat the additional fast dynamics accurately one can use a more general form of Floquet's theorem, where the evolution operator is partitioned as [106]

$$\begin{eqnarray}&&\hat{U}({t}_{f},{t}_{{i}})={{\rm{e}}}^{-{\rm{i}}\hat{K}({t}_{f})}{{\rm{e}}}^{-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}{\hat{H}}_{F}({t}_{f}-{t}_{{i}})}{{\rm{e}}}^{{\rm{i}}\hat{K}({t}_{{i}})}.\end{eqnarray} \tag{ 25 }$$

Here t_i,f denote the initial and final time of the evolution and we require that ${\rm{\exp }}({\rm{i}}\hat{K}(t))={\rm{\exp }}({\rm{i}}\hat{K}(t+T))$ is a time-periodic unitary operator. It contains information about the micro-motion and was given the intuitive name initial and final kick-operator in Ref. [100]. As mentioned above, the Hamiltonian ${\hat{H}}_{F}$ and the operator $\hat{K}(t)$ can be calculated perturbatively in powers of 1/ω

$$\begin{eqnarray}&&{\hat{H}}_{F}=\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{{\omega }^{n}}{\hat{H}}_{F}^{(n)},\qquad \hat{K}=\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{1}{{\omega }^{n}}{\hat{K}}^{(n)}.\end{eqnarray} \tag{ 26 }$$

The strategy is to compute ${\hat{H}}_{F}^{(n)}$ as a function of ${\hat{K}}^{(1)},\,\ldots ,\,{\hat{K}}^{(n+1)}$ and to choose ${\hat{K}}^{(n+1)}$ iteratively such that ${\hat{H}}_{F}^{(n)}$ is time-independent. This assures that the effective Hamiltonian is time-independent in all orders n [100, 106].

It indeed works efficiently for Hamiltonians where the driving frequency is off-resonant compared to any intrinsic energy difference of the static Hamiltonian [40, 90, 91, 107–112]. This regime is fulfilled in lattice-shaking experiments. In particular lattice shaking was employed to simulate classical frustrated magnetism on a triangular lattice [91] based on the generation of staggered artificial magnetic fields. Later these ideas were extended to realize tunable Peierls phases in 1D lattices [90]. More recently off-resonant shaking was used to implement the topological Haldane model [40] and interaction-induced gauge fields [112]. For the experimental techniques discussed in this tutorial, however, the frequency ℏω is resonant with an energy scale in the static Hamiltonian [41–43, 73, 92, 113–116]. Thus, it diverges in the limit $\omega \to \infty$, which, however, is the basic requirement for our perturbative expansion in equation (26). In the following two subsections, we will introduce two complementary methods to treat this problem: the first one is based on a formalism introduced by Rahav et al (section 3.1.1) and the second one on the Magnus expansion (section 3.1.2).

3.1.1. Perturbative analysis: resonant driving

We start by writing the static Hamiltonian as

$$\begin{eqnarray}&&{\hat{H}}_{0}=\displaystyle \sum _{\beta \gamma }{\hat{P}}^{\beta }{\hat{H}}_{\beta \gamma }^{(0)}{\hat{P}}^{\gamma }+{\rm{\hslash }}\omega \displaystyle \sum _{\beta }\beta {\hat{P}}^{\beta },\quad \beta ,\gamma \in {\mathbb{Z}},\end{eqnarray} \tag{ 27 }$$

where ${\hat{P}}^{\beta }$ is a projection operator. The integer β labels orthogonal sectors, ${\hat{P}}^{\beta }{\hat{P}}^{\gamma }={\delta }_{\beta \gamma }{\hat{P}}^{\beta }$ and ${\sum }_{\beta }{\hat{P}}^{\beta }={\mathbb{1}}$ [102]. This notation is useful for the treatment of distinct static lattice potentials and the number of sectors β will depend on it. For instance, figure 4(a) illustrates a staggered superlattice potential, where the potential energy of every other site is increased by ℏω and hence, β = {0, 1}. We will provide more specific examples in section 3.3. In general, terms of the form ${\hat{H}}_{\beta \beta }$ describe additional potential energy offsets or equivalent diagonal terms, while ${\hat{H}}_{\beta \gamma }$ (β $\ne$ γ) couple different sectors. We assume that all terms that diverge with ω are contained in the diagonal components of the static Hamiltonian and all remaining ones do not diverge. Similar to the static terms we write the time-periodic part

$$\begin{eqnarray}\hat{V}(t) & = & \displaystyle \sum _{j\gt 0}[{\hat{V}}^{(j)}{{\rm{e}}}^{{\rm{i}}{j}\omega t}+{\hat{V}}^{(-j)}{{\rm{e}}}^{-{\rm{i}}{j}\omega t}],\\ {\hat{V}}^{(j)} & = & \displaystyle \sum _{\beta \gamma }{\hat{P}}^{\beta }{\hat{H}}_{\beta \gamma }^{(j)}{\hat{P}}^{\gamma }.\end{eqnarray} \tag{ 28 }$$

Zoom In Zoom Out Reset image size

To treat the diverging terms we apply a time-dependent unitary transformation to the total Hamiltonian $\hat{H}(t)\,={\hat{H}}_{0}+\hat{V}(t)$ [102]:

$$\begin{eqnarray}| \psi (t)\rangle \to | {\psi }^{{\prime} }(t)\rangle & = & \hat{R}(t)| \psi (t)\rangle ,\\ \mathrm{with}\quad \hat{R}(t) & = & \exp \left[{\rm{i}}\displaystyle \sum _{\alpha }\alpha \omega t{\hat{P}}^{\alpha }\right].\end{eqnarray} \tag{ 29 }$$

This results in the transformed Hamiltonian

$$\begin{eqnarray}&&\hat{{ \mathcal H }}(t)=\hat{R}\hat{H}(t){\hat{R}}^{\dagger }-{\rm{i}}{\rm{\hslash }}\hat{R}{{\rm{d}}}_{t}{\hat{R}}^{\dagger }=\displaystyle \sum _{j}\hat{{ \mathcal H }}(j){{\rm{e}}}^{{\rm{i}}{j}\omega t},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\hat{{ \mathcal H }}(j)=\displaystyle \sum _{\beta \gamma }{\hat{P}}^{\beta }{\hat{H}}_{\beta \gamma }^{(j-\beta +\gamma )}{\hat{P}}^{\gamma },\end{eqnarray} \tag{ 31 }$$

which is indeed time-periodic and does not contain any divergent terms that scale with ω. Moreover we can directly apply the formalism [100] introduced above. In accordance with equation (25) the full-time evolution operator of the transformed systems is

$$\begin{eqnarray}\hat{U}({t}_{f},{t}_{{\rm{i}}}) & = & {{\rm{e}}}^{-{\rm{i}}\hat{{ \mathcal M }}({t}_{f})}{{\rm{e}}}^{-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}{\hat{{ \mathcal H }}}_{F}({t}_{f}-{t}_{{\rm{i}}})}{{\rm{e}}}^{{\rm{i}}\hat{{ \mathcal M }}({t}_{{\rm{i}}})},\\ {{\rm{e}}}^{{\rm{i}}\hat{{ \mathcal M }}(t)} & \equiv & {{\rm{e}}}^{{\rm{i}}\hat{{ \mathcal K }}(t)}\hat{R}(t),\end{eqnarray} \tag{ 32 }$$

where from now on we refer to $\hat{{ \mathcal M }}(t)$ as the micro-motion operator, which is a combination of the kick-operator defined above and the divergent terms of the static Hamiltonian. In the transformed system the Floquet Hamiltonian ${\hat{{ \mathcal H }}}_{F}$ and the kick-operator $\hat{{ \mathcal K }}(t)$ can be expanded in orders of 1/ω and one obtains [100]

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{F}={\hat{{ \mathcal H }}}^{(0)}+\displaystyle \frac{1}{{\rm{\hslash }}\omega }\displaystyle \sum _{j\gt 0}\displaystyle \frac{1}{j}[{\hat{{ \mathcal H }}}^{(+j)},{\hat{{ \mathcal H }}}^{(-j)}],\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&\hat{{ \mathcal K }}(t)=\displaystyle \frac{1}{{\rm{i}}{\rm{\hslash }}\omega }\displaystyle \sum _{j\gt 0}[{\hat{{ \mathcal H }}}^{(+j)}{{\rm{e}}}^{{\rm{i}}{j}\omega t}-{\hat{{ \mathcal H }}}^{(-j)}{{\rm{e}}}^{-{\rm{i}}{j}\omega t}],\end{eqnarray} \tag{ 34 }$$

where terms ${ \mathcal O }({\omega }^{-2})$ have been neglected for simplicity. This is a powerful approach for computing the full-time evolution of periodic Hamiltonians.

3.1.2. Magnus expansion approach

Alternatively the effective Hamiltonian can be calculated based on the Magnus expansion [117–119], which is particularly useful for describing the stroboscopic long-time evolution at integer multiples of the modulation period T, t = nT with $n\in {\mathbb{N}}$. An instructive comparison of the two methods for a simple double-well potential can be found in [51]. Similar to the discussion above, the Magnus expansion is a perturbative treatment, that works well in the high-frequency limit. The Floquet Hamiltonian is written as an infinite sum ${\hat{H}}_{F}={\sum }_{n=0}^{\infty }{\hat{H}}_{F}^{(n)}$, where the lowest order is given by

$$\begin{eqnarray}&&{\hat{H}}_{F}^{(0)}=\displaystyle \frac{1}{T}{\int }_{0}^{T}\hat{H}(t){\rm{d}}t\end{eqnarray} \tag{ 35 }$$

and higher orders of the expansion scale as (1/ω)ⁿ. For $\omega \to \infty$ the time-independent Floquet Hamiltonian is well approximated by the lowest order ${\hat{H}}_{F}^{(0)}\simeq {\hat{H}}_{F}$. In this description it is easy to see that there is a family of Floquet Hamiltonians with similar properties, which can be obtained by changing the time interval chosen for the integration, e.g. t ∈ [t₀, t₀ + T]. All these results are related by gauge transformations [119] and share the same spectral and topological properties.

For Hamiltonians of the form (27), (28) we again have to treat the diverging terms separately. We perform a slightly modified unitary transformation that includes explicitly the time-dependent term $\hat{V}(t)$ [119, 120] using

$$\begin{eqnarray}&&{\hat{R}}_{M}(t)=\exp \left[{\rm{i}}\displaystyle \sum _{\beta }\beta \omega t{\hat{P}}^{\beta }+{\rm{i}}\displaystyle \frac{1}{{\rm{\hslash }}}{\int }^{t}\hat{V}({t}^{{\prime} }){\rm{d}}{t}^{{\prime} }\right].\end{eqnarray} \tag{ 36 }$$

The transformed Hamiltonian is denoted as ${\hat{{ \mathcal H }}}_{M}(t)$, which is again time-periodic and according to equation (35) the Floquet Hamiltonian can be approximated by

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{F}^{M}=\displaystyle \frac{1}{T}{\int }_{0}^{T}{\hat{{ \mathcal H }}}^{M}(t){\rm{d}}t.\end{eqnarray} \tag{ 37 }$$

Let us start by illustrating the basic operating principle of laser-assisted tunneling using a two-level system with states $| 0\rangle$ and $| 1\rangle$. The two levels are separated by a large energy Δ and the system can be described by the static Hamiltonian

$$\begin{eqnarray}&&\hat{H}=J(| 0\rangle \langle 1| +| 1\rangle \langle 0| )+{\rm{\Delta }}{\hat{P}}^{1},\quad \quad {\hat{P}}^{1}=| 1\rangle \langle 1| ,\end{eqnarray} \tag{ 38 }$$

where J denotes the coupling strength between the levels, which is small compared to Δ. Due to the non-zero coupling J there is finite mixing of the two levels and the eigenenergies are separated by

$$\begin{eqnarray}&&{E}_{\mathrm{gap}}=\sqrt{{{\rm{\Delta }}}^{2}+4{J}^{2}},\end{eqnarray} \tag{ 39 }$$

which sets the exact value for the resonance condition ℏω = E_gap, where the two eigenstates are resonantly coupled. In actual experiments this is usually automatically ensured by the calibration measurements. If this is not the case, the dynamics can be largely modified if the system is modulated at ℏω = Δ [51, 101].

In the limit of a large energy offset Δ ≫ J the equation simplifies and the resonance condition is approximately given by E_gap ≃ Δ = ℏω. The dynamics of the resonantly-driven system is then described by the time-periodic Hamiltonian

$$\begin{eqnarray}\hat{H}(t) & = & J(| 0\rangle \langle 1| +| 1\rangle \langle 0| )+{\rm{\hslash }}\omega {\hat{P}}^{1}\\ & & +{V}_{0}\cos (\omega t+\varphi ){\hat{P}}^{0},\end{eqnarray} \tag{ 40 }$$

with ${\hat{P}}^{\beta }=| \beta \rangle \langle \beta |$ the projectors on the two sectors β = {0, 1} and V₀ the modulation amplitude. Following the theoretical description presented in section 3.1.2 we transform the system into the rotating frame [119, 120] according to equation (36) with

$$\begin{eqnarray}&&{\hat{R}}_{M}(t)=\exp \left[{\rm{i}}\omega t{\hat{P}}^{1}+{\rm{i}}\displaystyle \frac{{V}_{0}}{{\rm{\hslash }}\omega }\sin (\omega t+\varphi ){\hat{P}}^{0}\right].\end{eqnarray} \tag{ 41 }$$

This transformation effectively introduces the time-dependence into the coupling term and we obtain

$$\begin{eqnarray}{\hat{{ \mathcal H }}}_{M}(t) & = & J| 0\rangle \langle 1| {{\rm{e}}}^{{\rm{i}}\eta (t)}+J| 1\rangle \langle 0| {{\rm{e}}}^{-{\rm{i}}\eta (t)},\\ \eta (t) & = & -\left[\omega t-\displaystyle \frac{{V}_{0}}{{\rm{\hslash }}\omega }\sin (\omega t+\varphi )\right].\end{eqnarray} \tag{ 42 }$$

The transformed Hamiltonian ${\hat{{ \mathcal H }}}_{M}(t)$ is time-periodic without any static components. To lowest order (equation (37)) the Floquet Hamiltonian can be written as

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{F}^{M}=\mathop{\underbrace{J{{ \mathcal J }}_{1}\left(\displaystyle \frac{{V}_{0}}{{\rm{\hslash }}\omega }\right)}}\limits_{=: {J}_{\mathrm{eff}}}{{\rm{e}}}^{{\rm{i}}\varphi }| 0\rangle \langle 1| +{\rm{h.c.}},\end{eqnarray} \tag{ 43 }$$

where ${{ \mathcal J }}_{1}(x)$ is the first-order Bessel function of the first kind. The coupling between the two levels $| 0\rangle$ and $| 1\rangle$ is resonantly restored, renormalized and accompanied by an additional term ${{\rm{e}}}^{{\rm{i}}\varphi }$. It can be interpreted as a Peierls phase and is well suited to realize artificial magnetic fields (section 3.4). One can further show that for systems described by Hamiltonians of the form (38), the micro-motion operator in the lab-frame is to lowest order determined by the operator ${\hat{R}}_{M}(t)$ in equation (41). Note, however, that this is not true in general.

The resonantly coupled two-level example can be straightforwardly extended to a more general situation, where coupling is restored for energy offsets that are equal to an integer multiple of the modulation frequency Δ = νℏω, with $\nu \in {\mathbb{Z}}$ [121, 122]. These situations are also termed multi-photon processes, since multiple quanta of the driving field are required to bridge the energy gap Δ. The associated Floquet Hamiltonian in the rotating frame is

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{F}^{M}=J{{ \mathcal J }}_{\nu }\left(\displaystyle \frac{{V}_{0}}{{\rm{\hslash }}\omega }\right)[{{\rm{e}}}^{{\rm{i}}\nu \varphi }| 0\rangle \langle 1| +{{\rm{e}}}^{-{\rm{i}}\nu \varphi }| 1\rangle \langle 0| ],\end{eqnarray} \tag{ 44 }$$

where ${{ \mathcal J }}_{\nu }(x)$ denotes the ν-order Bessel function of the first kind. While for zero-order processes (ν = 0) the couplings remain real and are renormalized by ${{ \mathcal J }}_{0}({V}_{0}/[{\rm{\hslash }}\omega ])$ [107], first-order processes (ν = 1) lead to complex tunneling matrix elements and are typically termed photon- or laser-assisted tunneling [113]. Note that the phase factors depend on the resonance conditions, as summarized in the table below.

Resonance cond.	Eff. coupling strength	Phase factor
ℏω = Δ	$J{{ \mathcal J }}_{1}({V}_{0}/[{\rm{\hslash }}\omega ])$	φ
ℏω = −Δ	$J{{ \mathcal J }}_{1}({V}_{0}/[{\rm{\hslash }}\omega ])\,$	−φ + π

In the previous section we have introduced the problem of a resonantly coupled two-level system from a stroboscopic point of view based on the effective Hamiltonian ${\hat{{ \mathcal H }}}_{F}$. In an actual experiment the micro-motion $\hat{{ \mathcal M }}(t)$ may, however, have a large impact on experimental observables. Here, we will study its properties by means of a few specific examples. Common laser-assisted tunneling schemes rely on potential energy offsets Δ_m along one of the two lattice axes to inhibit tunneling, which can then be restored resonantly with the modulation.

We start our discussion with resonantly-driven 1D lattices and then extend it to the 2D situation in the next section. In particular we present two examples illustrated in figure 4, which are described by the tight-binding Hamiltonian

$$\begin{eqnarray}&&\hat{H}(t)=-{J}_{x}\displaystyle \sum _{m}({\hat{a}}_{m+1}^{\dagger }{\hat{a}}_{m}+{\rm{h.c.}})\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&+\,\displaystyle \sum _{m}({V}_{0}\sin (\omega t+{\varphi }_{m})+{{\rm{\Delta }}}_{m}){\hat{n}}_{m},\end{eqnarray} \tag{ 46 }$$

where φ_m is the site-dependent phase of the time-dependent potential and J_x is the tunnel coupling along x. In the following we consider only potentials with Δ_m+1 − Δ_m = ±Δ = ±ℏω but the model can be straightforwardly extended to more general energy landscapes [102]. Hamiltonian (46) is indeed time-periodic $\hat{H}(t+T)=\hat{H}(t)$ and of the form defined in equations (27), (28). The static divergent term can be written as

$$\begin{eqnarray}&&\hat{S}=\displaystyle \sum _{m}{{\rm{\Delta }}}_{m}{\hat{n}}_{m}={\rm{\hslash }}\omega \displaystyle \sum _{\beta }\beta {\hat{N}}_{\beta },{\hat{N}}_{\beta }=\displaystyle \sum _{m\in \beta }{\hat{n}}_{m},\end{eqnarray} \tag{ 47 }$$

where ${\hat{N}}_{\beta }$ is the number operator for the sector β and $\beta \in {\mathbb{Z}}$. Each sector is composed of several sites m, which are defined by the specific form of the on-site potential Δ_m. The staggered superlattice potential (figure 4(a)) separates into two different sectors β = {0, 1}, where every other lattice site belongs to the same sector. The Wannier–Stark ladder (figure 4(b)) instead consists of an infinite number of sectors, each of them containing only a single site. Independent of β the number operator ${\hat{N}}_{\beta }$ can be expanded in terms of projectors

$$\begin{eqnarray}&&{\hat{N}}_{\beta }=\displaystyle \sum _{{n}_{\beta }}{n}_{\beta }{\hat{P}}_{{n}_{\beta }},\quad \quad {\hat{N}}_{\beta }| {n}_{\beta }\rangle ={n}_{\beta }| {n}_{\beta }\rangle .\end{eqnarray} \tag{ 48 }$$

Using this equation and equation (47) one obtains

$$\begin{eqnarray}&&\hat{S}={\rm{\hslash }}\omega \displaystyle \sum _{\beta ,{n}_{\beta }}(\beta {n}_{\beta }){\hat{P}}_{{n}_{\beta }}.\end{eqnarray} \tag{ 49 }$$

Following the theoretical description presented in section 3.1 we find the micro-motion operator

$$\begin{eqnarray}\hat{{ \mathcal M }}(t) & \simeq & \displaystyle \sum _{m}\left(-\displaystyle \frac{{V}_{0}}{{\rm{\hslash }}\omega }\cos (\omega t+{\varphi }_{m})+\displaystyle \frac{{{\rm{\Delta }}}_{m}t}{{\rm{\hslash }}}\right){\hat{n}}_{m}\\ & =: & \displaystyle \sum _{m}{\chi }_{m}(t)\,{\hat{n}}_{m},\end{eqnarray} \tag{ 50 }$$

where we have applied the unitary transformation defined in equation (36) with the approximation $\exp [{\rm{i}}{ \mathcal M }(t)]\simeq \hat{R}(t)$. The transformed Hamiltonian reads

$$\begin{eqnarray}&&\hat{{ \mathcal H }}(t)=-{J}_{x}\displaystyle \sum _{m}({{\rm{e}}}^{{\rm{i}}{\eta }_{m}(t)}{\hat{a}}_{m+1}^{\dagger }{\hat{a}}_{m}+{\rm{h.c.}}),\end{eqnarray} \tag{ 51 }$$

with η_m(t) = χ_m+1(t) − χ_m(t) given by

$$\begin{eqnarray}{\eta }_{m}(t)\,=: & - & {\eta }_{0}\sin \left(\omega t+\displaystyle \frac{{\varphi }_{m+1}+{\varphi }_{m}}{2}\right)\\ & + & ({{\rm{\Delta }}}_{m+1}-{{\rm{\Delta }}}_{m})t/{\rm{\hslash }}.\end{eqnarray} \tag{ 52 }$$

Here ${\eta }_{0}=2{V}_{0}/({\rm{\hslash }}\omega )\sin [\delta \varphi /2]$ is the differential modulation amplitude and δφ = φ_m+1 − φ_m.

Looking at equation (50) we find that $\hat{{ \mathcal M }}(t)$ is proportional to the number operator ${\hat{n}}_{m}$. As a consequence it commutes with any terms of the effective Hamiltonian proportional to ${\hat{n}}_{m}$ and hence the density distribution of the eigenstates remains unaffected during one cycle in contrast to the quasi-momentum distribution [116], which is a typical observable probed with Bose–Einstein condensates (BEC) after standard time-of-flight (TOF) imaging [29]. The micro-motion operator consists of two distinct parts, one associated with the time-dependent potential proportional to γ = V₀/(ℏω) and a second one introduced by the static potential Δ_m. The first one gives rise to additional momentum components at positions related to the phase distribution φ_m. The amplitudes of these components oscillate periodically in time but their position remains unchanged. Their contribution to the signal vanishes in the high-frequency limit $\gamma \to 0$. The second term depends highly on the shape of the static potential. For the staggered superlattice ${{\rm{\Delta }}}_{m}={(-1)}^{m}{\rm{\Delta }}/2$ (figure 4(a)) additional momentum components appear at ±π/a, which oscillate in amplitude with period T = h/Δ but no additional components are populated during one cycle. This dynamics persists even in the limit $\gamma \to 0$. For the Wannier–Stark ladder Δ_m = mΔ (figure 4(b)) the corresponding term is associated with a constant drift in momentum space, δk(t) = ωt/a with constant amplitudes [116]. Consequently, all laser-assisted tunneling techniques based on Wanner-Stark ladders exhibit such a drift, also known as Bloch oscillations [123].

The complex couplings in 1D are mere gauge transformations and the system remains trivial. Extending those settings with a second perpendicular lattice potential enables the engineering of more interesting situations. We consider the following time-dependent Hamiltonian

$$\begin{eqnarray}\hat{H}(t) & = & \displaystyle \sum _{m,n}(-{J}_{x}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}-{J}_{y}{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}})\\ & & +\displaystyle \sum _{m,n}({V}_{0}\sin (\omega t+{\varphi }_{m,n})+{{\rm{\Delta }}}_{m}){\hat{n}}_{m,n},\end{eqnarray} \tag{ 53 }$$

where J_y is the tunnel coupling along the perpendicular axis and tunneling is inhibited along the x-axis for Δ ≫ J_x and V₀ = 0. Note, that the modulation potential involves a phase distribution φ_m,n that varies along both axes. Analog to the 1D case we find the micro-motion operator

$$\begin{eqnarray}\hat{{ \mathcal M }}(t) & = & \displaystyle \sum _{m,n}\left(-\displaystyle \frac{{V}_{0}}{{\rm{\hslash }}\omega }\cos (\omega t+{\varphi }_{m,n})+\displaystyle \frac{{{\rm{\Delta }}}_{m}t}{{\rm{\hslash }}}\right){\hat{n}}_{m,n}\\ & =: & \displaystyle \sum _{m,n}{\chi }_{m,n}(t)\ {\hat{n}}_{m,n}\end{eqnarray} \tag{ 54 }$$

and the transformed Hamiltonian

$$\begin{eqnarray}\hat{{ \mathcal H }}(t) & = & -{J}_{x}\displaystyle \sum _{m,n}({{\rm{e}}}^{{\rm{i}}{\eta }_{m,n}^{x}(t)}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}})\\ & & -{J}_{y}\displaystyle \sum _{m,n}({{\rm{e}}}^{{\rm{i}}{\eta }_{m,n}^{y}(t)}{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}}).\end{eqnarray} \tag{ 55 }$$

Indeed the periodic modulation generally affects tunneling along both directions. The functions ${\eta }_{m,n}^{x}(t)={\chi }_{m+1,n}(t)\,-{\chi }_{m,n}(t)$ and ${\eta }_{m,n}^{y}(t)={\chi }_{m,n+1}(t)-{\chi }_{m,n}(t)$ are determined by

$$\begin{eqnarray}{\eta }_{m,n}^{x}(t) & = & -{\eta }_{0x}\sin \left(\omega t+\displaystyle \frac{{\varphi }_{m+1,n}+{\varphi }_{m,n}}{2}\right)\\ & & +({{\rm{\Delta }}}_{m+1}-{{\rm{\Delta }}}_{m})\ t/{\rm{\hslash }},\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&{\eta }_{m,n}^{y}(t)=-{\eta }_{0y}\sin \left(\omega t+\displaystyle \frac{{\varphi }_{m,n+1}+{\varphi }_{m,n}}{2}\right),\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}&&{\eta }_{0i}=2\gamma \sin (\delta {\varphi }_{i}/2).\end{eqnarray} \tag{ 58 }$$

For $| {{\rm{\Delta }}}_{m+1}-{{\rm{\Delta }}}_{m}| ={\rm{\Delta }}$ tunneling can be resonantly restored with ℏω = Δ (laser-assisted tunneling) and one obtains the effective time-independent Hamiltonian

$$\begin{eqnarray}{\hat{{ \mathcal H }}}_{F} & = & -\mathop{\underbrace{{J}_{x}{{ \mathcal J }}_{1}({\eta }_{0x})}}\limits_{=: K}\displaystyle \sum _{m,n}{{\rm{e}}}^{{\rm{i}}{\phi }_{m,n}}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}\\ & & -\mathop{\underbrace{{J}_{y}{{ \mathcal J }}_{0}({\eta }_{0y})}}\limits_{=: J}\displaystyle \sum _{m,n}{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}},\end{eqnarray} \tag{ 59 }$$

where the Peierls phases ϕ_m,n ∝ (φ_m+1,n + φ_m,n)/2 depend on the distribution φ_m,n and on the on-site potential Δ_m (section 3.2). Along y tunneling is renormalized and in the high-frequency limit ($\gamma \to 0$) one finds

$$\begin{eqnarray}&&\displaystyle \frac{J}{{J}_{y}}=1-\displaystyle \frac{1}{4}{\eta }_{0y}^{2}+{ \mathcal O }({\eta }_{0y}^{4})\simeq 1-{\gamma }^{2}{\sin }^{2}(\delta {\varphi }_{y}/2),\end{eqnarray} \tag{ 60 }$$

which means it decreases with increasing γ. How strongly it decreases depends on the phase difference of the modulation between neighboring sites δφ_y = φ_m,n+1 − φ_m,n (figure 5(a)). For in-phase modulation the relative phase difference vanishes (δφ_y = 0) and the bare tunneling remains unchanged.

Zoom In Zoom Out Reset image size

In the perpendicular direction (x-axis) tunneling is resonantly restored, resulting in an effective tunnel coupling $K{{\rm{e}}}^{{\rm{i}}{\phi }_{m,n}}$. Its amplitude again depends on the relative phase difference δφ_x = φ_m+1,n − φ_m,n (figure 5(b)) and it is proportional to the driving amplitude γ. In the limit of large ω ($\gamma \to 0$) it can be approximated by

$$\begin{eqnarray}&&\displaystyle \frac{K}{{J}_{x}}=\displaystyle \frac{1}{2}{\eta }_{0x}+{ \mathcal O }({\eta }_{0x}^{3})\simeq \gamma \sin (\delta {\varphi }_{x}/2).\end{eqnarray} \tag{ 61 }$$

Again for in-phase modulation (δφ_x = 0) only a global energy shift is modulated in time, which cannot be used to restore tunneling. However, this is precisely how one can engineer more complex modulation geometries, which enable a local control of laser-assisted tunneling in the lattice; further details are discussed in [51].

The Peierls phases ϕ_m,n in Hamiltonian (59) appear only in the coupling terms with restored tunneling (x-axis). Note, that they give rise to a flux Φ_m,n = ϕ_m,n − ϕ_m,n+1 per plaquette (defined along +${\hat{{\bf{e}}}}_{z}$), whose value depends on the phase evolution in the perpendicular direction δφ_y, while the amplitude of the coupling K depends on δφ_x. The laser-assisted-tunneling scheme presented here naturally generates artificial fluxes on the order of ${\rm{\Phi }}={ \mathcal O }(2\pi )$ and therefore constitutes a good candidate to enter the strong-field regime required to study the Hofstadter model. Related techniques have been further identified, e.g. for ions or superconducting circuits [28].

The 2D lattice potential consists of a simple lattice along y generated at wavelength λ_s and a staggered superlattice potential along x with a staggered potential energy offset ±Δ/2 as illustrated in (figures 6(b), (c)). In the following we use the convention that high-energy sites are even lattice sites. For large energy offsets Δ ≫ J_x the dynamics along x is frozen, while atoms tunnel freely along y. The time-dependence is introduced with an additional pair of far-detuned running-wave beams with wavevectors $| {{\bf{k}}}_{1}| \simeq | {{\bf{k}}}_{2}| \equiv {k}_{R}$ and frequencies ω_1,2. As shown in figure 7(a), each of them is aligned along one of the main lattice axes. Tunneling along x is resonantly induced for ω = ω₁ − ω₂ = ±Δ/ℏ. The two beams interfere and create a local time-dependent potential of the form

$$\begin{eqnarray}&&{V}_{K}({\bf{r}},t)={V}_{0}\cos (\delta {\bf{k}}\cdot {\bf{r}}+\omega t+{\phi }_{0}),\end{eqnarray} \tag{ 64 }$$

with δk = k₁ − k₂ and ϕ₀ = ϕ₁ − ϕ₂. Strictly speaking there is an additional static term, which we have neglected because it only leads to a global energy shift. The relevant terms in V_K(r, t) are the site-dependent phases given by δk · R = k_Ra (m − n), where ${\bf{R}}={ma}{\hat{{\bf{e}}}}_{x}+{na}{\hat{{\bf{e}}}}_{y}$ defines the position in the lattice and ${\hat{{\bf{e}}}}_{x,y}$ are the unit vectors. For the implementation reported in [92] k_R = π/(2a) and the phases are given by

$$\begin{eqnarray}&&\delta {\bf{k}}\cdot {\bf{R}}=\displaystyle \frac{\pi }{2}(m-n).\end{eqnarray} \tag{ 65 }$$

Zoom In Zoom Out Reset image size

Along both lattice axes the phase advances by δφ_x = −δφ_y = π/2 as shown in figure 7(b). This value can be easily modified by changing the geometry of the running waves or a different ratio k_R/a, such that any value of the flux becomes accessible. In equation (64) there appears an additional phase offset ϕ₀, which is the relative phase between the two beams. In the experiment it is not actively stabilized and varies randomly for each realization but it has not influence on the distribution of artificial fluxes.

The dynamics of the system is described by a time-periodic Hamiltonian $\hat{H}(t+T)=\hat{H}(t)$ of the form (53) with Δ_m = (−1)^m Δ/2 and ${\varphi }_{m,n}=\tfrac{\pi }{2}(m-n+1)+{\phi }_{0}$ being the site-dependent phases generated by the running-wave beams (figures 7(a), (b)). As discussed in section 3.1 in the high-frequency limit ℏω ≫ J_x, J_y the stroboscopic evolution of the system can be described by an effective time-independent Hamiltonian of the form

$$\begin{eqnarray}{\hat{{ \mathcal H }}}_{\mathrm{eff}}\,= & & -K\displaystyle \sum _{m,n}{{\rm{e}}}^{{\rm{i}}{\phi }_{m,n}}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}\\ & & -J\displaystyle \sum _{m,n}{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}},\end{eqnarray} \tag{ 66 }$$

with $K={J}_{x}{{ \mathcal J }}_{1}({\eta }_{0})$ and $J={J}_{y}{{ \mathcal J }}_{0}({\eta }_{0})$. Due to the staggered energy offset the phase distribution ϕ_m,n is non-uniform (section 3.2) and varies according to

$$\begin{eqnarray}{\phi }_{m,n}=\left\{\begin{array}{ll}-({\varphi }_{m+1,n}+{\varphi }_{m,n})/2, & \mathrm{for}\,m\,\mathrm{odd}\\ +({\varphi }_{m+1,n}+{\varphi }_{m,n})/2+\pi , & \mathrm{for}\,m\,\mathrm{even}.\end{array}\right.\end{eqnarray} \tag{ 67 }$$

The distribution depends on the chosen laser-assisted-tunneling geometry and for the setting depicted in figure 7(a) the on-site phases are ${\varphi }_{m,n}=\tfrac{\pi }{2}(m-n+1)+{\phi }_{0}$ and one finds the corresponding Peierls phases

$$\begin{eqnarray}{\phi }_{m,n}=\left\{\begin{array}{ll}-\displaystyle \frac{\pi }{2}(m-n+3/2)-{\phi }_{0}, & \mathrm{for}\,m\,\mathrm{odd}\\ +\displaystyle \frac{\pi }{2}(m-n+7/2)+{\phi }_{0}, & \mathrm{for}\,m\,\mathrm{even}.\end{array}\right.\end{eqnarray} \tag{ 68 }$$

For simplicity we set ϕ₀ = −3π/4 and rewrite the expression for the phases as

$$\begin{eqnarray}&&{\phi }_{m,n}=\displaystyle \frac{\pi }{2}{(m+n)(-1)}^{m+1}.\end{eqnarray} \tag{ 69 }$$

This distribution gives rise to a staggered flux distribution ${{\rm{\Phi }}}_{m,n}=\tfrac{\pi }{2}{(-1)}^{m+1}$ with zero mean. It is uniform along y and staggered along x [92, 114].

The micro-motion operator required to understand the full-time evolution of the system is defined in equation (54). It leads to a very similar behavior as in the case of a periodically modulated 1D staggered superlattice potential (section 3.3). The part related to the staggered on-site potential Δ_m is trivial because it only leads to oscillations in the amplitudes of the individual momentum components. The modulation-induced term, however, gives rise to additional momentum components shifted by multiples of δk = (π/(2a), π/(2a)) and influences the measured momentum distributions presented in the following section.

In order to study the ground-state properties of the effective Hamiltonian defined in equation (66) we prepared a BEC in the lowest band of a 2D optical lattice with J_x/h = J_y/h = 31(2) Hz, which results in a 2D array of coupled 1D Bose gases, as reported in [92]. Then we ramped up the staggered potential energy offset Δ/h = 4.4(1) kHz within 0.7 ms. Subsequently, the lattice along x was decreased and the modulation was turned on to induce resonant tunneling along that axis with an effective coupling strength K/h = 30(2) Hz.

This is a highly non-adiabatic experimental protocol which was optimized experimentally. We found that after about 10 ms the phase coherence was restored, most likely because of entropy redistribution along the longitudinal modes. The quasi-momentum distribution was then observed after suddenly releasing the atoms and 20 ms TOF (figure 8(a)).

Zoom In Zoom Out Reset image size

Performing the same measurements with a simple cubic lattice and Φ = 0 would result in a single momentum component in the FBZ at k_x = k_y = 0 and additional ones at multiples of the reciprocal lattice vectors (2k_s, 0) and (0, 2k_s), where k_s = π/a. For the staggered flux lattice instead four momentum components per cubic lattice BZ were found. The observed distribution can be understood by solving the SE using the following ansatz for the wavefunction

$$\begin{eqnarray}{{\rm{\Psi }}}_{m,n}=\displaystyle \sum _{m,n}{{\rm{e}}}^{{{\rm{i}}{k}}_{x}{ma}}{{\rm{e}}}^{{{\rm{i}}{k}}_{y}{na}}\left\{\begin{array}{ll}{\psi }_{o}\ {{\rm{e}}}^{i\displaystyle \frac{\pi }{2}(m+n)}, & m\,\mathrm{odd}\\ {\psi }_{{\rm{e}}}, & m\,\mathrm{even}.\end{array}\right.\end{eqnarray} \tag{ 70 }$$

Note that, due to the magnetic translation symmetries, the magnetic unit cell is two times larger than the usual lattice unit cell with dimensions 2a × a. Using the above ansatz we can reduce the problem to a 2D eigenvalue equation ${\hat{H}}_{s}({\psi }_{e},{\psi }_{0})=E({\bf{k}})({\psi }_{e},{\psi }_{0})$ with

$$\begin{eqnarray}{\hat{H}}_{s}=\left(\begin{array}{cc}-2J{\rm{\cos }}({k}_{y}a) & -K\ ({\rm{i}}{{\rm{e}}}^{{{\rm{i}}{k}}_{x}a}+{{\rm{e}}}^{-{{\rm{i}}{k}}_{x}a})\\ -K\ (-{\rm{i}}{{\rm{e}}}^{-{{\rm{i}}{k}}_{x}a}+{{\rm{e}}}^{{{\rm{i}}{k}}_{x}a}) & 2J\sin ({k}_{y}a)\end{array}\right),\end{eqnarray} \tag{ 71 }$$

where k_x, k_y are defined within the FBZ, −k_s/2 ≤ k_x < k_s/2 and −k_s ≤ k_y < k_s, and ψ_e,o are the amplitudes on even and odd lattice sites respectively.

For equal couplings J/K = 1 a single energy minimum in the dispersion relation is found at k_e = (−k_s/4, −k_s/4). The corresponding ground-state wavefunction has equal amplitudes on all sites, $| {\psi }_{e}| =| {\psi }_{o}|$. According to the intrinsic structure of the eigenstates defined in equation (70) two momentum components are expected: one at k_e and an additional one at k_o = k_e + δk, with δk = (k_s/2, k_s/2). The second one is due to the additional phase factor of the wavefunction in equation (71) for odd lattice sites. The measured TOF images further show momentum components separated by the reciprocal magnetic lattice vectors (k_s, 0) and (0, 2k_s) as displayed in figure 8(a). The observed momenta are in agreement with theory.

The remaining question concerns the effect of the micro-motion operator on the measured distributions, since so far we have only discussed the results in comparison with the Floquet Hamiltonian (66) in the rotating frame. The term that may lead to additional momentum components (section 3.3) is proportional to V₀/(ℏω), which for our experimental parameters is not negligible V₀/(ℏω) ≃ 0.48. Indeed as can be seen in figure 8(a) there are two weak additional quasi-momenta in the FBZ located at (−k_s/4, 3k_s/4) and (k_s/4, − 3k_s/4) in agreement the micro-motion operator defined in equation (54). Looking at the full-time dynamics one finds that there are small oscillations in the amplitudes of the individual quasi-momenta but their position remains constant. Note that it is important to express all the operators using the gauge realized in the experiment because the TOF images depend on that [128, 134, 135]. For a different choice of the vector potential the momentum components would appear at different locations in the FBZ [51]. Intuitively this can be understood as follows: in TOF measurements all potentials are rapidly switched off and the system evolves according to the free-particle Hamiltonian. After long enough expansion times the images provide the Fourier transform of the condensate wavefunction. Phases that are imprinted on the condensate wavefunction, which is done in order to realize a particular vector potential, will directly show up as Fourier components in the TOF images, for instance, a constant phase difference between neighboring sites translates into a momentum shift of the distribution measured after TOF and hence, the particular gauge realized in the experiment will be visible in the TOF distributions.

Varying the coupling ratio J/K we find that there is a single minimum in the dispersion up to a critical value ${(J/K)}_{c}=\sqrt{2}$, above which the system exhibits a two-fold degenerate ground state (figure 8(b)). In future experiments it would be interesting to study the ground state of weakly-interacting bosonic atoms in this lattice and try to detect the incommensurate charge density wave that is expected to form [128].

Superlattices enable powerful local measurements even without single-site resolution [136–141]. In a 2D superlattice (figure 6) tunneling can be inhibited on every other bond along both lattice axes (figure 9(a)). This restricts the dynamics to a minimal lattice that consists of only four sites ${ \mathcal R }=\{A,B,C,D\}$. Using an additional vertical lattice at λ_z a 3D array of isolated four-site square plaquettes is realized. By preparing a single atom inside one such plaquette potential one can study the influence of the artificial flux on its dynamics and extract the amplitude of the artificial field.

Zoom In Zoom Out Reset image size

The measurements are performed in a plaquette potential which has an energy offset Δ along x and no offset along y, such that tunneling is inhibited along the x-axis by Δ. It is then restored resonantly with the running-wave setup described in section 5.1. The sign of the flux in each plaquette depends on the sign of the energy offset Δ and since only every other plaquette contributes to the measured signal each of them exhibits the same flux because each of them has a potential energy difference of the same sign.

The experimental sequence started by preparing a single atom in the ground state of the tilted plaquette $| {\psi }_{1}\rangle \,=(| A\rangle +| D\rangle )/\sqrt{2}$ with J_x/h = 1.2(1) kHz, J_y/h = 0.5(2) kHz and Δ/h = 4.9(1) kHz, as reported in [92]. To initiate the dynamics we suddenly turned on the modulation with V₀/(ℏω) ≃ 0.34, K/h = 0.28(1) kHz and J ≃ J_y. The following evolution is governed by the effective Hamiltonian, which can be expressed as a 4 × 4 matrix ${\hat{H}}_{P}$. We choose the basis $\{| A\rangle ,| B\rangle ,| C\rangle ,| D\rangle \}$ and obtain

$$\begin{eqnarray}{\hat{H}}_{P}=-\left(\begin{array}{cccc}0 & K & 0 & J\\ K & 0 & J & 0\\ 0 & J & 0 & K{{\rm{e}}}^{\mp {\rm{i}}{\rm{\Phi }}}\\ J & 0 & K{{\rm{e}}}^{\pm {\rm{i}}{\rm{\Phi }}} & 0\end{array}\right),\end{eqnarray} \tag{ 72 }$$

with Φ = π/2.

Experimentally we determine the real-space dynamics by measuring the site occupations in the plaquette ${N}_{{ \mathcal R }}$ as a function of time using a site-resolved band-mapping technique [125, 126]. From this we compute the mean position of the atom in the plaquette

$$\begin{eqnarray}&&\langle X\rangle =(-{N}_{A}+{N}_{B}+{N}_{C}-{N}_{D})\displaystyle \frac{a}{2N},\end{eqnarray} \tag{ 73 }$$

$$\begin{eqnarray}&&\langle Y\rangle =(-{N}_{A}-{N}_{B}+{N}_{C}+{N}_{D})\displaystyle \frac{a}{2N},\end{eqnarray} \tag{ 74 }$$

where N is the total atom number (figure 9(b)). There is an important difference between these measurements compared to conventional high-resolution detection techniques [136–141]. The measurements are performed simultaneously in a 3D array of plaquettes and one obtains averaged observables without true local information. The evolution starts at $\langle X\rangle (t=0)=-0.5$ and $\langle Y\rangle (t=0)=0$ with all atoms delocalized over the left bond. Shortly after tunneling is restored the atoms tunnel to the right (B- and C-sites) and without flux one would expect simple Rabi oscillations between the two bonds. Instead the atoms are deflected in the perpendicular direction (y-axis) compared to their initial tunneling direction, something that is characteristic of Lorentz-force-type dynamics. Moreover, the complete evolution appears to follow a path that is reminiscent of a classical cyclotron orbit (figure 9(b)). Using these results we can extract the value of the flux by solving the time-dependent SE numerically based on the Hamiltonian (72) and fitting it to the data. The parameters J and K have been calibrated independently in separate measurements such that the flux Φ was the only free fit parameter. We obtained the experimental value Φ_exp = 0.73(5) × π/2, which is slightly smaller compared to the ideal prediction Φ = π/2. The deviation is due to a reduced distance between lattice sites when partitioning the lattice into plaquettes and does not affect the value of the flux in the fully-coupled 2D lattice [92].

There are a few additional things that should be taken into account when studying the dynamics on isolated plaquette potentials [51], which we briefly mention. For our experimental parameters Δ/J_x ≃ 4, which means that the high-frequency limit is not well fulfilled and higher-order corrections that scale with J_x/Δ (section 3.2) might be relevant already on the level of the resonance condition, equation (39). Furthermore the dynamics (figure 9(b)) was measured non-stroboscopically, which requires additional justification for finite Δ/J_x, and the role of the initial phase of the driving was not discussed. Both effects are described by the micro-motion operator (section 3.2) determined by

$$\begin{eqnarray}&&{\hat{{ \mathcal M }}}_{P}(t)=\displaystyle \sum _{{ \mathcal R }}\left[\displaystyle \frac{{V}_{0}}{{\rm{\hslash }}\omega }\sin (\omega t+{\varphi }_{{ \mathcal R }})+\displaystyle \frac{{{\rm{\Delta }}}_{{ \mathcal R }}}{{\rm{\hslash }}}t\right]{\hat{n}}_{{ \mathcal R }}.\end{eqnarray} \tag{ 75 }$$

Here ${\varphi }_{{ \mathcal R }}$ is the phase of the modulation on site ${ \mathcal R }$, ${{\rm{\Delta }}}_{{ \mathcal R }}$ is the on-site energy and ${\hat{n}}_{{ \mathcal R }}$ is the number operator. The micro-motion gives rise to an initial kick ${\hat{{ \mathcal M }}}_{P}(0)$, which effectively transforms the prepared initial state $| {\psi }_{1}\rangle$ $\to {| {\psi }_{1}\rangle }^{{\prime} }=\exp [{\rm{i}}{\hat{{ \mathcal M }}}_{P}(0)]| {\psi }_{1}\rangle$, which then evolves according to the effective Hamiltonian. This can alter the measured evolution drastically depending on the initial phase ϕ₀. A few examples are calculated in [51]. In fact the measurement in figure 9(b) is only successful because we average over two specific phase configurations that are present in the 3D array of plaquettes, which makes the observable independent of ϕ₀ [51]. There are no additional dynamics during one cycle because the operator in equation (75) is proportional to ${\hat{n}}_{{ \mathcal R }}$ and commutes with ${N}_{{ \mathcal R }}$. If the high-frequency limit is not well fulfilled higher-order corrections become relevant and one finds additional oscillations in the real-space dynamics, whose amplitude scales with J_x/Δ. Figure 9(c) displays the expected oscillations for our experimental parameters. The small oscillations cannot be resolved in our measurements because we average over different realizations of ϕ₀.

Resonant tunneling is restored in the static lattice using two running-wave beams as discussed in the previous section (figure 10(a)). The two interfering beams create a local optical potential as defined in equation (64). In the high-frequency limit, the time-dependent problem can be mapped onto a time-independent one as discussed in section 3.4. The corresponding Hamiltonian can be written as

$$\begin{eqnarray}{\hat{{ \mathcal H }}}_{F} & = & -K\displaystyle \sum _{m,n}{{\rm{e}}}^{{\rm{i}}{\phi }_{m,n}}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}\\ & & -J\displaystyle \sum _{m,n}{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}},\end{eqnarray} \tag{ 76 }$$

with ${\phi }_{m,n}=\tfrac{\pi }{2}(n-m)$, $K={J}_{x}{{ \mathcal J }}_{1}({\eta }_{0})$, $J={J}_{y}{{ \mathcal J }}_{0}({\eta }_{0})$, ${\eta }_{0}\,=\sqrt{2}{V}_{0}/({\rm{\hslash }}\omega )$. This Hamiltonian indeed mimics the physics of the Hofstadter model for α = 1/4.

By exploiting an additional pseudo-spin degree of freedom the experimental configuration can be naturally extended to a TR symmetric topological Hamiltonian that underlies the quantum spin-Hall effect [12, 13]. The simplest version of a quantum spin-Hall insulator can be viewed as two independent copies of a QH insulator, where each spin state ($| \uparrow \rangle$ and $| \downarrow \rangle$) exhibits a magnetic field of the same strength but pointing in opposite directions depending on the state.

With cold atoms it can be realized by using two different states with opposite magnetic moments [41]. Here, we use two hyperfine states of ⁸⁷Rb to encode the pseudo-spin, which are denoted as $| \uparrow \rangle \equiv | F=1,{m}_{F}=-1\rangle$ and $| \downarrow \rangle \equiv | F=2,{m}_{F}=-1\rangle$.

The magnetic field gradient automatically generates a spin-dependent linear potential Δ_m = ±Δm, where the plus-sign corresponds to $| \uparrow \rangle$ and the minus-sign to $| \downarrow \rangle$ particles. In section 3.1 we have shown that the Peierls phases ±ϕ_m,n depend on the sign of the energy differences ±Δ between neighboring sites and we find the effective Hamiltonian

$$\begin{eqnarray}{\hat{{ \mathcal H }}}_{\mathrm{spin}} & = & \displaystyle \sum _{m,n}(-K{{\rm{e}}}^{{\rm{i}}{\hat{\phi }}_{m,n}}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}\\ & & -J{\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}},\end{eqnarray} \tag{ 77 }$$

where ${\hat{\phi }}_{m,n}={\phi }_{m,n}{\hat{\sigma }}_{z}$ and ${\hat{\sigma }}_{z}$ is the z-component of the Pauli spin matrix $\vec{\sigma }$. Without loss of generality we have neglected an additional π phase shift between the two spin-components, which has not effect on the properties of the system, since the two components are not coupled. It is easy to see that Hamiltonian (77) is TR symmetric and indeed describes two independent Harper–Hofstadter models, where Φ_↑ = π/2 and Φ_↓ = −π/2. Including additional terms in the Hamiltonian to couple the two individual components may indeed enable the study of non-Abelian gauge fields [143, 144].

In section 5.3 we have seen how superlattice potentials can be used to study the dynamics of single atoms in isolated plaquette potentials and how to extract experimentally the value of the artificially generated magnetic flux [92, 114]. In particular, the motion of the particle was found to be reminiscent of classical cyclotron orbits of charged particles in external magnetic fields. We have studied the micro-motion and its effect on experimental observables. Surprisingly, the results can be interpreted by considering the effective Floquet Hamiltonian in the rotating frame based on equation (72) since the micro-motion and the initial phase can be safely neglected for these measurements. In this section we are going to discuss how the same observables can be used to study the flux distribution across the lattice and the spin-dependence of the field.

As before the experiments were performed in a 3D optical lattice that consists of an array of isolated plaquette potentials (figure 9(a)). The final lattice depths were chosen such that the couplings within the plaquette were J_x/h = 1.2(1) kHz and J_y/h = 0.50(2) kHz and negligible coupling between neighboring plaquettes. Tunneling along x was suppressed by a magnetic field gradient of strength Δ/h ≃ 4.5 kHz.

The measurements started by preparing single atoms in the ground-state of the tilted plaquettes: $| {{\rm{\Psi }}}_{\uparrow }\rangle =(| A\rangle +| D\rangle )/\sqrt{2}$ or $| {{\rm{\Psi }}}_{\downarrow }\rangle =(| B\rangle +| C\rangle )/\sqrt{2}$ depending on the spin state. After suddenly switching on resonant coupling with K/h ≃ 0.28 kHz the mean atom position was measured as a function of time. All plaquette parameters Δ, J and K were calibrated independently for each measurement. A detailed description of these calibrations and further details about the experimental sequence can be found in [41, 51]. Using the site-resolved detection technique mentioned earlier [125, 126] we evaluated the mean atom positions $\langle X\rangle =({N}_{\mathrm{right}}-{N}_{\mathrm{left}})a/(2N)$ and $\langle Y\rangle =({N}_{\mathrm{up}}\,-{N}_{\mathrm{down}})a/(2N)$ within the plaquettes. Here we measured N_left = N_A + N_D, N_right = N_B + N_C, N_up = N_C + N_D and N_down = N_A + N_B, since it is technically more demanding to measure all ${N}_{{ \mathcal R }}$ at once.

We observe that once tunneling is restored the atoms start to tunnel to the opposite bond. Due to the gauge field, however, they get deflected along y, perpendicular to their initial velocity (figure 11). We observe that both spin components are deflected in the same direction, hence, the chirality of their cyclotron-like motion is opposite for the two spins as expected for a magnetic flux, whose direction depends on the spin Φ = ±π/2.

Zoom In Zoom Out Reset image size

These measurements alone cannot distinguish a staggered flux distribution Φ = (−1)^m π/2 from a uniform one because only plaquettes with equal sign of the flux contribute to the signal (section 5). However, by shifting the superlattice potential by one lattice constant (figure 11(b)), for instance, along the x-axis, one can select a different set of plaquette potentials, which have not participated in the previous measurements. We performed equivalent measurements in this configuration and observed exactly the same chiralities for the two spin components as before, which is a signature of the uniformity of the generated flux distribution.

The measured mean atom position $\langle X(t)\rangle$ along x, was fitted according to ${f}_{X}(t)={X}_{0}+{A}_{\langle X\rangle }X(t+\tau )$, where X(t) was computed using the effective time-independent Hamiltonian (72), X₀ and ${A}_{\langle X\rangle }$ were free fit parameters and we have introduced a small experimental time delay τ, which is caused by finite ramp durations. It is not know a priori but it was calibrated and held constant in the fitting algorithms. For the calculation we used the previously determined value of the magnetic flux per plaquette Φ = 0.73(5) × π/2 (section 5.3), whose reduced value stems from the smaller distance between lattice sites along y in the superlattice potential [92].

The spin-Hall effect is closely related to the usual Hall effect, with the difference that in this case a transverse spin current develops if a potential difference is applied to the material. There is no net charge current because spin-up and spin-down electrons are deflected in opposite directions. As a result a transverse spin imbalance builds up, which is analog to the transverse Hall voltage in conventional Hall systems. Experimentally the spin-Hall effect has been studied wit thin-film semiconductor devices [145, 146] and was simulated with harmonically trapped neutral atoms using Raman dressing [147].

Here we briefly review the first experimental study of spin-Hall physics with cold atoms trapped in optical lattices [41]. Based on the previous measurements we study the dynamics of single atoms in plaquette potentials as a function of the spin imbalance ${n}_{\uparrow }-{n}_{\downarrow }$, where n_i = N_i/N, i = $\{\uparrow ,\downarrow \}$ is the fraction of atoms in one spin component. We focus on the analysis of the dynamics in the non-trivial direction, i.e. perpendicular to the initial direction of motion. The experimental sequence was similar to the one described in section 6.3, but this time the atoms were initially prepared in a spin-superposition state ${| \psi \rangle }_{\mathrm{spin}}=\alpha | \uparrow \rangle +\beta | \downarrow \rangle$ before they were prepared in a superposition state on the lower bond of the plaquette potential, with equal weight on A and B sites independent of the spin ${| \psi \rangle }_{\mathrm{spin}}$ of the particle (insets in figure 12). This was achieved by setting the tunneling along y to a small value during state preparation. Increasing it suddenly to a larger value then initiated the dynamics and we measured the mean atom position in the non-trivial x-direction $\langle X\rangle /a$. Details of the experimental sequence are described in [41]. We found almost perfectly mirrored oscillations for the two spin components and extracted the oscillation amplitudes by fitting the numerical evolution as explained above. We found ${A}_{\langle X\rangle }=-0.28(2)$ for ${| \psi \rangle }_{\mathrm{spin}}=| \downarrow \rangle$ (blue) and ${A}_{\langle X\rangle }=0.26(4)$ for ${| \psi \rangle }_{\mathrm{spin}}=| \uparrow \rangle$ (green). In figure 12 we display the measured spin-dependent amplitudes as a function of spin imbalance ${n}_{\uparrow }-{n}_{\downarrow }$. We observe a linear dependence with a sign reversal when the spin is flipped, as expected for the spin-Hall effect.

Zoom In Zoom Out Reset image size

The experimental setup was modified compared to the one described in section 6 in order to increase the experimental flexibility in designing an adiabatic preparation sequence. We have developed an all-optical scheme as illustrated in figure 13(a). The key difference is that one of the running-wave beams is now retro-reflected to create a standing wave, which enables a spatial addressing of every other bond in the lattice. In order to address the remaining ones, the setup needs to be doubled, facilitating flux rectification in a staggered superlattice, which otherwise would result in a staggered flux distribution.

Zoom In Zoom Out Reset image size

The basic operating principle is as follows: along the x-axis tunneling is inhibited by a staggered superlattice potential of the form ${V}_{{\rm{SL}}}(x)={V}_{l}\,{{\rm{\cos }}}^{2}({k}_{s}x/2)+{V}_{s}\,{{\rm{\cos }}}^{2}({k}_{s}x+\pi /2)$, with alternating energy offset ±Δ and hopping J_x between neighboring sites whose strength is determined by V_l and V_s. Tunneling is then resonantly induced using two pairs of laser beams (red and blue arrows in figure 13(a)), with frequency differences $| {\omega }_{i}| \,={\rm{\Delta }}/{\rm{\hslash }}$, i = {r, b}. Both pairs create time-dependent potentials of the form ${V}_{i}(x,y,t)={V}_{0}\cos ({k}_{R}x\,\pm \xi )\cos (-{k}_{R}y+{\omega }_{i}t)$, with k_R = k_s/2 = π/(2a). Independent control over tunneling on bonds with energy offsets of the same sign (red and blue bonds in figure 13(a)) is obtained for ξ = π/4. In this configuration the red pair of beams induces tunneling only on the red bonds and the blue pair on the remaining ones. Setting $\omega := {\omega }_{r}=-{\omega }_{b}$ realizes a uniform flux Φ = ±π/2 per plaquette (figure 13(b)).

In the high-frequency limit ℏω ≫ J_x,y the system can be described by the Hofstadter model for α = 1/4. A detailed derivation of the effective Hamiltonian in the experimental gauge including the effect of the initial phase, which introduces additional corrections, can be found in [51]. The lattice parameters were set to have equal effective coupling strengths along both axes.

Following the theoretical description presented in section. 2 we find that the magnetic unit cell is four times larger than the lattice unit cell A_MU = 2a × 2a. It contains four non-equivalent sites and the lowest tight-binding band splits in into four subbands, where the two middle ones touch at four Dirac cones (figure 13(c)). Hence we speak of three well-separated Hofstadter bands. The lowest band exhibits topological properties equivalent to the lowest Landau level, with Chern number ν₁ = +1. Moreover it is characterized by a rather large flatness ratio of E_gap/ E_bw ≃ 7 making this system an excellent candidate to study fractional Chern insulators [150].

Having generated the desired bandstructure leaves one with the non-trivial task of preparing the atoms in the lowest Hofstadter band. Let us consider a sequence where we start with a BEC in the lowest band of the static trivial square lattice. Turning on the artificial flux immediately leads to a splitting of the band and the atoms might distribute between the higher Hofstadter bands in an uncontrolled way. To overcome this challenge we designed an alternative sequence, where we use additional staggered potentials along both directions to increase the unit cell of the static lattice. It is of the form

$$\begin{eqnarray}&&{\hat{H}}_{d}=\displaystyle \frac{\delta }{2}\displaystyle \sum _{m,n}[{(-1)}^{m}+{(-1)}^{n}]{\hat{n}}_{m,n}.\end{eqnarray} \tag{ 78 }$$

Turning on the modulation for large enough δ does not induce tunneling and the system remains trivial. Tuning the parameter δ allows us to connect these trivial bands to the Hofstadter bands in a controlled manner, since the number of subbands remains constant. The Hofstadter model with α = 1/4 is reached for δ = 0. Here we benefit from the increased experimental control facilitated by the all-optical setup [73].

The sequence started by preparing a BEC in the static lattice, which has a unit cell equivalent to the magnetic unit cell of the final configuration as depicted in figure 14(a). This is achieved with the additional staggered potentials, which introduce a detuning such that the laser-assisted tunneling is off-resonant. In the x-direction the detuning increases the already existing energy offset away from the resonance condition to Δ + δ, while along y it inhibits normal tunneling for δ ≫ J_y. In this regime the bandstructure is topologically trivial but the unit cell has the same size as the final magnetic unit cell. By slowly decreasing δ to zero we connect the spectrum to the desired Hofstadter spectrum (figure 14(b)). As a consequence, the number of bands is preserved during the loading sequence. It starts with large detuning (δ > 2J) with trivial bands, whose Chern numbers vanish. By decreasing δ one crosses a topological phase transition at δ = 2J, where the gaps in the spectrum close. Beyond this point the system is topologically non-trivial with non-zero Chern numbers and at the end of the sequence (δ = 0) the Harper–Hofstadter Hamiltonian for Φ = π/2 is realized, whose lowest band is characterized by ν₁ = +1. When crossing the topological phase transition the gaps in the spectrum close at one point in the FBZ (figure 14(b)). However, if we prepare a BEC in the lowest band it should be only minimally affected by this gap-closing point.

Zoom In Zoom Out Reset image size

Moreover this protocol can be reversed, which allows us to map the population of the different Hofstadter bands onto four topologically trivial subbands in the static lattice and perform band-mapping measurements to determine the populations experimentally. Assuming that heating effects and scattering processes are negligible during this sequence we further obtain access to the momentum distribution of the Hofstadter bands. Using this technique we found that we achieve an initial atom population in the lowest Hofstadter band of about 60%.

In section 2.4 we have seen how in electronic systems the quantized Hall conductance relates to the Chern numbers of the occupied energy bands. A complementary way to describe the properties of QH insulators can be established by looking at boundary modes. A quantum Hall system is insulating in the bulk but conducting at the edges and the conductivity is given by the number of metallic edge states, which is in turn related to the sum of Chern numbers below the Fermi energy, known as bulk-edge correspondence [8, 9].

In many cold-atom experiments there are no steep potential boundaries in contrast to solid-state measurements, hence, it offers a natural setting to look for observables that study the bulk topological properties. Indeed it was shown that a measurement of the anomalous Hall velocity [72], that emerges transverse to an applied potential gradient and is proportional to the Berry curvature, could reveal the Chern number of the corresponding energy band [148, 149].

Using a semi-classical approach one can calculate the velocity of a wavepacket in a state $| {u}_{\mu }({\bf{k}})\rangle$ of the μth energy band centered around momentum k in the presence of an external force ${\bf{F}}=F{\hat{{\bf{e}}}}_{y}$, according to [72]

$$\begin{eqnarray}{v}_{\mu }^{y}({\bf{k}}) & = & \displaystyle \frac{1}{{\rm{\hslash }}}{\partial }_{{k}_{y}}{E}_{\mu },\\ {v}_{\mu }^{x}({\bf{k}}) & = & \displaystyle \frac{1}{{\rm{\hslash }}}{\partial }_{{k}_{x}}{E}_{\mu }\mathop{\underbrace{-\displaystyle \frac{F}{{\rm{\hslash }}}{{\rm{\Omega }}}_{\mu }({\bf{k}})}}\limits_{\mathrm{anomalous}\,\mathrm{velocity}}.\end{eqnarray} \tag{ 79 }$$

As captured by the band velocity ${{\bf{v}}}_{\mu }^{\mathrm{band}}={\partial }_{{\bf{k}}}{E}_{\mu }/{\rm{\hslash }}$, the wavepacket performs Bloch oscillations [123] but in addition it exhibits a transverse velocity component proportional to the Berry curvature Ω_μ(k) defined in equation (20), which is the origin of the Hall velocity (figure 15(a)). Recently, the anomalous Hall velocity was observed with ultracold atoms in a shaken honeycomb lattice [40] and in an electronic graphene superlattice [151].

Zoom In Zoom Out Reset image size

This anomalous component can be employed to experimentally determine the Chern number of an energy band of interest. To isolate the non-trivial dynamics from the band velocity components, one may consider a uniformly populated energy band, where the particle density ρ_μ(k) ≡ ρ_μ at a certain momentum k and hence per state $| {u}_{\mu }({\bf{k}})\rangle$ is independent of k. Due to the symmetry of the dispersion relation the contributions from the band velocity identically vanish. For an infinite system the mean band velocity for the μth band is computed as

$$\begin{eqnarray}&&{\langle {v}_{\mu }^{\mathrm{band}}\rangle }_{i}=\displaystyle \frac{{\rho }_{\mu }}{{\rm{\hslash }}}{\int }_{\mathrm{FBZ}}{\partial }_{{k}_{i}}{E}_{\mu }{{\rm{d}}}^{2}k=0,\end{eqnarray} \tag{ 80 }$$

where the integration is performed over all momenta in the FBZ. For fermions this is achieved, if the Fermi energy lies within a spectral gap [149]. Here we consider bosonic atoms, which naturally condense in the lowest energy eigenstate. However, if we consider a finite-temperature system, whose temperature is large compared to the bandwidth of the lowest band but small compared to the energy gap, one can realize an incoherent distribution of atoms which populate the band approximately uniformly. In both cases the mean anomalous velocity of the μth band of an infinite system is calculated as

$$\begin{eqnarray}&&\langle {v}_{\mu }^{x}({\bf{k}})\rangle =-\displaystyle \frac{{\rho }_{\mu }F}{{\rm{\hslash }}}\mathop{\underbrace{{\int }_{\mathrm{FBZ}}{{\rm{\Omega }}}_{\mu }({\bf{k}}){{\rm{d}}}^{2}k}}\limits_{\to 2\pi {\nu }_{\mu }},\end{eqnarray} \tag{ 81 }$$

where ν_μ is the Chern number of the μth band. For well-separated non-degenerate energy bands the particle density in each band is found according to

$$\begin{eqnarray}&&{\rho }_{\mu }=\displaystyle \frac{{N}_{\mu }}{N}\displaystyle \frac{1}{{A}_{\mathrm{FBZ}}}={\eta }_{\mu }\displaystyle \frac{{a}^{2}}{{\pi }^{2}},\qquad {\eta }_{\mu }:= \displaystyle \frac{{N}_{\mu }}{N}\end{eqnarray} \tag{ 82 }$$

where N_μ is the particle number in each band, $N={\sum }_{\mu }{N}_{\mu }$ is the total atom number and A_FBZ = π²/a² is the area of the FBZ. Under the assumption that no additional dynamics occurs that would change the band populations over time (η_μ ≡ η_μ⁰), one can compute the contribution of the μth band to the center-of-mass (COM) displacement of the atoms transverse to the applied force as follows

$$\begin{eqnarray}&&{x}_{\mu }(t)={\eta }_{\mu }^{0}\displaystyle \frac{4{a}^{2}F}{h}{\nu }_{\mu }\,t=-4a{\eta }_{\mu }^{0}{\nu }_{\mu }\displaystyle \frac{t}{{\tau }_{{\rm{B}}}},\end{eqnarray} \tag{ 83 }$$

where τ_B = h/(Fa) is the Bloch oscillations period. For a filled lowest band η₁ = 1 (η_μ>1 = 0) we obtain x₁(t) = −4aν₁t/τ_B, which corresponds to a deflection of 4aν₁ per Bloch oscillation period and the Chern number ν₁ can be simply extracted from the slope of the linear displacement. It is at least on the order of four lattice constants and therefore can be detected with conventional imaging systems.

The initial state preparation is likely not perfect, hence there will be higher-band occupancies, which can be easily incorporated in the equation of motion. The Hofstadter model for Φ = π/2 consist of three energy bands. Note that the middle one consists of two subbands (figure 13(c)) and contains twice as many states. Thus, the particle density is determined by ρ₂ = η₂/(2A_FBZ); for each momentum there are two available quantum states. The combined COM-displacement $x(t)={\sum }_{\mu }{x}_{\mu }(t)$ is given by

$$\begin{eqnarray}&&x(t)=-4a\displaystyle \frac{t}{{\tau }_{{\rm{B}}}}\left({\eta }_{1}^{0}\,{\nu }_{1}+\displaystyle \frac{{\eta }_{2}^{0}}{2}\,{\nu }_{2}+{\eta }_{3}^{0}\,{\nu }_{3}\right).\end{eqnarray} \tag{ 84 }$$

Using the particle-hole symmetry of the Hofstadter model and ${\sum }_{\mu }{\nu }_{\mu }=0$ as discussed in section 2.4 this relation can be further simplified. Since the Chern number distribution is symmetric around zero energy one finds ν₁ = ν₃ = −ν₁ − ν₂, which can be used to rewrite equation (84) in terms of the lowest band Chern number:

$$\begin{eqnarray}&&x(t)=-4a\displaystyle \frac{t}{{\tau }_{{\rm{B}}}}\,{\gamma }_{0}\,{\nu }_{1},\qquad {\gamma }_{0}={\eta }_{1}^{0}-{\eta }_{2}^{0}+{\eta }_{3}^{0},\end{eqnarray} \tag{ 85 }$$

where we have introduced the initial band-filling factor γ₀, which was assumed to remain constant during the dynamics. The largest deflection occurs for γ₀ = 1, which corresponds to a filled lowest band, while it continuously decreases if higher-band populations are taken into account. For an infinite temperature state, γ₀ = 0 and no transverse motion is expected.

Finally the question on the importance of the micro-motion for the COM-measurement arises. In the high-frequency limit the dominant contribution to the micro-motion operator is defined in equation (54). We find that it is diagonal in the density operator ${\hat{n}}_{m,n}$ and we expect no influence of the micro-motion on our measurement. This lowest-order approximation is valid in the limit ℏω ≫ J_x, which is well fulfilled for the measurements discussed in the following section, where ℏω / J_x ≃ 13.5. For parameter regimes, where this limit is not well fulfilled additional contributions are expected that scale with an amplitude J_x/(ℏω) [102] as already discussed in section 5.3. For comparison, the residual oscillations shown in figure 9(c) for the cyclotron dynamics would be reduced by a factor of ∼3 for the experimental parameters used for the Chern number measurement.

The experimental sequence started by preparing the atoms in the Hofstadter bands using the loading scheme presented in section 7.1. Subsequently an external force F was applied. It was created using an additional laser beam, which was aligned along y, so that the atoms are positioned at the maximum slope of the Gaussian beam profile to induce the potential gradient, while along x the cloud was centered on the profile. We then tracked the COM-motion in situ using absorption imaging. In order to be less susceptible to systematic errors, that might be induced by a long-term drift of the absolute position of the atomic cloud, we measured the positions for opposite directions of the flux Φ = ±π/2 and evaluated the differential shift x(t, π/2) − x(t, −π/2) = 2x(t). In the following one dataset is defined as a set of ten images for each direction of the magnetic flux, which were taken in an alternating manner, averaged and then subtracted (figure 15(b) right column).

For short times there is an almost linear evolution of the differential COM-displacement in agreement with equation (85) and numerical studies [148, 149]. For longer times 2τ_B > 50 ms, however, we observe a rapid saturation of the transverse motion (figure 15(b)). Numerical simulations did indeed reveal small deviations from a perfectly linear evolution, which is due to Landau–Zener transitions. The observed saturation, however, happens on a much faster timescale. This suggests that additional effects associated with heating, most likely caused by the periodic modulation and its interplay with interactions between the particles, lead to repopulation between the bands [152, 153]. Moreover, we see a symmetric evolution with respect to the direction of the force ${\bf{F}}=\pm F{\hat{{\bf{e}}}}_{y}$ (black and gray data points) and no significant motion for a lattice without flux Φ = 0 (light blue data points) or staggered flux (dark blue data points), whose bands are topologically trivial with ν = 0 (section 5).

Making use of independent band population measurements we were able to confirm that atoms are transferred to higher Hofstadter bands during the dynamics, which induces a time-dependence for the filling factor γ(t) (figure 15(c)). This could be incorporated into our theoretical model by extending equation (85) and we obtain

$$\begin{eqnarray}&&x(t)=-\displaystyle \frac{4a}{{\tau }_{{\rm{B}}}}{\nu }_{1}{\int }_{0}^{t}\gamma ({t}^{{\prime} }){\rm{d}}{t}^{{\prime} },\end{eqnarray} \tag{ 86 }$$

with $\gamma (t)=[{\eta }_{1}(t)-{\eta }_{2}(t)+{\eta }_{3}(t)]$. The timescale associated with the decay of the filling factor $\gamma \to 0$ is consistent with the saturation of the observed COM-motion. Incorporating the additional population data in our model and fitting equation (86) to the COM-data enables a determination of ν₁, where the Chern number is the only free fit parameter. With this we found ν_exp = 0.99(5) (averaged over four independent measurements), which is in agreement with the theoretically expected value. Indeed this demonstrates that we have a good quantitative understanding of the dynamics and properties of the Floquet-engineered bandstructure. However, in order to improve temperatures by suppressing the heating mechanisms more detailed studies on modulation-induced heating processes are required [152, 153].