# Scalable, ab initio protocol for quantum simulating SU($N$)$\times$U(1)
Lattice Gauge Theories

Federica Maria Surace<EMAIL_ADDRESS>Department of Physics and Institute
for Quantum Information and Matter, California Institute of Technology,
Pasadena, California 91125, USA  Pierre Fromholz  The Abdus Salam
International Centre for Theoretical Physics (ICTP), strada Costiera 11, 34151
Trieste, Italy  International School for Advanced Studies (SISSA), via Bonomea
265, 34136 Trieste, Italy  Francesco Scazza  Department of Physics, University
of Trieste, 34127 Trieste, Italy  Istituto Nazionale di Ottica del Consiglio
Nazionale delle Ricerche (CNR-INO), 34149 Basovizza-Trieste, Italy  Marcello
Dalmonte  The Abdus Salam International Centre for Theoretical Physics (ICTP),
strada Costiera 11, 34151 Trieste, Italy  International School for Advanced
Studies (SISSA), via Bonomea 265, 34136 Trieste, Italy

###### Abstract

We propose a protocol for the scalable quantum simulation of
SU($N$)$\times$U(1) lattice gauge theories with alkaline-earth like atoms in
optical lattices in both one- and two-dimensional systems. The protocol
exploits the combination of naturally occurring SU($N$) pseudo-spin symmetry
and strong inter-orbital interactions that is unique to such atomic species. A
detailed ab initio study of the microscopic dynamics shows how gauge
invariance emerges in an accessible parameter regime, and allows us to
identify the main challenges in the simulation of such theories. We provide
quantitative results about the requirements in terms of experimental stability
in relation to observing gauge invariant dynamics, a key element for a deeper
analysis on the functioning of such class of theories in both quantum
simulators and computers.

##  1 Introduction

Figure 1: (a) Scheme of the 1+1d QLM: fermionic matter (in blue) resides on
the sites of a one-dimensional lattice; the gauge field is representd by
fermionic operators (rishons) that act on the two ends of a link (in orange).
The generators of the gauge symmetry are operators acting on a single block
(matter site and its neighboring gauge sites). (b) The matter/gauge sites
correspond to the minima of an optical lattice for alkaline-earth-like atoms
in the $g/e$ state respectively. (c) Scheme of QLM in 2+1D on a square
lattice. (d) The $g/e$ states correspond to the electronic ground state
(${}^{1}S_{0}$) and metastable clock state (${}^{3}P_{0}$) of the atoms. There
are $N_{I}$ nuclear spin states decoupled from the electronic degrees of
freedom.

Over the last ten years, there has been growing interest in utilizing quantum
simulators to tackle the physics of gauge theories [1, 2, 3, 4, 5, 6] \- the
back-bone of the standard model of particle physics. The main motivation
behind this research line is the fact that gauge theories represent formidable
computational challenges [7]: while classical simulations have been
tremendously successful in qualitatively and quantitatively clarifying key
aspects of quantum chromodynamics (including, e.g., its phase diagram at zero
chemical potential [8, 9] and its low-lying energy spectrum [10]), several
pivotal scenarios, such as real-time dynamics and finite baryon density,
remain inaccessible to any known classical algorithm.

Starting from early theoretical proposals [5] and from the pioneering trapped
ion experiments of Ref. [11], the quantum simulation of Abelian lattice gauge
theories (LGTs) [12] has already been shown to scale beyond simple building
blocks [13, 14, 15, 16, 17] to regimes where classical simulations are
challenging, if not prohibitive [18, 19, 20]. Oppositely, while digital
schemes for non-Abelian LGTs have been already demonstrated in small scale
experiments in Ref. [21, 22], no scalable non-Abelian LGTs quantum simulators
exist to date. Such dearth is particularly severe in view of the motivations
above, that are mostly dealing with SU$(N)$ gauge theories. The main roadblock
for quantum simulating non-Abelian models stems from the fundamentally more
complex nature of non-Abelian Gauss law in terms of quantum engineering: in
particular, it is extremely challenging to constrain the dynamics of a quantum
simulator in such a way that only gauge invariant states are populated.
Presently, there exist only a few proposals that have addressed this [23, 24,
25, 26, 27], focusing on conceptual aspects of such challenge.

Here, we pursue a different approach. Building on impressive experimental
developments over the last decade in harnessing quantum gases of fermionic
alkaline-earth-like atoms in optical lattices [28, 29, 30, 31, 32, 33, 34, 35,
36, 37], we present a proposal for the quantum simulation of
SU($N$)$\times$U(1) that combines a new conceptual framework to engineer gauge
invariant dynamics, with an ab initio description of the physical setup. This
combined approach has two key advantages, already demonstrated along a similar
route for Abelian theories [38]: firstly, it allows us to immediately identify
potential challenges in terms of microscopic constraints related to the
quantum hardware (e.g., inability to realize the desired background
potentials, etc.), that are known to be particularly relevant in the context
of quantum engineering of LGTs; and secondly, it enables us to make precise
predictions in terms of energy scales of the quantum simulator, and to
properly frame the theoretical issues of imperfections in realizing gauge
invariant dynamics.

##  2 Setup and model dynamics

Here, we are specifically interested in realizing quantum link models (QLMs)
[39, 40, 41, 42]: this formulation of lattice gauge theories is particularly
advantageous for both quantum computing and simulations [43] due to its finite
dimensional Hilbert spaces.

The main ingredients of the setup we propose are schematically depicted in
Fig. 1. We use cold alkaline-earth fermionic atoms in their stable
${}^{1}S_{0}\equiv g$ and metastable ${}^{3}P_{0}\equiv e$ states, with
hyperfine spin sublevels satisfying a SU($N_{\text{I}}$) symmetry [44, 28],
where $N_{I}=2I+1$ and $I$ is the nuclear spin of the atomic species. We
select $N\leq N_{\text{I}}$ of these $N_{\text{I}}$ hyperfine states when
initializing the system. Only these $N$ spin states are involved in the
dynamics due to the global SU($N_{\text{I}}$) symmetry leading to an effective
SU($N$) symmetry. The atoms are trapped via a combination of optical lattice
potentials: within this setting, $g$ atoms correspond to matter fields, while
double wells of $e$ atoms serve as gauge fields. In particular, the gauge
field is represented by fermionic operators (rishons), and the two sites of a
double well are associated with the operators acting on the left/right end of
the link.

We define a block as the ensemble of one $g$ site and its closest $e$ sites
(i.e. the closest side of a double well on each adjacent link, see Fig.
1-(a)). With this definition, gauge invariance translates to the conservation
of the number of atoms per block (the U(1) local charge) and total spin per
block (the SU($N$) local charge). As we will detail below, the strategy to
build the simulator consists in making all gauge-invariant particle
configurations in each block resonant with each others, while non-gauge
invariant configurations are off-resonant (or decoupled). Gauge-invariant
dynamics is then realized perturbatively in this manifold, in a similar spirit
as in other proposals based on linear or quadratic gauge protection [4, 45,
46, 47].

###  2.1 Gauss’ laws and QLM Hamiltonian

As for any gauge theory, the key element of QLMs is the presence of local
symmetries, that constrain the physical Hilbert space that shall be spanned by
the engineered dynamics. Within the QLM formulation, the generators of the
local symmetry can be written in terms of the left/right operators

$L^{a}_{x,x+\hat{k}}=c_{x,+\hat{k}}^{i\dagger}\lambda_{ij}^{a}c_{x,+\hat{k}}^{j}$
(1)
$R^{a}_{x,x+\hat{k}}=c_{x+\hat{k},-\hat{k}}^{i\dagger}\lambda_{ij}^{a}c_{x+\hat{k},-\hat{k}}^{j}.$
(2)

where $c_{x,+\hat{k}}^{i\dagger}$ is the rishon operator with color $i$ acting
on the left edge of the link $x,x+\hat{k}$ (see Fig. 1) and where we assume
summation on the color indices. The matrices $\lambda_{ij}^{a}$,
$a=1,\dots,N$, are the Pauli matrices for $N=2$ or their generalization for
SU($N$) such that $\Tr\lambda^{a}\lambda^{b}=2\delta^{ab}$.

We can also define the U$(1)$ electric field as

$E_{x,y}=\frac{1}{2}\left(c_{y,-\hat{k}}^{i\dagger}c_{y,-\hat{k}}^{i}-c_{x,+\hat{k}}^{i\dagger}c_{x,+\hat{k}}^{i}\right).$
(3)

With these definitions, the generators of the SU$(N)$ and U$(1)$ symmetries
have the form

$\displaystyle G^{a}_{x}$
$\displaystyle=\psi_{x}^{i\dagger}\lambda_{ij}^{a}\psi_{x}^{j}+\sum_{\hat{k}}(L^{a}_{x,x+\hat{k}}+R^{a}_{x-\hat{k},x})$
$\displaystyle=\psi_{x}^{i\dagger}\lambda_{ij}^{a}\psi_{x}^{j}+\sum_{\hat{k}}\sum_{r=\pm\hat{k}}c_{x,r}^{i\dagger}\lambda_{ij}^{a}c_{x,r}^{j}$
(4a) $\displaystyle G_{x}$
$\displaystyle=\psi_{x}^{i\dagger}\psi_{x}^{i}-\sum_{\hat{k}}(E_{x,x+\hat{k}}-E_{x-\hat{k},x})-q_{0}$
$\displaystyle=\psi_{x}^{i\dagger}\psi_{x}^{i}+\sum_{\hat{k}}\sum_{r=\pm\hat{k}}\left(c_{x,r}^{i\dagger}c_{x,r}^{i}-\frac{\mathcal{N}}{2}\right)-q_{0}.$
(4b)

The quantity $q_{0}$ is an additive constant that depends on the choice of the
vacuum sector. We used $\psi_{x}^{i\dagger}$ (resp. $\psi_{x}^{i}$) as the
creation (resp. annihilation) operator of fermionic matter with color index
$i$ on site $x$.

In our implementation, they correspond to creation/annihilation operators of a
$g$ particle of SU($N$) spin $i=1,...,N$ in the $g$ well at the center of
block $x$. The operators $c_{x,r}^{i\dagger}$, with $r=\pm\hat{k}$, are the
creation operator of an $e$ particle on the $e$ well of block $x$ located on
the link with block $x\pm\hat{k}$ (cf Fig. 1). All states satisfying
$G_{x}=0$, and $G_{x}^{a}=0$ for all $a$ (called Gauss’ laws) build up the
gauge invariant subspace.

The target Hamiltonian of the system in the resonant subspace is the following
SU($N$)$\times$U($1$) QLM Hamiltonian for massive staggered fermions (see [1]
for a review)

$\displaystyle H_{\text{QLM}}=$ $\displaystyle-\tau\sum_{\langle
x,y\rangle}\left(\psi^{i\dagger}_{x}U_{x,y}^{ij}\psi^{j}_{y}+\text{H.c.}\right)+m\sum_{x}s_{x}\psi^{i\dagger}_{x}\psi^{i}_{x}+\frac{g^{2}}{2}\sum_{\langle
x,y\rangle}(L_{x,y}^{2}+R_{x,y}^{2})+\frac{{g^{\prime}}^{2}}{2}\sum_{\langle
x,y\rangle}E_{x,y}^{2}$
$\displaystyle\qquad-\frac{1}{4g^{2}}\sum_{\Box}\Tr(U_{\Box}+U_{\Box}^{\dagger})$
(5)

with $s_{x}=+1$ $(-1)$ for even (odd) sites, and summation over repeated
indices. The matter/anti-matter particles have mass $m$; The gauge assisted-
hopping (matter-gauge coupling), of amplitude $\tau$, is mediated by the
SU($N$) parallel transporter
$U_{x,x+\hat{k}}^{ij}=c_{x,\hat{k}}^{i}c_{x+\hat{k},-\hat{k}}^{\dagger j}$,
according to the rishon formulation of QLM [1]. The terms proportional to
$g^{2}$ and $g^{\prime 2}$ are the electric terms for the non-Abelian and the
Abelian gauge field respectively. We used the shorthand notation where
$\langle x,y\rangle$ corresponds to the sum over the pairs of neighboring
sites $x,y=x+\hat{k}$. The magnetic term acts on the plaquettes of the
lattice: For each plaquette spanning sites $w,x,y,z$, we defined
$U_{\Box}=U_{w,x}^{ij}U_{x,y}^{jk}U_{y,z}^{kl}U_{z,w}^{li}$, where the SU$(N)$
indices $i,j,k,l$ are summed over. We work in the representation with
$\mathcal{N}\equiv
c_{x,\hat{k}}^{i\dagger}c_{x,\hat{k}}^{i}+c_{x+\hat{k},-\hat{k}}^{\dagger
i}c_{x+\hat{k},-\hat{k}}^{i}=1$ rishon per link: the corresponding
experimental setup has exactly one $e$ particle per double well, avoiding
2-body $e-e$ losses [28, 48]. In this representation, the electric terms are
constant and have no effect on the dynamics. We remark, however, that in our
approach we will focus on implementing Gauss’ law, and we will obtain an
effective gauge-invariant Hamiltonian which can, in general, contain
additional terms. This approach has been pursued often in quantum simulation
of gauge theories [23], and emphasizes the importance of gauge invariance over
the exact system dynamics: it is motivated by the fact that, on the long-term
goal of emulating continuum field theories, the specifics of interactions on
the lattice shall play a minor role (in a similar manner as improved actions
return the same continuum limit), once all local and global symmetries are
faithfully realized. While we do not have continuum limit ambitions here, we
still prefer to emphasize the role of symmetry over microscopic processes.

##  3 Microscopic model

We now elaborate on the derivation of a SU$(N)\times$U$(1)$ lattice gauge
theory when starting from the lattice discretization of the atomic
Hamiltonian. The continuous optical lattice Hamiltonian for alkaline-earth
atoms is [49, 44]

$\begin{split}\mathcal{H}=&\sum_{\alpha
i}\int\mathrm{d}^{3}\mathbf{r}\Psi^{\dagger}_{\alpha
i}\left(\mathbf{r}\right)\left(-\frac{\hbar^{2}}{2M}\nabla^{2}+V_{\alpha}\left(\mathbf{r}\right)\right)\Psi_{\alpha
i}\left(\mathbf{r}\right)+\hbar\omega_{0}\int\mathrm{d}^{3}\mathbf{r}\left(\rho_{e}\left(\mathbf{r}\right)-\rho_{g}\left(\mathbf{r}\right)\right)\\\
&\quad+\frac{g_{eg}^{+}+g_{eg}^{-}}{2}\int\mathrm{d}^{3}\mathbf{r}\rho_{e}\left(\mathbf{r}\right)\rho_{g}\left(\mathbf{r}\right)+\sum_{\alpha,i<j}g_{\alpha\alpha}\int\mathrm{d}^{3}\mathbf{r}\rho_{\alpha
i}\left(\mathbf{r}\right)\rho_{\alpha j}\left(\mathbf{r}\right)\\\
&\quad+\frac{g_{eg}^{+}-g_{eg}^{-}}{2}\sum_{ij}\int\mathrm{d}^{3}\mathbf{r}\Psi_{gi}^{\dagger}\left(\mathbf{r}\right)\Psi_{ej}^{\dagger}\left(\mathbf{r}\right)\Psi_{gj}\left(\mathbf{r}\right)\Psi_{ei}\left(\mathbf{r}\right),\end{split}$
(6)

