Realization of density-dependent Peierls phases to engineer quantized gauge fields coupled to ultracold matter

TL;DR

This paper demonstrates a method to engineer density-dependent gauge fields in ultracold fermions using Floquet driving, enabling simulation of complex gauge theories with controllable Peierls phases.

Contribution

It introduces a Floquet scheme to realize density-dependent Peierls phases in ultracold atoms, advancing quantum simulation of lattice gauge theories.

Findings

Successfully implemented density-assisted tunnelling with controllable phases.

Measured Peierls phases and identified regimes with different topological invariants.

Mapped the winding structure of the gauge field around a Dirac point.

Abstract

Gauge fields that appear in models of high-energy and condensed matter physics are dynamical quantum degrees of freedom due to their coupling to matter fields. Since the dynamics of these strongly correlated systems is hard to compute, it was proposed to implement this basic coupling mechanism in quantum simulation platforms with the ultimate goal to emulate lattice gauge theories. Here, we realize the fundamental ingredient for a density-dependent gauge field acting on ultracold fermions in an optical lattice by engineering non-trivial Peierls phases that depend on the site occupations. We propose and implement a Floquet scheme that relies on breaking time-reversal symmetry (TRS) by driving the lattice simultaneously at two frequencies which are resonant with the onsite interactions. This induces density-assisted tunnelling processes that are controllable in amplitude and phase. We demonstrate techniques in a Hubbard dimer to quantify the amplitude and to directly measure the Peierls phase with respect to the single-particle hopping. The tunnel coupling features two distinct regimes as a function of the modulation amplitudes, which can be characterised by a $\mathbb{Z}_2$-invariant. Moreover, we provide a full tomography of the winding structure of the Peierls phase around a Dirac point that appears in the driving parameter space.

Figures

Equations

$${\hat{H}}_{\text{eff}}^{(l)}=-\sum _{⟨\text{i},\text{j}⟩,\sigma }\left({\hat{t}}_{⟨\text{i},\text{j}⟩,\bar{\sigma }}^{(l)}{e}^{i{\hat{\mathcal{A}}}_{⟨\text{i},\text{j}⟩,\bar{\sigma }}^{(l)}}{\hat{c}}_{\text{j}\sigma }^{†}{\hat{c}}_{\text{i}\sigma }+\text{h.c.}\right).$$

$${\hat{H}}_{\text{eff}}^{(l)}={\vec{h}}^{(l)}\cdot \vec{\sigma }+\frac{{\delta }^{(l)}}{2}\text{mathds}{1}^{2\times 2}.$$

$${\vec{h}}^{(l)}=\left(-\sqrt{2}|t_{\text{eff}}^{(l)}|\cos ({\psi }^{(l)}),\sqrt{2}|t_{\text{eff}}^{(l)}|\sin ({\psi }^{(l)}),{\delta }^{(l)}/2\right),$$

$$\Delta (\tau )=\hslash \omega K_1\cos (\omega \tau +{\varphi }_{\text{c}})+2\hslash \omega K_2\cos (2\omega \tau +2{\varphi }_{\text{c}}+{\varphi }_{\text{r}}).$$

$$t_{\text{eff}}^{(l)}=t{e}^{-il{\varphi }_{\text{c}}}\sum _m{\mathcal{J}}_{-2m+l}(K_1){\mathcal{J}}_m(K_2)e^{-im{\varphi }_{\text{r}}}$$

$$H_{\text{Dirac}}=-t{c}_K{\tilde{K}}_1{\sigma }_x+t{c}_{\varphi }{\tilde{\varphi }}_{\text{r}}{\sigma }_y.$$

$${\hat{H}}_0=-t\sum _{j,\sigma }({\hat{c}}_{j+1,\sigma }^{†}{\hat{c}}_{j\sigma }+\text{h.c.})+U\sum _j{\hat{n}}_{j\uparrow }{\hat{n}}_{j\downarrow }$$

$$\hat{V}(\tau )=-f(\tau )\sum _{j,\sigma }j{\hat{n}}_{j\sigma }$$

$$f(\tau )$$

$${\hat{H}}_{\text{eff}}^{(l)}=$$

$$t_{\text{eff},\text{R/L}}^{(l)}=t{e}^{-il{\varphi }_{\text{c}}}\sum _m{\mathcal{J}}_{-2m\pm l}(K_1){\mathcal{J}}_m(K_2)e^{\mp im{\varphi }_{\text{r}}}$$

$${\hat{H}}_{\text{eff}}^{(l)}=-\sum _{j,\sigma }\left({\hat{t}}_{⟨j,j+1⟩,\bar{\sigma }}^{(l)}{e}^{i{\hat{\mathcal{A}}}_{⟨j,j+1⟩,\bar{\sigma }}^{(l)}}{\hat{c}}_{j+1,\sigma }^{†}{\hat{c}}_{j\sigma }+\text{h.c.}\right),$$

$${\hat{t}}_{⟨i,j⟩,\sigma }^{(l)}$$

$${\hat{\mathcal{A}}}_{⟨i,j⟩,\sigma }^{(l)}$$

$$V(x,y,z)$$

$$\hat{H}(\tau )={\hat{H}}_0+\hat{V}(\tau ).$$

$${\hat{H}}_0=-t\sum _{j,\sigma }({\hat{c}}_{j+1,\sigma }^{†}{\hat{c}}_{j\sigma }+\text{h.c.})+U\sum _j{\hat{n}}_{j\uparrow }{\hat{n}}_{j\downarrow },$$

$$\hat{V}(\tau )=-f(\tau )\sum _{j,\sigma }j{\hat{n}}_{j\sigma },$$

$$f(\tau )$$

$${\hat{R}}^{(l)}(\tau )=\exp \left[-i\left(l\omega \tau \sum _j{\hat{n}}_{j\uparrow }{\hat{n}}_{j\downarrow }-F(\tau )\sum _{j,\sigma }j{\hat{n}}_{j\sigma }\right)\right]$$

$$F(\tau )$$

$${\hat{H}}_{\text{rot}}^{(l)}$$

$$=$$

$${\hat{H}}_{\text{rot}}^{(l)}(\tau )=$$

$${\hat{a}}_{⟨i,j⟩,\sigma }$$

$${\hat{b}}_{⟨i,j⟩,\sigma \text{R}}$$

$${\hat{b}}_{⟨i,j⟩,\sigma \text{L}}$$

$${\hat{H}}_{\text{rot}}^{(l)}(\tau )=$$

$${\hat{H}}_{\text{eff}}^{(l)}=$$

$$t_{\text{eff},\text{R/L}}^{(l)}$$

$$t_{\text{eff},\text{L}}^{(l)}({\varphi }_{\text{c}},{\varphi }_{\text{r}})=t_{\text{eff},\text{R}}^{(l)}({\varphi }_{\text{c}}+\pi ,{\varphi }_{\text{r}}+\pi ).$$

