Exploring 4D Quantum Hall Physics with a 2D Topological Charge Pump

TL;DR

This paper demonstrates the first experimental observation of a 4D quantum Hall effect using a 2D topological charge pump with ultracold atoms, revealing a non-linear, quantized response linked to a 4D topological invariant.

Contribution

It introduces a method to realize and measure a 4D quantum Hall effect in a 2D system, enabling exploration of higher-dimensional topological phases.

Findings

First measurement of a bulk 4D topological response

Implementation of a dynamical 4D quantum Hall system

Full characterization of non-linear topological response

Abstract

The discovery of topological states of matter has profoundly augmented our understanding of phase transitions in physical systems. Instead of local order parameters, topological phases are described by global topological invariants and are therefore robust against perturbations. A prominent example thereof is the two-dimensional integer quantum Hall effect. It is characterized by the first Chern number which manifests in the quantized Hall response induced by an external electric field. Generalizing the quantum Hall effect to four-dimensional systems leads to the appearance of a novel non-linear Hall response that is quantized as well, but described by a 4D topological invariant - the second Chern number. Here, we report on the first observation of a bulk response with intrinsic 4D topology and the measurement of the associated second Chern number. By implementing a 2D topological charge pump with ultracold bosonic atoms in an angled optical superlattice, we realize a dynamical version of the 4D integer quantum Hall effect. Using a small atom cloud as a local probe, we fully characterize the non-linear response of the system by in-situ imaging and site-resolved band mapping. Our findings pave the way to experimentally probe higher-dimensional quantum Hall systems, where new topological phases with exotic excitations are predicted.

Figures

a, b is chosen such that the 4D Levi-Civita symbol is $\varepsilon_{wxyz}=+1$. The velocity of the wave packet $\mathbf{v}=\dot{\mathbf{r}}$ has two contributions: the group velocity arising from the dispersion of the band and the anomalous velocity due to the non-zero Berry curvature $\Omega^{\nu\mu}(\mathbf{k})=i\left(\left<\partial_{k_{\nu}}u|\partial_{k_{\mu}}u\right>-\left<\partial_{k_{\mu}}u|\partial_{k_{\nu}}u\right>\right)$. For a filled or homogeneously populated band, the group velocity term vanishes and with $\mathbf{E}=E_{z}\mathbf{e}_{z}$ and $\mathbf{B}=0$, the linear Hall respon

c leads to a systematic shift of $\nu_{2}^{\mathrm{exp}}$ by +0.05.

Equations

$${\nu }_1^x{d}_{\mathrm{l},\mathrm{x}}{\mathbf{e}}_x+{\nu }_2\theta d_{\mathrm{l},\mathrm{x}}{\mathbf{e}}_y$$

$${\hat{H}}_{2D}({\phi }_x,{\phi }_y)=$$

$$-\sum _{m_x,m_y}\left[J_x({\phi }_x)+(-1{)}^{m_x}\delta J_x({\phi }_x)/2\right]{\hat{a}}_{m_x+1,m_y}^{†}{\hat{a}}_{m_x,m_y}+h.c.$$

$$-\sum _{m_x,m_y}\left[J_y({\phi }_y)+(-1{)}^{m_y}\delta J_y({\phi }_y)/2\right]{\hat{a}}_{m_x,m_y+1}^{†}{\hat{a}}_{m_x,m_y}+h.c.$$

$$+\sum _{m_x,m_y}\frac{1}{2}\left[(-1{)}^{m_x}{\Delta }_x({\phi }_x)+(-1{)}^{m_y}{\Delta }_y({\phi }_y)\right]{\hat{a}}_{m_x,m_y}^{†}{\hat{a}}_{m_x,m_y}$$

$${\mathcal{I}}_y({\phi }_x)=n_{\mathrm{gs}}{\mathcal{I}}_y^{\mathrm{gs}}({\phi }_y)+n_{\mathrm{exc}}{\mathcal{I}}_y^{\mathrm{exc}}({\phi }_y)+n_2{\mathcal{I}}_y^{2,\mathrm{gs}}({\phi }_y)$$

$$\frac{\partial {\mathcal{I}}_y^{\mathrm{gs}}({\phi }_y^{\exp })}{\partial {\phi }_x}=\frac{\partial {\mathcal{I}}_y^{\mathrm{gs}}({\phi }_y^{\exp })}{\partial {\phi }_y^{\exp }}\frac{\partial {\phi }_y^{\exp }}{\partial {\phi }_x}=\alpha \frac{\partial {\mathcal{I}}_y^{\mathrm{gs}}({\phi }_y)}{\partial {\phi }_x}$$

$$\delta\mathcal{I}_{y}^{\mathrm{gs}}=\alpha\left[\mathcal{I}_{y}^{\mathrm{gs}}(\varphi_{y})\Big{|}_{\varphi_{x}=2\pi}-\mathcal{I}_{y}^{\mathrm{gs}}(\varphi_{y})\Big{|}_{\varphi_{x}=0}\right]$$

$${\dot{r}}^{\mu }=\frac{1}{\hslash }\frac{\partial \mathcal{E}(\mathbf{k})}{\partial k_{\mu }}+{\dot{k}}_{\nu }{\Omega }^{\nu \mu }(\mathbf{k})$$

$$\hslash {\dot{k}}_{\mu }=q{E}_{\mu }+q{\dot{r}}^{\nu }{B}_{\mu \nu }$$

$${\mathbf{v}}_{\mathrm{COM}}^{(0)}=\frac{q}{h}{A}_{\mathrm{M}}^{xz}{E}_z{\nu }_1^{zx}{\mathbf{e}}_x$$

$${\mathbf{v}}_{\mathrm{COM}}=\frac{q}{h}{A}_{\mathrm{M}}^{xz}{E}_z{\nu }_1^{zx}{\mathbf{e}}_x-{\left(\frac{q}{h}\right)}^2{A}_M{E}_z{B}_{xw}{\nu }_2{\mathbf{e}}_y$$

$${\hat{H}}_{2D}({\phi }_x,{\phi }_y)=$$

$$-\sum _{m_x,m_y}\left[J_x({\phi }_x)+\delta J_x^{m_x}({\phi }_x)\right]{\hat{a}}_{m_x+1,m_y}^{†}{\hat{a}}_{m_x,m_y}+h.c.$$

$$-\sum _{m_x,m_y}\left[J_y({\phi }_y)+\delta J_y^{m_y}({\phi }_y)\right]{\hat{a}}_{m_x,m_y+1}^{†}{\hat{a}}_{m_x,m_y}+h.c.$$

$$+\sum _{m_x,m_y}\left[{\Delta }_x^{m_x}({\phi }_x)+{\Delta }_y^{m_y}({\phi }_y)\right]{\hat{a}}_{m_x,m_y}^{†}{\hat{a}}_{m_x,m_y}$$

$$\delta J_x^{m_x}({\phi }_x)=-\frac{\delta J_x^{(0)}}{2}\sin \left({\tilde{\Phi }}_{xz}{m}_x-{\phi }_x\right)$$

$$\delta J_y^{m_y}({\phi }_y)=-\frac{\delta J_y^{(0)}}{2}\sin \left({\tilde{\Phi }}_{yw}{m}_y-{\phi }_y\right)$$

$${\Delta }_x^{m_x}({\phi }_x)=\frac{{\Delta }_x^{(0)}}{2}\sin \left({\tilde{\Phi }}_{xz}(m_x-1/2)-{\phi }_x\right)$$

$${\Delta }_y^{m_y}({\phi }_y)=\frac{{\Delta }_y^{(0)}}{2}\sin \left({\tilde{\Phi }}_{yw}(m_y-1/2)-{\phi }_y\right)$$

$${\hat{H}}_{4D}=\frac{1}{4{\pi }^2}{\int }_0^{2\pi }{\hat{H}}_{2D}({\phi }_x,{\phi }_y)\mathrm{d}{\phi }_x\mathrm{d}{\phi }_y^{(0)}$$

$${\hat{a}}_{m_x,m_y}^{†}=\sum _{m_z,m_w}{e}^{i({\phi }_x{m}_z+{\phi }_y^{(0)}{m}_w)}{\hat{a}}_{\mathbf{m}}^{†}$$

$${\hat{a}}_{m_x,m_y}=\sum _{m_z,m_w}{e}^{-i({\phi }_x{m}_z+{\phi }_y^{(0)}{m}_w)}{\hat{a}}_{\mathbf{m}}$$

$${\hat{H}}_{4D}={\hat{H}}_{xz}+{\hat{H}}_{yw}+{\hat{H}}_{\delta J}$$

$${\hat{H}}_{xz}=$$