with $\Psi^{\dagger}_{\alpha i}$ the field creation operator at $\mathbf{r}$
in the state $\alpha=g,e$ with SU($N$) spin $i$ and
$\rho_{\alpha}\left(\mathbf{r}\right)=\sum_{i}\rho_{\alpha,i}\left(\mathbf{r}\right)=\sum_{i}\Psi^{\dagger}_{\alpha
i}\Psi_{\alpha i}$. $M$ is the mass of the atomic species, $V_{\alpha}$ is the
lattice potential, $\hbar\omega_{0}$ is the transition frequency between $g$
and $e$, $g_{\alpha\alpha}$, and $g_{eg}^{\pm}$ are the two-bodies contact
interaction scattering lengths. The lattice potential $V_{\alpha
i}(\mathbf{r})$ is realized using counter-propagating laser beams and is
designed such that the minima of the lattice potential correspond to the
lattice sites in Fig. 1.

The lattice version of the Hamiltonian is derived by decomposing $\Psi_{\alpha
i}$ into series of maximally localized Wannier functions [50, 51] from the
populated bands of the lattices only, separated by the gap $\Delta_{g}$ and
$\Delta_{e}$ from the unpopulated bands. Neglecting those higher bands, we are
left with one Wannier function per block for the $g$ atoms and one Wannier
function on each end of a link for the $e$ atoms (Figs. 5, 6). The
decomposition in Wannier functions of the lowest bands reads:

$\Psi_{gi}(\mathbf{r})=\sum_{x}w_{g}(\mathbf{r}-\mathbf{R}_{x})\psi_{x}^{i},$
(7)

$\Psi_{ei}(\mathbf{r})=\sum_{x}\sum_{s}w_{e,s}(\mathbf{r}-\mathbf{R}_{x,s})c_{x,s}^{i},$
(8)

where $\mathbf{R}_{x}$ is the center of the block $x$, and the sum over $s$
runs over all the $e$ sites in a block (for example, $s=\pm$ in 1d,
$s=\pm\hat{k}_{1},\pm\hat{k}_{2}$ in 2d). The lattice Hamiltonian is obtained
by substituting Eqs. (7) and (8) in the Hamiltonian $\mathcal{H}$ in Eq. (6).

##  4 Gauge-invariant effective Hamiltonian

We can write the lattice Hamiltonian as $H=H_{0}+H_{1}+\dots$ where
$H_{0}=\sum_{x}h_{x}^{0}$ contains all the terms acting on a single block,
$H_{1}=\sum_{\langle xy\rangle}h^{1}_{x,y}$ contains all the terms connecting
two neighbouring blocks, and the dots refer to all the other terms with longer
range, which (as we check for the chosen realistic set of parameters in
Appendix A) can be neglected. In our setup, the dominant energy scale is set
by the single-block Hamiltonian $H_{0}$, and $H_{1}$ can be regarded as a
perturbation. This property, together with the intrinsic global
SU$(N)\times$U$(1)$ symmetry of the alkaline-earth atoms, is the crucial
ingredient for gauge-invariance: the global symmetries imply
$[H_{0},\sum_{x}G_{x}]=0$ and $[H_{0},\sum_{x}G_{x}^{a}]=0$ for every $a$;
since $h_{y}^{0}$ and the local generators $G_{x}$ ($G_{x}^{a}$) trivially
commute for $x\neq y$, we get $[h_{x}^{0},G_{x}]=0$ and
$[h_{x}^{0},G_{x}^{a}]=0$ for every $a$ and for every block $x$. In other
words, $H_{0}$ is manifestly gauge-invariant, but is non-interacting, because
it does not couple different blocks. The perturbation $H_{1}$, on the other
hand, is not gauge invariant (it does not commute with the local generators,
only with the global ones), but couples different blocks. To obtain gauge-
invariant dynamics in perturbation theory, we require that all the gauge-
invariant eigenstates of $H_{0}$ are degenerate (resonant) with respect to
$H_{0}$: the perturbation $H_{1}$ can then induce non-trivial dynamics in the
resonant gauge-invariant sector. Conversely, we want all gauge-variant states
to be off-resonant. If some gauge-variant states are resonant, we can still
obtain approximate gauge-invariant dynamics as long as the perturbation does
not couple these states with the gauge-invariant sector.

###  4.1 Single block Hamiltonian

Figure 2: Example of correlated hopping obtained as a second-order process in
perturbation theory. The transparent intermediate states are outside of the
resonant subspace and are responsible of the second-order perturbation.
Antisymmetric nuclear spin singlet (S)
---
Graphical representation | Electronic state (sym.) | Energy
| $\ket{g,g}$ | $E_{gg}=2\mu_{g}+U_{gg}$
| $\frac{1}{\sqrt{2}}\left(\ket{g,e_{r_{1}}}+\ket{e_{r_{1}},g}\right)$ | $E_{eg}^{S}=\mu_{g}+\mu_{e}+U_{eg}^{+}$
| $\frac{1}{\sqrt{2}}\left(\ket{e_{r_{1}},e_{r_{2}}}+\ket{e_{r_{2}},e_{r_{1}}}\right)$ | $E_{ee}^{S}=2\mu_{e}$
Symmetric nuclear spin multiplet (M)
Graphical representation | Electronic state (antisym.) | Energy
| $\frac{1}{\sqrt{2}}\left(\ket{g,e_{r_{1}}}-\ket{e_{r_{1}},g}\right)$ | $E_{eg}^{M}=\mu_{g}+\mu_{e}+U_{eg}^{-}$
| $\frac{1}{\sqrt{2}}\left(\ket{e_{r_{1}},e_{r_{2}}}-\ket{e_{r_{2}},e_{r_{1}}}\right)$ | $E_{ee}^{M}=2\mu_{e}$
Table 1: The two particles eigenstates and eigenenergies of a single block $x$
for $N=2$. “$e_{r_{1}}$” designate a state with one $e$ particle in one of the
rishons in block $x$, and $r_{1},r_{2}$ are two different rishons. Only the
singlet states satisfy gauge invariance ($G^{a}_{x}=0$ for all $a$).

The lattice formulation of Eq. 6 involving only block $x$ is

$\begin{split}h_{x}^{0}&=\mu_{g}n_{x}^{g}+\frac{U_{gg}}{2}n_{x}^{g}(n_{x}^{g}-1)+\sum_{r}\left[\mu_{e}n_{x,r}^{e}+\frac{U_{eg}^{+}+U_{eg}^{-}}{2}n_{x}^{g}n_{x,r}^{e}+\frac{U_{eg}^{+}-U_{eg}^{-}}{2}\psi_{x}^{i\dagger}c_{x,r}^{j\dagger}\psi_{x}^{j}c_{x,r}^{i}\right]\\\
&\qquad+\sum_{r\neq
s}t_{rs}(c_{x,r}^{i\dagger}c_{x,s}^{i}+\text{H.c.})+\dots,\end{split}$ (9)

with $n_{x}^{g}=\psi_{x}^{i\dagger}\psi_{x}^{i}$ and
$n_{x,r}^{e}=c_{x,r}^{i\dagger}c_{x,r}^{i}$. The sums over $r$ and $s$ run
over all the rishons within the block (i.e., the closer of the two rishons for
each link connected to $x$), and the repeated spin indices $i,j$ are summed
over. We assume that the $e$ double wells are separated by a large barrier,
such that all hoppings $t_{rs}$ and the omitted density-assisted hoppings that
change the number $\mathcal{N}$ of rishons per link are suppressed (see
Appendix A for an estimate of these values for a set of realistic experimental
parameters). The Hamiltonian $h_{x}^{0}$ has a U$(N)$ gauge symmetry, as shown
by the commutation with the local generators
$[h_{x}^{0},G_{x}^{a}]=[h_{x}^{0},G_{x}]=0$.

###  4.2 Resonance conditions

We say that a process is resonant when it connects degenerate eigenstates of
$H_{0}=\sum_{x}h^{0}_{x}$. To induce the desired gauge-invariant dynamics, we
require the gauge-assisted hopping to be resonant between
compatible111“Compatible” means that the configurations considered must all
verify Gauss’ laws simultaneously. This includes having exactly one $e$
particle on each link at all times. gauge-invariant states.

For example for $N=2$ in one spatial dimension, the correlated hopping in Fig.
2 is resonant if

$E_{ee}^{S}+E_{gg}=2E_{eg}^{S},$ (10)

with the notations of Table 1. This condition imposes $U_{gg}=2U_{eg}^{+}$.

We now demonstrate the resonance conditions for all $N$ and all dimensions. We
want to prove that all states with $G_{x}=0$, $G_{x}^{a}=0$ are degenerate
eigenstates of $H_{0}=\sum_{x}h_{x}^{0}$, provided that $U_{gg}=2U_{eg}^{+}$.
As shown in Appendix C, we can write the Hamiltonian $H_{0}$ as

$\displaystyle H_{0}=$
$\displaystyle\sum_{x}\left[U_{eg}^{+}G_{x}n_{x}-\frac{U_{eg}^{+}-U_{eg}^{-}}{4}M^{a}_{x}G^{a}_{x}\right]$
$\displaystyle+\frac{U_{gg}-2U_{eg}^{+}}{2}\sum_{x}n_{x}(n_{x}-1)+\text{const.},$
(11)

with $M^{a}_{x}=\psi^{i\dagger}_{x}\lambda_{ij}\psi^{j}_{x}$. From this
expression we see that, when $U_{gg}=2U_{eg}^{+}$, $H_{0}$ is proportional to
the identity on the gauge invariant subspace, and all gauge-invariant states
are degenerate.

Together with the resonance condition, we require that gauge-variant states
are made off-resonant. We note that states with $G_{x}=0$, $M_{x}^{a}=0$ and
$G_{x}^{a}\neq 0$ are resonant but non-gauge-invariant. An example is the
gauge-breaking $e-e$ multiplet state in Table 1, which has the same energy as
the $e-e$ singlet state. In general, when $g$-wells are allowed to be empty,
such as in 1D for $N=2$, the gauge-breaking $e-e$ multiplet states have the
same energy $E_{ee}^{M}$ as the $e-e$ gauge-invariant singlet state
($E_{ee}^{S}$) such that the two belong to the same degenerate subspace.
However, the hybridization between the two types of states is suppressed as no
on-block or nearest-block term resonantly connect the two ensembles222Because
of the global SU$(N)$ symmetry, it is not possible to break Gauss’ law only on
a single block. If we imagine breaking Gauss’ law on two neighboring blocks
through a resonant second-order process, the final states of the two blocks
have to be the $e-e$ multiplets. This final state has two rishon on the link
between the two blocks: since the hopping between different links is
negligible, these states are not connected to gauge-invariant states, which
have one rishon per link. . For $N>2$ in 1D or $N>3$ in 2D hexagonal, the
available multiplet states always contain at least one $g$ particle, such that
they are separated from the singlet states by an energy difference
$\propto(U_{eg}^{+}-U_{eg}^{-})/4$. A crucial requirement for this mechanism
to work is to have large energy separation
$U_{eg}^{+}-U_{eg}^{-}=U_{eg}^{+}\cdot(1-g_{eg}^{-}/g_{eg}^{+})$: this
requirement is easily realized in fermionic alkaline-earth-like species, since
the ratio of scattering lengths $g_{eg}^{-}/g_{eg}^{+}$ is very different from
$1$. For example, $g_{eg}^{-}/g_{eg}^{+}=0.1170,1.616,0.428$ for
${}^{173}\mathrm{Yb}$, ${}^{171}\mathrm{Yb}$, and ${}^{87}\mathrm{Sr}$,
respectively [28, 30, 52, 53, 35].

###  4.3 Effective Hamiltonian

The effective Hamiltonian can be obtained using standard perturbative methods
(such as, e.g., the Schrieffer-Wolff method). Before showing an example for
the case $N=2$ in 1 spatial dimension, we are now going to discuss the general
features of the effective Hamiltonian.

We assume that the dominant terms generating the dynamics in perturbation
theory are single particle hoppings of $g$ particle between neighbouring
blocks or $e$ particles between the two ends of a link:

$\displaystyle H_{1}^{\text{hop}}=$ $\displaystyle-t_{g}\sum_{\langle
x,y\rangle}(\psi^{i\dagger}_{x}\psi^{i}_{y}+\text{H.c.})$ $\displaystyle-
t_{e}\sum_{x,\hat{k}}(c^{i\dagger}_{x,+\hat{k}}c^{i}_{x+\hat{k},-\hat{k}}+\text{H.c.}).$
(12)

The induced dynamics to second order in perturbation theory contains terms
that act on pairs of neighboring blocks while conserving the number of rishons
per link ($\mathcal{N}=1$) and preserving the U$(1)\times$SU$(N)$ gauge
symmetry. While the effective Hamiltonian will only contain gauge-invariant
terms, it will generally not coincide with the QLM Hamiltonian in Eq. (2.1).

Matter-gauge coupling, for example, can be obtained at second order in
perturbation theory, through a correlated hopping of the $g$ and the $e$
particle (Fig. 2). The amplitude for this process, however, depends on the
configuration of the other particles in the blocks involved (see Appendix D):
the corresponding Hamiltonian term is not simply given by
$(\psi_{x}^{i\dagger}U_{x,y}^{ij}\psi_{y}^{j}+\text{H.c.})$, but will contain
some dependence on the particle densities on neighboring $e$ sites.

The mass term, which breaks space translation symmetry, can be realized by a
modulation of the optical lattice potential with an additional mode, with a
spatial periodicity of two blocks. The $g$ and $e$ amplitudes for this
modulation are chosen to be much smaller than the ones of the main lattice
potential shown in Fig. 1-(b), such that the Wannier functions are almost
unaffected. This mode contributes to the effective lattice Hamiltonian as a
site-dependent potential for the $g$ and $e$ atoms and can be chosen to be of
the same order of the other terms in the effective second-order Hamiltonian.
Since this mode allows to tune the mass term independently from the other
parameters, we will focus, from now on, on the case $m=0$. The case $m\neq 0$
can then be straightforwardly derived, as discussed in Appendix A.

Together with the matter-gauge coupling, diagonal terms (in the electric
basis) are generated through second order processes. Some of these terms can
be easily compensated at this perturbative order: for example an on-site
interaction of the $g$ particles, of the form $n_{x}^{2}$ can be removed with
a small detuning from the resonance condition $U_{gg}=2U_{eg}^{+}$. In order
to observe some non-trivial dynamics, it is important that the diagonal terms
are not too large compared to the off-diagonal correlated hopping.

Plaquette terms are the other off-diagonal terms in $H_{\text{QLM}}$ [Eq.
(2.1)], but they only arise through forth-order processes in perturbations
theory (on a square lattice), and we therefore expect them to be very small.
Electric terms, on the other hand, are diagonal but they are trivially
constant for this representation of the gauge field with one rishon per link.

As an example, we compute the effective Hamiltonian to second order in
perturbation theory for the case $N=2$ in 1D:

