Effective multi-body SU($N$)-symmetric interactions of ultracold fermionic atoms on a 3-D lattice

Rapid advancements in the experimental capabilities with ultracold alkaline-earth-like atoms (AEAs) bring to a surprisingly near term the prospect of performing quantum simulations of spin models and lattice field theories exhibiting SU($N$) symmetry. Motivated in particular by recent experiments preparing high density samples of strongly interacting ${}^{87}$Sr atoms in a three-dimensional optical lattice, we develop a low-energy effective theory of fermionic AEAs which exhibits emergent multi-body SU($N$)-symmetric interactions, where $N$ is the number of atomic nuclear spin levels. Our theory is limited to the experimental regime of (i) a deep lattice, with (ii) at most one atom occupying each nuclear spin state on any lattice site. The latter restriction is a consequence of initial ground-state preparation. We fully characterize the low-lying excitations in our effective theory, and compare predictions of many-body interaction energies with direct measurements of many-body excitation spectra in an optical lattice clock. Our work makes the first step in enabling a controlled, bottom-up experimental investigation of multi-body SU($N$) physics.


I. INTRODUCTION
Fermionic alkaline-earth atoms (AEAs), in addition to other atoms such as ytterbium (Yb) sharing similar electronic structure, are currently the building blocks of the most precise atomic clocks in the world [1][2][3]. These atoms have a unique, ultra-narrow optical transition between metastable 1 S 0 and 3 P 0 electronic orbital states, i.e. the "clock states", that allows for coherence times which can exceed 100 seconds [4,5]. Furthermore, AEAs can be trapped in fully controllable optical lattice potentials and interrogated with ultrastable lasers that can resolve and probe their rich hyperfine spectra, consisting of N different nuclear spin levels with N as large as 10 in strontium ( 87 Sr) and 6 in ytterbium ( 173 Yb).
Last year (2017), a new generation of OLCs became operational at JILA, interrogating a Fermi degenerate gas of 87 Sr atoms in a 3-D lattice at nanokelvin temperatures [8]. All of these atoms' degrees of freedom, including the electronic orbital, nuclear spin, and motional states, can be fully controlled with high fidelity in a 3-D lattice [9][10][11][12]. With frequency measurements reaching the 10 −19 fractional uncertainty level, the new OLCs are thus enabling an exciting opportunity to probe, for the first time, quantum dynamics with sub-millihertz spectral resolution [8].
A wonderful consequence of the efforts to build better clocks is the development of highly controllable quantum simulators of many-body systems in the strongly-interacting regime, where inter-particle interactions set the largest energy scale relevant for system dynamics [8,13,14]. The marriage between precision clock spectroscopy and quantum many-body physics [15][16][17][18][19][20] has an enormous potential to enable novel explorations of physics for the same reason that makes AEAs such remarkable time-keepers. Specifically, due to the lack of electronic orbital angular momentum in the 1 S 0 and 3 P 0 states, AEAs exhibit decoupled orbital and nuclear spin degrees of freedom. For atoms with N nuclear spin levels, this decoupling leads to nuclear-spin-conserving SU(N )-symmetric interactions governed entirely by orbital-state parameters [15,19,21].
The presence of this exotic SU(N ) symmetry in a highly controllable experimental platform opens the door to experimental studies of e.g. the SU(N ) Heisenberg model, whose phase diagram is believed to exhibit features such as a chiral spin liquid (CSL) phase with topological order and fractional statistics [22][23][24]. In addition to illuminating open questions in our understanding of the fractional quantum Hall effect and unconventional superconductivity [25][26][27], the CSL can support non-Abelian excitations which allow for universal topological quantum computation [23,28]. Harnessing the SU(N )-symmetric interactions of AEAs might also enable the simulation of various lattice gauge theories [29,30], some of which share important qualitative features with quantum chromodynamics such as fewbody bound states and confinement [31,32]. These direct, quantum simulations have an extraordinary potential to provide novel insights by circumventing e.g. severe sign problems which plague classical simulations of strongly interacting fermionic systems [29,33].
In this work, we investigate the first experimental capabilities with ultracold fermionic AEAs to prepare high-density samples in a 3-D lattice with multiple occupation of individual lattice sites [34]. Specifically, we consider the preparation of isolated few-body systems in the deep-lattice limit, and carry out a bottom-up investigation of emergent multi-body interactions on multiply-occupied lattice sites. These multi-body interactions appear in a low-energy effective theory of the atoms, and inherit the SU(N ) symmetry of their bare, pair-wise interactions, thereby enabling experimental studies of multi-body SU(N ) physics through the exquisite capabilities with OLCs.
Though effective multi-body interactions have previously been studied in the context of harmonically [35,36] and lattice-confined [37] neutral bosons prepared in a single hyperfine state, our work deals for the first time with fermions that have internal degrees of freedom and multiple collisional parameters. Nonetheless, due to the SU(N ) symmetry of these collisions, we are able to find a simple way to express effective multi-body Hamiltonians and fully characterize their low-lying many-body energy eigenstates. While some past work has detected experimental signatures of multi-body interactions in the form of quantum phase revivals [38], we perform a direct comparison of the many-body interaction energies predicted by our low-energy effective theory with experimental measurements of densitydependent atomic energy spectra performed in ref. [34], similarly to the measurements with bosons performed in ref. [39].
The remainder of this paper is structured as follows. In section II we provide an overview of the relevant one-and two-body physics of ultracold atoms on a lattice, and summarize our main result of effective multi-body interactions between these atoms. In section III we discuss our method for deriving a low-energy effective theory, provide a perturbative expansion for the net effective Hamiltonian, and compute all M -body Hamiltonians through third order in the low-energy effective theory. We then analyze the low-lying excitations of the effective theory in section IV, comparing spectral predictions with experimental measurements, and study the orbital-state dynamics of nuclear spin mixtures interrogated via Rabi spectroscopy.
Finally, we summarize and conclude our findings in section V, and provide some discussion of future outlooks.