$$-\sum _{\mathbf{m}}\frac{{\Delta }_x^{(0)}}{4}{e}^{i\left[{\tilde{\Phi }}_{xz}(m_x-1/2)+\pi /2\right]}{\hat{a}}_{\mathbf{m}+{\mathbf{e}}_z}^{†}{\hat{a}}_{\mathbf{m}}+h.c.$$

$${\hat{H}}_{yw}$$

$$-$$

$$B_{xw}=-\frac{{\Phi }_0}{d_{\mathrm{s},\mathrm{w}}{d}_{\mathrm{l},\mathrm{y}}}\theta$$

$${\Omega }^x(k_x,{\phi }_x,k_y,{\phi }_y)=i\left(⟨{\partial }_{{\phi }_x}u|{\partial }_{k_x}u⟩-⟨{\partial }_{k_x}u|{\partial }_{{\phi }_x}u⟩\right)$$

$${\nu }_1^x=\frac{1}{2\pi }\oint {\Omega }^x\mathrm{d}{k}_x\mathrm{d}{\phi }_x$$

$$v_y(k_x,{\phi }_x,k_y,{\phi }_y)={\Omega }^y{\partial }_t{\phi }_y=\frac{2\pi }{d_{\mathrm{l},\mathrm{y}}}\theta {\Omega }^x{\Omega }^y{\partial }_t{\phi }_x$$

$$\delta y_{\mathrm{COM}}=\underset{\overline{\Omega }({\phi }_y)}{\underbrace{\frac{1}{2\pi }\oint {\Omega }^x{\Omega }^y\mathrm{d}{k}_x\mathrm{d}{k}_y\mathrm{d}{\phi }_x}}\theta d_{\mathrm{l},\mathrm{x}}$$

$$\delta y_{\mathrm{COM}}=\frac{1}{2\pi }\oint \overline{\Omega }({\phi }_y)\theta d_{\mathrm{l},\mathrm{x}}\mathrm{d}{\phi }_y={\nu }_2\theta d_{\mathrm{l},\mathrm{x}}$$

$${\nu }_2=1/(4{\pi }^2){\oint }_{\mathrm{BZ}}{\Omega }^x{\Omega }^y\mathrm{d}{k}_x\mathrm{d}{k}_y\mathrm{d}{\phi }_x\mathrm{d}{\phi }_y$$

$$\left(\begin{array}[]{c}r_{1}\\ r_{2}\\ r_{3}\\ r_{4}\end{array}\right)=\frac{1}{4}\left(\begin{array}[]{cccc}1&1&-1&-1\\ 1&1&1&1\\ 1&-1&-1&1\\ 1&-1&1&-1\end{array}\right)\cdot\left(\begin{array}[]{c}\delta J_{x}/\delta J_{x}^{(0)}\\ \Delta_{x}/\Delta_{x}^{(0)}\\ \delta J_{y}/\delta J_{y}^{(0)}\\ \Delta_{y}/\Delta_{y}^{(0)}\end{array}\right)$$

Full text

Exploring 4D Quantum Hall Physics with a 2D Topological Charge Pump

Michael Lohse1,2, Christian Schweizer1,2, Hannah M. Price3, Oded Zilberberg4 & Immanuel Bloch1,2

1 Fakultät für Physik, Ludwig-Maximilians-Universität, Schellingstraße 4, 80799 München, Germany

2 Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Straße 1, 85748 Garching, Germany

3 INO-CNR BEC Center & Dipartimento di Fisica, Università di Trento, Via Sommarive 14, 38123 Povo, Italy

4 Institut für Theoretische Physik, ETH Zürich, Wolfgang-Pauli-Straße 27, 8093 Zürich, Switzerland

The discovery of topological states of matter has profoundly augmented our understanding of phase transitions in physical systems. Instead of local order parameters, topological phases are described by global topological invariants and are therefore robust against perturbations. A prominent example thereof is the two-dimensional integer quantum Hall effect Klitzing et al. (1980). It is characterized by the first Chern number which manifests in the quantized Hall response induced by an external electric field Thouless et al. (1982). Generalizing the quantum Hall effect to four-dimensional systems leads to the appearance of a novel non-linear Hall response that is quantized as well, but described by a 4D topological invariant – the second Chern number Yang (1978); Zhang and Hu (2001). Here, we report on the first observation of a bulk response with intrinsic 4D topology and the measurement of the associated second Chern number. By implementing a 2D topological charge pump with ultracold bosonic atoms in an angled optical superlattice, we realize a dynamical version of the 4D integer quantum Hall effect Thouless (1983); Kraus et al. (2013). Using a small atom cloud as a local probe, we fully characterize the non-linear response of the system by in-situ imaging and site-resolved band mapping. Our findings pave the way to experimentally probe higher-dimensional quantum Hall systems, where new topological phases with exotic excitations are predicted Zhang and Hu (2001).

Topology, originally a branch of mathematics, has become an important concept in different fields of physics, ranging from particle physics Yang and Mills (1954); Yang (1978) to solid state physics Qi and Zhang (2011) and quantum computation Nayak et al. (2008). In this context, a hallmark achievement was the discovery of the 2D integer quantum Hall (QH) effect Klitzing et al. (1980). It demonstrated that the Hall conductance, in response to an electric field $\mathbf{E}$, is quantized for charged particles moving in 2D in the presence of a perpendicular magnetic field. In a cylindrical geometry, following Laughlin’s gedankenexperiment, $\mathbf{E}$ can be generated by varying the magnetic flux $\phi_{x}(t)$ along the axis of the cylinder Laughlin (1981) (Fig. 1a). Due to the interplay of the perpendicular magnetic field and the induced electric field $E_{z}$, the particles start moving in the $x$-direction. This gives rise to the quantized Hall response through an integer number of particles that is transported from one edge to the other per flux quantum threaded through the cylinder. This quantized response is characterized by an integer topological invariant, the first Chern number Thouless et al. (1982).

Dimensionality plays a crucial role for topological phases and many intriguing states have recently been discovered in 3D systems, e.g. Weyl semimetals Lu et al. (2015); Xu et al. (2015) and 3D topological insulators Hsieh et al. (2008). Ascending further in dimensions, a generalization of the QH effect to 4D was proposed Zhang and Hu (2001). It has received much attention in theoretical studies Qi et al. (2008); Edge et al. (2012); Li et al. (2013), in particular since it might exhibit novel strongly correlated QH phases Zhang and Hu (2001). This interest was renewed recently due to the unprecedented control and flexibility offered by engineered systems like ultracold atoms and photonics. Indeed, such systems have already been used to study various topological effects Goldman et al. (2016); Lu et al. (2014), including in synthetic dimensions Mancini et al. (2015); Stuhl et al. (2015), and offer a direct route for realizing 4D physics Price et al. (2015); Ozawa et al. (2016).

In the simplest case, a 4D QH system can be composed of two 2D QH systems in orthogonal subspaces (Fig. 1a, b). In addition to the quantized linear response underlying the 2D QH effect, it exhibits a quantized non-linear 4D Hall response Kraus et al. (2013). The latter arises when – simultaneously with the perturbing electric field $\mathbf{E}$ – a magnetic perturbation $\mathbf{B}$ is added, coupling the motion in the two 2D QH systems. Due to the 4D symmetry, there are multiple possibilities for the orientation of $\mathbf{E}$ and $\mathbf{B}$. The non-linear response, however, is always characterized by the same 4D topological invariant, the second Chern number. For simplicity, we therefore focus on the geometry depicted in Fig. 1a, b. In this case, the non-linear response can be understood semiclassically as originating from a Lorentz force created by $\mathbf{B}$ Price et al. (2016). The direction of this response is in an orthogonal subspace compared to both perturbing fields and it can thus only occur in four or more dimensions. Therefore, such a response has never been observed in any physical system. Note that a measurement of the second Chern number has recently been reported in an artificially generated parameter space of a four-level system Sugawa et al. (2016).

Topological charge pumps exhibit topological transport properties similar to higher-dimensional QH systems and thus provide a way to probe 4D QH physics in lower-dimensional dynamical systems. The first example for such a topological charge pump is the 1D Thouless pump Thouless (1983), where a quantized particle transport is generated by varying the system’s parameters in an adiabatic and periodic way. This modulation can be parametrized by a pump parameter and at each point in the cycle, the 1D system constitutes a Fourier component of a 2D QH system Qi et al. (2008); Kraus and Zilberberg (2012). The induced motion is thus equivalent to the linear Hall response and characterized by the same 2D topological invariant, the first Chern number. Indeed, the QH effect on a cylinder can be mapped to a 1D charge pump with the threaded flux $\phi_{x}$ acting as the pump parameter Laughlin (1981) (Fig. 1a). Building on pioneering experiments in condensed matter systems Pothier et al. (1992); Switkes et al. (1999), topological charge pumps have recently been realized in photonic waveguides Kraus et al. (2012); Verbin et al. (2015) and with ultracold atoms Lohse et al. (2016); Nakajima et al. (2016); Lu et al. (2016); Schweizer et al. (2016).