$\displaystyle H_{\text{eff}}=$
$\displaystyle-\sum_{x}[\tau_{1}+\tau_{2}E_{x-1,x}E_{x+1,x+2}+\tau_{3}(E_{x-1,x}-E_{x+1,x+2})](\psi_{x}^{i\dagger}U_{x,x+1}^{ij}\psi_{x+1}^{j}+\text{H.c.})$
$\displaystyle+\frac{u}{2}\sum_{x}\psi_{x}^{j\dagger}\psi_{x}^{j}(\psi_{x}^{i\dagger}\psi_{x}^{i}-1)+w\sum_{x}\psi_{x}^{j\dagger}\psi_{x}^{j}\psi_{x+1}^{i\dagger}\psi_{x+1}^{i},$
(13)

where

$\displaystyle\tau_{1}$
$\displaystyle=\frac{t_{g}t_{e}}{6(U^{+}_{eg}-U^{-}_{eg})}-\frac{t_{g}t_{e}}{2(U^{+}_{eg}+3U^{-}_{eg})}+\frac{t_{g}t_{e}}{U_{eg}^{+}-3U_{eg}^{-}},$
(14) $\displaystyle\tau_{2}$
$\displaystyle=\frac{2t_{g}t_{e}}{3(U^{+}_{eg}-U^{-}_{eg})}-\frac{2t_{g}t_{e}}{U^{+}_{eg}+3U^{-}_{eg}}-\frac{4t_{g}t_{e}}{U_{eg}^{+}-3U_{eg}^{-}},\hskip
25.6073pt\tau_{3}=0,$ (15) $\displaystyle u$
$\displaystyle=\frac{8t_{g}^{2}+5t_{e}^{2}}{(U_{eg}^{+}-3U_{eg}^{-})}-\frac{3t_{e}^{2}}{(U_{eg}^{+}+U_{eg}^{-})}+\frac{8t_{g}^{2}}{U_{eg}^{+}+3U_{eg}^{-}}-\frac{20t_{e}^{2}}{3(U_{eg}^{+}-U_{eg}^{-})},$
(16) $\displaystyle w$
$\displaystyle=\frac{8t_{g}^{2}+5t_{e}^{2}}{2(U_{eg}^{+}-3U_{eg}^{-})}-\frac{3t_{e}^{2}}{2(U_{eg}^{+}+U_{eg}^{-})}+\frac{2(t_{g}^{2}+t_{e}^{2})}{U_{eg}^{+}+3U_{eg}^{-}}-\frac{4t_{g}^{2}+5t_{e}^{2}}{3(U_{eg}^{+}-U_{eg}^{-})}.$
(17)

While this effective Hamiltonian is obtained assuming that the hoppings
$t_{g}$ and $t_{e}$ are the dominant terms in the perturbation, interactions
between neighboring blocks can in general have a comparable amplitude: in
fact, since we work in a regime with $U_{gg},U_{eg}^{+}\gg t_{g},t_{e}$, we
need a strong transverse confinement, which will enhance all the interaction
terms. A perturbative calculation that includes some of these terms is
presented in Appendix D.

In general, the effective Hamiltonian obtained from perturbation theory is a
useful guidance to understand the possible phases that can arise for various
parameters regimes, and to easily show the interpretation of the atomic model
as a non-Abelian lattice gauge theory. It stands as an open question whether
the range of parameters that can be probed in the experimental setup can
realize the full phase diagram of the QLM in Eq. (2.1), and whether other
phases can arise from the additional terms. Nevertheless, we expect that many
of the physical properties of non-Abelian QLMs will be present in the
experimental model we propose.

##  5 Experimental setup

We here discuss more in detail the experimental setup. In particular, we focus
on the implementation of the optical lattice and discuss how the amplitudes of
the lasers can be tuned to realize the desired parameter regime for the
effective model.

###  5.1 Optical lattice

The 1D optical lattice potential displayed in Fig. 1-(b) can be realized by
interfering laser beams at two different wavelengths $\lambda_{0}$,
$\lambda_{1}$, and at angles $\theta_{0},\theta_{1}$, creating sinusoidal
potentials along the $x$-direction with amplitude $V^{0}_{\alpha}$ and
$V^{1}_{\alpha}$, and periodicity $a_{0}\equiv\lambda_{0}/2\sin(\theta_{0}/2)$
and $a_{1}\equiv\lambda_{1}/2\sin(\theta_{1}/2)=a_{0}/2$, respectively. In
addition, a tight transverse confinement can be obtained by creating an
optical lattice with spacing $d$ and potential strengths $V_{\alpha}^{\perp}$
along the transverse $y$ and $z$ directions. Finally, a small potential
$V_{\alpha}^{S}$ at angle $\theta_{S}$ with periodicity
$a_{S}\equiv\lambda_{S}/2\sin(\theta_{S}/2)=2a_{0}$ is needed to produce the
mass term. These lattices add together to give the optical potential
$V_{\alpha}$ in Eq. (6):

$\begin{split}V_{\alpha}(\textbf{r})=\>&V_{\alpha}^{0}\sin^{2}\left(\frac{\pi
x}{a_{0}}\right)+V^{1}_{\alpha}\sin^{2}\left(\frac{2\pi
x}{a_{0}}+\,\varphi\right)\\\ &+V^{S}_{\alpha}\sin^{2}\left(\frac{\pi
x}{2a_{0}}\right)\\\ &+V_{\alpha}^{\perp}\left[\sin^{2}\left(\frac{\pi
y}{d}\right)+\sin^{2}\left(\frac{\pi z}{d}\right)\right].\end{split}$ (18)

The wavelengths $\lambda_{0},\lambda_{1},\lambda_{S}$ of the lasers are chosen
such that, for each component, the ratio between the desired amplitudes of the
$g$ and $e$ potentials coincides with the ratio of the dynamical
polarizabilities. We set $\theta_{1}=\pi$, such that the lattice spacing is
$a_{0}=\lambda_{1}$. The phase $\varphi$ can be conveniently modulated to
prepare some initial states of interest (see Appendix E), but will be
otherwise set to $\varphi=0$.

The amplitudes
$V_{\alpha}^{0},V_{\alpha}^{1},V_{\alpha}^{S},V_{\alpha}^{\perp}$ and the
width $d$ of the transverse confinement are the tunable parameters in the
experimental setup (see Appendix B). In particular, $V_{\alpha}^{S}$
determines the mass $m$ of the (staggered) fermions in the lattice gauge
theory. Note that the single-particle terms (such as the hoppings
$t_{g},t_{e}$) do not depend on the transverse confinement, while the
interaction terms are enhanced when $V_{\alpha}^{\perp}$ is large: by tuning
the amplitude and the width of the transverse potential with respect to the
longitudinal terms $V_{\alpha}^{0},V_{\alpha}^{1}$ we can achieve the desired
regime $U_{gg},U_{eg}^{+}\gg t_{g},t_{e}$ (see Appendix A). Moreover, the
relative ratio between $V_{g}^{\perp}$ and $V_{e}^{\perp}$ can be adjusted to
satisfy the resonance condition.

Similar arguments can be applied in a bidimensional case. In Fig. 3-(a,b) we
show the profile of an optical lattice that can be used to simulate a QLM in a
bidimensional square lattice:

$\displaystyle V_{\alpha}(\mathbf{r})$
$\displaystyle=V_{\alpha}^{1}\left[\sin^{2}\left(\frac{\pi(x+y)}{a_{0}}\right)+\sin^{2}\left(\frac{\pi(x-y)}{a_{0}}\right)\right]$
$\displaystyle+V_{\alpha}^{0}\left[\sin^{2}\left(\frac{3\pi
x}{a_{0}}\right)+\sin^{2}\left(\frac{3\pi y}{a_{0}}\right)\right]$
$\displaystyle+V_{\alpha}^{S}\left[\sin^{2}\left(\frac{\pi
x}{2a_{0}}\right)+\sin^{2}\left(\frac{\pi y}{2a_{0}}\right)\right]$
$\displaystyle+V_{\alpha}^{\perp}\sin^{2}\left(\frac{\pi z}{d}\right).$ (19)

Figure 3: Optical lattice for the two-dimensional QLM: the $g$ potential (a)
has minima (dark blue) in the vertices of the lattice, the $e$ potential (b)
has two minima (dark orange) on each link. (c) The brickwall lattice is
realized with the same potential of the square lattice ($g$ sites (blue)
correspond to matter, the $e$ sites (orange) to rishons): the shaded-out
orange boxes represent double wells that are initialized with no $e$ atoms,
such that no gauge variable is present on the corresponding link. The
staggering lattice shifts even and odd blocks (in green/purple). (d) The
brickwall lattice is equivalent to the honeycomb lattice.

The superlattice potential $V_{\alpha}^{S}$ ensures the staggering of the mass
term [Fig. 3-(c)]. Similarly to the $1D$ case, the amplitude of the transverse
potential $V_{\alpha}^{\perp}$ is tuned such that the resonance condition is
satisfied.

The same optical lattice can be used to realize a brickwall lattice
[equivalent to a honeycomb lattice, Fig. 3-(d)] QLM: it is sufficient to
initialize the system in a state with “empty” links as in Fig. 3-(c).

###  5.2 Ab initio derivation of the model

We now show an example of realistic parameters, obtained through an ab initio
calculation: the single particle Hamiltonians for the atoms in the $g/e$
optical lattices are solved using Bloch theorem; maximally localized Wannier
functions are computed from the lowest band for the $g$ lattice (which is
separated from the second band by an energy gap $\Delta_{g}$) and from the two
lowest bands of the $e$ lattice (whose gap with the third band is
$\Delta_{e}$). From the Wannier functions, we can compute the parameters in
Eqs. (9) and (4.3) using Eqs. (6), (7) and (8). An example of realistic values
of the experimental parameters chosen for $N=2$ in 173Yb, and the
corresponding computed amplitudes are summarized in Tab. 2 and considered
below. The large $g_{eg}^{+}/g_{eg}^{-}$ ratio makes 173Yb a particularly
convenient atomic species for our proposal. Note that another suitable choice
is 87Sr, while 171Yb cannot be used because the negative sign of the ratio
$g_{gg}/g_{eg}^{+}<0$ implies that the resonance condition cannot be
satisfied.

$M$ |  $172.93$ u
---|---
$g_{gg}$ | $7.748\,\mathrm{Hz}\,\mu\mathrm{m}^{3}$
$g_{eg}^{+}$ | $72.97\,\mathrm{Hz}\,\mu\mathrm{m}^{3}$
$g_{eg}^{-}$ | $8.536\,\mathrm{Hz}\,\mu\mathrm{m}^{3}$

$a_{0}$ | $\sim 0.45\,\mathrm{\mu m}$
---|---
$V^{0}_{g}$ | $13.531\,h\cdot\mathrm{kHz}$
$V^{1}_{g}$ | $-4.273\,h\cdot\mathrm{kHz}$
$V^{0}_{e}$ | $-131.40\,h\cdot\mathrm{kHz}$
$V^{1}_{e}$ | $-149.56\,h\cdot\mathrm{kHz}$
$d_{\perp}$ | $\sim 0.80\,\mathrm{\mu m}$
$V_{g}^{\perp}$ | $2168.87\,h\cdot\mathrm{kHz}$
$V_{e}^{\perp}$ | $23.29\,h\cdot\mathrm{kHz}$

Table 2: First column: parameters of the 173Yb in Eq. (6) [28, 29, 30]. Second
column: realistic parameters for the optical lattices in Eq. (18). Third
column: corresponding values of the lattice Hamiltonian, Eqs. (9) and (4.3).
Fourth column: parameters of the effective Hamiltonian in Eq. (4.3).

The parameters $V_{g}^{\perp}$ and $d$ of the transverse potential are chosen
to satisfy $U_{gg}\gg t_{g},t_{e}$, such that our perturbative treatment
successfully applies. We then tune $V_{e}^{\perp}$ to realize the resonance
condition $U_{gg}\simeq 2U_{eg}^{+}$ (the interactions $U_{eg}^{\pm}$ grow
with $V_{e}^{\perp}$, while $t_{g},t_{e},U_{gg}$ are unaffected). For the
values in Tab. 2, in particular, we chose $U_{eg}^{+}$ such that
$U_{gg}-2U_{eg}^{+}+u=0$, i.e., such that the effective on-site matter-matter
interaction vanishes up to second order in perturbation theory. We find that,
since the $e$ atoms are already strongly localized in the longitudinal
direction, to satisfy the resonance condition we need a stronger transverse
confinement on the $g$ atoms than on the $e$ atoms.

The desired values of the potentials $V_{\alpha}^{1}$, $V_{\alpha}^{0}$,
$V_{\alpha}^{S}$ and $V_{\alpha}^{\perp}$ can be obtained by an appropriate
choice of the wavelengths of the lasers, as discussed in Appendix B.

For the values considered in Tab. 2, we found that the $g$ hopping mediated by
$e$ particles with/without spin exchange has amplitude $(A_{eg}^{-}\pm
A_{eg}^{+})/2$ comparable to $t_{g}$ and $t_{e}$. We therefore performed a
calculation of the effective Hamiltonian $H_{\text{eff}}$ including these
terms in $H_{1}$. The expressions for $u,w,\tau_{1},\tau_{2},\tau_{3}$ in this
case are presented in Appendix D. Their values, computed for the parameters
considered here, are also reported in Tab. 2.

###  5.3 Preparation of the initial state and readout

In order to experimentally prepare initial states within the gauge invariant
subspace, where each block is occupied by a spin singlet state, one can
exploit the spectroscopic resolution of the clock transition to resolve two-
particle states with different electronic wavefunctions. In particular, $N=2$
gauge-invariant (spin singlet) states can be created starting from a band
insulator of $g$-state atoms in the lattice potential with periodicity
$a_{0}$, and transferring each on-site $\ket{g,g}$ (singlet) state into one of
the states given in the top part of Table 1. For homogeneous or two-block
translation invariant states, this can be achieved by only employing global
clock laser pulses after turning on the short lattice with periodicity
$a_{0}/2$ (see Appendix E for details). To prepare more complex gauge-
invariant states such as the one depicted in Fig. 4, local $g$-site addressing
is required in order to manipulate site-dependent clock light shifts [54, 55],
allowing to prepare an arbitrary singlet state starting from $\ket{g,g}$ on
any given block. Time-resolved clock spectroscopic could be exploited to
measure the real-time populations of the singlet and triplet sectors, allowing
for experimentally tracking the violation of SU($N$) gauge invariance.
Further, U$(1)$ gauge invariance would be validated spectroscopically within
the same measurement, since any single-block occupation number differing from
$N$ would produce a clearly resolvable spectroscopic signal [35]. The coherent
evolution of each rishon-site pair, and specifically the mutually exclusive
occupation of neighboring rishon sites, could be monitored using standard
band-mapping techniques exploited in double-well lattices [56]. Finally, site-
resolved in situ imaging of $e$-state atoms will grant access to the electric
field evolution.

##  6 Robustness of the gauge symmetries

In order to evaluate the violation of gauge invariance, we simulated the real-
time dynamics of the lattice Hamiltonian in one spatial dimension for $N=2$
using exact diagonalisation on a system of $4$ blocks with periodic boundary
conditions. This method allows us to check the validity of our proposal,
beyond perturbative arguments.