II. OVERVIEW
Consider a collection of fermionic atoms with total nuclear spin I and two atomic orbital states loaded into a translationally invariant 3-D optical lattice. While an external trapping potential will generally break discrete translational symmetry of the lattice, any background inhomogeneity can be made negligible by spectroscopically addressing a sufficiently small region of the lattice [34]. Throughout this paper, we work strictly in the deep-lattice regime with negligible tunneling between lattice sites. We also assume a socalled "magic-wavelength" lattice for which both atomic orbital states experience identical lattice potentials [40]. The single-particle Hamiltonian of the atoms can then be written in the form whereĉ inµs is a fermionic operator which annihilates a single atom on lattice site i ∈ Z 3 in motional state n ∈ N 3 0 with nuclear spin µ ∈ {−I, I + 1, · · · , I} (i.e. projected onto a quantization axis) and orbital state s ∈ {g, e}; and E n is the energy of a single atom in motional state n. In a harmonic trap approximation we would have E n = (3/2 + n x + n y + n z ) ω for an on-site angular trap frequency ω (with the reduced Planck constant = 1 throughout this paper), but in general the aharmonicity of the lattice potential will cause a non-negligible shift in motional state energies.
In the absence of hyperfine coupling, as when addressing the spinless 1 S 0 (g) and 3 P 0 (e) orbital states, interactions between any two atoms are governed by their orbital states alone, and are therefore characterized by four scattering lengths a X with X ∈ {gg, eg − , eg + , ee}, where the + (−) superscript denotes symmetrization (anti-symmetrization) of a two-body orbital state under particle exchange. In the low-energy limit, we can write the bare two-body interaction Hamiltonian in the form [41] whereψ µs is a fermionic field operator for atoms with nuclear spin µ and orbital state s; ρ µs ≡ψ † µsψ µs is an atomic density field operator; and the coupling constants G X are defined in terms of the scattering lengths by where m A is the mass of a single atom. Defining for brevity where q, r, s, t ∈ {g, e} are orbital state indices, we can alternately write the bare two-body interaction Hamiltonian in the more compact form For nuclear spins µ, ν, the symbol G qr st gives the coupling constant between the two-atom states (µ, q) + (ν, r) ↔ (µ, s) + (ν, t).
We now expand the field operatorsψ µs in the Wannier basis for a 3-D lattice, such thatψ µs (x) = i,n φ in (x)ĉ inµs with spatial wavefunctions φ in and fermionic annihilation operatorsĉ inµs indexed by lattice sites i and motional states n. Invoking the tight-binding approximation, we assume that the spatial overlap integral in (5) is negligible unless all wavefunctions are localized at the same lattice site; we discuss the breakdown of this approximation and its consequences for our low-energy effective theory in Appendix E. The relevant spatial overlap integral is then which for a lattice with discrete translational invariance is independent of the lattice site i.
The two-body interaction Hamiltonian can be written in terms of this overlap integral as where for brevity we will henceforth suppress the identical site index (i) on all operators, and implicitly sum over all free indices in a summand (i.e. indices which do not have a fixed value). We may also at times suppress motional state indices on the overlap integral K k mn , in which case the suppressed indices are implicitly zero (corresponding to a motional ground state); i.e.
For simplicity, we will also generally work in a gauge for which all two-body overlap integrals are real, such that K k mn = K k mn * = K mn k . The existence of such a gauge is guaranteed by the analytic properties of the Wannier orbitals φ in [42].
We are interested in the ground and low-lying orbital eigenstates of several interacting atoms on a single lattice site, as such states can be readily prepared and coherently addressed in current experiments with ultracold atoms [34]. Although these experiments operate at temperatures sufficiently low for negligible thermal occupation of excited motional states, simply neglecting excited single-particle motional states fails to reproduce the observed spectra of multiply-occupied lattice sites. The reason for this failure lies in the fact that the bare two-body Hamiltonian H int in (7) is not diagonal in the basis of single-particle motional states. Even at zero temperature, therefore, the presence of excited motional states modifies many-body atomic energy spectra when the magnitude of off-diagonal elements of H int approach the scale of low-lying single-particle motional state energy gaps. To address this problem, we develop a low-energy effective theory for ultracold two-level fermionic atoms on a lattice. As our work is motivated by progress in experiments which initialize all atoms in identical (i.e. ground) orbital and motional states, to simplify our theory we assume that any N atoms on a single lattice site occupy N nuclear spin states. Multiple occupation of a single nuclear spin state on the same lattice site is initially forbidden by fermionic statistics, and subsequently cannot be violated in the absence of tunneling and hyperfine coupling.
Our low-energy effective theory exhibits SU(N )-symmetric multi-body interactions, such that the effective interaction Hamiltonian can be written in the form where H M is an M -body Hamiltonian, I is the total nuclear spin of each atom (e.g. I = 9/2 for 87 Sr), and the sum terminates at 2I +1 because this is the largest number of atoms which may initially occupy a single lattice site. Focusing on the low-lying orbital state excitations of this theory, we find that SU(N ) symmetry allows us to express multi-body Hamiltonians in the simple form whereĉ µs ≡ĉ 0,µs is a fermionic annihilation operator for an atom in the motional ground state with nuclear spin µ and orbital state s;n µs ≡ĉ † µsĉ µs is a number operator; and the sum is performed over all choices of nuclear spins µ j with j = 1, 2, · · · , M for which all µ j are distinct, or equivalently all choices of µ j for which the set {µ j } contains M elements, for a total of 2I+1 M × (M !) nuclear spin combinations. The coefficients U X can be expressed in terms of the coupling constants of the bare two-body interaction Hamiltonian (see Appendix F). The key feature of the M -body interactions in (10) is that they take a form identical to the two-body interactions between motional ground-state atoms, but with the addition of M − 2 spectator atoms. This form is a direct consequence of the SU(N ) symmetry of underlying two-body interactions.