$${\left[t_{\text{eff},\text{R}}^{(l)}({\varphi }_{\text{c}},{\varphi }_{\text{r}})\right]}^{\ast }=t_{\text{eff},\text{R}}^{(l)}(-{\varphi }_{\text{c}},-{\varphi }_{\text{r}})$$

Full text

The fundamental manifestation of a gauge field in electromagnetism is the Lorentz force acting on charged particles. In ultracold Bose and Fermi gases, the charge neutrality of the atoms requires to engineer synthetic magnetic fields Goldman2014a ; Cooper2019 . This has been achieved for bulk systems by a rotation of the gas or a suitable coupling of momentum states via Raman lasers Cooper2008 ; Lin2009 . For a tight-binding model on a lattice, the equivalent of an Aharonov-Bohm phase can be synthesized with Peierls phases resulting from a complex-valued tunnelling matrix element. Such phases can be engineered in a Floquet approach by a suitable driving scheme Bukov2015 ; Eckardt2017 , which has been used in cold atom experiments to generate static gauge fields Aidelsburger2011 ; Miyake2013 ; Struck2013 ; Jotzu2014 . So far, these synthetic fields for atoms in optical lattices were intrinsically classical, as they did not feature a back-action from the atoms. In contrast, gauge fields appearing in nature are dynamical in the sense that they are influenced by the spatial configuration and motion of the matter field Cheng1991 ; Levin2005 ; Savary2017 . Therefore, as a first step towards the simulation of lattice gauge theories Wiese2013 ; Zohar2015 ; Dalmonte2016 ; Martinez2016 , it is necessary to implement a coupling mechanism between the gauge and matter fields. One possibility is to engineer density-dependent gauge fields by making use of interactions Edmonds2013a . Such a scheme has recently been implemented experimentally by adding a directional mean-field shift in momentum space to a Bose-Einstein condensate Clark2018 . For tight-binding models a back-action mechanism encoded in Peierls phases that depend on the occupation of the lattice sites has been suggested theoretically Keilmann2011 ; Greschner2014 ; Greschner2015 ; Bermudez2015 ; Cardarelli2016 ; Strater2016 ; Barbiero2018 .

We propose and experimentally realize a scheme based on a Floquet-engineering approach, which allows us to control both the amplitude and Peierls phase of density-dependent tunnelling matrix elements. In order to obtain non-trivial hopping phases, we explicitly break TRS by modulating the position of an optical lattice simultaneously at two frequencies $\omega/(2\pi)$ and $2\omega/(2\pi)$ Struck2012 . By tuning the Hubbard on-site interactions close to a resonance $U=l\hbar\omega$ ($l\in\mathbb{Z}$), we induce complex density-assisted tunnelling processes $t^{(l)}_{\text{eff}}=|t^{(l)}_{\text{eff}}|\exp(i\psi^{(l)})$ by exchanging photons with the drive (for realizations of real-valued density-assisted hoppings induced by a periodic modulation see Refs. Stoferle2004 ; Jordens2008 ; Greif2011 ; Chen2011 ; Ma2011 ; Meinert2016 ; Desbuquois2017 ; Gorg2018 ; Messer2018 ; Xu2018 ; Sandholzer2018 ; Schweizer2019 ). Using Floquet theory, we derive an effective static Hamiltonian describing the long-term dynamics of the system Bukov2015 (see Methods and supplementary information (SI)), which is given by

[TABLE]

Here, the operators $\hat{c}^{\dagger}_{\textbf{j}\sigma}$ and $\hat{c}_{\textbf{j}\sigma}$ create and annihilate a fermion on site j in spin state $\sigma\in\{\uparrow,\downarrow\}$, respectively. In general, both the tunneling amplitude $\hat{t}_{\langle\textbf{i},\textbf{j}\rangle,\bar{\sigma}}^{(l)}$ and phase $\hat{\mathcal{A}}_{\langle\textbf{i},\textbf{j}\rangle,\bar{\sigma}}^{(l)}$ for atoms in state $\sigma$ are operators, since they depend on the configuration of particles in the opposite state $\bar{\sigma}$, which therefore act as a link variable on the nearest-neighbor bond $\langle\textbf{i},\textbf{j}\rangle$. Using a spin-1 representation in the eigenbasis of $\hat{\tau}_{\langle\textbf{i},\textbf{j}\rangle,\sigma}^{z}=\hat{n}_{\textbf{i}\sigma}-\hat{n}_{\textbf{j}\sigma}$ where $\hat{n}_{\textbf{i}\sigma}=\hat{c}^{\dagger}_{\textbf{i}\sigma}\hat{c}_{\textbf{i}\sigma}$ is the number operator, the tunneling operators can be expressed as $\hat{t}_{\langle\textbf{i},\textbf{j}\rangle,\bar{\sigma}}^{(l)}=\text{diag}(|t_{\text{eff},\text{L}}^{(l)}|,|t_{\text{eff}}^{(0)}|,|t_{\text{eff},\text{R}}^{(l)}|)$ and $\hat{\mathcal{A}}_{\langle\textbf{i},\textbf{j}\rangle,\bar{\sigma}}^{(l)}=\text{diag}(-\psi_{\text{L}}^{(l)},\psi^{(0)},\psi_{\text{R}}^{(l)})$, respectively (see SI). Here, $l=0$ corresponds to the single-particle hopping and R and L denote density-assisted tunnelings ($|l|>0$) that involve an atom in the opposite spin state on the right and left site, respectively. Importantly, all tunneling amplitudes $|t^{(l)}_{\text{eff}}|$ and Peierls phases $\psi^{(l)}$ can be tuned independently in our scheme via the driving strengths and phases. Since $\hat{\mathcal{A}}_{\langle\textbf{i},\textbf{j}\rangle,\bar{\sigma}}^{(l)}$ is an operator, it acts as a dynamical gauge field for atoms in spin state $\sigma$. The Hamiltonian (1) is reminiscent of a lattice gauge theory and we show in the SI how to engineer dynamical $\mathbb{Z}_{2}$ and $\mathbb{Z}_{3}$ gauge fields with a two-frequency driving scheme.