Figure 4: (a) Initial state. (b) Abelian Gauss’ law: Expectation value of the
generator $G_{x}$. (c) Non-Abelian Gauss’ law: the total gauge violation of a
block $x$ is normalized such that $\sum_{\alpha}(G_{x}^{\alpha})^{2}/8=1$ for
a spin multiplet state. In (b) and (c), the lines corresponding to $x=0$ and
$x=2$ perfectly overlap because of the symmetry of the initial state, so only
the latter is shown. (d) Evolution of the Abelian electric field $E_{x-1,x}$.
Dashed lines: effective Hamiltonian, second-order perturbation theory.

In Fig. 4-(b-d) we plot the results for the time evolution starting from the
initial state in Fig. 4-(a), assuming periodic boundary conditions. The
simulated Hamiltonian contains $H_{0}$, $H_{1}^{\text{hop}}$, as well as the
aforementioned density-mediated hoppings (see Appendix D). We find that both
the Abelian [Fig.4-(b)] and the Non-Abelian [Fig.4-(c)] gauge generators show
small violations of Gauss’ laws, while, on the same time scale, the electric
field [Fig.4-(d)] has a non-trivial evolution, performing large oscillations.
The evolution of the electric field is compared with the one generated by the
effective Hamiltonian $H_{\text{eff}}$: the deviations are small at short
times, and become more significant at longer times, suggesting that higher
order perturbative corrections come into play (see, e.g., Ref. [38] for a
similar effect). A full calculation of higher order terms is beyond the scope
of our work.

While the finite system size of our simulation may have an effect on the long-
time dynamics, we expect that the robustness will not depend significantly on
the system size [57]. We note that the Gauss’ law violations have a fast
growth at short times (shorter than the resolution of Fig. 4), but then
saturate or oscillate in a bounded interval. This saturation can be understood
from perturbation theory: the effective second-order Hamiltonian can be
obtained from the original lattice Hamiltonian through a Schrieffer-Wolff
transformation [58, 59, 60]; the same Schrieffer-Wolff transformation
“dresses” the generators of the gauge symmetry. The “bare” generators differ
from the dressed generators, which are quasi-conserved in the dynamics, by
terms of order $\lambda^{2}$, where $\lambda$ is the perturbative parameter,
leading to a saturation of the Gauss’ law violations to order $\lambda^{2}$.
The violations are then expected to slowly grow, departing from the saturation
value at longer times, when terms beyond the second-order approximation become
significant.

##  7 Discussion and conclusion

We proposed an ab initio scalable experimental setup of an SU($N$)$\times$U(1)
non-Abelian matter-gauge lattice gauge theory quantum simulator in 1D and 2D.
The experiment uses ultra-cold 173Yb in an optical superlattice. The system is
restricted to evolve within the resonant gauge-invariant subspace set by the
initial state, while gauge-breaking states are off-resonant. For the regime of
parameters we chose, the gauge-breaking dynamics is bounded and negligible.
Our scheme thus allows for an experimental verification of gauge-invariant
dynamics from first principles.

The setup we propose can be used to investigate the real-time evolution of
non-Abelian quantum link models with fermionic matter. Dynamical phenomena of
interest include, for example, chiral symmetry restoration, ergodicity
breaking, particle collisions, and the evolution of local defects.

The first principle character of our investigation reveals possibilities, but
also challenges along the realization of non-Abelian lattice gauge theory
dynamics, that are typically not easily captured in more conceptually-minded
proposals, and were thus not fully clear so far (at least, to us). In
particular, similar to Abelian theories [38], the spatial engineering of
lattice potentials introduces rather complicated Wannier functions, that
definitely go beyond what it is typically assumed in tight-binding
approximations. Moreover, the role of exchange interactions shall be
considered with great care: while those are in principle quite large energy
scales, the fact that they need to be employed off-site (combined with the
exotic structure of Wannier functions) makes them not as effective as expected
in terms of engineering gauge invariant dynamics (this aspect was implicitly
hinted to in Ref. [61]). This observation suggests that more detailed
microscopic computations of effective Hamiltonians in systems involving spin
exchange interactions are pivotal to quantitatively understand timescales in
those promising settings [61, 62, 63], as well as in those where only local
interactions have been used [23, 24].

At the conceptual level, an interesting perspective is to investigate
alternative setups, where mixtures of alkaline-earth and alkali atoms are
trapped simultaneously [64]. It might be possible that, thanks to the much
wider flexibility in trapping potentials offered by the distinct atomic
structures, such systems can help circumvent some of the challenges single-
species platforms face. Given recent experimental progresses in the field,
investigations along this line are now, in our opinion, a promising way
forward, that shall be approached with a similar first principles spirit as
highlighted here.

We acknowledge useful discussions with Monika Aidelsburger. F. M. S.
acknowledges support provided by the U.S. Department of Energy Office of
Science, Office of Advanced Scientific Computing Research, (DE-SC0020290), by
Amazon Web Services, AWS Quantum Program, and by the DOE QuantISED program
through the theory consortium “Intersections of QIS and Theoretical Particle
Physics” at Fermilab. The work of M. D. was partly supported by the MUR
Programme FARE (MEPH), by QUANTERA DYNAMITE PCI2022-132919, by the PNRR MUR
project PE0000023-NQSTI, and by the EU-Flagship programme Pasquans2. M.D. also
thanks the Yukawa Institute for Theoretical Physics for hospitality during the
QIMG2023 workshop (YITP-T-23-01). F. S. acknowledges funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (Grant agreement No. 949438), from the Italian MUR under
the FARE programme (project FastOrbit) and the PRIN programme (project CoQuS).

## References

  * [1] Uwe Jens Wiese.  “Ultracold quantum gases and lattice systems: Quantum simulation of lattice gauge theories”.  Annalen der Physik 525, 777–796 (2013).  arXiv:1305.1602.
  * [2] John Preskill.  “Simulating quantum field theory with a quantum computer” (2018).  url: arxiv.org/abs/1811.10085.
  * [3] Erez Zohar, J Ignacio Cirac, and Benni Reznik.  “Quantum simulations of lattice gauge theories using ultracold atoms in optical lattices”.  Reports on Progress in Physics 79, 014401 (2015).
  * [4] M. Dalmonte and S. Montangero.  “Lattice gauge theories simulations in the quantum information era”.  Contemp. Phys. 57, 388 (2016).
  * [5] Mari Carmen Bañuls, Rainer Blatt, Jacopo Catani, Alessio Celi, Juan Ignacio Cirac, Marcello Dalmonte, Leonardo Fallani, Karl Jansen, Maciej Lewenstein, Simone Montangero, and et al.  “Simulating lattice gauge theories within quantum technologies”.  The European Physical Journal D74 (2020).
  * [6] Christian W. Bauer, Zohreh Davoudi, A. Baha Balantekin, Tanmoy Bhattacharya, Marcela Carena, Wibe A. de Jong, Patrick Draper, Aida El-Khadra, Nate Gemelke, Masanori Hanada, Dmitri Kharzeev, Henry Lamm, Ying-Ying Li, Junyu Liu, Mikhail Lukin, Yannick Meurice, Christopher Monroe, Benjamin Nachman, Guido Pagano, John Preskill, Enrico Rinaldi, Alessandro Roggero, David I. Santiago, Martin J. Savage, Irfan Siddiqi, George Siopsis, David Van Zanten, Nathan Wiebe, Yukari Yamauchi, Kübra Yeter-Aydeniz, and Silvia Zorzetti.  “Quantum simulation for high-energy physics”.  PRX Quantum 4, 027001 (2023).
  * [7] I. Montvay and G. Muenster.  “Quantum Fields on a lattice”.  Cambridge Univ. Press, Cambridge. (1994).
  * [8] Kenji Fukushima and Tetsuo Hatsuda.  “The phase diagram of dense QCD”.  Rep. Progr. Phys. 74, 014001 (2010).
  * [9] R.A. Soltz, C. DeTar, F. Karsch, Swagato Mukherjee, and P. Vranas.  “Lattice QCD thermodynamics with physical quark masses”.  Annual Review of Nuclear and Particle Science 65, 379–402 (2015).
  * [10] Zoltan Fodor and Christian Hoelbling.  “Light hadron masses from lattice QCD”.  Rev. Mod. Phys. 84, 449–495 (2012).
  * [11] E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz, P. Zoller, and R. Blatt.  “Real-time dynamics of lattice gauge theories with a few-qubit quantum computer”.  Nature 534, 516 (2016).
  * [12] Kenneth G. Wilson.  “Confinement of quarks”.  Phys. Rev. D 10, 2445 (1974).
  * [13] Marco Anderlini, Patricia J. Lee, Benjamin L. Brown, Jennifer Sebby-Strabley, William D. Phillips, and J. V. Porto.  “Controlled exchange interaction between pairs of neutral atoms in an optical lattice”.  Nature 448, 452–456 (2007).
  * [14] S. Trotzky, P. Cheinet, S. Folling, M. Feld, U. Schnorrberger, A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. Lukin, and I. Bloch.  “Time-resolved observation and control of superexchange interactions with ultracold atoms in optical lattices”.  Science 319, 295–299 (2008).
  * [15] Christian Schweizer, Fabian Grusdt, Moritz Berngruber, Luca Barbiero, Eugene Demler, Nathan Goldman, Immanuel Bloch, and Monika Aidelsburger.  “Floquet approach to Z2 lattice gauge theories with ultracold atoms in optical lattices”.  Nature Physics 15, 1168–1173 (2019).
  * [16] Alexander Mil, Torsten V. Zache, Apoorva Hegde, Andy Xia, Rohit P. Bhatt, Markus K. Oberthaler, Philipp Hauke, Jürgen Berges, and Fred Jendrzejewski.  “A scalable realization of local U(1) gauge invariance in cold atomic mixtures”.  Science 367, 1128–1130 (2020).
  * [17] Bing Yang, Hui Sun, Robert Ott, Han-Yi Wang, Torsten V. Zache, Jad C. Halimeh, Zhen-Sheng Yuan, Philipp Hauke, and Jian-Wei Pan.  “Observation of gauge invariance in a 71-site Bose–Hubbard quantum simulator”.  Nature 587, 392–396 (2020).
  * [18] Hannes Bernien, Sylvain Schwartz, Alexander Keesling, Harry Levine, Ahmed Omran, Hannes Pichler, Soonwon Choi, Alexander S Zibrov, Manuel Endres, Markus Greiner, et al.  “Probing many-body dynamics on a 51-atom quantum simulator”.  Nature 551, 579 (2017).
  * [19] Zhao-Yu Zhou, Guo-Xian Su, Jad C Halimeh, Robert Ott, Hui Sun, Philipp Hauke, Bing Yang, Zhen-Sheng Yuan, Jürgen Berges, and Jian-Wei Pan.  “Thermalization dynamics of a gauge theory on a quantum simulator”.  Science 377, 311–314 (2022).
  * [20] Roland C Farrell, Marc Illa, Anthony N Ciavarella, and Martin J Savage.  “Scalable circuits for preparing ground states on digital quantum computers: The Schwinger model vacuum on 100 qubits” (2023).  arXiv:2308.04481.
  * [21] Natalie Klco, Martin J. Savage, and Jesse R. Stryker.  “SU(2) non-Abelian gauge field theory in one dimension on digital quantum computers”.  Phys. Rev. D 101, 074512 (2020).
  * [22] Yasar Atas, Jinglei Zhang, Randy Lewis, Amin Jahanpour, Jan F. Haase, and Christine A. Muschik.  “SU(2) hadrons on a quantum computer” (2021).  arXiv:2102.08920.
  * [23] D. Banerjee, M. Bögli, M. Dalmonte, E. Rico, P. Stebler, U.-J. Wiese, and P. Zoller.  “Atomic Quantum Simulation of U(N) and SU(N) Non-Abelian Lattice Gauge Theories”.  Physical Review Letters 110, 125303 (2013).
  * [24] Erez Zohar, J. Ignacio Cirac, and Benni Reznik.  “Cold-atom quantum simulator for SU(2) Yang-Mills lattice gauge theory”.  Phys. Rev. Lett. 110, 125304 (2013).
  * [25] L. Tagliacozzo, A. Celi, P. Orland, M. W. Mitchell, and M. Lewenstein.  “Simulation of non-Abelian gauge theories with optical lattices”.  Nature Communications 4 (2013).
  * [26] E. Rico, M. Dalmonte, P. Zoller, D. Banerjee, M. Bögli, P. Stebler, and U.-J. Wiese.  “SO(3) “nuclear physics” with ultracold gases”.  Annals of Physics 393, 466–483 (2018).
  * [27] Valentin Kasper, Torsten V. Zache, Fred Jendrzejewski, Maciej Lewenstein, and Erez Zohar.  “Non-Abelian gauge invariance from dynamical decoupling”.  Phys. Rev. D 107, 014506 (2023).
  * [28] F Scazza, C Hofrichter, M Höfer, P. C. De Groot, I Bloch, and S Fölling.  “Observation of two-orbital spin-exchange interactions with ultracold SU(N)-symmetric fermions”.  Nature Phys. 10, 779–784 (2014).
  * [29] G. Cappellini, M. Mancini, G. Pagano, P. Lombardi, L. Livi, M. Siciliani de Cumis, P. Cancio, M. Pizzocaro, D. Calonico, F. Levi, C. Sias, J. Catani, M. Inguscio, and L. Fallani.  “Direct observation of coherent interorbital spin-exchange dynamics”.  Phys. Rev. Lett. 113, 120402 (2014).
  * [30] M. Höfer, L. Riegger, F. Scazza, C. Hofrichter, D. R. Fernandes, M. M. Parish, J. Levinsen, I. Bloch, and S. Fölling.  “Observation of an orbital interaction-induced feshbach resonance in ${}^{173}\mathrm{Yb}$”.  Phys. Rev. Lett. 115, 265302 (2015).
  * [31] Christian Hofrichter, Luis Riegger, Francesco Scazza, Moritz Höfer, Diogo Rio Fernandes, Immanuel Bloch, and Simon Fölling.  “Direct probing of the Mott crossover in the $\mathrm{SU}(N)$ Fermi-Hubbard model”.  Phys. Rev. X 6, 021030 (2016).
  * [32] L. F. Livi, G. Cappellini, M. Diem, L. Franchi, C. Clivati, M. Frittelli, F. Levi, D. Calonico, J. Catani, M. Inguscio, and L. Fallani.  “Synthetic dimensions and spin-orbit coupling with an optical clock transition”.  Phys. Rev. Lett. 117, 220401 (2016).
  * [33] S. L. Campbell, R. B. Hutson, G. E. Marti, A. Goban, N. Darkwah Oppong, R. L. McNally, L. Sonderhouse, J. M. Robinson, W. Zhang, B. J. Bloom, and J. Ye.  “A Fermi-degenerate three-dimensional optical lattice clock”.  Science 358, 90–94 (2017).
  * [34] L. Riegger, N. Darkwah Oppong, M. Höfer, D. R. Fernandes, I. Bloch, and S. Fölling.  “Localized magnetic moments with tunable spin exchange in a gas of ultracold fermions”.  Phys. Rev. Lett. 120, 143601 (2018).
  * [35] A. Goban, R. B. Hutson, G. E. Marti, S. L. Campbell, M. A. Perlin, P. S. Julienne, J. P. D’Incao, A. M. Rey, and J. Ye.  “Emergence of multi-body interactions in a fermionic lattice clock”.  Nature 563, 369–373 (2018).
  * [36] Hideki Ozawa, Shintaro Taie, Yosuke Takasu, and Yoshiro Takahashi.  “Antiferromagnetic spin correlation of $\mathrm{SU}(\mathcal{N})$ Fermi gas in an optical superlattice”.  Phys. Rev. Lett. 121, 225303 (2018).
  * [37] Ren Zhang, Yanting Cheng, Peng Zhang, and Hui Zhai.  “Controlling the interaction of ultracold alkaline-earth atoms”.  Nature Reviews Physics 2, 213–220 (2020).
  * [38] Federica Maria Surace, Pierre Fromholz, Nelson Darkwah Oppong, Marcello Dalmonte, and Monika Aidelsburger.  “Ab initio derivation of lattice-gauge-theory dynamics for cold gases in optical lattices”.  PRX Quantum 4, 020330 (2023).
  * [39] D. Horn.  “Finite matrix models with continuous local gauge invariance”.  Physics Letters B 100, 149–151 (1981).
  * [40] Peter Orland and Daniel Rohrlich.  “Lattice gauge magnets: Local isospin from spin”.  Nuclear Physics B 338, 647–672 (1990).
  * [41] R. Brower, S. Chandrasekharan, and U.-J. Wiese.  “QCD as a quantum link model”.  Physical Review D 60, 094502 (1999).  arXiv:9704106.
  * [42] S Chandrasekharan and U.-J Wiese.  “Quantum link models: A discrete approach to gauge theories”.  Nuclear Physics B 492, 455–471 (1997).
  * [43] Uwe-Jens Wiese.  “From quantum link models to D-theory: a resource efficient framework for the quantum simulation and computation of gauge theories”.  Philosophical Transactions of the Royal Society A 380, 20210068 (2022).
  * [44] Alexey Vyacheslavovich Gorshkov, M Hermele, V Gurarie, C Xu, Paul S Julienne, J Ye, Peter Zoller, Eugene Demler, Mikhail D Lukin, and AM Rey.  “Two-orbital SU(N) magnetism with ultracold alkaline-earth atoms”.  Nature Physics 6, 289–295 (2010).
  * [45] Jad C Halimeh and Philipp Hauke.  “Stabilizing gauge theories in quantum simulators: a brief review” (2022).  arXiv:2204.13709.
  * [46] Jad C. Halimeh, Lukas Homeier, Christian Schweizer, Monika Aidelsburger, Philipp Hauke, and Fabian Grusdt.  “Stabilizing lattice gauge theories through simplified local pseudogenerators”.  Phys. Rev. Res. 4, 033120 (2022).
  * [47] Lukas Homeier, Annabelle Bohrdt, Simon Linsel, Eugene Demler, Jad C Halimeh, and Fabian Grusdt.  “Realistic scheme for quantum simulation of $Z_{2}$ lattice gauge theories with dynamical matter in (2+1) D”.  Communications Physics 6, 127 (2023).
  * [48] X. Zhang, M. Bishof, S. L. Bromley, C. V. Kraus, M. S. Safronova, P. Zoller, A. M. Rey, and J. Ye.  “Spectroscopic observation of SU(N)-symmetric interactions in Sr orbital magnetism”.  Science 345, 1467–1473 (2014).
  * [49] M A Cazalilla, A F Ho, and M Ueda.  “Ultracold gases of ytterbium: ferromagnetism and Mott states in an SU(6) Fermi system”.  New Journal of Physics 11, 103033 (2009).
  * [50] S. Kivelson.  “Wannier functions in one-dimensional disordered systems: Application to fractionally charged solitons”.  Phys. Rev. B 26, 4269–4277 (1982).
  * [51] Nicola Marzari and David Vanderbilt.  “Maximally localized generalized Wannier functions for composite energy bands”.  Physical Review B - Condensed Matter and Materials Physics 56 (1997).
  * [52] Koki Ono, Jun Kobayashi, Yoshiki Amano, Koji Sato, and Yoshiro Takahashi.  “Antiferromagnetic interorbital spin-exchange interaction of Yb171”.  Physical Review A 99, 032707 (2019).
  * [53] Oscar Bettermann, Nelson Darkwah Oppong, Giulio Pasqualetti, Luis Riegger, Immanuel Bloch, and Simon Fölling.  “Clock-line photoassociation of strongly bound dimers in a magic-wavelength lattice” (2020).  arXiv:2003.10599.
  * [54] Christof Weitenberg, Manuel Endres, Jacob F. Sherson, Marc Cheneau, Peter Schauß, Takeshi Fukuhara, Immanuel Bloch, and Stefan Kuhr.  “Single-spin addressing in an atomic mott insulator”.  Nature 471, 319–324 (2011).
  * [55] Christian Gross and Waseem S. Bakr.  “Quantum gas microscopy for single atom and spin detection”.  Nature Physics 17, 1316–1323 (2021).
  * [56] S. Trotzky, P. Cheinet, S. Fölling, M. Feld, U. Schnorrberger, A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. Lukin, and I. Bloch.  “Time-resolved observation and control of superexchange interactions with ultracold atoms in optical lattices”.  Science 319, 295–299 (2008).
  * [57] Jad C. Halimeh and Philipp Hauke.  “Reliability of lattice gauge theories”.  Phys. Rev. Lett. 125, 030503 (2020).
  * [58] Nilanjana Datta, Jürg Fröhlich, Luc Rey-Bellet, and Roberto Fernández.  “Low-temperature phase diagrams of quantum lattice systems. ii. convergent perturbation expansions and stability in systems with infinite degeneracy”.  Helvetica Physica Acta 69, 752–820 (1996).
  * [59] Sergey Bravyi, David P. DiVincenzo, and Daniel Loss.  “Schrieffer–wolff transformation for quantum many-body systems”.  Annals of Physics 326, 2793–2826 (2011).
  * [60] A. H. MacDonald, S. M. Girvin, and D. Yoshioka.  “$\frac{t}{U}$ expansion for the hubbard model”.  Phys. Rev. B 37, 9753–9756 (1988).
  * [61] K. Stannigel, P. Hauke, D. Marcos, M. Hafezi, S. Diehl, M. Dalmonte, and P. Zoller.  “Constrained dynamics via the Zeno effect in quantum simulation: Implementing non-Abelian lattice gauge theories with cold atoms”.  Phys. Rev. Lett. 112, 120406 (2014).
  * [62] Jad C Halimeh, Lukas Homeier, Annabelle Bohrdt, and Fabian Grusdt.  “Spin exchange-enabled quantum simulator for large-scale non-Abelian gauge theories” (2023).  arXiv:2305.06373.
  * [63] Pierpaolo Fontana, Joao C. Pinto Barros, and Andrea Trombettoni.  “Quantum simulator of link models using spinor dipolar ultracold atoms”.  Phys. Rev. A 107, 043312 (2023).
  * [64] Kevin Singh, Shraddha Anand, Andrew Pocklington, Jordan T. Kemp, and Hannes Bernien.  “Dual-element, two-dimensional atom array with continuous-mode operation”.  Phys. Rev. X 12, 011040 (2022).
  * [65] Thomas Uehlinger, Gregor Jotzu, Michael Messer, Daniel Greif, Walter Hofstetter, Ulf Bissbort, and Tilman Esslinger.  “Artificial graphene with tunable interactions”.  Phys. Rev. Lett. 111, 185307 (2013).