III. LOW-ENERGY EFFECTIVE THEORY
As the bare two-body Hamiltonian H int is not diagonal with respect to single-particle motional states, the problem of determining many-body spectra and eigenstates nominally involves all motional degrees of freedom of the atoms involved. If we are interested in addressing internal atomic states at zero temperature, however, in principle we need only to consider the motional ground state associated with each many-body internal atomic state.
We can then ignore all excited many-body motional states, as they will be neither thermally occupied nor externally interrogated. This procedure makes the many-body problem more tractable by greatly reducing the dimensionality of our the relevant Hilbert space.
We denote the space of all internal states of atoms in single-particle, non-interacting (multi-particle, interacting) motional ground states by H single ground (H multi ground ). For sufficiently weak interactions, one can construct a unitary transformation U acting on the full manybody Hilbert space which rotates H multi ground into H single ground [43]. The unitary U can then be used to construct an effective Hamiltonian H eff int on H single ground that reproduces the spectrum of H int on H multi ground , namely H eff int = U H int U † . This method for constructing an effective theory is commonly known as the Schrieffer-Wolff transformation, named after the authors of its celebrated application in relating the Anderson and Kondo models of magnetic impurities in metals [44].
Using the machinery developed in ref. [43] for performing a rotation between low-energy subspaces of a perturbed (i.e. interacting) and unperturbed (i.e. non-interacting) Hamiltonian, we derive an expansion for the effective interaction Hamiltonian H eff int in terms two-body interaction Hamiltonian H int (see Appendix A). This expansion takes the form where H denote an eigenbasis of E 0 with respect to the single-particle Hamiltonian H 0 , and E α denote the motional energy (with respect to H 0 ) of a state |α ∈ E 0 relative to the corresponding motional ground-state energy, we define the operator which sums over projections onto excited states with corresponding energetic suppression factors. The operator I together with the projector P 0 onto H single ground allows us concisely write the first few terms in (11) as where [X, Y ] + ≡ XY + Y X. The naive idea of neglecting all excited states and using H int directly to describe the spectrum of single-particle motional ground states is thus equivalent to truncating our expansion for H eff int at first order. In addition to this first order term, the expansion involves effective corrections to the action of H int on the single-particle motional ground states (i.e. on H single ground ) in the form of higher-order terms with intermediate or virtual occupation of excited states, manifest in I.
Substituting the definition of I into (13)- (14) yields expressions that are highly reminiscent of standard non-degenerate perturbation theory in quantum mechanics, but which nonetheless exhibit crucial differences. The first, and most obvious difference is that these expressions are operator equations, and that the sums over virtual states are performed over a basis for the orthogonal complement of the subspace H single ground , rather than a basis for the orthogonal complement of a single state, as in non-degenerate perturbation theory. Second, the non-degeneracy condition in standard perturbation theory is here elevated to a restriction on the magnitude of the perturbation H int relative to the spectral gap ∆ of the non-interacting Hamiltonian H 0 between H single ground and the excited subspace E 0 . Specifically, the validity of (13)- (14) is conditional only on H int ≤ ∆/2, where X ≡ max |ψ ∈H ψ|X † X|ψ is the operator norm, with no restrictions on spectral gaps or degeneracies within H single ground [43,45]. Finally, the effective theory involves no corrections to the many-body energy eigenstates. Including such corrections would invalidate the effective theory, which preserves the eigenstates of the non-interacting Hamiltonian H 0 on H single ground . As a last comment, we note that our chosen method for constructing an effective Hamiltonian is distinct from adiabatic elimination methods which are commonly used in the atomic physics and quantum optics communities to develop effective theories for e.g. the low-lying levels of a Lambda system [46][47][48][49]. Unlike the perturbative, but exact Schieffer-Wolff transformation, adiabatic elimination methods use approximations which rely on the fast dynamics of excited states. While generally reasonable, these approximations must be made carefully to avoid potential problems with self-consistency (see section 3 of ref. [46]), and yield no obvious or straightforward means to compute effective corrections beyond second order in the couplings between low-and high-energy sectors of a Hilbert space [47][48][49]. While at least one attempt at systematically computing higher-order corrections in the framework of adiabatic elimination has recently been made [49], the resulting expressions do not lend themselves as nicely to analytical treatment, and were in any case found by the authors to be equivalent to a Schrieffer-Wolff transformation.