A dynamical version of the 4D QH effect can be realized by extending the concept of topological charge pumping to 2D Kraus et al. (2013). Using dimensional reduction Qi et al. (2008); Kraus and Zilberberg (2012), the Fourier components of a 4D QH system can be mapped onto a 2D system. For the geometry in Fig. 1a, b, the corresponding 2D model is a square superlattice (Fig. 1c and Supplementary Information). It consists of two 1D superlattices along $x$ and $y$, each formed by superimposing two lattices $V_{\mathrm{s,\mu}}\sin^{2}\left(\pi\mu/d_{\mathrm{s,\mu}}\right)+V_{\mathrm{l,\mu}}\sin^{2}\left(\pi\mu/d_{\mathrm{l,\mu}}-\varphi_{\mu}/2\right)$, $\mu\in\{x,y\}$. Here, $d_{\mathrm{s,\mu}}$ and $d_{\mathrm{l,\mu}}>d_{\mathrm{s,\mu}}$ denote the lattice periods and $V_{\mathrm{s,\mu}}$ ($V_{\mathrm{l,\mu}}$) the depth of the short (long) lattice potential, respectively. The position of the long lattices relative to the short ones is determined by the corresponding superlattice phases $\varphi_{\mu}$.

The phase $\varphi_{x}$ is chosen as the pump parameter, i.e. pumping is performed by moving the long lattice along $x$. This is equivalent to threading $\phi_{x}$ in the 4D model and therefore leads to a quantized motion along $x$ – the linear response (Fig. 1c, d). The magnetic perturbation $B_{xw}$ corresponds to a transverse superlattice phase $\varphi_{y}$ that depends linearly on $x$. This can be realized by tilting the long $y$-lattice relative to the short one by an angle $\theta\ll 1$ in the $xy$-plane (Fig. 1c). Then $\varphi_{y}(x)=\varphi_{y}^{(0)}+2\pi\theta\,x/d_{\mathrm{l,y}}$ to first order in $\theta$, thereby coupling the two orthogonal axes. When $\varphi_{x}$ is varied, the motion along $x$ thus changes $\varphi_{y}$ and – analogous to the Lorentz force in 4D – induces a quantized non-linear response along $y$ Kraus et al. (2013) (Fig. 1d).

For a uniformly populated band in an infinite system, the change in the centre-of-mass (COM) position during one cycle $\varphi_{x}=0\rightarrow 2\pi$ is

[TABLE]

with $\mathbf{e}_{x}$ ($\mathbf{e}_{y}$) the unit vector along $x$ ($y$) (see Supplementary Information). The first term describes the quantized linear response in the $x$-direction, proportional to the pump’s first Chern number $\nu_{1}^{x}$. Here, $\nu_{1}^{x}$ is defined using a generalized 2D Brillouin zone spanned by the quasimomentum $k_{x}$ and $\varphi_{x}$. It is obtained by integrating the Berry curvature $\Omega^{x}(k_{x},\varphi_{x})=i\left(\left<\partial_{\varphi_{x}}u|\partial_{k_{x}}u\right>-\left<\partial_{k_{x}}u|\partial_{\varphi_{x}}u\right>\right)$ over the entire Brillouin zone, where $\left|u(k_{x},\varphi_{x})\right>$ denotes the eigenstate of a given non-degenerate band at $k_{x}$ and $\varphi_{x}$. As $\nu_{1}^{x}$ can only take integer values, the motion is quantized Lohse et al. (2016). The second term is the non-linear response in the $y$-direction. It is quantified by a 4D integer topological invariant, the pump’s second Chern number $\nu_{2}=1/(4\pi^{2})\oint_{\mathrm{BZ}}\Omega^{x}\Omega^{y}\mathrm{d}k_{x}\mathrm{d}k_{y}\mathrm{d}\varphi_{x}\mathrm{d}\varphi_{y}$, where $\mathrm{BZ}$ indicates the generalized 4D Brillouin zone (Fig. 1e). Therefore, the non-linear response is quantized as well and has intrinsic 4D symmetries resulting from the higher-dimensional non-commutative geometry.

We implement a 2D topological charge pump with bosonic 87Rb atoms forming a Mott insulator in isolated planes of a 3D optical lattice with superlattices along $x$ and $y$ with $d_{\mathrm{s}}\equiv d_{\mathrm{s,x}}=d_{\mathrm{s,y}}$ and $d_{\mathrm{l}}\equiv d_{\mathrm{l,x}}=d_{\mathrm{l,y}}=2d_{\mathrm{s}}$ (see Methods). This creates double well potentials along both $x$ and $y$ (Fig. 1c). In the tight-binding limit, this realizes a 2D Rice-Mele model Rice and Mele (1982) in each plane, where the on-site energies of the lattice sites are modulated alternatingly by $(-1)^{m_{x}}\Delta_{x}(\varphi_{x})/2+(-1)^{m_{y}}\Delta_{y}(\varphi_{y})/2$. Here, $m_{x}$ ($m_{y}$) denotes the $m_{x}$-th ($m_{y}$-th) lattice site along $x$ ($y$). The long lattices also lead to dimerized tunnel couplings between neighbouring sites, $J_{x}(\varphi_{x})+(-1)^{m_{x}}\delta J_{x}(\varphi_{x})/2$ and $J_{y}(\varphi_{y})+(-1)^{m_{y}}\delta J_{y}(\varphi_{y})/2$. The corresponding unit cell is a four-site plaquette and the lowest band in the short lattices splits into four subbands.

In the experiment, we study the non-linear bulk response of the lowest subband, for which $\nu_{2}=+1$ is expected for $d_{\mathrm{l}}=2d_{\mathrm{s}}$. As the initial state, a Mott insulator with quarter filling is prepared at $\varphi_{x}=0$, where one atom is localized in the ground state of each unit cell, leading to a uniform occupation of the lowest subband (see Methods). The pumping is performed along $x$ by adiabatically varying $\varphi_{x}$ and we examine the resulting motion of the atoms. We locally probe the system’s transport properties by using a small atom cloud extending over approximately 20 sites in the $x$-direction. In this case, the variation of $\Omega^{y}(\varphi_{y})$ over the cloud size is negligible and the $y$-displacement per cycle is given by $\overline{\Omega}(\varphi_{y}^{(0)})\,\theta\,d_{\mathrm{l}}\,\mathbf{e}_{y}$ with $\overline{\Omega}=1/(2\pi)\oint\Omega^{x}\Omega^{y}\mathrm{d}k_{x}\mathrm{d}k_{y}\mathrm{d}\varphi_{x}$ (see Supplementary Information). From this local response, the global quantized non-linear response of an infinite system can be reconstructed by sampling all possible values of $\varphi_{y}^{(0)}$, thereby integrating over the entire 4D generalized Brillouin zone.

To probe the motion of the cloud, we measure its COM position as a function of $\varphi_{x}$. As the non-linear response results from two weak perturbations, the displacement per cycle is typically only a fraction of $d_{\mathrm{l}}$. It is thus too small to be resolved experimentally, where the number of cycles is limited due to heating. For suitable lattice parameters, however, where $\overline{\Omega}$ is strongly peaked at $\varphi_{y}^{(0)}=(l+1/2)\pi,l\in\mathbb{Z}$ (c.f. Fig. 1e), signatures of the 4D-like non-linear drift can be seen at $\varphi_{y}^{(0)}=\pi/2$ (Fig. 2).

To quantify this non-linear response, we instead use site-resolved band mapping, which measures the atom number on even ($N_{e}$) and odd sites ($N_{o}$) along $y$. This allows for an accurate determination of the average population imbalance of the double wells in the $y$-direction, $\mathcal{I}_{y}=(N_{o}-N_{e})/(N_{o}+N_{e})$. If no transitions between neighbouring unit cells along $y$ occur, $\mathcal{I}_{y}$ is directly related to the COM motion (see Methods). An example for a measurement of $\mathcal{I}_{y}(\varphi_{x})$ is shown in Fig. 3a. The measured non-linear response is smaller than expected for an ideal system. The deviation can be attributed to the appearance of doubly-occupied plaquettes and band excitations along $y$ during the pumping as well as a finite pumping efficiency along $x$ (see Methods). The latter is determined by the fraction of atoms remaining in the lowest subband during half a cycle and amounts to 98.6(2)%. Taking these imperfections into account, we find excellent agreement between the experimental data and the expected imbalance (Fig. 3a).