In order to directly measure the matrix elements of the tunneling amplitude $\hat{t}_{\langle\textbf{i},\textbf{j}\rangle,\bar{\sigma}}^{(l)}$ and the gauge field operator $\hat{\mathcal{A}}_{\langle\textbf{i},\textbf{j}\rangle,\bar{\sigma}}^{(l)}$ in the experiment, we project the Hamiltonian (1) onto a link $\langle\textbf{i},\textbf{j}\rangle$ by realizing individual Hubbard dimers (Fig. 1a). We introduce an asymmetry between the two sites of the dimer with a static energy bias $\Delta_{0}$, such that we can selectively address the density-assisted tunneling processes $t^{(l)}_{\text{eff,R}}$ and $t^{(l)}_{\text{eff,L}}$ with the drive. If $U$ is much larger than both $\Delta_{0}$ and the static tunneling $t$, the ground state for static double wells occupied by two atoms in states $\uparrow$ and $\downarrow$, respectively, is given by the singlet $\left|\text{s}\right\rangle=\left(\left|\uparrow,\downarrow\right\rangle-\left|\downarrow,\uparrow\right\rangle\right)/\sqrt{2}$. When driving the system resonantly such that $U\approx l\hbar\omega+\Delta_{0}$, the singlet is coupled to the double occupancy state $\left|\text{d}\right\rangle=\left|0,\uparrow\downarrow\right\rangle$ via $t^{(l)}_{\text{eff}}\equiv t^{(l)}_{\text{eff,R}}$ by absorbing $l$ photons from the drive (see Fig. 1b and Supp. Fig. 1). The effective Hamiltonian in this two-level system can be written as

[TABLE]

Here, $\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the vector of the Pauli spin matrices and

[TABLE]

where $\delta^{(l)}=U-l\hbar\omega-\Delta_{0}$ is the detuning from the $l$-th order resonance (for the measurement of $t^{(l)}_{\text{eff,L}}$ at $U\approx l\hbar\omega-\Delta_{0}$ see Supp. Fig. 7). In order to directly measure $\psi^{(l)}$ we perform an interference measurement, in which the single-particle phase $\psi^{(0)}$ acts as a reference (for our parameters $\psi^{(0)}\approx 0$, see Supp. Fig. 3). Instead of comparing $\psi^{(l)}$ to the phase acquired by single atoms in spatially separated dimers (see Fig. 1a), we can conduct an interference measurement within each doubly occupied dimer by switching the state $\downarrow$ to a third internal state labeled $\rightarrow$ via a radio frequency (RF) pulse (see Fig. 1b). The interaction between the spin states $\uparrow$ and $\rightarrow$ can be set to $U=\Delta_{0}\ll\hbar\omega$, such that the atoms experience the effective Hamiltonian in Eq. (2) with $l=0$, which contains the single-particle tunnelling $t_{\text{eff}}^{(0)}$. On the Bloch sphere, the vector $\vec{h}^{(0)}$ representing the Hamiltonian $H_{\text{eff}}^{(0)}$ is pointing along the $x$-axis for $\delta^{(0)}=0$, while the resonant Hamiltonian $\vec{h}^{(l)}$ is rotated around the $z$-axis by an angle $\psi^{(l)}$ (see Fig. 1c). To characterise any quantum state $\left|\varphi\right\rangle$, we can measure for both combinations of spins the fraction of double occupancies $\mathcal{D}=\left|\langle d\,|\varphi\rangle\right|^{2}$ and singlets $\mathcal{S}=\left|\langle s\,|\varphi\rangle\right|^{2}$. Here, $\left\langle...\right\rangle$ denotes the average over the inhomogeneous distribution of $\Delta_{0}$ in different dimers resulting from the underlying harmonic trapping potential.

In the experiment, we use a harmonically confined cloud of $N=41(4)\times 10^{3}$ ultracold fermionic ${}^{40}\text{K}$ atoms in a balanced mixture of two initial internal states $\uparrow$ and $\downarrow$, which are loaded into a dimerized, three-dimensional optical lattice (see Fig. 1a). The sites constituting the dimers are connected with a static tunnelling amplitude of $t/h=260(30)\>\text{Hz}$ and are offset in energy by $\Delta_{0}/h=660(20)\>\text{Hz}$. Using a suitable loading procedure, $56(2)\%$ of the atoms occupy dimers that are populated by two opposite spins (see Methods). The interaction $U/h$ between atoms in states $\uparrow$ and $\downarrow$ can be tuned in a range between $5$ and $10\>\text{kHz}$ using a magnetic Feshbach resonance. The drive consists of a time-periodic modulation of the lattice position at two frequencies $\omega/(2\pi)=2.75\>\text{kHz}$ and $2\omega/(2\pi)=5.5\>\text{kHz}$, which in a co-moving frame correponds to a modulation of the energy offset $\Delta_{\text{tot}}(\tau)=\Delta_{0}+\Delta(\tau)$ within the dimers Desbuquois2017 with the time-dependent part

[TABLE]

Here, $K_{1}$ and $K_{2}$ are the dimensionless driving amplitudes and $\phi_{\text{c}}$ is a common phase which shifts the waveform in time without changing its shape (see Fig. 1d). It can be set to zero by choosing an appropriate origin of time. In contrast, the relative phase $\phi_{\text{r}}$ explicitly breaks TRS for $\phi_{\text{r}}\neq 0,\pi$. Therefore, it will both affect the absolute value $|t_{\text{eff}}^{(l)}|$ and, crucially, lead to a non-trivial phase $\psi^{(l)}$ that cannot be eliminated by a suitable gauge choice.

To derive the effective tunnelling matrix element for our driving scheme, we perform a high-frequency expansion in a rotating frame (see Bukov2015 and SI) and find

[TABLE]

to lowest order, where $\mathcal{J}_{m}$ is the $m$-th order Bessel function. The effective tunnel coupling is given by the interference of all multi-photon processes in which $m$ photons are absorbed from the $2\omega$ drive and $2m-l$ photons are re-emitted into the $\omega$ drive, such that the total energy added to the system is $l\hbar\omega=U-\Delta_{0}$. In the experiment, we investigate the case $l=2$, for which the leading terms of the sum can be written as $t_{\text{eff}}^{(2)}=t(\alpha^{(2)}+\beta^{(2)}e^{-i\phi_{\text{r}}}$), where $\alpha^{(2)},\beta^{(2)}>0$ depend on $K_{1}$ and $K_{2}$ and we fixed the gauge such that $\phi_{\text{c}}=0$. It can be seen that if $\phi_{\text{r}}=0$ or $\pi$ such that TRS is not broken, the tunnelling matrix element is real. Furthermore, if $\alpha^{(2)}=\beta^{(2)}$ and $\phi_{\text{r}}=\pi$, the tunnelling amplitude vanishes. Away from this singular point, $|t_{\text{eff}}^{(2)}|$ increases linearly with $\left|\alpha^{(2)}-\beta^{(2)}\right|$ and $\phi_{\text{r}}$ and is therefore forming a Dirac point in this generalized parameter space. At the same time, the Peierls phase $\psi^{(2)}$ has a vortex structure around the singularity.