A. Diagrammatic representation of effective Hamiltonians
The form of the bare two-body interaction Hamiltonian H int in (7) motivates a diagrammatic representation of terms in the effective Hamiltonian H eff in (11), similarly to the diagrams used to represent elements of the scattering matrix in field theory. Effective M -body interaction terms at order p in H int can be represented by directed, weakly connected graphs 1. An example second-order diagram and the corresponding three-body interaction term in H (2) int , with n > 0 andĉ µs ≡ĉ 0,µs . For the sake of presentation, this diagram has colors associated with nuclear spin and orbital states, an arrow on each line to emphasize that they are directed left-to-right, and an explicit coupling constant written next to each vertex; we will generally not include these features, as they are not necessary to uniquely identify the term represented by a diagram. We will also drop explicit appearances of the ground-state projector P 0 in our expressions, with the understanding that the low-energy effective theory implicitly addresses only single-atom motional ground states.  (7) in the definition of a diagram, in all but the two-body case these factors of 1/2 will be cancelled out by corresponding symmetry factors, i.e. the appearance of duplicate diagrams which are equal up to a relabeling of indices (see Appendix B). The explicit signs and factor of 1/2 which appear in the effective Hamiltonians in (13) and (14) are not included in the definition of a diagram, and must be kept track of manually.

B. Effective two-body interactions and renormalization
The effective two-body Hamiltonian in (11) has contributions at all orders in the coupling constants, and can be expanded in the form where the blob on the right schematically represents the net effective two-body interaction.
In general, the sums over excited states in the loop diagrams of (15) will diverge [35], which implies the need to renormalize the coupling constants in our theory. Even if these diagrams do not diverge, we would like to express the effective two-body interactions in terms of net, physically observable two-body coupling constants as on the right-hand side of (15), rather than in terms of long sums at all orders of the bare coupling constants. We therefore renormalize our coupling constants by introducing counter-termsG X into the effective interaction Hamiltonians via G X → G X +G X . The values of these counter-terms are fixed by enforcing that G X is equal to the net two-body interaction strength. Diagrammatically, this renormalization condition is where we represent counter-terms by a crossed dot (i.e. ⊗) at a vertex. To leading order in the coupling constants, it follows that the counter-terms are second order in the coupling constants, i.e.G ∼ G 2 . This finding implies that counter-terms only come into play at nextto-leading order (NLO) in the calculation of any effective M -body interaction Hamiltonian.
By construction, these counter-terms exactly cancel all terms beyond leading order in (15), which implies that the effective two-body interaction Hamiltonian is simply where for consistency with existing literature [35] we define α (1) 2 ≡ K as the overlap integral between two atoms occupying single-particle motional ground states.
Before moving on to consider effective three-body interactions, there are two comments we must make concerning the result in (17). First, this effective two-body interaction takes the same form as the bare two-body interaction in (7), but without excited motional states, and with renormalized coupling constants in place of bare ones. Our choice of renormalization scheme, while convenient for the analytical development of a low-energy effective theory, no longer allows us to use the coupling constants G X as defined by the free-space scattering lengths a X in (3) to compute interaction energies. Instead, we must first compute the effective coupling constants G lattice X (U) in a lattice, which now depend on the lattice depth U, and in turn use the effective coupling constants to compute interaction energies. We discuss the calculation of effective coupling constants in Appendix C.
Second, the overview we provided in section II, and consequently our result in (17), neglects the effect of momentum-dependent two-body scattering. When the effective range of momentum-dependent interactions is comparable to the scattering lengths a X , as is the case for ultracold 87 Sr, these interactions are third order in the coupling constants G X . Just as the O (G 2 ) counter-terms do not affect M -body interactions until next-to-leading order (NLO), the O (G 3 ) momentum-dependent interactions do not come into play until nextnext-leading order (NNLO). Given that we develop our low-energy effective theory through third order in the coupling constants, these interactions will not appear in any of our threeand four-body calculations, but they do need to be considered in the calculation of pair-wise interaction energies. The primary interest of our work, however, concerns effective M -body interactions for M ≥ 3; we therefore defer the calculation of momentum-dependent two-body interactions to Appendix D.

C. Effective three-body interactions at second order
At second order in the coupling constants, the expansion for H (2) int in (13) gives rise to effective three-body interactions. The virtual state of three-body terms in H (2) int cannot have two atoms in excited motional states, as otherwise the second application of H int in H (2) int would have to address both of theses atoms to bring them back down to the ground state, resulting in a two-body process as in the second diagram of (15). All second-order three-body terms must therefore have only one excited atom in the virtual state, and take the form µr νs µr ρt ρt Unlike for the two-body diagram in (17), the explicit factors of 1/2 which appear in the bare two-body Hamiltonian H int in (7) are now cancelled out by symmetry factors which account for duplicate diagrams; this cancellation will generally occur for all connected Mbody diagrams with M > 2 (see Appendix B). The net effective three-body interaction Hamiltonian at second order is then given by the sum over all diagrams of the form in (18), where α 3 ≡ n>0 K 2 n /E n , and the preceding minus sign is as prescribed by H int in (13).