##  Appendix A Ab-initio derivation of lattice parameters

In this Section, we will show how to explicitly derive the parameters of the
lattice Hamiltonian from the amplitude of the optical potential. For a very
similar discussion, see also Ref. [38]. We focus here on the case of one
spatial dimension, but the procedure can also be applied to the case of two
spatial dimensions, with minor changes. As explained in the main text, we can
first focus on the case $m=0$ and compute the effective Hamiltonian for
$V^{S}_{\alpha}=0$. The effect of $V_{S}^{\alpha}\neq 0$ will then be
discussed below as a small perturbation.

As the first step of the procedure, we have to solve the single particle
Schroedinger equation in the periodic potential of the optical lattice. Since
the potential in Eq. (18) consists of a sum of $x$, $y$ and $z$ components,
the solution can be factorised in the three direction. We find each of these
factors separately. We first solve the Schroedinger equation for the
longitudinal potential $V_{\alpha}^{x}(x)$ for $\alpha=g,e$ [shown in grey in
Figs. 5-(a), 6-(a) respectively]:

$V_{\alpha}^{x}(\textbf{r})=V_{\alpha}^{0}\sin^{2}\left(\frac{\pi
x}{a_{0}}\right)+V^{1}_{\alpha}\sin^{2}\left(\frac{2\pi x}{a_{0}}\right).$
(20)

The Bloch bands obtained as solutions for the $g$, $e$ atoms are shown in
Figs. 5-(b) and 6-(b) respectively.

Figure 5: (a) Longitudinal lattice potential $V_{g}^{x}$ (in grey). Maximally
localized Wannier function from the first Bloch band (in blue). (b) Bloch
bands.  Figure 6: (a) Longitudinal lattice potential $V_{e}^{x}$ (in grey).
Maximally localized Wannier functions from a combination of the first two
Bloch bands (in orange). (b) Bloch bands. The inset shows the small gap
between the two lowest bands.

For the $g$ lattice, the lowest band is separated from the second-lowest band
by a large gap $\Delta_{g}$. We can focus on the lowest band to get a single
localized Wannier function $w_{g}(x)$ per unit cell [Fig. 5-(a)]. For the $e$
lattice we have to consider the two lowest bands, which are separated by a
small energy gap, in order to get the two maximally localized Wannier
functions $w^{x}_{e,-}(x)$ and $w^{x}_{e,+}(x)$ for the $e,-$ and $e,+$ sites
in a single unit cell [Fig. 6-(a)]. To find the maximally localized Wannier
functions we follow Refs. [50, 65] and look for the eigenstates of the
position operator projected on the bands of interest. The same procedure can
be applied to compute the Wannier functions in the transverse directions. In
this case, the lattice is a simple sinusoidal potential and it suffices to
consider the Wannier functions $w_{g}^{\perp}$ and $w_{e}^{\perp}$ from the
lowest $g$ and $e$ bands (for the values of our simulations, the potential is
very deep for both species). The profile of $w_{g}^{\perp}$ and
$w_{e}^{\perp}$ thus exclusively depends on the parameter $d$, $V_{g}^{\perp}$
and $V_{e}^{\perp}$

The lattice Hamiltonian can then be obtained from Eqs. (6), (7) and (8), with
the following definitions

$w_{g}(x,y,z)=w_{g}^{x}(x)w_{g}^{\perp}(y)w_{g}^{\perp}(z),\qquad
w_{e\pm}(x,y,z)=w_{e\pm}^{x}(x)w_{e}^{\perp}(y)w_{e}^{\perp}(z).$ (21)

We give here the definition of some of the most relevant parameters:

$t_{g}=-\int\mathrm{d}x\;\left(w^{x}_{g}(x)\right)^{*}\left(-\frac{\hbar^{2}}{2M}\partial_{x}^{2}+V_{g}^{x}(x)\right)w^{x}_{g}(x-a_{0}),$
(22)
$t_{e}=-\int\mathrm{d}x\;\left(w^{x}_{e,+}(x)\right)^{*}\left(-\frac{\hbar^{2}}{2M}\partial_{x}^{2}+V_{e}^{x}(x)\right)w^{x}_{e,-}(x-a_{0}),$
(23)
$U_{gg}=g_{gg}\left(\int\mathrm{d}z\;|w_{g}^{\perp}(z)|^{4}\right)^{2}\int\mathrm{d}x\;|w_{g}^{x}(x)|^{4},$
(24)

$U_{eg}^{\pm}=g_{eg}^{\pm}\left(\int\mathrm{d}z\;|w_{g}^{\perp}(z)|^{2}|w_{e,+}^{\perp}(z)|^{2}\right)^{2}\int\mathrm{d}x\;|w_{g}^{x}(x)|^{2}|w_{e,+}^{x}(x)|^{2},$
(25)

$A_{eg}^{\pm}=g_{eg}^{\pm}\left(\int\mathrm{d}z\;|w_{g}^{\perp}(z)|^{2}|w_{e,+}^{\perp}(z)|^{2}\right)^{2}\int\mathrm{d}x\;\left(w_{g}^{x}(x)\right)^{*}w_{g}^{x}(x-a_{0})|w_{e,+}^{x}(x)|^{2}.$
(26)

Examples of realistic values computed using this procedure are reported in
Tab. 2. Note however that $U_{gg}$ is of the same order of $\Delta_{g}$;
therefore, we expect our prediction to overestimate the interaction strength,
since the proximity to the higher bands effectively lowers the interactions
[28, 29, 34].

We also compute other terms, such as next-nearest neighbor hoppings, and other
density-mediated terms, to show that these are negligible:

$t_{g}^{(1)}=-\int\mathrm{d}x\;\left(w^{x}_{g}(x)\right)^{*}\left(-\frac{\hbar^{2}}{2M}\partial_{x}^{2}+V_{g}^{x}(x)\right)w^{x}_{g}(x-2a_{0}),$
(27)
$t_{e}^{(1)}=-\int\mathrm{d}x\;\left(w^{x}_{e,+}(x)\right)^{*}\left(-\frac{\hbar^{2}}{2M}\partial_{x}^{2}+V_{e}^{x}(x)\right)w^{x}_{e,-}(x),$
(28)

$B_{eg}^{\pm}=g_{eg}^{\pm}\left(\int\mathrm{d}z\;|w_{g}^{\perp}(z)|^{2}|w_{e,+}^{\perp}(z)|^{2}\right)^{2}\int\mathrm{d}x\;\left(w_{e,+}^{x}(x)\right)^{*}w_{e,-}^{x}(x-a_{0})|w_{g}^{x}(x)|^{2},$
(29)

$C_{eg}^{\pm}=g_{eg}^{\pm}\left(\int\mathrm{d}z\;|w_{g}^{\perp}(z)|^{2}|w_{e,+}^{\perp}(z)|^{2}\right)^{2}\int\mathrm{d}x\;\left(w_{e,+}^{x}(x)\right)^{*}w_{e,-}^{x}(x-a_{0})\left(w_{g}^{x}(x-a_{0})\right)^{*}w_{g}^{x}(x),$
(30)

$D_{eg}^{\pm}=g_{eg}^{\pm}\left(\int\mathrm{d}z\;|w_{g}^{\perp}(z)|^{2}|w_{e,+}^{\perp}(z)|^{2}\right)^{2}\int\mathrm{d}x\;|w_{g}^{x}(x)|^{2}|w_{e,-}^{x}(x-a_{0})|^{2},$
(31)