In our experiment, we measure both the absolute value of the effective tunnelling on the resonance $l=2$ (see Eq. (5)) and its phase compared to the single-particle tunnelling $l=0$. In order to quantify $|t_{\text{eff}}^{(2)}|$ without a bias resulting from the inhomogeneity of the harmonic trap, we perform a Landau-Zener type measurement. Starting from a singlet state in the static system at $U/h=5.41(7)\>\text{kHz}$, we first ramp up the modulation in $5.45\>\text{ms}$ while being detuned from the resonance and subsequently sweep the interactions over the avoided crossing to $U/h=7.9(1)\>\text{kHz}$ in $20\>\text{ms}$ (see Fig. 2a). If the size of the gap at the resonance given by $2\sqrt{2}|t_{\text{eff}}^{(2)}|$ is large enough, we adiabatically follow the Floquet eigenstate and convert $\left|\text{s}\right\rangle$ to $\left|0,\uparrow\downarrow\right\rangle$ Desbuquois2017 . According to the Landau-Zener formula, the measured double occupancy fraction after the interaction sweep will be given by $\mathcal{D}=\mathcal{D}_{\text{max}}[1-\exp{(-\Gamma^{2})}]$ with $\Gamma=|t_{\text{eff}}^{(2)}|/(\kappa\cdot t)$, where $\mathcal{D}_{\text{max}}=0.56(2)$ is the maximum value of $\mathcal{D}$ given by the initial preparation. The sensitivity of the measurement is characterised by $\kappa$, which is given by $\kappa=0.15(2)$ for our interaction ramp speed. To confirm the dependence of $\mathcal{D}$ on $|t_{\text{eff}}^{(2)}|$, we benchmark our gap measurement by driving only at a single frequency $\omega/(2\pi)$ or $2\omega/(2\pi)$. For $K_{2}=0$ or $K_{1}=0$, $|t_{\text{eff}}^{(2)}|$ reduces to $\mathcal{J}_{2}(K_{1})$ or $\mathcal{J}_{1}(K_{2})$, respectively. Fig. 2b shows that the transfer fraction to the double occupancy state first increases with the magnitude of the effective tunnelling before it saturates for a gap size which corresponds to $|t_{\text{eff}}^{(2)}|\approx 0.2t$. This gives us a high sensitivity for small absolute values of the effective tunnelling.

In order to directly measure the Peierls phase $\psi^{(2)}$, we implement a scheme similar to a Ramsey experiment. We design our protocol for the two-level system depicted in Fig. 1c, in which we use the near- and off-resonant Hamiltonians represented by the vectors $\vec{h}^{(2)}$ and $\vec{h}^{(0)}$ as distinct rotation axes. Instead of scanning the evolution time of the state as in a typical Ramsey sequence, we vary the initial phase between $\vec{h}^{(2)}$ and $\vec{h}^{(0)}$ by changing the common phase $\phi_{\text{c}}$. The angle between the two rotation axes, which determines the interference fringes for the populations of $\left|\text{d}\right\rangle$ and $\left|\text{s}\right\rangle$, is given by $\psi^{(2)}=-2\phi_{\text{c}}+\psi_{\text{r}}^{(2)}(K_{1},K_{2},\phi_{\text{r}})$. In addition to $\phi_{\text{c}}$, it contains the non-trivial part of the Peierls phase $\psi_{\text{r}}^{(2)}$ given by the amplitudes and relative phase of the two-frequency modulation (see Eq. (5)). By varying $\phi_{\text{c}}$, we can extract the Peierls phase $\psi_{\text{r}}^{(2)}$ from the phase of the resulting fringes.

More precisely we first prepare an eigenstate of the Hamiltonian $H_{\text{eff}}^{(2)}$ (see Eq. (2)) for $\delta^{(2)}=0$, which is given by $(-e^{i\psi^{(2)}},1)/\sqrt{2}$. This is achieved by ramping up the drive within $5.45\>\text{ms}$ away from the resonance followed by a sweep of the interactions on resonance to $U/h=6.23(8)\>\text{kHz}$ within $100\>\text{ms}$ (see Fig. 2c). After that, we project the system onto the off-resonant Hamiltonian $H_{\text{eff}}^{(0)}$ at $\delta^{(0)}=0$. The quench is achieved by applying an RF pulse which lasts $9.5\;\mu\text{s}$ and converts $95(4)\%$ of the $\downarrow$ atoms to spin $\rightarrow$. If $H_{\text{eff}}^{(0)}$ is not (anti-)parallel to $H_{\text{eff}}^{(2)}$, the state will start to rotate around the new Hamiltonian, leading to oscillations of the singlet and double occupancy fractions (see Fig. 2c). When fixing the evolution time to the point where the Bloch vector has rotated by an angle of $\pi/2$ around $H_{\text{eff}}^{(0)}$, we observe Ramsey fringes for the observables as a function of $\phi_{\text{c}}$ given by $\mathcal{D}(\phi_{\text{c}}),\mathcal{S}(\phi_{\text{c}})=\left[1\pm\sin{\left(-2\phi_{\text{c}}+\psi_{\text{r}}^{(2)}\right)}\right]/2$ from which we can directly extract $\psi_{\text{r}}^{(2)}$ (see Fig. 2d). Due to the evolution of the state during the initial preparation of the eigenstate of $H_{\text{eff}}^{(2)}$, we measure an overall phase offset of the Ramsey fringes, which we determine in independent measurements to be $-0.15(4)\pi$ (see Methods). For all data shown in the following, the tunnelling phase was extracted from the Ramsey fringes for $\mathcal{D}$.

We begin our investigation of the effective tunnel coupling induced by the two-frequency drive by mapping out the transition for which $|t_{\text{eff}}^{(2)}|=0$ in the $K_{1}$-$K_{2}$ parameter space. As discussed above, this occurs at the TR-symmetric point $\phi_{\text{r}}=\pi$ and to lowest order for $\alpha^{(2)}=\beta^{(2)}$, i.e. $\mathcal{J}_{2}(K_{1})\mathcal{J}_{0}(K_{2})=\mathcal{J}_{0}(K_{1})\mathcal{J}_{1}(K_{2})$. Fig. 3a shows the result of the gap measurement in the $K_{1}$-$K_{2}$ parameter space following the experimental protocol in Fig. 2a,b. The gap shows a clear minimum along the diagonal, separating two distinct regions with large values of $|t_{\text{eff}}^{(2)}|$. The gap closing nicely follows the theoretical prediction derived from Eq. (5) without free parameters (see also Supp. Fig. 2). While the double occupancy fraction goes to almost zero for small values of both $K_{1}$ and $K_{2}$, the minimum is less pronounced if both amplitudes are high. In this region, the two-level approximation in Eq. (2) breaks down and the singlet is transformed into the other double occupancy state $\left|\uparrow\downarrow,0\right\rangle$ during the interaction sweep. This is demonstrated by a full numerical simulation of the gap measurement protocol (see Supp. Fig. 5) and results from the proximity of the resonances at $U=l\hbar\omega\pm\Delta_{0}$, which can be avoided by choosing higher modulation frequencies.