D. Effective three-body interactions at third order
The third-order effective interaction Hamiltonian H int in (14) contains both three-and four-body terms. To compactly enumerate and evaluate all three-body diagrams at third order, we introduce an expanded coupling symbol for more general (µ, q) + (ν, r) ↔ (ρ, s) + (σ, t) coupling induced by terms proportional tô c † µqĉ † νrĉ σtĉρs . The minus sign in (20) accounts for fermionic statistics: if (µ, ν) = (σ, ρ), then we are considering a term of the form At the cost of introducing an additional sum over new nuclear spin indices, the expanded coupling symbol allows us to collect together diagrams which have the same graph topology, but represent different matrix elements of the effective Hamiltonian due to the exchange of nuclear spins at a vertex. The third order three-body diagrams in H µr νs ρt µr νs ρt µr mνs nνs and the mirror image of (23). As prescribed by H int in (14), these diagrams have an associated minus sign if they contain only one excited virtual state, and a factor of 1/2 if they contain a virtual ground state. Remembering that counter-terms are O (G 2 ), there are additionally two third-order three-body diagrams in H and its mirror image, whereG qr st is equal to the counter-term associated with G qr st . To leading order in the coupling constants, the renormalization condition in (16) which means that the sum over all counter-term diagrams in (24) is The net contribution of the these counter-term diagrams together with the two-particle-loop diagrams in (23) is therefore An identical cancellation occurs for the mirror images of the diagrams in (23) and (24), which results in a contribution to the effective Hamiltonian that is equal to (27). The contribution from the three-particle-loop diagrams in (22) is The third-order three-body interaction Hamiltonian is therefore where, letting H 3,2 ≡ |{µ,ν,ρ}|=3 G rs r s G r s r s G s t s t ĉ † µr ĉ † νs ĉ † ρt ĉ ρtĉνsĉµr .
mµq nνr µq νr ρs σt ρs µq νr σt and the mirror image of (34). As we are computing the leading-order contribution to effective four-body interactions, there are no counter-terms contributions. In principle, there is now also the possibility to make the disconnected diagrams of the form , , and .
As prescribed by H int in (14), however, the second and third of these diagrams pick up a factor of −1/2, so the sum over disconnected diagrams vanishes.
Observing that the diagrams of the form in (34) are equal to their mirror image, the positive contributions to the third-order four-body interaction Hamiltonian are while the negative contributions are such that the net third-order four-body interaction Hamiltonian is where wheren µs ≡ĉ † µsĉ µs is a number operator, and the coefficients U X can be determined from the coupling constants G Y and prefactors α As these states are fully symmetric in their orbital degrees of freedom, they are antisymmetric in their nuclear spin degrees of freedom, forming an SU(N ) singlet. Furthermore, the symmetric state is particularly interesting as its orbital degrees of freedom form an Nbody entangled W state, which belongs to a special class of multi-partite entangled states that are robust against disposal or loss of particles. This state thus constitutes an important resource for many quantum information processing and quantum communication tasks [50].
In addition to the states in (43), the multi-body Hamiltonian H M in (42) has an (N − 1)fold degenerate excited-state eigenspace which is asymmetric in the orbital degrees of freedom, spanned by the states for j = 2, · · · , N . If N > 2, the asymmetric states are not separable in their orbital and nuclear spin degrees of freedom. An important feature of the excited states in (43) and (44) is that they are entirely independent of M , which implies that the effect of multi-body interactions is simply to modify the many-body atomic energy spectra without affecting the energy eigenstates. The eigenvalues E  Eigenstate