By fitting our model for the double well imbalance to the data (see Methods), we can extract the change of the population imbalance for ground state atoms during one cycle, $\delta\mathcal{I}_{y}^{\mathrm{gs}}=\mathcal{I}_{y}^{\mathrm{gs}}(\varphi_{x}=2\pi)-\mathcal{I}_{y}^{\mathrm{gs}}(\varphi_{x}=0)$. This slope is determined by $\overline{\Omega}$ and thus characterizes the transport properties of the lowest band. To reconstruct the quantized response of an infinite system and thereby obtain $\nu_{2}$, the measurement of $\mathcal{I}_{y}(\varphi_{x})$ is repeated for different values of $\varphi_{y}^{(0)}$. This corresponds to using the small atom cloud as a local probe at different positions along $x$ (Fig. 3b). The experimentally determined slope of the non-linear response for the lowest band agrees very well with the one expected in an ideal system (Fig. 3c).

An experimental value for the second Chern number of the lowest subband, $\nu_{2}^{\mathrm{exp}}$, can be obtained by averaging $\delta\mathcal{I}_{y}^{\mathrm{gs}}$ over $\varphi_{y}^{(0)}\in[0,2\pi[$. For this, the ideal slope is fitted to the measured one by scaling it with a global amplitude, $\left(\nu_{2}^{\mathrm{exp}}/\nu_{2}\right)\delta\mathcal{I}_{y}^{\mathrm{gs}}(\varphi_{y}^{(0)})$. For symmetry reasons, it is sufficient to restrict $\varphi_{y}^{(0)}$ to $[0,\pi[$ for $d_{\mathrm{l}}=2d_{\mathrm{s}}$ Marra et al. (2015). In this interval, the non-linear response has significant contributions only in the vicinity of $\varphi_{y}^{(0)}=\pi/2$. Fitting the data of Fig. 3c yields $\nu_{2}^{\mathrm{exp}}=1.07(8)$, in agreement with the expected value $\nu_{2}=+1$ for $d_{\mathrm{l}}=2d_{\mathrm{s}}$. The error takes into account the fit error and the uncertainties in the lattice depths and $\theta$. Note that while this approach requires prior knowledge of the response of an ideal system, $\nu_{2}^{\mathrm{exp}}$ can also be obtained directly from the measurements of $\mathcal{I}_{y}(\varphi_{x})$ without any additional input, which gives $\nu_{2}^{\mathrm{exp}}=0.9(2)$ (see Supplementary Information).

In the 4D QH system, the non-linear response scales linearly with the magnetic perturbation. The same dependence is thus expected for the non-linear response of our 2D charge pump with respect to $\theta$. We verify this by measuring the peak slope $\delta\mathcal{I}_{y}^{\mathrm{gs}}$ at $\varphi_{y}^{(0)}=\pi/2$ versus $\theta$ (Fig. 4). This provides another way to obtain the second Chern number by determining the slope of $\delta\mathcal{I}_{y}^{\mathrm{gs}}(\theta)$ (see Methods). A linear fit gives $\nu_{2}^{\mathrm{exp}}=1.01(8)$, where the error is determined as described above. Furthermore, we confirm that the peak slope at fixed $\theta$ scales with the depth of the short $y$-lattice $V_{\mathrm{s,y}}$ as expected (see Supplementary Information). In particular, the direction of the non-linear response is independent of $V_{\mathrm{s,y}}$, indicating its robustness against perturbations of the system.

In conclusion, we presented the first measurement of the 4D second Chern number in a bulk system. This first realization of a 2D topological charge pump opens the route to studying higher-dimensional QH physics in experimentally accessible setups. Extending this work, additional density-type non-linear Hall responses can be measured which are implied by the intrinsic 4D symmetry of a 2D charge pump Kraus et al. (2013). They manifest differently in the COM drift, but are characterized by the same topological invariant Price et al. (2016). Including interactions might allow for the detection of intriguing fractional phases originating in the 4D fractional QH effect Zhang and Hu (2001), similar to previous proposals for 1D charge pumps Zhu et al. (2013). Furthermore, a QH system with four extended dimensions might be realized with cold atoms Price et al. (2015) using recently demonstrated techniques for creating synthetic dimensions Mancini et al. (2015); Stuhl et al. (2015).

Note: Simultaneously with this work, complementary results on topological edge states in 2D photonic pumps have been obtained [O. Zilberberg et al., to be published].

We acknowledge insightful discussions with M. Aidelsburger and I. Carusotto. This work was funded by the European Commission (UQUAM, SIQS), the Deutsche Forschungsgemeinschaft (DIP, FOR2414) and the Nanosystems Initiative Munich. M. L. was additionally supported by the Elitenetzwerk Bayern (ExQM), H. M. P. by the European Commission (FET Proactive grant AQuS and Marie Skłodowska–Curie Action, Grant No. 656093 “SynOptic”) and the Autonomous Province of Trento (SiQuro) and O. Z. by the Swiss National Science Foundation.

Methods

Lattice configuration. All experiments are performed in a mutually orthogonal retro-reflected 3D optical lattice consisting of superlattices along $x$ and $y$ and a simple lattice in the $z$-direction. Each superlattice is created by superimposing two standing waves, a short lattice with wavelength $\lambda_{\mathrm{s}}=767$ nm and a long lattice with $\lambda_{\mathrm{l}}=2\lambda_{\mathrm{s}}$. The vertical lattice along $z$ is formed by a standing wave with $\lambda_{z}=844$ nm.

Tight-binding Hamiltonian of the 2D superlattice. In the tight-binding limit, a 2D superlattice with $d_{\mathrm{l}}=2d_{\mathrm{s}}$ is described by the 2D Rice-Mele Hamiltonian

[TABLE]

with $\hat{a}_{m_{x},m_{y}}^{\dagger}$ $(\hat{a}_{m_{x},m_{y}})$ being the creation (annihilation) operator acting on the ($m_{x}$,$m_{y}$)-th site in the $xy$-plane. The first (second) term describes the hopping between neigbouring sites along the $x$-axis ($y$-axis) and the last term contains the on-site potential of each lattice site. The long lattices lead to a modulation of the on-site energies, $(-1)^{m_{\mu}}\Delta_{\mu}/2$, and the tunnelling matrix elements, $J_{\mu}+(-1)^{m_{\mu}}\delta J_{\mu}/2$, with $\mu\in\{x,y\}$, which depend on the respective superlattice phase $\varphi_{\mu}$.

Initial state preparation for band-mapping measurements. For all sequences, a Mott insulator with quarter filling consisting of about 5000 87Rb atoms is prepared in the lowest subband of the 2D superlattice. To this end, a Bose-Einstein condensate is loaded from a crossed dipole trap into the lattice by first ramping up the blue-detuned short lattices along $x$ and $y$ to $3.0(1)E_{\mathrm{r,s}}$ during 50 ms to lower the initial density of the atom cloud. Then these lattices are switched off again within 50 ms while at the same time the vertical lattice as well as both long lattices are increased to $30(1)E_{\mathrm{r,z}}$ and $30(1)E_{\mathrm{r,l}}$, respectively, with $\varphi_{x}=0.000(5)\pi$ and $\varphi_{y}=\varphi_{y}^{(0)}$. Subsequently, doubly-occupied lattice sites are converted to singly-occupied ones (see Supplementary Information), creating a Mott insulator with unit filling and a negligible fraction of doublons. Afterwards, each lattice site is split to a four-site plaquette by ramping up the short lattices along $x$ and $y$ to their final depth of $7.0(2)E_{\mathrm{r,s}}$ and decreasing the long lattices to $20.0(6)E_{\mathrm{r,l}}$ within 5 ms.

Sequence for pumping. The superlattice phase can be controlled by slightly changing the frequency of the lasers used for generating the long lattices and thereby moving the relative position between the short and long lattice at the position of the atoms. The pumping along $x$ is performed by slowly changing $\varphi_{x}$, starting from the staggered configuration at $\varphi_{x}=0.000(5)\pi$, where the energy difference between neighbouring sites $|\Delta_{x}|$ is largest and the tunnel couplings are equal, $\delta J_{x}=0$. To minimize non-adiabatic transitions to higher bands, each pump cycle consists of three s-shaped ramps $\varphi_{x}\in\left[0,0.5\pi\right],\left[0.5\pi,1.5\pi\right]$ and $\left[1.5\pi,2\pi\right]$. This reduces the ramp speed in the vicinity of the symmetric double well configuration ($\Delta_{x}=0$) at $\varphi_{x}=(l+1/2)\pi$, $l\in\mathbb{Z}$, where the gap to the first excited band is smallest. The duration of the $\pi/2$ ramps is 7 ms and 14 ms for the ramp by $\pi$. Due to the limited tuning range of a single laser, a second laser is required for implementing multiple pump cycles, which is set to a constant phase of $\varphi_{x}=0.000(5)\pi$. At the end of each cycle, an instantaneous switch from the primary laser to the second one is made and within 5 ms the frequency of the former is ramped back to its initial value — corresponding to an identical lattice configuration. After switching back to the first laser, the next cycle continues as described above. We checked experimentally that this handover between the two lasers does not create any measurable band excitations.

Measurement of the in-situ position. To determine the non-linear COM displacement along $y$, a double-differential measurement is conducted to minimize the effect of shot-to-shot fluctuations of the atom position. In order to do this, the COM position is measured before ($y_{\mathrm{i}}$) and after the pumping ($y_{\mathrm{f}}$) and compared to a reference sequence ($y_{\mathrm{i}}^{(0)}$, $y_{\mathrm{f}}^{(0)}$). In the latter, the pumping is performed with only the short lattice along $y$ (at $V_{\mathrm{s,y}}=40(1)E_{\mathrm{r,s}}$) and therefore the non-linear response is zero. The initial position is obtained during the doublon removal sequence, where the atoms are initially prepared in the $(F=1,m_{F}=0)$ hyperfine state and one atom from each doubly-occupied site is transferred to $(F=2,m_{F}=-2)$ using microwave-dressed spin-changing collisions (see Supplementary Information). Here, $F$ denotes the total angular momentum of the atoms. In addition, we transfer 50% of the atoms on singly-occupied sites to the $F=2$ manifold as well by applying a microwave $\pi$-pulse resonant on the $(F=1,m_{F}=0)\rightarrow(F=2,m_{F}=0)$ transition. The $F=2$ atoms thus have the same density distribution as the remaining $F=1$ atoms and are imaged prior to the push-out pulse, which removes them from the lattice. The motion of the atoms due to the non-linear response is then given by $\Delta y=(y_{\mathrm{f}}-y_{\mathrm{i}})-(y_{\mathrm{f}}^{(0)}-y_{\mathrm{i}}^{(0)})$. The difference of the COM displacement along $y$ between $\theta_{1}$ and $\theta_{2}$ is defined as $\Delta r_{y}=\Delta y(\theta_{1})-\Delta y(\theta_{2})$. For the $x$-direction, it is obtained directly from $\Delta x=(x_{\mathrm{f}}-x_{\mathrm{i}})-\delta\overline{x}$ without comparing it to the reference sequence. Here, $\delta\overline{x}$ is the average displacement of all data points for a given angle, accounting for a small constant offset between the measured initial and final positions.

Relation between centre-of-mass position and double well imbalance. If there are no inter-double-well transitions along $y$, the change in the double well imbalance $\delta\mathcal{I}_{y}=\mathcal{I}_{y}(\varphi_{x})-\mathcal{I}_{y}(\varphi_{x}=0)$ can be directly related to the COM motion along $y$. The COM position in the $y$-direction is given by $y_{\mathrm{COM}}=d_{\mathrm{l}}/N\,\sum_{ij}\left[(j-1/4)N_{e,ij}+(j+1/4)N_{o,ij}\right]$, where the sum is over all unit cells, $N_{e,ij}$ ($N_{o,ij}$) is the occupation of the even (odd) sites along $y$ in the $(i,j)$-th unit cell and $N$ is the total atom number. Expressing this in terms of the total number of atoms on even and odd sites, $N_{e}=\sum_{ij}N_{e,ij}$ and $N_{o}=\sum_{ij}N_{o,ij}$, and assuming that there are no transitions between neighbouring unit cells along $y$, i.e. $\sum_{i}N_{e,ij}+N_{o,ij}$ remains constant, the change in the COM position can be written as $\delta y_{\mathrm{COM}}=y_{\mathrm{COM}}(\varphi_{x})-y_{\mathrm{COM}}(\varphi_{x}=0)=d_{\mathrm{l}}\,\delta\mathcal{I}_{y}/4$. Note that this derivation implicitly assumes that the COM of the maximally-localized Wannier functions on the lattice sites along $y$ is independent of $\varphi_{y}$, which is a valid approximation deep in the tight-binding regime. Otherwise, the proportionality factor $d_{\mathrm{l}}/2$ has to be replaced by the distance between the COM of the Wannier functions on the even and odd site of a double well.

Model for double well imbalance including experimental imperfections. To determine the slope of the non-linear response from the band mapping data, we use a simple model that takes into account band excitations and double occupation of plaquettes and the experimental pumping efficiency of the linear response. The average double well imbalance $\mathcal{I}_{y}(\varphi_{x})$ can be written as

[TABLE]

where $n_{\mathrm{gs}}$ ($n_{\mathrm{exc}}$) is the fraction of atoms on singly-occupied plaquettes in the ground (first excited) state along $y$ and $n_{\mathrm{2}}$ is the fraction of atoms on doubly-occupied plaquettes, which we assume to be in the ground state. These quantities can be determined experimentally at each point in the pumping sequence. $\mathcal{I}_{y}^{\mathrm{gs}}$, $\mathcal{I}_{y}^{\mathrm{exc}}$ and $\mathcal{I}_{y}^{\mathrm{2,gs}}$ denote the imbalance of the corresponding state, which is determined by the local phase of the $y$-superlattice at the position of the cloud along $x$, $\varphi_{y}(x_{\mathrm{COM}})$, and can be calculated using the respective double well Hamiltonian (see Supplementary Information). The COM position in turn depends on the pump parameter $\varphi_{x}$ and includes corrections for the finite pumping efficiency, $x_{\mathrm{COM}}(\varphi_{x})=\mathrm{sgn}(\varphi_{x})\sum_{i=1}^{|\varphi_{x}|/\pi}\left(2\beta_{0}\beta^{i}-\beta\right)$ for $\varphi_{x}/\pi\in\mathbb{Z}$. Here, $\beta_{0}=0.980(4)$ is the initial ground state occupation along $x$ and $\beta=0.986(2)$ is the pumping efficiency, given by the fraction of atoms that are transferred by one lattice site along $x$ during each half of a pump cycle. The main contributions limiting the pumping efficiency are band excitations in the pumping direction as well as non-adiabatic transitions between neighbouring double wells induced by the external harmonic confinement.

Fit function for non-linear response. Based on the above model, the experimental data is fitted with the function $\mathcal{I}_{y}(\varphi_{x})+\mathcal{I}_{0}$ with $\varphi_{y}\rightarrow\varphi_{y}^{\mathrm{exp}}=\varphi_{y}^{(0)}+\alpha\,(\varphi_{y}-\varphi_{y}^{(0)})$. The two fit parameters are the prefactor $\alpha$, which describes the change of the superlattice phase along $y$ with $\varphi_{x}$ compared to the ideal case $\varphi_{y}^{\mathrm{exp}}=\varphi_{y}$, and an overall offset $\mathcal{I}_{0}$. The transport properties of the lowest band are encoded in the slope of the ground state imbalance at $\varphi_{x}=0$. Knowing $\alpha$, it can be related to the ideal slope via

[TABLE]

Per cycle, this gives a change of the population imbalance for ground state atoms of

[TABLE]

Determination of the second Chern number from scaling of the non-linear response with $\theta$. The COM displacement per cycle along $y$ for an infinite system, $\delta y_{\mathrm{COM}}=\nu_{2}\theta d_{\mathrm{l,x}}$, scales linearly with the perturbing angle $\theta$. The second Chern number can thus be extracted from the slope of $\delta y_{\mathrm{COM}}(\theta)$. Having confirmed that the measured shape of $\delta\mathcal{I}_{y}^{\mathrm{gs}}(\varphi_{y}^{(0)})$ is the same as expected theoretically, the response of an infinite system at a given angle $\theta$ can be inferred from a single measurement of $\delta\mathcal{I}_{y}^{\mathrm{gs}}$ at a fixed $\varphi_{y}^{(0)}$. This holds for all angles since the shape of $\overline{\Omega}(\varphi_{y}^{(0)})$ is independent of $\theta$. To obtain $\nu_{2}$, it is therefore sufficient to determine the slope of $\delta\mathcal{I}_{y}^{\mathrm{gs}}(\theta)$ at a constant $\varphi_{y}^{(0)}$.

Supplementary Information for:

Exploring 4D Quantum Hall Physics with a 2D Topological Charge Pump

Michael Lohse1,2, Christian Schweizer1,2, Hannah M. Price3, Oded Zilberberg4 & Immanuel Bloch1,2

1**Fakultät für Physik, Ludwig-Maximilians-Universität, Schellingstraße 4, 80799 München, Germany

2**Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Straße 1, 85748 Garching, Germany

3**INO-CNR BEC Center & Dipartimento di Fisica, Università di Trento, Via Sommarive 14, 38123 Povo, Italy

4**Institut für Theoretische Physik, ETH Zürich, Wolfgang-Pauli-Straße 27, 8093 Zürich, Switzerland

S.XI Hall response of the 4D quantum Hall system

Assuming perfect adiabaticity, the Hall response of the 4D system shown in Fig. 1a, b can be evaluated from the semiclassical equations of motion for a wave packet centred at position $\mathbf{r}$ and quasimomentum $\mathbf{k}$ Xiao et al. (2010):

[TABLE]

Here, $\mathcal{E}(\mathbf{k})$ is the energy of the respective eigenstate at $\mathbf{k}$, $q$ the charge of the particle and the Einstein notation is used for the spatial indices ${\mu,\nu}\in\{w,x,y,z\}$111Note that the orientation of the axes in Fig. 1a, b is chosen such that the 4D Levi-Civita symbol is $\varepsilon_{wxyz}=+1$. The velocity of the wave packet $\mathbf{v}=\dot{\mathbf{r}}$ has two contributions: the group velocity arising from the dispersion of the band and the anomalous velocity due to the non-zero Berry curvature $\Omega^{\nu\mu}(\mathbf{k})=i\left(\left<\partial_{k_{\nu}}u|\partial_{k_{\mu}}u\right>-\left<\partial_{k_{\mu}}u|\partial_{k_{\nu}}u\right>\right)$. For a filled or homogeneously populated band, the group velocity term vanishes and with $\mathbf{E}=E_{z}\mathbf{e}_{z}$ and $\mathbf{B}=0$, the linear Hall response is given by the COM velocity

[TABLE]

where $A_{\mathrm{M}}^{xz}$ denotes the size of the magnetic unit cell and $\nu_{1}^{zx}=1/(2\pi)\oint_{\mathrm{BZ}}\Omega^{zx}\,\mathrm{d}^{2}k$ the first Chern number of the 2D QH system in the $xz$-plane. The integration is performed over the 2D Brillouin zone spanned by $k_{x}$ and $k_{z}$.

Adding the perturbing magnetic field $B_{xw}$ generates a Lorentz force acting on the moving cloud, $\hbar\dot{\mathbf{k}}=qE_{z}\,\mathbf{e}_{z}-qv_{x}^{(0)}B_{xw}\,\mathbf{e}_{w}$ Price et al. (2015). Note that this additional force can alternatively be interpreted as arising from a Hall voltage in the $w$-direction that is created by the current along $x$ in the presence of $B_{xw}$. This force in turn induces an additional anomalous velocity along $y$, giving rise to the non-linear Hall response. The resulting average velocity is then

[TABLE]

with $A_{M}$ being the size of the 4D magnetic unit cell. The second Chern number is given by $\nu_{2}=1/(4\pi^{2})\oint_{\mathrm{BZ}}\Omega^{xw}\Omega^{zy}+\Omega^{xy}\Omega^{wz}+\Omega^{zx}\Omega^{wy}\mathrm{d}^{4}k$, where $\mathrm{BZ}$ denotes the 4D Brillouin zone.

S.XII Mapping of a 2D Topological Charge Pump to a 4D Quantum Hall System

The Hamiltonian of a 2D topological charge pump for a given set of parameters $\{\varphi_{x},\varphi_{y}\}$ can be interpreted as a Fourier component of a higher-dimensional quantum Hall system. Using the approach of dimensional extension Kraus and Zilberberg (2012), a 2D charge pump can be mapped onto a 4D QH system, whose Fourier components are sequentially sampled during a pump cycle. This is demonstrated in the following for the deep tight-binding limit $V_{\mathrm{s,\mu}}\gg V_{\mathrm{l,\mu}}^{2}/(4E_{\mathrm{r,s}})$, $\mu\in\{x,y\}$, where the corresponding 4D system consists of two 2D Harper-Hofstadter-Hatsugai models Harper (1955); Azbel (1964); Hofstadter (1976); Hatsugai and Kohmoto (1990) in the $xz$- and $yw$-plane. A similar analogy can be made in the opposite limit of a vanishing short lattice, $V_{\mathrm{s,x}}\rightarrow 0$ and $V_{\mathrm{s,y}}\rightarrow 0$. In this case, each axis of the 2D lattice maps onto the Landau levels of a free particle in an external magnetic field in 2D Lohse et al. (2016). For the lowest band, these two limiting cases are topologically equivalent, i.e. they are connected by a smooth crossover without closing the gap to the first excited band. The topological invariants governing the linear and non-linear response are thus independent of the depth of the short lattices.

For non-interacting atoms in the tight-binding limit, the motion in a 2D superlattice is captured by the following Hamiltonian:

[TABLE]

Here, $\hat{a}_{m_{x},m_{y}}^{\dagger}$ $(\hat{a}_{m_{x},m_{y}})$ are the creation (annihilation) operator acting on the ($m_{x}$,$m_{y}$)-th site in the $xy$-plane, $J_{\mu}+\delta J_{\mu}^{m_{\mu}}$ with $\mu\in\{x,y\}$ are the tunneling matrix elements between neighbouring sites along $\mu$ and $\Delta_{x}^{m_{x}}+\Delta_{y}^{m_{y}}$ is the on-site potential of a given site. In the presence of the long lattices, the tunnel couplings as well as the on-site energies are modulated periodically by $\delta J^{m_{\mu}}_{\mu}$ and $\Delta_{x}^{m_{x}}+\Delta_{y}^{m_{y}}$, respectively. For the lattice configuration used in the experiment, where $d_{\mathrm{l,\mu}}=2d_{\mathrm{s,\mu}}$, these modifications can be expressed as $(-1)^{m_{\mu}}\delta J_{\mu}(\varphi_{\mu})/2$ and $(-1)^{m_{\mu}}\Delta_{\mu}(\varphi_{\mu})/2$.

In the deep tight-binding regime, $J_{x}$ and $J_{y}$ become independent of the superlattice phases and the modulations can be approximated as

[TABLE]

with $\tilde{\Phi}_{xz}=2\pi d_{\mathrm{s,x}}/d_{\mathrm{l,x}}$ and $\tilde{\Phi}_{yw}=2\pi d_{\mathrm{s,y}}/d_{\mathrm{l,y}}$. In this case, $\hat{H}_{2D}$ is equivalent to the generalized 2D Harper model Harper (1955); Roux et al. (2008) which describes the Fourier components of a 4D lattice model with two uniform magnetic fields. The 4D parent Hamiltonian can be obtained by performing an inverse Fourier transform Kraus et al. (2013)

[TABLE]

with

[TABLE]

where $\mathbf{m}=\{m_{x},m_{y},m_{z},m_{w}\}$ indicates the position in the 4D lattice. This yields

[TABLE]

The first term $\hat{H}_{xz}$ describes a 2D Harper-Hofstadter model Harper (1955); Azbel (1964); Hofstadter (1976) in the $xz$-plane with a uniform magnetic flux per unit cell $\Phi_{xz}=\Phi_{0}\tilde{\Phi}_{xz}/(2\pi)=d_{\mathrm{s,x}}/d_{\mathrm{l,x}}\,\Phi_{0}$ with $\Phi_{0}$ denoting the magnetic flux quantum:

[TABLE]

Correspondingly, the second term $\hat{H}_{yw}$ is an independent 2D Harper-Hofstadter model in the $yw$-plane with $\Phi_{yw}=d_{\mathrm{s,y}}/d_{\mathrm{l,y}}\,\Phi_{0}$. Due to the position dependence of the transverse superlattice phase $\varphi_{y}$, it also contains the magnetic perturbation, i.e. a weak homogeneous magnetic field in the $xw$-plane:

[TABLE]

with $\tilde{\Phi}_{xw}=-2\pi\theta d_{\mathrm{s,x}}/d_{\mathrm{l,y}}$. The strength of the perturbing magnetic field is thus given by

[TABLE]

where $d_{\mathrm{s,w}}$ is the lattice spacing along $w$. For $\delta J_{\mu}^{(0)}\neq 0$, the third contribution $\hat{H}_{\delta J}$ leads to the appearance of additional diagonal tunnel coupling elements in the $xz$- and $yw$-plane with an amplitude of $\delta J_{x}^{(0)}/4$ and $\delta J_{y}^{(0)}/4$, respectively. The individual 2D models without the magnetic perturbation $B_{xw}$ then correspond to the Harper-Hofstadter-Hatsugai model Hatsugai and Kohmoto (1990) with a uniform magnetic flux $\Phi_{xz}$ and $\Phi_{yw}$, respectively, i.e. the same flux as for $\delta J_{\mu}^{(0)}=0$.

S.XIII Transport properties of a 2D topological charge pump

When the pump parameter $\varphi_{x}$ is changed slowly, a particle that is initially in an eigenstate $\left|u(k_{x},\varphi_{x}(t=0),k_{y},\varphi_{y})\right>$ of the 2D superlattice Hamiltonian $\hat{H}_{2D}$ [Eq. (S.5)] will adiabatically follow the corresponding instantaneous eigenstate $\left|u(k_{x},\varphi_{x}(t),k_{y},\varphi_{y})\right>$. In absence of a tilt, $\theta=0$, the particle acquires an anomalous velocity $\Omega^{x}\partial_{t}\varphi_{x}\mathbf{e}_{x}$ during this evolution, analogous to the linear Hall response in a QH system. In this case, the Berry curvature $\Omega^{x}$ is defined in a 4D generalized Brillouin zone $(k_{x},\varphi_{x},k_{y},\varphi_{y})$:

[TABLE]

For a homogeneously populated band, the COM displacement along $x$ during one cycle, obtained by integrating the average anomalous velocity over one period, can be expressed as an integral of the Berry curvature over the 2D generalized Brillouin zone spanned by $k_{x}$ and $\varphi_{x}$. It is thus determined by the pump’s first Chern number

[TABLE]

When a tilt is present, $\theta\neq 0$, this motion along $x$ leads to a change in $\varphi_{y}$. This induces an additional anomalous velocity in the $y$-direction, giving rise to the non-linear response. Neglecting the contribution from the group velocity (which averages to zero for a homogeneously populated band), we obtain for a given eigenstate:

[TABLE]

The distribution of $\Omega^{x}\Omega^{y}$ in the 4D generalized Brillouin zone is shown in Fig. 1e for the lattice parameters used for the measurements in Fig. 3 and 4 of the main text. It exhibits a pronounced peak around $\varphi_{x}\in\{\pi/2,3\pi/2\}$ and $\varphi_{y}\in\{\pi/2,3\pi/2\}$.

For a small cloud that homogeneously populates a single band as in the experiment, the variation of $\Omega^{x}\Omega^{y}$ over the size of the cloud due to the position dependence of $\varphi_{y}$ is negligible for $\theta\ll 1$. The average velocity for the non-linear response can then be calculated by averaging Eq. (S.19) over both quasimomenta $k_{x}$ and $k_{y}$. The COM displacement after a complete cycle can be determined by integrating the velocity over one period. We can thus express the change in the COM position per cycle as

[TABLE]

If the number of pump cycles is small, the change of $\varphi_{y}$ as a result of the linear pumping response can be neglected and the non-linear displacement per cycle is very well approximated by $\delta y_{\mathrm{COM}}\approx\overline{\Omega}(\varphi_{y}^{(0)})\,\theta\,d_{\mathrm{l,x}}$.

The response of a large system with size $L_{x}\gg d_{\mathrm{l,y}}/\theta$ can be obtained by averaging Eq. (S.20) over $\varphi_{y}(x)\in[0,2\pi[$, yielding

[TABLE]

where the second Chern number $\nu_{2}$ is calculated by integrating $\Omega^{x}\Omega^{y}$ over the entire 4D generalized Brillouin zone:

[TABLE]

S.XIV Pump Path

Varying the pump parameter $\varphi_{x}$ periodically modulates the tight-binding parameters $\delta J_{x}(\varphi_{x})$ and $\Delta_{x}(\varphi_{x})$ describing the superlattice along $x$ [Eq. (S.5)]. For $d_{\mathrm{l}}=2d_{\mathrm{s}}$, the modulation of $\delta J_{x}$ and $\Delta_{x}$ is out of phase and the system therefore evolves along a closed trajectory in the $\delta J_{x}$–$\Delta_{x}$ parameter space (Fig. S1a). This pump path encircles the degeneracy point ($\delta J_{x}=0$, $\Delta_{x}=0$), where the two lowest subbands of the Rice-Mele model touch. This singularity can be interpreted as the source of the non-zero Berry curvature $\Omega^{x}$ in the generalized Brillouin zone, which gives rise to the linear pumping response. All pump paths that encircle the degeneracy can be continuously transformed into one another without closing the gap to the first excited subband and are thus topologically equivalent with respect to the linear response, i.e. the value of $\nu_{1}^{x}$ does not change.

Similarly, the tight-binding parameters $\delta J_{y}$ and $\Delta_{y}$ depend on the phase of the transverse superlattice $\varphi_{y}$. For a large cloud, all possible values of $\varphi_{y}$ and thus $\delta J_{y}$ and $\Delta_{y}$ are sampled simultaneously (Fig. S1b). During a pump cycle, the system therefore traces out a closed surface in the 4D parameter space of $\delta J_{x}$, $\Delta_{x}$, $\delta J_{y}$ and $\Delta_{y}$ (Fig. S1c). In this parameter space, the two lowest subbands touch in the two planes ($\delta J_{x}=0$, $\Delta_{x}=0$) and ($\delta J_{y}=0$, $\Delta_{y}=0$), which intersect in a single point at the origin (Fig. S1d). Analogous to the linear response, this degeneracy generates the non-zero Berry curvatures $\Omega^{x}$ and $\Omega^{y}$, which cause the non-linear motion in the $y$-direction. Due to the 4D character of the parameter space, the 4D pump path can enclose the degeneracy (Fig. S1e). Whenever this is the case, the topology of the cycle does not change and the value of $\nu_{2}$ remains the same.

To visualize the pump path in the 4D parameter space in Fig. S1, we apply the following transformation:

[TABLE]

where the tight-binding parameters are normalized by their respective maximum values. The degeneracy planes are then given by $r_{1}=-r_{2}$, $r_{3}=-r_{4}$ and $r_{1}=r_{2}$, $r_{3}=r_{4}$, respectively, i.e. they become perpendicular planes in the $(r_{1},r_{2},r_{3})$-space.

S.XV Calculation of the double well imbalance along $y$

The measurement of the population imbalance in the $y$-direction as a function of $\varphi_{x}$ for Fig. 3 and 4 of the main text is performed after an integer or half-integer number of pump cycles, i.e. $\varphi_{x}=l\pi$, $l\in\mathbb{Z}$. At these points, the superlattice along $x$ is in the staggered configuration with the maximum energy offset $|\Delta_{x}|\gg J_{x}$ and $\delta J_{x}=0$. The atoms are thus fully localized on either even or odd sites along $x$ for $\varphi_{x}=2l\pi$ or $\varphi_{x}=(2l+1)\pi$, respectively. The four-site unit cell of the 2D superlattice therefore effectively reduces to a double well along $y$.

For singly-occupied double wells, the expected imbalance in the $y$-direction for atoms in the ground ($\mathcal{I}_{y}^{\mathrm{gs}}$) and first excited state ($\mathcal{I}_{y}^{\mathrm{exc}}$) can then be calculated from the single-particle double well Hamiltonian:

[TABLE]

with $J_{y}^{0}(\varphi_{y})=J_{y}(\varphi_{y})+\delta J_{y}(\varphi_{y})/2$ and using the Fock basis for the atom on the even and odd site, $\left|1,0\right>$ and $\left|0,1\right>$, respectively.

Correspondingly, the imbalance for the ground state of a doubly-occupied double well ($\mathcal{I}_{y}^{\mathrm{2,gs}}$) can be determined using the two-particle double well Hamiltonian:

[TABLE]

in the Fock basis $\{\left|2,0\right>,\left|1,1\right>,\left|0,2\right>\}$. Here, $U$ denotes the on-site interaction energy for two atoms localized on the same lattice site.

S.XVI Removal of doubly-occupied sites

After preparing the Mott insulator with unit filling in the long lattices, sites containing two atoms are converted to singly-occupied ones using microwave-dressed spin-changing collisions Widera et al. (2005) and a resonant optical push-out pulse. For this, the lattice depths are increased to $V_{\mathrm{s,x}}=70(2)E_{\mathrm{r,s}}$, $V_{\mathrm{l,x}}=30(1)E_{\mathrm{r,s}}$, $V_{\mathrm{l,y}}=70(2)E_{\mathrm{r,l}}$ and $V_{z}=100(3)E_{\mathrm{z}}$ in 5 ms to maximize the on-site interaction energy. The atoms, which are initially in the $(F=1,m_{F}=-1)$ hyperfine state, are converted to $(F=1,m_{F}=0)$ with an adiabatic radio-frequency transfer. By ramping a magnetic offset field in the presence of a microwave field, a Landau-Zener sweep is performed that adiabatically converts pairs of $m_{F}=0$ atoms on the same lattice site to an $m_{F}=+1$ and an $m_{F}=-1$ atom via coherent spin-changing collisions. Subsequently, the $m_{F}=-1$ atoms are removed by an adiabatic microwave transfer to $(F=2,m_{F}=-2)$ followed by a resonant optical pulse after lowering the lattices to $V_{\mathrm{s,x}}=0E_{\mathrm{r,s}}$, $V_{\mathrm{l,x}}=30(1)E_{\mathrm{r,l}}$, $V_{\mathrm{l,y}}=40(1)E_{\mathrm{r,l}}$ and $V_{z}=40(1)E_{\mathrm{z}}$.

S.XVII Measurement of band excitations

Band excitations in the $y$-direction are measured by adiabatically ramping the superlattice phase $\varphi_{y}^{(0)}$ from its initial value to $\pi/2\pm 0.156(5)\pi$ and subsequently increasing the short lattice depth to $V_{\mathrm{s,y}}=40(1)E_{\mathrm{r,s}}$. In this lattice configuration, ground state atoms on both singly- and doubly-occupied plaquettes are fully localized on the lower-lying site along $y$ due to the large double well tilt $\Delta_{y}$ and the suppression of tunnelling $J_{y}^{0}\rightarrow 0$. Atoms in the excited band along $y$, on the other hand, localize on the higher-lying site and can be detected directly by measuring the resulting double well imbalance.

S.XVIII Detection of doubly-occupied plaquettes

The doublon fraction can be determined by taking advantage of the fact that two atoms in the same double well localize on the lower-lying site only at much larger double well tilts than a single atom due to the repulsive on-site interaction. For this, the double wells along $y$ are first merged to a single site by removing the short lattice and increasing the long lattice to $V_{\mathrm{l,y}}=30(1)E_{\mathrm{r,l}}$ within 5 ms. At the same time, the orthogonal lattice depths are ramped up to $V_{\mathrm{s,x}}=70(2)E_{\mathrm{r,s}}$ and $V_{z}=100(3)E_{\mathrm{r,z}}$ to increase the interaction energy. After that, $\varphi_{y}^{(0)}$ is shifted adiabatically to either $0.474(5)\pi$ or $0.431(5)\pi$ and the sites are split into double wells again by ramping up the short lattice to $V_{\mathrm{s,y}}=40(1)E_{\mathrm{r,s}}$. At $\varphi_{y}^{(0)}=0.431\pi$, both single atoms and doublons are fully localized on the lower-lying site. At $\varphi_{y}^{(0)}=0.474\pi$, on the other hand, single atoms are still very well localized on the lower site, but two atoms in the same double well localize on different sites due to the large interaction energy of $U/h=5.4$ kHz. By determining the site occupations for both phases, one can thus infer the doublon fraction from the difference in the even-odd imbalance between the two measurements.

S.XIX Alignment of the tilted superlattice

Each optical lattice is created by retroreflecting a laser beam, which is focussed onto the atoms by a lens on either side of the cloud. For the superlattices, the incoming beams of the short and long lattice are overlapped with a dichroic mirror in front of the first lens. In order to control the tilt angle $\theta$ of the long lattice along $y$, a glass block is placed in the beam path prior to the overlapping. By rotating this glass block, a parallel displacement of the incoming beam can be induced, which is then converted into an angle $\theta$ relative to the short lattice beam at the first lens. The two beams intersect at the focus point of the lens, which corresponds to the position of the atom cloud. After passing through the second lens behind the cloud, both beams are retroreflected by the same mirror. The counterpropagating beams travel along the paths of the incoming beams, thereby creating the lattice potentials with the same relative angle $\theta$.

S.XX Determination of the angle $\theta$

When the long lattice in the $y$-direction is tilted by an angle $\theta$ with respect to the short lattice, the phase of the superlattice along $y$ depends on the position along $x$. This leads to a modification of the on-site potential, which for small angles can be approximated as a linear gradient along the $x$-axis, pointing in opposite directions on even and odd sites in $y$: $\Delta_{y}^{m_{y}}(\varphi_{y})\approx\Delta_{y}^{m_{y}}(\varphi_{y}^{(0)})+(-1)^{m_{y}}\delta\,m_{x}$. The strength of the gradient is given by \delta=\pi d_{\mathrm{s}}/d_{\mathrm{l}}\,\partial\Delta_{y}/\partial\varphi_{y}\big{|}_{x=0}\,\theta for a given superlattice phase $\varphi_{y}^{(0)}$ and can thus be used to determine $\theta$. In order to do this in the experiment, a superfluid is prepared at $\mathbf{k}=0$ in a 2D lattice with $V_{\mathrm{s,x}}=13.0(4)E_{\mathrm{r,s}}$, $V_{\mathrm{s,y}}=20.0(6)E_{\mathrm{r,s}}$ and $V_{\mathrm{l,y}}=70(2)E_{\mathrm{r,l}}$. The superlattice phase $\varphi_{y}^{(0)}$ is set to either $0.344(5)\pi$ or $0.656(5)\pi$ such that the atoms are fully localized on even or odd sites along $y$, respectively. The Bloch oscillations induced by the gradient are probed by measuring the momentum distribution of the atoms after a variable hold time. The angle $\theta$ is then calculated from the average Bloch oscillation period of both phases to minimize the influence of additional residual gradients.

S.XXI Non-linear response versus lattice depth

The technique for detecting the non-linear response with site-resolved band mapping introduced in the main text allows to accurately determine the slope over a wide range of lattice parameters. To demonstrate this, we measure the slope of the non-linear response at $\varphi_{y}^{(0)}=0.500(5)\pi$ and $\theta=0.54(3)\,\mathrm{mrad}$ for various values of the transverse short lattice depth $V_{\mathrm{s,y}}$ (Fig. S2). As expected, the slope increases with larger depths as the band gap decreases and the Berry curvature $\Omega^{y}$ becomes more and more localized around $\varphi_{y}^{(0)}=(l+1/2)\pi$ with $l\in\mathbb{Z}$.

At $V_{\mathrm{s,y}}=6.25E_{\mathrm{r,s}}$, the first and second excited subband along $y$ touch for $\varphi_{y}^{(0)}=l\pi$, leading to a topological transition where the signs of the first and second Chern number of the first excited subband change from $+1$ for $V_{\mathrm{s,y}}<6.25E_{\mathrm{r,s}}$ to $-1$ for $V_{\mathrm{s,y}}>6.25E_{\mathrm{r,s}}$. This corresponds to a transition between the Landau and Hofstadter regimes Lohse et al. (2016). For the lowest band, the two regimes are topologically equivalent and the atoms thus move in the same direction. In both limits, the experimentally determined slope matches very well with the one expected in an ideal system. This nicely illustrates that the transport properties of the lowest band can be extracted correctly in both regimes, even in the presence of atoms in the first excited band.

S.XXII Direct Determination of the second Chern number

The method for determining the second Chern number from the measurement of the non-linear response versus $\varphi_{y}^{(0)}$ presented in the main text relies on prior knowledge of the response of an ideal system, both for the ground state and the excited states. While this significantly improves the accuracy for $\nu_{2}^{\mathrm{exp}}$, it can also be determined directly from the measured double well imbalance $\mathcal{I}_{y}(\varphi_{x})$ without any additional information about the system. To this end, the average change of the imbalance per cycle for the entire cloud, $\delta\mathcal{I}_{y}(\varphi_{y}^{(0)})$, is obtained from a linear fit of the differential imbalance $\mathcal{I}_{y}(\varphi_{x})-\mathcal{I}_{y}(-\varphi_{x})$ for each value of $\varphi_{y}^{(0)}$. The influence of the excitations can be reduced by restricting the fitting region to a small number of pump cycles. The response of an infinite system is reconstructed by averaging $\delta\mathcal{I}_{y}(\varphi_{y}^{(0)})$ over $\varphi_{y}^{(0)}$ using linear interpolation between the data points. When taking into account all points with $\varphi_{x}/(2\pi)\leq 3$ as well as the finite pumping efficiency along $x$, this gives $\nu_{2}^{\mathrm{exp}}=0.94(19)$ for the data from Fig. 3. Note that the linear interpolation for the discrete sampling used in Fig. 3c leads to a systematic shift of $\nu_{2}^{\mathrm{exp}}$ by +0.05.