The line along which the gap closes separates two distinct regions in parameter space, which are characterized by a $\mathbb{Z}_{2}$-invariant (see SI). When going from top left to bottom right in Fig. 3a, the parameter $\alpha^{(2)}-\beta^{(2)}$ and hence the tunneling matrix element $t_{\text{eff}}^{(2)}$ change from negative to positive values. At the phase transition, $\alpha^{(2)}-\beta^{(2)}=0$ such that $|t_{\text{eff}}^{(2)}|=0$ (see Fig. 3b showing a cut along the black line indicated in Fig. 3a). The sign change of the effective tunnelling amplitude can be demonstrated by measuring its phase across the transition line, which exhibits a sharp jump by $\pi$ for the critical values of $K_{1}$ and $K_{2}$ where the gap closes (see Fig. 3c). This in turn proves that the gap fully closes, since the tunnelling amplitude is continuous in the modulation parameters.

After mapping out the gap closing in the parameter space of the driving amplitudes, we additionally investigate the influence of the relative modulation phase $\phi_{\text{r}}$. To this end, we always fix the parametrization in the $K_{1}$-$K_{2}$ space to be along the black line in Fig. 3a. If we expand the tunnel coupling around the point where the gap closes up to linear order in $\tilde{K}_{1}=K_{1}-K_{1,\text{crit}}$ and $\tilde{\phi}_{\text{r}}=\phi_{\text{r}}-\pi$, we find $t_{\text{eff}}^{(2)}=t[c_{K}\tilde{K}_{1}+ic_{\phi}\tilde{\phi}_{\text{r}}]/\sqrt{2}$ (see SI). Here, the numerical factors $K_{1,\text{crit}}=1.06(1)$, $c_{K}=0.537(1)$ and $c_{\phi}=0.123(1)$ depend on the $K_{1}$-$K_{2}$ parametrization. For $\delta^{(2)}=0$, the low-energy Hamiltonian around the gap closing point can therefore be written as

[TABLE]

This is a Dirac Hamiltonian in the driving parameters, which only affects the density-assisted tunnelling processes, while the single-particle hopping remains trivial.

Fig. 4a shows the measurement of the gap near the Dirac point located at $K_{1,\text{crit}}$ and $\phi_{\text{r}}=\pi$. It demonstrates that $|t_{\text{eff}}^{(2)}|$ has a clear minimum at the singularity and increases away from it. The gap closes at $\phi_{\text{r}}=1.03(6)\pi$ as expected from theory (see analytical results in Fig. 4b and Supp. Fig. 4 and a numerical simulation of the gap measurement in Supp. Fig. 6). In addition, Fig. 4c shows a full tomography of the tunnelling phase $\psi_{\text{r}}^{(2)}$ around the Dirac point. It has a vortex structure and the phase increases by $2\pi$ when going clockwise around the singularity. For high values of $K_{1}>K_{1,\text{crit}}$, $\psi_{\text{r}}^{(2)}$ only changes little as a function of the relative phase (see Fig. 4d). In this regime, the driving component at $\omega$ is dominant, which corresponds to the lower right corner in Fig. 3a. Around the critical point $K_{1}=K_{1,\text{crit}}$ and $\phi_{\text{r}}=\pi$, the tunnelling phase is very sensitive to the exact driving parameters and suddenly jumps from [math] to $\pi$ when lowering $K_{1}$ (see also Fig. 3c). For $K_{1}<K_{1,\text{crit}}$, we enter the upper left region in Fig. 3a and suddenly observe a running phase $\psi_{\text{r}}^{(2)}$, which means that the state vector is winding once around the Bloch sphere when $\phi_{\text{r}}$ is swept from [math] to $2\pi$ (see Fig. 4e). While the detailed shape of the phase vortex depends on the parametrization of the driving waveform in Eq. (4), the phase difference between two configurations at the TR-symmetric points is forced to be either [math] or $\pi$. This quantity is therefore a $\mathbb{Z}_{2}$-invariant which can be used to characterise the corresponding regimes (see Fig. 3a and SI).

In future experiments, the full control over both the amplitude and Peierls phase of the density-dependent tunnelling matrix element demonstrated in this work can be further extended by introducing a temporal or spatial dependence for the driving parameters, which maps the Dirac point into another parameter space. In addition, it is straightforward to couple the individual dimers by allowing for tunnelling in all directions in order to study the intriguing interplay between the interaction induced gauge field and the atomic density. Recent experiments have shown that driven many-body systems can be well understood in the effective Hamiltonian picture and that problems associated with interacting Floquet systems such as heating can be mitigated in certain lattice geometries Gorg2018 ; Messer2018 ; Sandholzer2018 . In higher dimensions, a variety of phenomena related to density-dependent gauge fields could be studied, such as anyonic statistics in one dimension Keilmann2011 ; Greschner2015 ; Cardarelli2016 ; Strater2016 or flux attachment Barbiero2018 . Finally, the symmetry between atoms in different internal states can be broken by using a spin-selective drive Jotzu2015 , such that distinct gauge and matter particles can be identified. As shown in the SI, using such a driving scheme involving two modulation frequencies enables to engineer both a dynamical $\mathbb{Z}_{2}$ gauge field $\hat{\mathcal{A}}_{\langle\textbf{i},\textbf{j}\rangle}=\text{diag}\,(0,\pi)$, which only requires real-valued tunneling matrix elements Barbiero2018 ; Schweizer2019 , as well as a $\mathbb{Z}_{3}$ gauge field $\hat{\mathcal{A}}_{\langle\textbf{i},\textbf{j}\rangle}=\text{diag}\,(0,2\pi/3,4\pi/3)$.

METHODS

.1 Driving scheme and effective Hamiltonian

In the following, we derive the effective Hamiltonian for a Fermi-Hubbard model, which is driven at two frequencies that are resonant with the onsite interaction $U$. We show that this scheme is suited to control both the amplitude and Peierls phase of density-assisted tunneling processes independently from the single-particle hopping. For simplicity, we discuss the case of a one-dimensional lattice. For a more details and extensions of the scheme to realize dynamical $\mathbb{Z}_{2}$ and $\mathbb{Z}_{3}$ gauge fields see SI.

The full time-dependent Hamiltonian can be written as $\hat{H}(\tau)=\hat{H}_{0}+\hat{V}(\tau)$, where the static part corresponds to the usual Fermi-Hubbard model

[TABLE]

and the drive is given by

[TABLE]

in a frame that is co-moving with the shaken lattice. For the two-frequency modulation scheme, the oscillating inertial force reads

[TABLE]

in analogy to the site offset in the case of a double well (see Eq. (4) of the main text). Using Floquet theory, we can derive an effective static Hamiltonian around the resonance $U=l\hbar\omega$ ($l\in\mathbb{Z}$) Bukov2015 , which yields to lowest order

[TABLE]