A. Many-body state spectroscopy
Spectroscopic interrogation is a powerful means to probe the internal structure and dynamics of a system under examination. Consequently, we consider Rabi spectroscopy of the low-lying energy eigenstates in multiply-occupied lattice sites. If we interrogate a lattice site by a laser red-detuned by ∆ from the single-atom orbital state excitation energy, we realize the Hamiltonian where E X is an eigenvalue of the effective interaction Hamiltonian H eff int , P X is a projector onto the corresponding eigenspace, and are single-atom pseudospin operators. The Rabi frequency Ω µ is proportional to the Clebsch-Gordan coefficient I, µ; 1, 0|I, µ ∝ µ for a photon-induced nuclear-spin-conserving orbital state transition of an atom with nuclear spin µ. We therefore define the "bare" Rabi frequency Ω 0 ≡ Ω µ /µ to explicitly factor out dependence on nuclear spins µ.
Consider now a single lattice site in the orbital ground state |N , 0 with nuclear spins N ≡ {µ j } for j = 1, 2, · · · , N . If we red-detune the interrogation laser by δ from a many-body orbital state excitation energy, i.e. set ∆ = ∆ N X − δ for ∆ N X ≡ E N X − E N,0 and X ∈ {+, −}, then in the subspace of the target states {|N , 0 , |N X } the Hamiltonian in (45) becomes where are many-body pseudospin operators, and the Rabi frequencies Ω N X are determined by While the symmetric excited state |N , + is given in (43), at this point we have not explicitly solved for the asymmetric excited state |N , − . The asymmetric state is implicitly defined by (50), and lies somewhere in the span of the N − 1 asymmetric states given in (44).
Determining the symmetric-state Rabi frequency Ω N ,+ is simply a matter of projecting the expression in (50) onto |N , + , which yields whereμ N ≡ µ∈N µ/N is the average nuclear spin µ ∈ N . In order to determine the asymmetric-state Rabi frequency Ω N ,− , we rearrange (50) to find Denoting the standard deviation of nuclear spins µ ∈ N by σ N , normalization of |N , − thus determines the asymmetric-state Rabi frequency which in turn implies that the asymmetric excited state |N , − is   (with s-wave scattering parameters retrieved from ref. [19]), averaged over all nuclear spin combinations of N ∈ {1, · · · , 5} atoms per lattice site for a fixed total atom number. (Right) Experimental measurements of 3 P 0 populations retrieved from ref. [34], with Lorentzian fits to each peak as a visual guide. Resonance peaks are identified by the many-body orbital states which are excited at the peak.
dressed by an external interrogation laser. Figure 3 shows the population of the excited 3 P 0 orbital state when these atoms are interrogated for a time t = π/Ω I by a laser with Rabi frequency Ω I = 50 × 2π Hz (i.e. for individual atoms with nuclear spin µ = I) and detuning ∆ from the single-atom 1 S 0 → 3 P 0 orbital excitation energy. The 3 P 0 population peaks when the laser detuning ∆ is equal to the many-body excitation energy Identifying peaks in excitation spectra such as in figure 3 constitutes a measurement of many-body excitation energies, which was performed in ref. [34] to detect signatures of effective multi-body interactions. Figure 4 shows a comparison between experimental where the Rabi frequencies Ω N ,+ , Ω N ,− and pseudo-spin operators S x N ,± are respectively given in (51), (53), and (49). Denoting the eigenstates of S x N ,± by and defining the identity operator projected to the relevant subspace, we can write the state ρ and from which it follows that the net excited-state population at time t is where X |N |=N ≡ |N |=N X/ Tr P N,0 is an average of X over all choices of N distinct nuclear spins. To understand the periodic collapse and revival of symmetric-state populations in figure   6, we observe from (51) that the symmetric-state phases in (59) and (61) take the form where for fermionic atoms with half-integer nuclear spin, 2µ is always an odd integer, which implies that the sum in (62) is an integer with the same parity (i.e. even/odd) as the occupation number N (i.e. the number of elements in N ). At reduced times τ N = nπ with integer n, therefore, if the occupation number N is even then all phases 2tΩ N ,+ are integer multiples of 2π, which leads to a collapse of the excited-state populations as ρ (±) If the occupation N is odd, meanwhile, then the phases 2tΩ N ,+ are all odd (even) integer multiples of π for odd (even) n. This alignment of phases implies a complete population transfer to the excited state ρ N,− ≡ P N,− / Tr P N,− for odd n, and a collapse back to the orbital ground state ρ N,0 for even n, precisely as observed in figure 6.

V. SUMMARY AND OUTLOOK
Current 3-D optical lattice experiments with fermionic AEAs are capable of operating in the low-temperature, high-density, strongly-interacting limit where inter-atomic interactions set the dominant energy scale governing system dynamics. For AEAs with total nuclear spin I and N = 2I + 1 nuclear spin states, these interactions exhibit an exotic SU(N ) symmetry which is of great interest for near-term quantum simulations of SU(N ) spin models and lattice field theories. Working in the deep-lattice limit, we have derived a low-energy effective theory of such atoms. This theory exhibits emergent multi-body interactions that inherit the SU(N ) symmetry of the atoms' bare, two-body interactions. The SU(N ) symmetry of M -body Hamiltonians in the effective theory has allowed us to express them in a simple form, and fully characterize their corresponding eigenstates and spectra. Capitalizing on the extreme precision of state-of-the-art clock spectroscopy, we have tested spectral predictions of our theory against direct experimental measurements of the many-body 87 Sr excitation spectrum. This comparison shows good agreement between theory and experiment, clearly demonstrating the need to consider multi-body effects for understanding the low-energy physics of high density AEA samples on a 3-D lattice. Finally, we analyzed the many-body orbital-state dynamics of multiply-occupied lattice sites prepared in a nuclear spin mixture and interrogated via Rabi spectroscopy. This analysis is useful for future experimental probes of many-body state structures, as well as for the preparation of long-lived states with multi-partite entanglement (i.e. |N , ± ) which may be used as a resource to perform quantum information processing tasks.
Despite the nominal success of our low-energy effective theory in reproducing experimental observations, there remains room for improvement in the form of controlled, systematic treatment of higher-order and tunneling processes. Nonetheless, our work makes a major step towards the experimental investigation of multi-body SU(N ) physics, providing the necessary framework for future studies going beyond the deep-lattice limit to realize multi-body super-exchange dynamics and orbital SU(N ) quantum magnetism with AEAs.
{|α 0 } for H 0 on G 0 and {|α } for H on G, this transformation is implemented by a unitary U for which |α 0 = U |α and U → 1 as V → 0. The effective Hamiltonian is then simply The prescription in (A1) for constructing an effective Hamiltonian is commonly known as a Schieffer-Wolff transformation [44]. Unitaries U which follow this prescription are not unique, and different choices of U amount to different realizations of the Schieffer-Wolff transformation. In ref. [43], the authors construct the unique operator S which generates a direct or minimal rotation U min = e S between G and G 0 , and use this construction to expand (A1) as a perturbative series in V . The rotation U min is minimal in the sense that it minimizes the distance of candiate unitaries U from the identity 1 with respect to the Euclidian operator norm 1 X E ≡ Tr (X † X). This rotation is determined uniquely by enforcing (i) that the generator S is strictly block-off-diagonal with respect to G 0 and E 0 , (ii) that the norm S E < π/2, and (iii) that the block-off-diagonal parts of (A1) are zero.
To summarize the solution in ref. [43], the effective Hamiltonian H eff induced by a direct rotation can be expanded as where H (p) eff is order p in V . Letting Q 0 ≡ 1 − P 0 denote the projector onto E 0 and X denote any operator on H, we define the superoperators which select out the diagonal (D) and off-diagonal (O) parts of X with respect to G 0 and E 0 , and The first few terms of the expansion in (A2) are then, as derived in ref. [43], Exploiting the fact that in our case ψ|H 0 |ψ = 0 for all |ψ ∈ G 0 , we let B 0 (E 0 ) denote an eigenbasis of H 0 for E 0 and define the operator which sums over projections onto excited states with corresponding energetic suppression factors. We then expand which simplifies the expression for H (2) eff as Working toward a similar expansion for H eff , we compute where [X, Y ] + ≡ XY + Y X. The expressions in (A5), (A9), and (A14) complete the derivation for our expansion of the effective interaction Hamiltonian H eff int in (11) through third order. In the case of ultracold atoms on a lattice, the motional ground-state subspace G 0 actually contains many internal atomic states with different energies. Nonetheless, the total Hilbert space is completely separable into uncoupled subspaces associated with each symmetrized many-body internal atomic state. One can therefore diagonalize the interaction Hamiltonian with respect to these internal states and derive an effective theory within each of the corresponding subspaces, in each case setting the appropriate ground-state energy to zero. This procedure is equivalent to simultaneously calculating the effective Hamiltonian H eff int for all internal states via the prescriptions we have provided, but letting E α denote only the motional excitation energy of states |α 0 ∈ B 0 (E 0 ).