$J_{gg}=g_{gg}\left(\int\mathrm{d}z\;|w_{g}^{\perp}(z)|^{4}\right)^{2}\int\mathrm{d}x\;|w_{g}^{x}(x)|^{2}|w_{g}^{x}(x-a_{0})|^{2},$
(32)
$K_{gg}=g_{gg}\left(\int\mathrm{d}z\;|w_{g}^{\perp}(z)|^{4}\right)^{2}\int\mathrm{d}x\;\left(w_{g}^{x}(x)\right)^{*}\left(w_{g}^{x}(x-a_{0})\right)^{*}w_{g}^{x}(x)^{2}.$
(33)

The values of these amplitudes for the same parameters as in Tab. 2 are given
in Tab. 3. Note that some of these terms have non-zero matrix elements within
the resonant sector of $H_{0}$ of our interest: $t_{e}^{(1)}$ for example,
connects the single-block states
$\frac{1}{\sqrt{2}}(\ket{g^{1}e^{2}_{+}}-\ket{g^{2}e^{1}_{+}})$ and
$\frac{1}{\sqrt{2}}(\ket{g^{1}e^{2}_{-}}-\ket{g^{2}e^{1}_{-}})$, and would
break the condition of having a single rishon ($e$ atom) per link. The error
induced by this term is expected to remain small for the times $\sim
250\;\mathrm{ms}$ considered here. The next-nearest neighbor $g$ hopping
$t_{g}^{(1)}$, as well as the density mediated hoppings $B_{eg}^{\pm}$,
$K_{gg}$ have no matrix elements between states in the resonant sector. Like
$t_{g}$, $t_{e}$ and $A_{eg}^{+}$, they would contribute to the effective
Hamiltonian to second order in perturbation theory, but since their amplitudes
are smaller we do not expect them do have a drastic effect on the dynamics.
The terms corresponding to $C_{eg}^{\pm}$ could directly contribute to desired
the matter-gauge coupling, without invoking second-order process, but their
amplitudes are so small that they have no visible effect on the time scales we
consider. The terms $D_{eg}^{\pm}$ and $J_{gg}$ have diagonal matrix elements
on gauge-invariant states and contribute as small correction to $u$ and $w$.

Table 3: Values of additional parameters in the lattice Hamiltonian, obtained
for the same values as in Tab. 2. For the case of $eg$ interactions, we only
report the "$+$" amplitudes, as the "$-$" can be easily computed from them by
multiplying by the ratio $g_{eg}^{-}/g_{eg}^{+}$.

Finally, let us discuss the case $m\neq 0$, obtained with a small but non-zero
$V_{\alpha}^{S}$. This term gives the following contribution to the effective
lattice Hamiltonian

$H_{S}=\delta_{g}\sum_{x}s_{x}n_{x}^{g}+\delta_{e}\sum_{x}s_{x}(n_{x,+}^{e}+n_{x,-}^{e}),$
(34)

where $s_{x}=(-1)^{x}$ and

$\delta_{g}=\int\mathrm{d}x\;V_{g}^{S}(x)|w_{g}(x)|^{2},\hskip
28.45274pt\delta_{e}=\int\mathrm{d}x\;V_{e}^{S}(x)|w_{e,+}(x)|^{2}.$ (35)

Using $G_{x}=n_{x}^{g}+\sum_{r}n^{e}_{x,r}-2$, we get

$H_{S}=(\delta_{g}-\delta_{e})\sum_{x}s_{x}n_{x}^{g}+\delta_{e}\sum_{x}s_{x}(G_{x}+2).$
(36)

In the gauge-invariant subspace, this corresponds to a mass
$m=\delta_{g}-\delta_{e}$ for the fermionic matter [see Eq. (2.1)].

##  Appendix B Tuning of the optical lattice potential

The amplitudes of the terms in the lattice Hamiltonian, obtained with our ab-
initio calculation, depend on six tunable parameters $V_{\alpha}^{0}$,
$V_{\alpha}^{1}$, $V_{\alpha}^{\perp}$, where $\alpha=g,e$. These parameters
can be independently controlled by selecting the appropriate optical intensity
and the wavelength for each laser: the amplitude of the potential (for both
the $g$ and $e$ atoms) is proportional to the optical intensity, while the
ratio between the $g$ and $e$ amplitudes (which is determined by the ratio
between the dynamical polarizability, i.e. the ratio of the state-dependent
AC-Stark shifts) depends only on the wavelength. In Fig. 7 we plot the AC-
Stark shift in units of the optical intensity for the $g\equiv{}^{1}S_{0}$ and
$e\equiv{}^{3}P_{0}$ states of ytterbium as a function of the wavelength.

The large (positive) ratio $V_{e}^{1}/V_{g}^{1}$, required for the desired
strong localization of the $e$ Wannier functions, can be achieved by choosing
$\lambda_{1}\simeq 0.45\,\mathrm{\mu m}$ (leftmost dashed line in Fig. 7),
where the AC-Stark shifts have the same sign and that of the $e$ state is
large, because of the proximity with a resonance. This choice, compared to
other possible wavelengths that yield a similar large ratio
$V_{e}^{1}/V_{g}^{1}$, has the advantage that it gives the smallest lattice
spacing $a_{0}=\lambda_{1}\simeq 0.45\,\mathrm{\mu m}$, and, consequently, the
largest recoil energy. A large (in absolute value) negative
$V_{e}^{0}/V_{g}^{0}$ is achieved by choosing $\lambda_{0}\simeq
0.64\;\mathrm{\mu m}$ (rightmost dashed line in Fig. 7), while the large
positive ratio of transverse confining potentials
$V_{g}^{\perp}/V_{e}^{\perp}$ is obtained for $\lambda_{\perp}\simeq
0.56\;\mathrm{\mu m}$ or $\lambda_{\perp}\simeq 0.59\;\mathrm{\mu m}$. The
condition $a_{0}=\lambda_{1}$ can then be met by using an angle
$\theta_{0}\simeq 0.518\cdot\pi$.

Figure 7: AC-Stark shift $V$ (normalized to the optical intensity, i.e. in
units of Hz W/cm2) of the ${}^{1}S_{0}$ and ${}^{3}P_{0}$ states of ytterbium
as a function of wavelength. Both curves have been calculated exploiting
energy level data from the NIST Atomic Spectra database.

##  Appendix C Single block Hamiltonian

We now demonstrate how to derive Eq. (4.2) from the single block Hamiltonian
in Eq. (9). We first note that, by using the identity
$\lambda^{a}_{ij}\lambda^{a}_{kl}=2\delta_{il}\delta_{jk}-\frac{2}{N}\delta_{ij}\delta_{kl}$,
we can write

$\displaystyle\begin{split}\psi_{x}^{i\dagger}c_{x,r}^{j^{\dagger}}\psi_{x}^{j}c_{x,r}^{i}&=-\frac{1}{2}\psi_{x}^{i\dagger}\lambda_{ij}\psi_{x}^{j}c_{x,r}^{k\dagger}\lambda_{kl}^{a}c_{x,r}^{l}-\frac{1}{N}\psi_{x}^{i^{\dagger}}\psi_{x}^{i}c_{x,r}^{j\dagger}c_{x,r}^{j}\\\
&=-\frac{1}{2}M_{x}^{a}(G_{x}^{a}-M_{x}^{a})-\frac{1}{N}n_{x}^{g}n_{x,r}^{e},\end{split}$
(37)

where $M_{x}^{a}=\psi_{x}^{i\dagger}\lambda_{ij}\psi_{x}^{j}$. The quadratic
Casimir operator in $M_{x}^{a}$ can be expressed as

$\displaystyle\begin{split}\frac{1}{2}\sum_{a}(M_{x}^{a})^{2}&=\psi_{x}^{i\dagger}\psi_{x}^{j}\psi_{x}^{j\dagger}\psi_{x}^{i}-\frac{1}{N}\psi_{x}^{i\dagger}\psi_{x}^{i}\psi_{x}^{j\dagger}\psi_{x}^{j}\\\
&=(N+1)n_{x}^{g}-\left(1+\frac{1}{N}\right)(n_{x}^{g})^{2}.\end{split}$ (38)

We get

$\displaystyle\begin{split}h_{x}^{0}&=\mu_{g}n_{x}^{g}+\frac{U_{gg}}{2}n_{x}^{g}(n_{x}^{g}-1)+\sum_{r}\left[\mu_{e}n_{x,r}^{e}+\frac{U_{eg}^{+}+U_{eg}^{-}}{2}n_{x}^{g}n_{x,r}^{e}+\frac{U_{eg}^{+}-U_{eg}^{-}}{2}\psi_{x}^{i\dagger}c_{x,r}^{j\dagger}\psi_{x}^{j}c_{x,r}^{i}\right]\\\
&=\left[\mu_{g}-\frac{U_{gg}}{2}+\frac{U_{eg}^{+}-U_{eg}^{-}}{2}(N+1)\right]n_{x}^{g}+\frac{(N-1)U_{eg}^{+}+(N+1)U_{eg}^{-}}{2N}\left(n_{x}^{g}+\sum_{r}n_{x,r}^{e}\right)n_{x}^{g}\\\
&\qquad+\left(\frac{U_{gg}}{2}+U_{eg}^{+}\right)(n_{x}^{g})^{2}-\frac{U_{eg}^{+}-U_{eg}^{-}}{4}M_{x}^{a}G_{x}^{a}.\end{split}$
(39)

We can write $G_{x}=n_{x}^{g}+\sum_{r}n_{x,r}^{e}-C$ where $C$ is an
appropriate constant. With this definition, the single block Hamiltonian
$h_{x}^{0}$ becomes

$h_{x}^{0}=\mu_{g}n_{x}^{g}+\mu_{e}(G_{x}-n_{x}^{g}+C)+\frac{U_{gg}}{2}n_{x}^{g}(n_{x}^{g}-1)+U_{eg}^{+}n_{x}^{g}(G_{x}-n_{x}^{g}+C)-\frac{U_{eg}^{+}-U_{eg}^{-}}{4}M^{a}_{x}G^{a}_{x}.$
(40)

We can then use that $\sum_{x}n_{x}^{g}=\text{const.}$ and
$\sum_{x}G_{x}=\text{const.}$ to write

$H_{0}=\sum_{x}h_{x}^{0}=\sum_{x}\left[\frac{U_{gg}-2U_{eg}^{+}}{2}n_{x}^{g}(n_{x}^{g}-1)+U_{eg}^{+}n_{x}^{g}G_{x}-\frac{U_{eg}^{+}-U_{eg}^{-}}{4}M^{a}_{x}G^{a}_{x}\right]+\text{const.}$
(41)

###  C.1 Example: $N=2$

As an example, here we compute the eigenstates of $h^{0}_{x}$ on a single
block for the simple case $N=2$ in one spatial dimension. For convenience of
notation, we define the following quantities

$U_{eg}^{n}=\frac{U_{eg}^{+}+3U_{eg}^{-}}{4},\qquad
U_{eg}^{s}=\frac{U_{eg}^{-}-U_{eg}^{+}}{2},$ (42)

$n^{e}_{x}=\sum_{r}n_{x,r}^{e},\qquad
s_{x}^{g,a}=\frac{1}{2}\psi_{x}^{i\dagger}\sigma^{a}_{ij}\psi_{x}^{j},\qquad
s_{x}^{e,a}=\frac{1}{2}\sum_{r}c_{x,r}^{i\dagger}\sigma^{a}_{ij}c_{x,r}^{j}.$
(43)

The single block Hamiltonian can be conveniently expressed as

$h^{0}_{x}=\mu_{g}n^{g}_{x}+\mu_{e}n_{x}^{e}+U_{gg}\frac{n^{g}_{x}(n^{g}_{x}-1)}{2}+U_{eg}^{n}n^{e}_{x}n^{g}_{x}+2U_{eg}^{s}s_{x}^{e,a}s_{x}^{g,a}.$
(44)

In Tab. 4 we report the eigenstates of $h^{0}_{x}$, their quantum numbers, and
their energies. In Appendix D we use these states to compute the effective
Hamiltonian in perturbation theory.