Here, single-particle tunnelings associated with the operator $\hat{a}_{\langle i,j\rangle,\sigma}=(1-\hat{n}_{i\sigma})(1-\hat{n}_{j\sigma})+\hat{n}_{i\sigma}\hat{n}_{j\sigma}$ are described by the matrix element $t_{\text{eff}}^{(0)}$, while density-assisted processes that create or annihilate a double occupancy via hopping to the right (R) or left (L) corresponding to the operators $\hat{b}_{\langle i,j\rangle,\sigma\text{R}}=(1-\hat{n}_{i\sigma})\hat{n}_{j\sigma}$ and $\hat{b}_{\langle i,j\rangle,\sigma\text{L}}=\hat{n}_{i\sigma}(1-\hat{n}_{j\sigma})$, respectively, are governed by $t_{\text{eff},\text{R/L}}^{(l)}$. The matrix elements are given by

[TABLE]

for $l\in\mathbb{Z}$ (see Eq.(5) of the main text).

For $U=l\hbar\omega$, the effective Hamiltonian in Eq. (M4) can be written as in Eq. (1) in the main text

[TABLE]

where the tunneling amplitudes and phases are described by the operators

[TABLE]

In Eq. (M6) it becomes apparent that atoms in the opposite spin state $\bar{\sigma}$ act as link variables and determine both the hopping amplitude and phase for atoms in state $\sigma$. Since there are three distinct tunneling processes corresponding to the operators $\hat{b}_{\langle i,j\rangle,\sigma\text{L}}$, $\hat{a}_{\langle i,j\rangle,\sigma}$ and $\hat{b}_{\langle i,j\rangle,\sigma\text{R}}$ for which the occupation difference between sites $i$ and $j$ equals $-1$, [math] or $1$, respectively, it is natural to express the operators in Eqs. (M7) and (M8) in the eigenbasis of the spin-1 operator $\hat{\tau}_{\langle i,j\rangle,\sigma}^{z}=\hat{n}_{i\sigma}-\hat{n}_{j\sigma}$, which has eigenvalues $m_{\langle i,j\rangle,\sigma}\in\{-1,0,1\}$. In this basis, the tunneling and gauge field operators are given by $\hat{t}_{\langle i,j\rangle,\sigma}^{(l)}=\text{diag}(|t_{\text{eff},\text{L}}^{(l)}|,|t_{\text{eff}}^{(0)}|,|t_{\text{eff},\text{R}}^{(l)}|)$ and $\hat{\mathcal{A}}_{\langle i,j\rangle,\sigma}^{(l)}=\text{diag}(-\psi_{\text{L}}^{(l)},\psi^{(0)},\psi_{\text{R}}^{(l)})$, respectively. All tunneling amplitudes and phases can be tuned independently via the driving amplitudes $K_{1}$ and $K_{2}$ and the relative phase $\phi_{\text{r}}$.

.2 Optical lattice

We perform our experiments in a three-dimensional, optical lattice which is formed by a combination of four orthogonal, retro-reflected laser beams at a wavelength of $\lambda=1064\>\text{nm}$ (for more details, see earlier work Tarruell2012 ). While the beams $\overline{\text{X}}$ and Y are effectively not interfering with any other beam due to a frequency detuning, the beams X and Z are interfering with each other and are actively phase-stabilized to $\varphi=0.00(3)\pi$. The resulting potential for the atoms is given by

[TABLE]

where $k=2\pi/\lambda$ and $V_{\overline{X},X,Y,Z}$ are the lattice depths in units of the recoil energy $E_{\text{R}}=h^{2}/2m\lambda^{2}$ of each laser beam in the three different directions $x,y,z$ ($h$ is the Planck constant and $m$ the mass of the atoms). Both the interference term $\alpha=1.01(1)$ and the individual lattice depths $V_{\overline{X},X,Y,Z}$ are calibrated via amplitude modulation of the lattice depth with a ${}^{87}\text{Rb}$ Bose-Einstein condensate. To account for systematic errors from the calibration and fluctuations of the lattice depths, we include a relative error of $4.9\%$ (for $\overline{\text{X}}$ and X) or $2.3\%$ (for Y and Z) on the lattice depths for the calculation of the tight-binding parameters. In addition, we take into account a relative uncertainty of the magnetic field of $10^{-4}$. In the final lattice configuration, the depths are given by $V_{\overline{X},X,Y,Z}=[22(1),1.00(5),38.8(9),29.3(7)]\;E_{\text{R}}$. The corresponding potential consists of an array of individual double wells aligned in the $x$-direction with an intra-dimer tunnel coupling $t/h=260(30)\>\text{Hz}$. All dynamics between different dimers are suppressed by adjusting the inter-dimer tunnelling amplitudes in all spatial directions to be below $t/h=2\>\text{Hz}$. In addition, the phase $\theta$ in Eq. (.2) is adjusted to $1.0050(1)\pi$, which introduces a static energy bias between the two sites of the double well of $\Delta_{0}/h=660(20)\>\text{Hz}$. Due to the underlying harmonic confinement with trapping frequencies of $\omega_{x,y,z}/(2\pi)=[121(2),107(1),151(2)]\>\text{Hz}$, an additional inhomogeneous site offset is introduced in the dimers. For a double well which is located at a typical distance of 20 lattice sites away from the center of the trap, the additional tilt is on the order of $\Delta_{0}/2$.

.3 Preparation of the atoms in the double wells

The preparation procedure of two distinguishable Fermions in the double wells is very similar to earlier work Desbuquois2017 . In brief, the starting point of our experiment is a balanced mixture of atoms in the $F=9/2,m_{F}=-9/2$ and $F=9/2,m_{F}=-7/2$ hyperfine states of ${}^{40}\text{K}$ (called $\uparrow$ and $\rightarrow$ in the main text), which are confined in an optical harmonic trap. After evaporatively cooling the atoms, we end up with $N=41(4)\times 10^{3}$ atoms at a temperature $T/T_{\text{F}}=0.09(2)$ ($T_{\text{F}}$ denotes the Fermi-temperature). After this, we tune the s-wave scattering length between the atoms to be very strongly attractive $a\rightarrow-\infty$ by means of a magnetic Feshbach resonance located at $202.1\>\text{G}$. Then, we first load the atoms into a checkerboard lattice with $V_{\overline{X},X,Y,Z}=[0,3,7,3]\;E_{\text{R}}$ within $200\>\text{ms}$ followed by a second ramp to a deep checkerboard lattice with $V_{\overline{X},X,Y,Z}=[0,30,40,30]\;E_{\text{R}}$ in $20\>\text{ms}$. Due to the strong attractive interactions, $56(2)\%$ of the atoms form a doubly occupied site, while no site is occupied by more than two atoms due to Pauli-blocking. The next step is to perform a Landau-Zener sweep with a radio frequency (RF) pulse to transfer the atoms in the $F=9/2,m_{F}=-7/2$ state to the $F=9/2,m_{F}=-5/2$ state (called $\downarrow$ in the main text). The interactions of the mixture in the $m_{F}=-9/2$ and $m_{F}=-5/2$ states can be tuned by a second magnetic Feshbach resonance at $224.2\>\text{G}$, which allows us to access strong repulsive interactions with $a>175\;a_{0}$ ($a_{0}$ is the Bohr radius). We adjust the scattering length to $a=237(1)\;a_{0}$ and subsequently split the wells of the checkerboard lattice into two sites by ramping to the final lattice configuration with $V_{\overline{X},X,Y,Z}=[22(1),1.00(5),38.8(9),29.3(7)]\;E_{\text{R}}$ in $20\>\text{ms}$. Here, the atoms interact with an on-site interaction energy $U/h=5.41(7)\>\text{kHz}$ and this is the starting point of the experiments.