The fact that we include factors of 1/2 from the bare two-body interaction Hamiltonian
H int in the definition of diagrams implies that p-vertex diagrams acquire a factor of 1/2 p . In practice, however, these factors are exactly cancelled out by corresponding symmetry factors in all M -body diagrams with M > 2. As an illustrative example, consider the second-order effective Hamiltonian H (2) int in (13), which expanded in full reads The three-body terms in this Hamiltonian have |{µ, ν, ρ, σ}| = 3 and only one virtually excited atom. The non-vanishing three-body terms must therefore either have ρ ∈ {µ, ν} and contain a factor of the formĉ † Xĉ † σt ĉ σtĉY , or have σ ∈ {µ, ν} with a factor of the form c † ρs ĉ † Xĉ Yĉρs , where the labels X, Y both address whichever nuclear spin (i.e. µ or ν) was excited in the corresponding term. Diagrammatically, we have terms of the form Observing thatĉ † ρs ĉ † Xĉ Yĉρs =ĉ † Xĉ † ρs ĉ ρsĉY , however, it is clear that both of the terms represented in (B2) are equal up to the re-indexing (σ, t, t ) ↔ (ρ, s, s ). There is therefore a symmetry factor of 2 associated with the second vertex of the diagrams in (B2), which cancels out with the explicit factor of 1/2 at that vertex, i.e. the first factor of 1/2 in (B1).
A symmetry factor of essentially identical origin appears at every vertex with an "incoming" These terms are equal up to the re-indexing (ν, r, r , r , n) ↔ (µ, q, q , q , m), which implies that there is a symmetry factor of 2 associated with the first vertex of the diagrams in (B6).
A symmetry factor of identical origin is associated with the first vertex of the diagrams in (B3), and the first two vertices of the diagram in (B4).
The final case we must consider is (ii), which occurs in the first vertex of (B5). In this case, the symmetry factor of 2 which appears in case (i) to account for the possibility of a nuclear spin exchange simply gets "pushed forward" to the vertex at which the two nuclear spins in question part ways, e.g. to account for the two possibilities After summing over all free indices, therefore, all two-body diagrams have a remaining factor of 1/2 from the first vertex. In all connected M -body diagrams with M > 2, meanwhile, every factor of 1/2 can be identified one-to-one with a corresponding symmetry factor of 2.