State | $n_{x}^{g}$ | $n_{x}^{e}$ | $s_{x}^{g}$ | $s_{x}^{e}$ | $s_{x}^{\mathrm{tot}}$ | $s_{x}^{e,a}s_{x}^{g,a}$ | Energy
---|---|---|---|---|---|---|---
$\ket{0}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | 0
$\ket{g^{1}},\;\ket{g^{2}}$ | $1$ | $0$ | $1/2$ | $0$ | $1/2$ | $0$ | $\mu^{g}$
|  $\ket{e_{+}^{1}},\;\ket{e_{+}^{2}}$,
---
$\ket{e_{-}^{1}},\;\ket{e_{-}^{2}}$
$0$ | $1$ | $0$ | $1/2$ | $1/2$ | $0$ | $\mu^{e}$
$\ket{g^{1}g^{2}}$ | $2$ | $0$ | $0$ | $0$ | $0$ | $0$ | $2\mu^{g}+U_{gg}$
|  $\frac{1}{\sqrt{2}}(\ket{g^{1}e^{2}_{+}}-\ket{g^{2}e^{1}_{+}})$,
---
$\frac{1}{\sqrt{2}}(\ket{g^{1}e^{2}_{-}}-\ket{g^{2}e^{1}_{-}})$
$1$ | $1$ | $1/2$ | $1/2$ | $0$ | $-3/4$ | $\mu^{g}+\mu^{e}+U_{eg}^{+}$
|  $\frac{1}{\sqrt{2}}(\ket{g^{1}e^{2}_{+}}+\ket{g^{2}e^{1}_{+}})$,
---
$\ket{g^{1}e^{1}_{+}},\;\ket{g^{2}e^{2}_{+}},$
$\frac{1}{\sqrt{2}}(\ket{g^{1}e^{2}_{-}}+\ket{g^{2}e^{1}_{-}})$,
$\ket{g^{1}e^{1}_{-}},\;\ket{g^{2}e^{2}_{-}},$
$1$ | $1$ | $1/2$ | $1/2$ | $1$ | $1/4$ | $\mu^{g}+\mu^{e}+U_{eg}^{-}$
| $\ket{e^{1}_{+}e^{2}_{+}},\;\ket{e^{1}_{-}e^{2}_{-}},$
---
$\frac{1}{\sqrt{2}}(\ket{e^{1}_{+}e^{2}_{-}}-\ket{e^{1}_{+}e^{2}_{-}})$
$0$ | $2$ | $0$ | $0$ | $0$ | $0$ | $2\mu_{e}$
| $\ket{e^{1}_{+}e^{1}_{-}},\;\ket{e^{2}_{+}e^{2}_{-}},$
---
$\frac{1}{\sqrt{2}}(\ket{e^{1}_{+}e^{2}_{-}}+\ket{e^{1}_{+}e^{2}_{-}})$
$0$ | $2$ | $0$ | $1$ | $1$ | $0$ | $2\mu_{e}$
| $\ket{g^{1}g^{2}e^{1}_{+}},\;\ket{g^{1}g^{2}e^{2}_{+}},$
---
$\ket{g^{1}g^{2}e^{1}_{-}},\;\ket{g^{1}g^{2}e^{2}_{+}}$
$2$ | $1$ | $0$ | $1/2$ | $1/2$ | $0$ | $2\mu^{g}+\mu^{e}+U_{gg}+2U_{eg}^{n}$
| $\ket{g^{1}e^{1}_{+}e^{2}_{+}},\;\ket{g^{2}e^{1}_{+}e^{2}_{+}},$
---
$\ket{g^{1}e^{1}_{-}e^{2}_{-}},\;\ket{g^{2}e^{1}_{-}e^{2}_{-}}$
$1$ | $2$ | $0$ | $1/2$ | $1/2$ | $0$ | $\mu^{g}+2\mu^{e}+2U_{eg}^{n}$
| $\ket{e^{1}_{+},e^{2}_{+},e^{1}_{-}},\;\ket{e^{1}_{+},e^{2}_{+},e^{2}_{-}},$
---
$\ket{e^{1}_{+},e^{1}_{-},e^{2}_{-}},\;\ket{e^{2}_{+},e^{1}_{-},e^{2}_{-}}$
$0$ | $3$ | $0$ | $1/2$ | $1/2$ | $0$ | $3\mu^{e}$
|
$\frac{1}{\sqrt{2}}(\ket{g^{1}e^{1}_{+}e^{2}_{-}}-\ket{g^{1}e^{2}_{+}e^{2}_{-}}),$
---
$\frac{1}{\sqrt{2}}(\ket{g^{2}e^{1}_{+}e^{2}_{-}}-\ket{g^{2}e^{2}_{+}e^{2}_{-}})$
$1$ | $2$ | $1/2$ | $0$ | $1/2$ | $0$ | $\mu^{g}+2\mu^{e}+2U_{eg}^{n}$
| $\ket{g^{1}e^{1}_{+}e^{1}_{-}},\;\ket{g^{2}e^{2}_{+}e^{2}_{-}},$
---
$\frac{1}{\sqrt{3}}(\ket{g^{1}e^{2}_{+}e^{2}_{-}}+\ket{g^{2}e^{1}_{+}e^{2}_{-}}+\ket{g^{2}e^{2}_{+}e^{1}_{-}}),$
$\frac{1}{\sqrt{3}}(\ket{g^{2}e^{1}_{+}e^{1}_{-}}+\ket{g^{1}e^{2}_{+}e^{1}_{-}}+\ket{g^{1}e^{1}_{+}e^{2}_{-}})$
$1$ | $2$ | $1/2$ | $1$ | $3/2$ | $1/2$ | $\mu^{g}+2\mu^{e}+2U_{eg}^{-}$
|
$\frac{1}{\sqrt{6}}(2\ket{g^{1}e^{2}_{+}e^{2}_{-}}-\ket{g^{2}e^{1}_{+}e^{2}_{-}}-\ket{g^{2}e^{2}_{+}e^{1}_{-}})$
---
$\frac{1}{\sqrt{6}}(2\ket{g^{2}e^{1}_{+}e^{1}_{-}}-\ket{g^{1}e^{2}_{+}e^{1}_{-}}-\ket{g^{1}e^{1}_{+}e^{2}_{-}})$
$1$ | $2$ | $1/2$ | $1$ | $1/2$ | $-1$ | $\mu^{g}+2\mu^{e}+2(U_{eg}^{n}-U_{eg}^{s})$
Table 4: Eigenstates of $h^{0}$ for number of particles
$n_{x}^{e}+n_{x}^{g}\leq 3$. The notation implies proper antisymmetrization of
the fermionic wavefunction (for example
$\ket{g^{1}g^{2}e_{+}^{2}}\equiv\psi_{x}^{1\dagger}\psi_{x}^{2\dagger}c_{x,+}^{2\dagger}\ket{0}$).

##  Appendix D Effective Hamiltonian

We here derive the effective Hamiltonian to second order in perturbation
theory for the case $N=2$ in one spatial dimension. The perturbation can be
written as a sum of terms that act over pairs of neighboring blocks:

$H_{1}=\sum_{x}h_{x,x+1}^{1},$ (45)
$h_{x,x+1}^{1}=-t_{g}(\psi^{i\dagger}_{x}\psi^{i}_{x+1}+\text{H.c.})-t_{e}(c^{i\dagger}_{x,+}c^{i}_{x+1,-}+\text{H.c.})+\frac{A_{eg}^{+}+A_{eg}^{-}}{2}(n_{x,+}^{e}+n_{x+1,-}^{e})(\psi_{x}^{i\dagger}\psi_{x+1}^{i}+\text{H.c.})\\\
+\frac{A_{eg}^{-}-A_{eg}^{+}}{2}(c_{x,+}^{i\dagger}c_{x,+}^{j}\psi_{x}^{j\dagger}\psi_{x+1}^{i}+c_{x+1,-}^{i\dagger}c_{x+1,-}^{j}\psi_{x}^{j\dagger}\psi_{x+1}^{i}+\text{H.c.}),$
(46)

where longer range terms and other interactions are neglected (we check the
validity of this approximation for realistic parameters in Appendix A). It is
easy to check that $h_{x,x+1}$ has no matrix elements between gauge-invariant
states. The second-order effective Hamiltonian on the gauge-invariant subspace
can then be obtained then as a sum of two body operators

$H_{\text{eff}}=\sum_{x}h_{x,x+1}^{\text{eff}},\qquad\bra{j}{h_{x,x+1}^{\text{eff}}}\ket{l}=\sum_{k}(\bra{k}{h_{x,x+1}^{1}}\ket{j})^{*}\bra{k}{h_{x,x+1}^{1}}\ket{l}\omega_{jk}^{-1},$
(47)

where $\ket{j}$ and $\ket{l}$ belong to the Hilbert space of two blocks and
can be represented as the tensor product of gauge-invariant eigenstates of
$h^{0}$, with the additional constraint that the link connecting them hosts a
single $e$ particle. The states $\ket{k}$ also belong to the Hilbert space of
two blocks, but are the tensor products of eigenstates of $h^{0}$ that are not
gauge-invariant. The states $\ket{j}$ and $\ket{k}$ are degenerate eigenstates
of $h_{x}^{0}+h_{x}^{1}$ (we assume the resonance condition
$U_{gg}=2U_{eg}^{+}$), with eigenvalue $\epsilon_{j}$, while the intermediate
state $\ket{k}$ has eigenvalue $\epsilon_{k}$. We defined
$\omega_{jk}=\epsilon_{j}-\epsilon_{k}$.

Since $h_{x,x+1}^{1}$ conserves the total number of $g$ atoms and the total
number of $e$ atoms on the two blocks, the effective Hamiltonian has a block
structure, where each block is labelled by these two numbers of atoms. We now
compute the effective Hamiltonian separately for each of these blocks.

###  D.1 Sector with 3 $g$ atoms and 1 $e$ atom

A basis for the gauge-invariant states on this block is defined by the two
states

$\ket{0_{\text{gi}}^{(3,1)}}=\frac{1}{\sqrt{2}}\epsilon_{ab}\ket{g^{1}g^{2}}_{x}\ket{g^{a}e_{-}^{b}}_{x+1},\qquad\ket{1_{\text{gi}}^{(3,1)}}=\frac{1}{\sqrt{2}}\epsilon_{ab}\ket{g^{a}e_{+}^{b}}_{x}\ket{g^{1}g^{2}}_{x+1},$
(48)

where we defined

$\displaystyle\epsilon_{ab}=\begin{cases}1\;&a=1,b=2\\\ -1\;&a=2,b=1\\\
0\;&a=b.\end{cases}$ (49)

The non-zero matrix elements of $h_{x,x+1}^{(1)}$ between these states and
non-gauge invariant states take the form

$\left(h_{x,x+1}^{1}\right)_{kj}^{(3,1)}=\begin{pmatrix}t_{g}+\dfrac{1}{2}(A_{eg}^{+}+A_{eg}^{-})&-t_{e}\\\
-t_{e}&t_{g}+\dfrac{1}{2}(A_{eg}^{+}+A_{eg}^{-})\end{pmatrix},$ (50)

where the column index $j=0,1$ labels the gauge-invariant basis states
$\ket{j_{\text{gi}}^{(3,1)}}$ defined above, while the row index $k=0,1$ is
for the non-gauge-invariant states $\ket{k_{\text{ngi}}^{(3,1)}}$

$\ket{0_{\text{ngi}}^{(3,1)}}=\frac{1}{\sqrt{2}}\epsilon_{ab}\ket{g^{a}}_{x}\ket{g^{1}g^{2}e^{b}_{-}}_{x+1},\qquad\ket{1_{\text{ngi}}^{(3,1)}}=\frac{1}{\sqrt{2}}\epsilon_{ab}\ket{e^{a}_{+}}_{x}\ket{g^{1}g^{2}e^{b}_{-}}_{x+1}.$
(51)

For both values of $k$, one finds
$\omega_{j,k}=U_{eg}^{+}-2U_{eg}^{n}=(U_{eg}^{+}-3U_{eg}^{-})/2$. The
effective Hamiltonian in this sector then reads

$\left(h_{x,x+1}^{\text{eff}}\right)^{(3,1)}_{j_{1},j_{2}}=\dfrac{2}{U_{eg}^{+}-3U_{eg}^{-}}\begin{pmatrix}\left(t_{g}+\frac{A_{eg}^{+}+A_{eg}^{-}}{2}\right)^{2}+t_{e}^{2}&-t_{e}(2t_{g}+A_{eg}^{+}+A_{eg}^{-})\\\
-t_{e}(2t_{g}+A_{eg}^{+}+A_{eg}^{-})&\left(t_{g}+\frac{A_{eg}^{+}+A_{eg}^{-}}{2}\right)^{2}+t_{e}^{2}\end{pmatrix}$
(52)

###  D.2 Sector with 2 $g$ atoms and 2 $e$ atoms

Similarly to the previous case, we define a basis
$\ket{j_{\text{gi}}^{(2,2)}}$ for the gauge-invariant states

$\ket{0_{\text{gi}}^{(2,2)}}=\frac{1}{\sqrt{2}}\epsilon_{ab}\ket{g^{1}g^{2}}_{x}\ket{e_{+}^{a}e_{-}^{b}}_{x+1},\qquad\ket{1_{\text{gi}}^{(2,2)}}=\frac{1}{2}\epsilon_{ab}\epsilon_{cd}\ket{g^{a}e_{+}^{b}}_{x}\ket{g^{c}e_{+}^{d}}_{x+1},$
(53)

and a basis $\ket{k_{\text{ngi}}^{(2,2)}}$ for the non-gauge invariant states:

$\ket{0_{\text{ngi}}^{(2,2)}}=\frac{1}{2}\epsilon_{ab}\epsilon_{cd}\ket{g^{a}}_{x}\ket{g^{b}e^{c}_{+}e^{d}_{-}}_{x+1},\qquad\ket{1_{\text{ngi}}^{(2,2)}}=\frac{1}{\sqrt{2}}\epsilon_{ab}\ket{g^{1}g^{2}e^{a}_{+}}_{x}\ket{e^{b}_{+}}_{x+1},$
$\ket{2_{\text{ngi}}^{(2,2)}}=\frac{1}{2\sqrt{3}}(1-\delta_{ab})\ket{g^{a}}_{x}\left(2\ket{g^{a}e^{b}_{+}e^{b}_{-}}_{x+1}-\ket{g^{b}e^{a}_{+}e^{b}_{-}}_{x+1}-\ket{g^{b}e^{b}_{+}e^{a}_{-}}_{x+1}\right),$
$\ket{3_{\text{ngi}}^{(2,2)}}=\frac{1}{\sqrt{2}}\epsilon_{ab}\ket{e^{a}_{+}}_{x}\ket{g^{1}g^{2}e^{b}_{-}}_{x+1}.$
(54)

Since $h_{x,x+1}^{1}$ does not change the occupations $n_{x,-}^{e}$ and
$n_{x+1,+}^{e}$, so we only consider here the two states that have
$n_{x,-}^{e}=0,n_{x+1,+}^{e}=1$ without loss of generality. Because of
inversion symmetry, the results can be directly generalised to the sector with
$n_{x,-}^{e}=1,n_{x+1,+}^{e}=0$. The matrix elements of $h_{x,x+1}^{1}$ in
this basis read

$\left(h_{x,x+1}^{1}\right)_{kj}^{(2,2)}=\begin{pmatrix}-\dfrac{1}{2\sqrt{2}}(4t_{g}+3A_{eg}^{-}+A_{eg}^{+})&\dfrac{1}{2}t_{e}\\\
-t_{e}&\dfrac{1}{2\sqrt{2}}(2t_{g}+A_{eg}^{-}+A_{eg}^{+})\\\
\dfrac{1}{2\sqrt{6}}(A_{eg}^{+}-A_{eg}^{-})&-\dfrac{\sqrt{3}}{2}t_{e}\\\
0&\dfrac{1}{\sqrt{2}}(t_{g}+A_{eg}^{-})\end{pmatrix},$ (55)

and the energy differences are

$\displaystyle\omega_{j,0}$
$\displaystyle=2U_{eg}^{+}-2U_{eg}^{n}=3(U_{eg}^{+}-U_{eg}^{-})/2,$ (56)
$\displaystyle\omega_{j,1}$
$\displaystyle=2U_{eg}^{+}-U_{gg}-2U_{eg}^{n}=-(U_{eg}^{+}+3U_{eg}^{-})/2,$
(57) $\displaystyle\omega_{j,2}$
$\displaystyle=2U_{eg}^{+}-2U_{eg}^{n}+2U_{eg}^{s}=(U_{eg}^{+}-U_{eg}^{-})/2,$
(58) $\displaystyle\omega_{j,3}$
$\displaystyle=2U_{eg}^{+}-U_{gg}-2U_{eg}^{n}=-(U_{eg}^{+}+3U_{eg}^{-})/2.$
(59)

We then get the effective Hamiltonian for this sector

$\left(h_{x,x+1}^{\text{eff}}\right)^{(2,2)}_{j_{1},j_{2}}=\begin{pmatrix}\frac{(A_{eg}^{+}-A_{eg}^{-})^{2}+(4t_{g}+3A_{eg}^{-}+A_{eg}^{+})^{2}}{12(U_{eg}^{+}-U_{eg}^{-})}-\frac{2t_{e}^{2}}{(U_{eg}^{+}+3U_{eg}^{-})}&-\frac{\sqrt{2}t_{e}(t_{g}+A_{eg}^{+})}{3(U^{+}_{eg}-U^{-}_{eg})}+\frac{t_{e}(A_{eg}^{+}+A_{eg}^{-}+2t_{g})}{\sqrt{2}(U^{+}_{eg}+3U^{-}_{eg})}\\\
-\frac{\sqrt{2}t_{e}(t_{g}+A_{eg}^{+})}{3(U^{+}_{eg}-U^{-}_{eg})}+\frac{t_{e}(A_{eg}^{+}+A_{eg}^{-}+2t_{g})}{\sqrt{2}(U^{+}_{eg}+3U^{-}_{eg})}&-\frac{4(t_{g}+A_{eg}^{-})^{2}+(2t_{g}+A_{eg}^{-}+A_{eg}^{+})^{2}}{4(U^{+}_{eg}+3U^{-}_{eg})}+\frac{5t_{e}^{2}}{3(U^{+}_{eg}-U^{-}_{eg})}.\end{pmatrix}$
(60)

###  D.3 Sector with 1 $g$ atom and 3 $e$ atoms