.4 Periodic driving

The driving is implemented as in previous work Desbuquois2017 ; Gorg2018 . We sinusoidally modulate the position of the retro-reflecting mirror with a piezoelectric actuator along the direction of the dimers. As a result, the entire lattice potential is moving in space and the time-dependent potential for the two-frequency drive is given by $V(x,y,z,\tau)\equiv V(x-A_{1}\cos(\omega\tau+\phi_{\text{c}})-A_{2}\cos(2\omega\tau+2\phi_{\text{c}}+\phi_{\text{r}}),y,z)$. Here, $A_{n}$ are the real-space amplitudes of the modulation at frequency $n\omega/(2\pi)=n\cdot 2.75\>\text{kHz}$ ($n=1,2$). They are related to the normalised drive amplitude by $K_{n}=m\,A_{n}\,\omega\,d/\hbar$, where $d$ is the distance between the two sites of the double well. In our lattice geometry, $d\neq\lambda/2$ and the distance has to be calculated for the specific lattice geometry used in the experiment. For that, we calculate the maximally localised Wannier functions as eigenstates of the band-projected position operator and extract their center-of-mass position. For our lattice configuration $V_{\overline{X},X,Y,Z}=[22(1),1.00(5),38.8(9),29.3(7)]\;E_{\text{R}}$ we find $d=0.792(7)\cdot\lambda/2$. Furthermore, $\phi_{\text{c}}$ and $\phi_{\text{r}}$ are the common and relative phase of the two-frequency drive. During the modulation, we make sure that the phase relation between the X- and Z-beams is stabilized to $\varphi=0.00(3)\pi$ by modulating the frequency of the respective incoming laser beams using acousto-optical modulators. This phase modulation of the incoming beams is also used as a calibration of the driving amplitude and phase of the retro-reflecting mirror. From this we infer an uncertainty on the relative phase $\phi_{\text{r}}$ of at most $0.01\pi$, while the relative error on the amplitudes $K_{n}$ is $0.5\%$. An additional uncertainty for the amplitudes results from the imprecise knowledge of the site distance $d$ coming from the uncertainty of the lattice calibration. Due to a residual experimental mismatch of the phase modulations of the incoming and retro-reflected beams, the interference amplitude of the lattice is periodically reduced by at most $0.4\%$. To enter the driven regime, we linearly ramp up the modulation within $5.454\>\text{ms}$ and subsequently keep fixed driving amplitudes $K_{1}$ and $K_{2}$. For the gap measurements, we ramp the interactions while modulating the lattice from the initial value $U/h=5.41(7)\>\text{kHz}$ across the resonance to the final interaction $U/h=7.9(1)\>\text{kHz}$ within $20\>\text{ms}$. Alternatively, for the measurement of the tunnelling phase, we prepare an eigenstate of the resonantly driven double well by ramping the interactions to the resonance at $U/h=6.23(8)\>\text{kHz}$ within $100\>\text{ms}$.

.5 Experimental measurement of the Peierls phase

For the measurement of the tunnelling phase, we first prepare an eigenstate on the resonance $U=2\hbar\omega+\Delta_{0}$ as described above, followed by a projection onto an off-resonantly driven double well by switching the internal state of the atoms with an RF pulse. To realize the near-resonant condition, we work with a $\{\uparrow,\downarrow\}$-pair of atoms in the $m_{F}=-9/2$ and $m_{F}=-5/2$ states which are strongly interacting. Afterwards, we switch to a $\{\uparrow,\rightarrow\}$-pair of atoms in the $m_{F}=-9/2$ and $m_{F}=-7/2$ states, which have a much weaker on-site interaction energy. Importantly, we have to simultaneously match the two resonance conditions for the interactions $U_{-5/2,-9/2}=2\hbar\omega+\Delta_{0}$ and $U_{-7/2,-9/2}=\Delta_{0}$ at the same strength of the magnetic offset field. For our lattice configuration and choices of $\omega/(2\pi)=2.75\>\text{kHz}$ and $\Delta_{0}/h=660(20)\>\text{Hz}$, this condition is fulfilled for a magnetic field of $B_{\text{res}}=210.82(2)\>\text{G}$. Here, the scattering lengths are $a_{-5/2,-9/2}=273(1)\;a_{0}$ and $a_{-7/2,-9/2}=25(2)\;a_{0}$ and the corresponding on-site interactions are $U_{-5/2,-9/2}/h=6.23(8)\>\text{kHz}$ and $U_{-7/2,-9/2}/h=0.56(5)\>\text{kHz}$. To switch between the two different regimes, we transfer the atoms from the $m_{F}=-5/2$ to the $m_{F}=-7/2$ state with a fidelity of $95(4)\%$ by applying an RF pulse with a duration of $9.5\;\mu\text{s}$ and a frequency of $48.692\>\text{MHz}$. After the interaction quench, the quantum state will start to rotate around the new off-resonant Hamiltonian on the Bloch sphere with a frequency of $2\sqrt{2}|t_{\text{eff}}^{(0)}|/h$. To measure the Ramsey fringes, we fix the evolution time to $\tau=h/(8\sqrt{2}|t_{\text{eff}}^{(0)}|)$ where the rotation angle is equal to $\pi/2$. Since $|t_{\text{eff}}^{(0)}|$ is changing as a function of our driving parameters $K_{1}$, $K_{2}$ and $\phi_{\text{r}}$ (see Supp. Fig. 3), we have to adjust the timing for each choice of parameters. We do this experimentally by projecting a pure singlet state with the RF pulse onto the off-resonant Hamiltonian, which results in coherent oscillations between the singlet and double occupancy states. From these oscillation, we extract the $\pi/2$ time for a certain set of driving parameters and interpolate between them.

.6 Fit of the Ramsey fringes