Appendix C: Effective coupling constants in a lattice
Due to our choice of renormalization scheme in section III B, the interaction energies prescribed by our low-energy effective theory for multiply-occupied lattice sites are not given directly by the coupling constants G X defined by the free-space scattering lengths a X in (3).
Instead, we must first compute effective coupling constants G lattice X (U) in a lattice with depth U, and in turn use the effective coupling constants to compute interaction energies. As the renormalization procedure G X → G lattice X is identical for all coupling constants, we henceforth drop the subscript X ∈ {gg, eg − , eg + , ee} on coupling constants G X in the remainder of this appendix. To further simplify notation, we will also neglect the explicit dependence of parameters on the lattice depth U, which we generally keep fixed.
Proper calculations of the interaction energy of two ultracold fermions in an optical lattice were performed in refs. [51] and [52] using a two-channel model of a Feshbach resonance, yielding prescriptions for computing effective coupling constants in a lattice from free-space interaction parameters. These calculations, however, are both analytically and numerically involved. We therefore instead opt to use a modified version of the considerably simpler single-channel calculation in ref. [53] of the interaction energy of two ultracold atoms in a harmonic trap. Our approach is equivalent to the calculation of Hubbard parameters performed in ref. [54], and has been demonstrated to reproduce correct results in the limit of a deep lattice (compared to the lattice photon recoil energy) and small positive scattering lengths (compared to the effective harmonic oscillator length) [51].
The exact result in Eq. 16 of ref. [53] for the interaction energy of two ultracold atoms in a harmonic oscillator with angular trap frequency ω can be written in the form where G HO is an effective coupling constant in the harmonic trap, φ HO 0 (x) is the corresponding single-particle ground-state wavefunction, and Γ is the gamma function. The expression in (C1) can be solved numerically as is, or expanded about G HO K HO /ω = 0 to get where the first few coefficients are The series in (C2) can in turn be inverted to solve for G HO with an expansion of the form where if we truncate the series in (C2) at n = 2, the first few coefficients of (C4) arẽ The coefficientsc n thus found are consistent with the coefficients c (n+1) 2 reported in table 1 of ref. [35], in which the authors compute the first few terms of the two-body Hamiltonian H 2 directly as expressed in (15) by using a renormalization scheme which subtracts off divergences term by term.
All of the above results are exact for two atoms in a harmonic oscillator interacting via swave scattering. In order to adapt these results for a lattice, we expand the lattice potential about a lattice site centered at x = (0, 0, 0) as where k L = (k x L , k y L , k z L ) is the lattice wavenumber, m A is the atomic mass, and ω eff is an effective angular harmonic trap frequency. We then use ω eff in place of ω in (C1), and use an overlap integral K computed with the ground-state wavefunctions φ 0 in a lattice rather than those in a harmonic oscillator. We retrieve free-space s-wave scattering lengths a free for 87 Sr from ref. [19] to determine the free-space coupling constants G free ≡ (4π/m A ) a free .
This procedure yields an effective coupling constant G lattice given by with a solution where the first few coefficients are provided in (C5).

Appendix D: Momentum-dependent s-wave interactions
In addition to the renormalization of coupling constants discussed in Appendix C, computing two-body interaction energies E (2) N X at third order in the low-energy effective theory requires accounting for the contribution of momentum-dependent s-wave interactions.
At next-to-leading order in the relative momentum k between two atoms, the effective momentum-dependent scattering length a eff is given in terms of the zero-momentum scattering length a by [55,56] which for r eff ak 2 1, implies that a eff ≈ a 1 + 1 2 r eff ak 2 = a + 1 2 r eff a 2 k 2 .
Here r eff is an effective range of O (k 2 ) interactions, determined in atomic units by the scattering length a and van der Waals C 6 coefficient by [57] r eff = 1 3 ξ −2 χ 1 − 2χ + 2χ 2 a, where ξ ≡ Γ (3/4) Γ (1/4) , χ ≡ √ 2 ξ (m A C 6 ) 1/4 a , (D3) states occupying highly excited motional levels, whose spatial wavefunctions can span multiple lattice sites. Nonetheless, we can place upper bounds on the magnitude of inter-site corrections to the effective on-site interaction Hamiltonians by treating tunneling and intersite interactions of virtual excited states perturbatively and assuming no energetic penalty for nearest-neighbor hopping. These bounds can be used to diagnose the breakdown of the on-site effective theory, and signal when a more careful consideration of inter-site effects is necessary to make precise predictions about many-body spectra and dynamics.
If we still assume negligible overlap between single-particle ground-state wavefunctions in different lattice sites but consider nearest-neighbor wavefunction overlaps of states with motional excitations, our one-body and bare two-body Hamiltonians become 3,1 G 2 , γ ∼ γ 3,2 G 2 , γ 3,2 ≡ n,m>0 where we have identified, up to an assignment of coupling constants G, the magnitude of all nonzero matrix elements of the diagrams with respect to an eigenbasis of the on-site in a primitive cubic lattice) is roughly bounded as 2 , (N − 1) γ 3,1 + γ where the factor of b accounts for the multiplicity of neighboring sites; the factor of N 2 accounts for the number of on-site pairs of atoms which are addressed by the diagrams in (E3)-(E5); and the factor of N −1 on γ 3,X accounts for the number of atoms in a neighboring site which are addressed by the corresponding processes. These factors count the number of matrix elements in the Hamiltonian with magnitude ∼ γ (2) X G 2 . The maximization in (E6) is performed because the relevant two-and three-body processes are mutually exclusive, requiring a different number of atoms on neighboring lattice sites. For a conservative bound of |δE N |, the coupling factor G 2 in (E6) can simply be maximized over its allowed values for a given state of atoms on a lattice site, e.g. G 2 g for a state with no orbital excitations, or max G 2 g , G 2 + , G 2 − for a state with one net orbital excitation (in both cases, assuming no orbital excitations in neighboring sites). In the latter case, the bound in (E6) can also be reduced by observing that to conserve energy, it must be the excited atom which moves to a 3,+ = α 4,3 + α 3,− = α 3 − α 4,3 + α In terms of the spatial overlap integrals defined in (6) and (8), the prefactors α 4,2 ≡ m,n>0