We now apply the same procedure to the sector with 1 $g$ atom and 3 $e$ atoms.
The basis states $\ket{j_{\text{gi}}^{(1,3)}}$ for the gauge-invariant
subspace in this sector are

$\ket{0_{\text{gi}}^{(1,3)}}=\frac{1}{2}\epsilon_{ab}\epsilon_{cd}\ket{g^{a}e_{-}^{b}}_{x}\ket{e_{+}^{c}e_{-}^{d}}_{x+1},\qquad\ket{1_{\text{gi}}^{(1,3)}}=\frac{1}{2}\epsilon_{ab}\epsilon_{cd}\ket{e_{+}^{a}e_{-}^{b}}_{x}\ket{g^{c}e_{+}^{d}}_{x+1}.$
(61)

We then define the non-gauge-invariant states $\ket{k_{\text{ngi}}^{(1,3)}}$
as

$\ket{0_{\text{ngi}}^{(1,3)}}=\frac{1}{2}\epsilon_{ab}\epsilon_{cd}\ket{g^{a}e^{c}_{+}e^{d}_{-}}_{x}\ket{e^{b}_{+}}_{x+1},\qquad\ket{1_{\text{ngi}}^{(1,3)}}=\frac{1}{2}\epsilon_{ab}\epsilon_{cd}\ket{e^{a}_{-}}_{x}\ket{g^{b}e^{c}_{+}e^{d}_{-}}_{x+1},$
$\ket{2_{\text{ngi}}^{(1,3)}}=\frac{1}{2\sqrt{3}}(1-\delta_{ab})\left(2\ket{g^{a}e^{b}_{+}e^{b}_{-}}_{x}-\ket{g^{b}e^{a}_{+}e^{b}_{-}}_{x+1}-\ket{g^{b}e^{b}_{+}e^{a}_{-}}_{x+1}\right)\ket{e^{a}_{+}}_{x+1},$
$\ket{3_{\text{ngi}}^{(1,3)}}=\frac{1}{2\sqrt{3}}(1-\delta_{ab})\ket{e^{a}_{-}}_{x}\left(2\ket{g^{a}e^{b}_{+}e^{b}_{-}}_{x+1}-\ket{g^{b}e^{a}_{+}e^{b}_{-}}_{x+1}-\ket{g^{b}e^{b}_{+}e^{a}_{-}}_{x+1}\right).$
(62)

We compute the matrix elements of $h_{x,x+1}^{1}$:

$\left(h_{x,x+1}^{1}\right)_{kj}^{(1,3)}=\begin{pmatrix}\dfrac{t_{e}}{2}&-t_{g}-\dfrac{1}{4}(A_{eg}^{+}+3A_{eg}^{-})\\\
-t_{g}-\dfrac{1}{4}(A_{eg}^{+}+3A_{eg}^{-})&\dfrac{t_{e}}{2}\\\
-\dfrac{\sqrt{3}}{2}t_{e}&\dfrac{1}{4\sqrt{3}}(A_{eg}^{+}-A_{eg}^{-})\\\
\dfrac{1}{4\sqrt{3}}(A_{eg}^{+}-A_{eg}^{-})&-\dfrac{\sqrt{3}}{2}t_{e}\end{pmatrix},$
(63)

and the energy differences

$\displaystyle\omega_{j,0}=\omega_{j,1}$
$\displaystyle=U_{eg}^{+}-2U_{eg}^{n}=(U_{eg}^{+}-3U_{eg}^{-})/2,$ (64)
$\displaystyle\omega_{j,2}=\omega_{j,3}$
$\displaystyle=U_{eg}^{+}-2U_{eg}^{n}+2U_{eg}^{s}=-(U_{eg}^{+}+U_{eg}^{-})/2.$
(65)

We finally get the effective Hamiltonian in this sector

$\left(h_{x,x+1}^{\text{eff}}\right)^{(1,3)}_{j_{1},j_{2}}=\begin{pmatrix}\frac{(4t_{g}+A_{eg}^{+}+3A_{eg}^{-})^{2}+4t_{e}^{2}}{8(U^{+}_{eg}-3U^{-}_{eg})}-\frac{36t_{e}^{2}+(A_{eg}^{+}-A_{eg}^{-})^{2}}{24(U^{+}_{eg}+U^{-}_{eg})}&\frac{t_{e}(A_{eg}^{+}-A_{eg}^{-})}{2(U_{eg}^{+}+U_{eg}^{-})}-\frac{t_{e}(4t_{g}+A_{eg}^{+}+3A_{eg}^{-})}{2(U^{+}_{eg}-3U^{-}_{eg})}\\\
\frac{t_{e}(A_{eg}^{+}-A_{eg}^{-})}{2(U_{eg}^{+}+U_{eg}^{-})}-\frac{t_{e}(4t_{g}+A_{eg}^{+}+3A_{eg}^{-})}{2(U^{+}_{eg}-3U^{-}_{eg})}&\frac{(4t_{g}+A_{eg}^{+}+3A_{eg}^{-})^{2}+4t_{e}^{2}}{8(U^{+}_{eg}-3U^{-}_{eg})}-\frac{36t_{e}^{2}+(A_{eg}^{+}-A_{eg}^{-})^{2}}{24(U^{+}_{eg}+U^{-}_{eg})}.\end{pmatrix}$
(66)

###  D.4 Parameters of the effective Hamiltonian

Once we have the effective 2-body Hamiltonian $h_{x,x+1}$ in its matrix forms,
with the blocks defined in Eqs. (52), (60) and (66), it is useful to express
it in terms of the local densities and creation/annihilation operators. To
this end, we split $h_{x,x+1}$ in its diagonal and off-diagonal parts

$h_{x,x+1}^{\text{eff}}=h_{x,x+1}^{\text{eff, diag}}+h_{x,x+1}^{\text{eff,
off-diag}},$ (67)

and we make the following ansatz for the diagonal part:

$h_{x,x+1}^{\text{eff,diag}}=wn_{x}^{g}n_{x+1}^{g}+\frac{1}{2}\sum_{y=x,x+1}\left(\tilde{\mu}_{g}n_{y}^{g}+\tilde{\mu}_{e}n_{y,+}^{e}+\tilde{\mu}_{e}n_{y,-}^{e}+\frac{u}{2}n_{y}^{g}(n_{y}^{g}-1)\right).$
(68)

The coefficients $\tilde{\mu}_{g},\tilde{\mu}_{e},u,w$ can then be obtained by
evaluating Eq. (68) for the gauge-invariant basis states defined in each
sector, and setting them to be equal to the diagonal matrix elements computed
above. We thus have to solve the linear system of equations

$\begin{dcases}\left(h_{x,x+1}^{\text{eff}}\right)^{(3,1)}_{0,0}=2w+\frac{3}{2}\tilde{\mu}_{g}+\frac{1}{2}\tilde{\mu}_{e}+\frac{u}{2}\\\
\left(h_{x,x+1}^{\text{eff}}\right)^{(2,2)}_{0,0}=\tilde{\mu}_{g}+\tilde{\mu}_{e}+\frac{u}{2}\\\
\left(h_{x,x+1}^{\text{eff}}\right)^{(2,2)}_{1,1}=w+\tilde{\mu}_{g}+\tilde{\mu}_{e}\\\
\left(h_{x,x+1}^{\text{eff}}\right)^{(1,3)}_{0,0}=\frac{1}{2}\tilde{\mu}_{g}+\frac{3}{2}\tilde{\mu}_{e}.\end{dcases}$
(69)

Since the total numbers of $g$ and $e$ atoms are conserved, the terms
$\tilde{\mu}_{g}$ and $\tilde{\mu}_{g}$ only contribute as inconsequential
constants to the effective Hamiltonian. The on-site and nearest-neighbor mater
interactions ($u$ and $w$ respectively) play a role in the dynamics. The
solution of the linear systems yields

$u=2\left(h_{x,x+1}^{\text{eff}}\right)^{(3,1)}_{0,0}-4\left(h_{x,x+1}^{\text{eff}}\right)^{(2,2)}_{1,1}+2\left(h_{x,x+1}^{\text{eff}}\right)^{(1,3)}_{0,0},$
(70)
$w=\left(h_{x,x+1}^{\text{eff}}\right)^{(3,1)}_{0,0}-\left(h_{x,x+1}^{\text{eff}}\right)^{(2,2)}_{0,0}-\left(h_{x,x+1}^{\text{eff}}\right)^{(2,2)}_{1,1}+\left(h_{x,x+1}^{\text{eff}}\right)^{(1,3)}_{0,0}.$
(71)

Similarly to the diagonal part, we formulate the following ansatz for the off-
diagonal terms:

$h_{x,x+1}^{\text{eff, off-
diag}}=-[\tau_{1}+\tau_{2}E_{x-1,x}E_{x+1,x+2}+\tau_{3}(E_{x-1,x}-E_{x+1,x+2})](\psi_{x}^{i\dagger}U_{x,x+1}^{ij}\psi_{x+1}^{j}+\text{H.c.}).$
(72)

The system of equations that we aim to solve is

$\begin{dcases}\left(h_{x,x+1}^{\text{eff}}\right)^{(3,1)}_{0,1}=-2\left(\tau_{1}-\frac{1}{4}\tau_{2}-\tau_{3}\right)\\\
\left(h_{x,x+1}^{\text{eff}}\right)^{(2,2)}_{0,1}=-\sqrt{2}\left(\tau_{1}+\frac{1}{4}\tau_{2}\right)\\\
\left(h_{x,x+1}^{\text{eff}}\right)^{(1,3)}_{0,1}=-\left(\tau_{1}-\frac{1}{4}\tau_{2}+\tau_{3}\right).\end{dcases}$
(73)

We obtain

$\displaystyle\tau_{1}$
$\displaystyle=-\frac{1}{8}\left(h_{x,x+1}^{\text{eff}}\right)^{(3,1)}_{0,1}-\frac{1}{2\sqrt{2}}\left(h_{x,x+1}^{\text{eff}}\right)^{(2,2)}_{0,1}-\frac{1}{4}\left(h_{x,x+1}^{\text{eff}}\right)^{(1,3)}_{0,1},$
(74) $\displaystyle\tau_{2}$
$\displaystyle=\frac{1}{2}\left(h_{x,x+1}^{\text{eff}}\right)^{(3,1)}_{0,1}-\sqrt{2}\left(h_{x,x+1}^{\text{eff}}\right)^{(2,2)}_{0,1}+\left(h_{x,x+1}^{\text{eff}}\right)^{(1,3)}_{0,1},$
(75) $\displaystyle\tau_{3}$
$\displaystyle=\frac{1}{4}\left(h_{x,x+1}^{\text{eff}}\right)^{(3,1)}_{0,1}-\frac{1}{2}\left(h_{x,x+1}^{\text{eff}}\right)^{(1,3)}_{0,1}.$
(76)

##  Appendix E Initial state preparation

In order to initialize the atomic system such that only gauge-invariant (spin
singlet) states are occupied, it is convenient to prepare first a one-
dimensional $g$-state band-insulating phase in the $a_{0}$-spacing lattice
$V_{g}^{0}(x)$. Thereby, each site of the lattice is occupied by $N$ atoms in
a spin antisymmetric state. For $N=2$, this corresponds to the $\ket{g,g}$
state. Let us now consider the available two-particle spin-singlet states for
a certain block in 1D, i.e.:

$\left\\{\ket{g,g},\frac{1}{\sqrt{2}}\left(\ket{g,e\pm}+\ket{e\pm,g}\right),\frac{1}{\sqrt{2}}\left(\ket{e+,e-}+\ket{e-,e+}\right)\right\\}.$
(77)

Among these, states with different electronic state populations can be coupled
together by laser light driving the 1S${}_{0}\rightarrow{}^{3}$P0 clock
transition, suitably polarized and tuned to be selectively resonant with the
transitions between different two-particle states [28], which are split owing
to the different interaction strenghts in the singlet and triplet channels. In
particular, the relevant Rabi couplings for $\pi$-polarized clock light are
given by:

$\displaystyle\frac{\hbar}{2}\Omega^{eg\pm}$
$\displaystyle=\left(\bra{g,e\pm}\pm\bra{e\pm,g}\right)H_{\pi}^{(2)}\ket{g,g}=\frac{\hbar}{\sqrt{2}}(\Omega_{\uparrow}\mp\Omega_{\downarrow})$
(78a) $\displaystyle\frac{\hbar}{2}\Omega^{ee}$
$\displaystyle=\left(\bra{g,e\pm}\pm\bra{e\pm,g}\right)H_{\pi}^{(2)}\left(\ket{e+,e-}+\ket{e-,e+}\right)=-\frac{\hbar}{\sqrt{2}}(\Omega_{\uparrow}\mp\Omega_{\downarrow})$
(78b)

where $H_{\pi}^{(2)}$ is the two-particle atom-light interaction Hamiltonian
with clock $\pi$-polarized light, and $\Omega_{i}$ is the single-particle
clock Rabi frequency of spin state $i=\uparrow,\downarrow$ with $\pi$
polarization. It is thus possible to selectively couple all onsite $\ket{g,g}$
pairs to the desired singlet state by means of a global clock light pulse. The
relative phase $\varphi$ between the long and the short lattice potentials
$V_{\alpha}^{0}(x)$ and $V_{\alpha}^{1}(x)$ can be further exploited to
selectively prepare the states
$\left(\bra{g,e\pm}+\bra{e\pm,g}\right)/\sqrt{2}$. The experimental protocols
to create the $\left(\ket{e+,e-}+\ket{e-,e+}\right)/\sqrt{2}$ on each odd or
even block (see Fig. 1(c) of the main text), or the
$\left(\bra{g,e-}+\bra{e-,g}\right)$ state on each block are sketched in Fig.
8.

Figure 8: Experimental protocols for gauge-invariant state preparation.
(a1-a2) A $SU(N=2)$ band insulator of $g$-state atoms is prepared first in the
long lattice $V^{0}_{\alpha}$, and the short lattice $V^{1}_{\alpha}$ is
subsequently added. (b1-b2) Preparation sequence for obtaining a gauge-
invariant staggered block configuration. A finite shift
$\delta_{e}-\delta_{g}\neq 0$ created by the staggering lattice
$V^{S}_{\alpha}$ allows for selecting only even or odd blocks, where the
onsite $\ket{g,g}$ singlet is transferred to
$\left(\ket{e+,e-}+\ket{e-,e+}\right)$, occupying both left and right rishon
sites of the block. (c) Preparation sequence for obtaining a homogeneous
gauge-invariant configuration. A relative phase $\varphi\neq 0$ between the
long $V^{0}_{\alpha}(x)$ and the short $V^{1}_{\alpha}(x)$ lattices allows for
transferring on each block the onsite $\ket{g,g}$ singlet to a $e-g$ singlet,
occupying the left rishon site and possessing an interaction energy
$U_{eg}^{+}$. Subsequently, $\varphi$ is ramped to 0, and the staggering
lattice $V^{S}_{\alpha}$ is turned on adiabatically. (d) The combination of
two 1D protocols illustrated above allows to initialize gauge-invariant states
in a 2D brick-wall lattice (greyed-out orange boxes highlight double wells
that are initialized with no $e$ atoms).