For the fringes, we perform 3 independent measurements of the final double occupancy fraction for 7 different values of the common phase $\phi_{\text{c}}$ between 0 and $\pi$ (see Fig. 2d). To extract the tunnelling phase $\psi_{\text{r}}$, we fit the resulting double occupancy fringe with a function $\mathcal{D}(\phi_{\text{c}})=A\sin{(2\phi_{\text{c}}+\psi_{\text{r}})}+o$, where the period is fixed to $\pi$. To estimate the error, we use a resampling method which assumes that the measurement results for each value of $\phi_{\text{c}}$ follow a normal distribution according to the measured values of the mean and standard deviation of $\mathcal{D}$. Afterwards, we randomly sample a value for the double occupancy fraction for each common phase and refit the dataset. We repeat this procedure $1500$ times while additionally varying the initialization values for the fit parameters $A$ and $o$ by $\pm 10\%$. The mean +($-$) standard deviation of the distribution of phases fitted on the resampled data is used as an upper (lower) bound for the fitted value of $\psi_{\text{r}}$ of the measured data, which is expressed in asymmetric error bars in Figs. 3c and 4d,e. The same resampling method is also employed to estimate the uncertainty on the center position of the Lorentzian fits that are used to determine the gap closing (see Figs. 3b, 4b and Supp. Fig. 7b).

.7 Phase offset of the Ramsey fringes

In the measurements of the tunnelling phase, we observe an overall offset, i.e. the phase of the Ramsey fringes is not vanishing for $\psi_{\text{r}}=0$. This can be explained by the evolution of the state during the adiabatic preparation of the eigenstate of $H_{\text{eff}}^{(2)}$. In particular, the relative phase between the singlet and double occupancy states is not only given by $-2\phi_{\text{c}}+\psi_{\text{r}}$, but it has an additional dynamical phase contribution in the lab frame given by $-2\omega\tau$. Therefore, in addition to the slow adiabatic following to the equator of the Bloch sphere (see Fig. 2c), the state vector rotates at a frequency of $\omega/(2\pi)$ around the $z$-axis. Even when fixing the total preparation time of the eigenstate to a multiple of the driving period, any residual detuning from the resonance will lead to a modified rotation frequency and therefore to a finite phase accumulation up to the point at which the RF pulse is applied. Since the preparation takes hundreds of driving cycles, this phase offset can be significant. Furthermore, finite frequency effects and dynamics that depend on the exact launching protocol of the drive lead to additional phase shifts. To calibrate the resulting phase offset in the experiment, we take 4 reference Ramsey fringes for a single frequency drive with $\omega/(2\pi)=5.5\>\text{kHz}$, both for positive and negative site offsets $\Delta_{0}$. For this single-frequency drive, the non-trivial contribution $\psi_{\text{r}}$ vanishes, which allows us to directly measure the phase offset. From these measurements, we obtain an offset of $-0.15(4)\pi$ (uncertainty denotes the standard error).

.8 Detection

The detection of the double occupancy and singlet fractions is similar to earlier work Jordens2008 ; Greif2013 . To characterise the state of the atoms, we first freeze all dynamics by quickly ramping up the tunnelling barrier in the double well within $100\;\mu\text{s}$ to a $V_{\overline{X},X,Y,Z}=[30,0,40,30]\;E_{\text{R}}$ cubic lattice. We can detect double occupancies both for a $\{\uparrow,\downarrow\}$-pair and a $\{\uparrow,\rightarrow\}$-pair of atoms. For this, we ramp down the magnetic field below the Feshbach resonance of the $m_{F}=-9/2$ and $m_{F}=-7/2$ atoms. We then selectively transfer one of the atoms forming the double occupancy from the $m_{F}=-5/2$ to the $m_{F}=-7/2$ state (or vice versa) with an RF sweep by making use of the interaction shift. We can count the number of atoms in each $m_{F}$-state by applying a Stern-Gerlach pulse during a time-of-flight expansion followed by absorption imaging. To detect singlets, we apply a magnetic field gradient after the lattice freeze, which leads to coherent oscillations between the singlet and triplet state. After properly adjusting the evolution time, we detect the singlet state by merging two adjacent sites by going to a $V_{\overline{X},X,Y,Z}=[0,30,40,30]\;E_{\text{R}}$ checkerboard lattice. In this process, the singlet state will be adiabatically tranformed to a double occupancy in the final lattice, which we can detect as outlined above.

.9 Theoretical treatment of the driven double well

We perform both analytical and numerical studies of a double well subject to a two-frequency drive (for an analytical derivation of the effective Hamiltonian and tunnelling matrix element see supplementary material, for details about the numerical simulation see Desbuquois2017 ). To calculate the numerical quasi-energy spectrum and the Floquet-eigenstates of the time-dependent problem (see Supp. Figs. 1a-c and 5e-h), we use a Trotter decomposition to compute the evolution operator over one modulation cycle. The content of double occupancy and singlet states for a given Floquet-eigenstate $\left|f_{i}\right\rangle$ ($i=1,...,4$) is then given by $\mathcal{D}_{i}=\left|\left\langle\uparrow\downarrow,0\right.\left|f_{i}\right\rangle\right|^{2}+\left|\left\langle 0,\uparrow\downarrow\right.\left|f_{i}\right\rangle\right|^{2}$ and $\mathcal{S}_{i}=\left|\left\langle\text{s}\right.\left|f_{i}\right\rangle\right|^{2}$, respectively. In addition, we perform a numerical simulation of the full gap measurement protocol described above (see Supp. Figs. 5a-c and 6). These calculations capture the full time-dependence of the system, i.e. the drive at frequencies $\omega/(2\pi)$ and $2\omega/(2\pi)$ as well as the ramps of the driving amplitudes $K_{1}$ and $K_{2}$ and the interaction $U$. In detail, we initiate the system in a singlet state at $U=5.41\>\text{kHz}$ and first increase the amplitudes $K_{1}$ and $K_{2}$ of the two-frequency drive within $5.454\>\text{ms}$ to their final values at a fixed relative phase $\phi_{\text{r}}$ and $\phi_{\text{c}}=0$. Then, $U$ is ramped to $U=7.9\>\text{kHz}$ within $20\>\text{ms}$. To determine the double occupancy in the final state $\left|\varphi_{\text{fin}}\right\rangle$, we compute both overlaps with the double occupancy states $\left|\left\langle\uparrow\downarrow,0\right.\left|\varphi_{\text{fin}}\right\rangle\right|^{2}$ and $\left|\left\langle 0,\uparrow\downarrow\right.\left|\varphi_{\text{fin}}\right\rangle\right|^{2}$, respectively. The results for the same driving parameters as in Fig. 3 (Fig. 4a) are shown in Supp. Fig. 5a-c (Supp. Fig. 6).

.10 Data availability

All data files are available from the corresponding author upon request. Source Data for Figs. 2-4 and Supp. Figs. 1d and 7 are provided with the online version of the paper.

.11 Code availability

The source code for the fit of the Ramsey fringes is available from the corresponding author upon request.