# Engineering infinite-range SU($n$) interactions with spin-orbit-coupled
fermions in an optical lattice

Michael A. Perlin<EMAIL_ADDRESS>Diego Barberena  Mikhail Mamaev
Bhuvanesh Sundar  JILA, National Institute of Standards and Technology and
University of Colorado, 440 UCB, Boulder, Colorado 80309, USA  Center for
Theory of Quantum Matter, University of Colorado, Boulder, CO, 80309, USA
Robert J. Lewis-Swan  Homer L. Dodge Department of Physics and Astronomy, The
University of Oklahoma, Norman, Oklahoma 73019, USA  Center for Quantum
Research and Technology, The University of Oklahoma, Norman, Oklahoma 73019,
USA  Ana Maria Rey  JILA, National Institute of Standards and Technology and
University of Colorado, 440 UCB, Boulder, Colorado 80309, USA  Center for
Theory of Quantum Matter, University of Colorado, Boulder, CO, 80309, USA

(22 September 2021)

###### Abstract

We study multilevel fermions in an optical lattice described by the Hubbard
model with on site SU($n$)-symmetric interactions. We show that in an
appropriate parameter regime this system can be mapped onto a spin model with
all-to-all SU($n$)-symmetric couplings. Raman pulses that address internal
spin states modify the atomic dispersion relation and induce spin-orbit
coupling, which can act as a synthetic inhomogeneous magnetic field that
competes with the SU($n$) exchange interactions. We investigate the mean-field
dynamical phase diagram of the resulting model as a function of $n$ and
different initial configurations that are accessible with Raman pulses.
Consistent with previous studies for $n=2$, we find that for some initial
states the spin model exhibits two distinct dynamical phases that obey simple
scaling relations with $n$. Moreover, for $n>2$ we find that dynamical
behavior can be highly sensitive to initial intra-spin coherences. Our
predictions are readily testable in current experiments with ultracold
alkaline-earth(-like) atoms.

##  I Introduction

SU($n$) symmetries play an important role in physics. Underpinning much of
high energy physics, the SU($n$) gauge theory known as Yang-Mills theory is
central to our understanding of the electroweak and strong forces. Extensions
of Yang-Mills and SU($n$) symmetry feature in the most well-studied examples
of holographic duality [1] and the connection between entanglement and gravity
[2] through the anti-de Sitter/conformal field theory (AdS/CFT)
correspondence. In a condensed matter setting, SU(2) appears ubiquitously as a
symmetry of the Hubbard model, with important consequences for the study of
quantum magnetism and high temperature superconductivity [3]. The extension of
SU(2) Hubbard and spin models to SU($n$) has led to predictions of exotic
phases of matter such as valence bond solids [4, 5, 6, 7] and chiral spin
liquids [8, 7, 9, 10], as well as the potential to perform universal
topological quantum computation [11, 12] and other phenomena [13, 14].
Furthermore, disordered SU($n$) spin models have opened analytically tractable
avenues for studying quantum chaos and information scrambling [15].

The tremendous theoretical significance of SU($n$) symmetries makes it all the
more exciting that they appear naturally in experimental atomic, molecular,
and optical (AMO) platforms with exquisite degrees of microscopic control.
This symmetry arises through the independence of atomic orbital and
interaction parameters on the $n$ nuclear spin states of alkaline-earth(-like)
atoms, with e.g. $n=10$ for 87Sr [16, 17, 18, 19]. As a result, AMO
experiments can directly probe the role of SU($n$) interactions in
controllable settings. Recent progress includes studies of the thermodynamic
properties of SU($n$) fermionic gases [20, 21, 22, 23, 24, 25, 26, 27],
SU($n$) Hubbard phases and phase transitions [28, 29, 30], single- [31] and
two-orbital [32, 33, 34, 35] SU($n$) magnetism, and multi-body
SU($n$)-symmetric interactions [36, 37].

In the spirit of quantum simulation, further investigations in controlled
settings will play an important role in understanding the consequence of
SU($n$) symmetries for fundamental questions in physics, as well as their
practical use in technological applications. For example, SU(2)-symmetric spin
interactions can be harnessed to develop quantum sensors that surpass
classical limits on measurement precision [38, 39]. The prospect of similarly
exploiting more general SU($n$) symmetries to achieve a technological
advantage is still an unexplored avenue of research with untapped potential.

In this work, we consider an experimentally relevant and theoretically
tractable regime of the SU($n$) Hubbard model, highlighting differences and
similarities with the more familiar case of SU(2). Working at ultracold
temperatures and unit spatial filling (one atom per lattice site), we begin by
mapping the SU($n$) Hubbard model onto a multilevel spin model with all-to-all
SU($n$)-symmetric interactions in Section II. In Section III we consider the
use of control fields to address nuclear spins, finding a simple three-laser
driving scheme that allows for the preparation of interesting states with
nontrivial intra-spin correlations when $n>2$. We consider the effect of spin-
orbit coupling (SOC) induced by control fields in Section IV, finding in
particular that the weak-SOC limit generally gives rise to a (synthetic)
inhomogeneous magnetic field, extending previously known results to $n>2$ [40,
41, 42, 43, 44, 38]. Finally, we combine these ingredients to examine mean-
field dynamical behaviors of the SU($n$) spin model in Section V, finding
that:  1. long-time-averaged observables obey simple scaling relations with
the spin dimension $n$, exhibiting (for spin-polarized initial states)
dynamical ferromagnetic and dynamical paramagnetic phases, as previously seen
for the case of $n=2$ [45, 46], and  2. for $n>2$ the long-time dynamics can
be highly sensitive to the intra-spin coherences of the initial state.  We
conclude and discuss future directions in Section VI.

##  II From lattice fermions to an SU($n$) spin model

Here we derive a collective SU($n$) spin model for a system of ultracold
alkaline-earth(-like) atoms trapped in an optical lattice. Without external
driving fields, the evolution of such atoms in their electronic ground state
is governed by the single-body kinetic and two-body interaction Hamiltonians

$\displaystyle\hat{H}_{\text{kin}}$
$\displaystyle=-J\sum_{\braket{j,j^{\prime}},\mu}\hat{c}_{j\mu}^{\dagger}\hat{c}_{j^{\prime}\mu}+\text{h.c.},$
(1) $\displaystyle\hat{H}_{\text{int}}$
$\displaystyle=\dfrac{U}{2}\sum_{j,\mu,\nu}\hat{c}_{j\mu}^{\dagger}\hat{c}_{j\mu}\hat{c}_{j\nu}^{\dagger}\hat{c}_{j\nu},$
(2)

where $\braket{j,j^{\prime}}$ denotes neighboring lattice sites $j$ and
$j^{\prime}$; $\mu,\nu\in\left\\{s,s-1,\cdots,-s\right\\}$ index orthogonal
spin states of a spin-$s$ nucleus, with $s=\frac{n-1}{2}$ (e.g.
$s=\frac{9}{2}$ in the case of 87Sr with $10$ nuclear spin states);
$\hat{c}_{j\mu}$ is a fermionic annihilation operator, $J$ is a tunneling
amplitude (for simplicity assumed to be the same in all directions); and $U$
is a two-body on-site interaction energy. In the present work, we neglect
inter-site interactions and interaction-assisted hopping, which is a good
approximation for a sufficiently deep lattice, namely when $J\lesssim
E_{\text{R}}$, where $E_{\text{R}}$ is the atom recoil energy. For simplicity,
we now assume a one-dimensional periodic lattice of $L$ sites, and expand the
on-site fermionic operators in terms of operators addressing (quasi-)momentum
modes $q$ (in units with lattice spacing $a=1$),
$\hat{c}_{j\mu}=\frac{1}{\sqrt{L}}\sum_{q}e^{-\mathrm{i}\mkern 1.0muq\cdot
j}\hat{c}_{q\mu}$, finding that

$\displaystyle\hat{H}_{\text{kin}}$
$\displaystyle=-2J\sum_{q,\mu}\cos\left(q\right)\hat{c}_{q\mu}^{\dagger}\hat{c}_{q\mu},$
(3) $\displaystyle\hat{H}_{\text{int}}$
$\displaystyle=\dfrac{u}{2N}\sum_{k,\ell,p,q,\mu,\nu}\hat{c}_{k\mu}^{\dagger}\hat{c}_{\ell\mu}\hat{c}_{p\nu}^{\dagger}\hat{c}_{q\nu}\times\delta_{k+p,\ell+q},$
(4)

where $N$ is the total number of atoms on the lattice, we define $u\equiv
U\times N/L$ for convenience, $\delta_{k+p,\ell+q}=1$ if $k+p=\ell+q$ and zero
otherwise (enforcing conservation of momentum).

If the interaction energy $U$ is small compared to the single-particle
bandwidth $4J$, then the mode-changing collisions in $\hat{H}_{\text{int}}$
become off-resonant, motivating the frozen-mode approximation
$\left\\{k,p\right\\}=\left\\{\ell,q\right\\}$ (i.e. either $k=\ell$ and
$p=q$, or $k=q$ and $p=\ell$)aaaNote that the frozen-mode approximation
neglects correlated momentum-hopping terms of the form
$\hat{c}_{\pi-p,\mu}^{\dagger}\hat{c}_{\pi-q,\mu}\hat{c}_{p\nu}^{\dagger}\hat{c}_{q\nu}$,
which conserve both momentum and energy. We defer a careful treatment of these
terms to future work, noting only that they vanish on the manifold of
permutationally symmetric spin states with one atom per lattice site, and that
the frozen-mode approximation is benchmarked in Refs. [38, 45] and Appendix
A.. The terms with $k=\ell$ and $p=q$ are
$\frac{u}{2N}\sum\hat{c}_{k\mu}^{\dagger}\hat{c}_{k\mu}\hat{c}_{p\nu}^{\dagger}\hat{c}_{p\nu}=\frac{1}{2}Nu$,
which is a constant energy shift that we can freely neglect. Defining the spin
operators $\hat{s}_{\mu\nu q}\equiv\hat{c}_{q\mu}^{\dagger}\hat{c}_{q\nu}$,
the remaining terms of the kinetic and interaction Hamiltonians are

$\displaystyle\hat{H}_{\text{kin}}$
$\displaystyle=-2J\sum_{q,\mu}\cos\left(q\right)\hat{s}_{\mu\mu q},$ (5)
$\displaystyle\hat{H}_{\text{int}}$
$\displaystyle=-\dfrac{u}{2N}\sum_{p,q,\mu,\nu}\hat{s}_{\mu\nu
p}\hat{s}_{\nu\mu q}.$ (6)

Throughout this work, we will assume that atomic modes are singly-occupied,
e.g. due to the initialization of a spin-polarized state with one atom per
lattice site, in which multiple occupation of an atomic mode is forbidden by
fermionic statistics (Pauli exclusion). In this case we can simply treat our
system as $N$ distinguishable $n$-level quantum spins at “lattice sites”
$p,q$. Note that the “kinetic” terms of this spin model
($\hat{H}_{\text{kin}}$) are proportional to the identity operator,
contributing an overall shift in energy that we can neglect at this point.
Nevertheless, these kinetic terms will become important in the presence of an
external drive, which we discuss in Section IV. The validity of approximating
the Hubbard model in Eqs. (1)–(2) by the spin model in Eqs. (5)–(6) has been
previously benchmarked for SU(2)-symmetric interactions [38, 45], and we
provide additional benchmarking for SU(4) and SU(6) in Appendix A.

To further simplify the interaction Hamiltonian $\hat{H}_{\text{int}}$ and
write it in a form reminiscent of more familiar SU(2) spin models, we now
construct the operator-valued spin matrix

$\displaystyle\bm{\hat{s}}_{q}\equiv\sum_{\mu,\nu}\hat{s}_{\mu\nu
q}\outerproduct{\mu}{\nu},$ (7)

and for any pair of such operator-valued matrices $\bm{\hat{A}},\bm{\hat{B}}$,
we define the inner product

$\displaystyle\bm{\hat{A}}\cdot\bm{\hat{B}}\equiv\sum_{\mu,\nu}\hat{A}_{\mu\nu}^{\dagger}\hat{B}_{\mu\nu}.$
(8)

These definitions allow us to write the spin Hamiltonian in Eq. (6) as

$\displaystyle\hat{H}_{\text{int}}=-\dfrac{u}{2N}\sum_{p,q}\bm{\hat{s}}_{p}\cdot\bm{\hat{s}}_{q}=-\dfrac{u}{2N}\bm{\hat{S}}\cdot\bm{\hat{S}},$
(9)

where $\bm{\hat{S}}\equiv\sum_{q}\bm{\hat{s}}_{q}$ is a collective spin
matrix, analogous to the collective spin vector
$\vec{S}=(\hat{S}_{\text{x}},\hat{S}_{\text{y}},\hat{S}_{\text{z}})$ in the
case of SU(2) [38], with
$\frac{1}{2}\bm{\hat{S}}\cdot\bm{\hat{S}}\simeq\vec{S}\cdot\vec{S}=\hat{S}_{\text{x}}^{2}+\hat{S}_{\text{y}}^{2}+\hat{S}_{\text{z}}^{2}$
when $n=2$ (here $\simeq$ denotes equality up to identity terms).

Figure 1:  (a) Ultracold atoms on a lattice of $L$ sites tunnel between
neighboring lattice sites at a rate $J$, and locally repel each other with
interaction energy $U$. (b,c) When the interaction energy $U$ is small
compared to the single-particle bandwidth $4J$, the frozen-mode approximation
enables the interaction Hamiltonian to be written as a spin model consisting
of exchange terms $\bm{\hat{s}}_{p}\cdot\bm{\hat{s}}_{q}$, which swap the
states of two spins pinned to modes $p,q$. (d) Interactions open an energy gap
$u=U\times N/L$ between the manifold of permutationally symmetric states of
$N$ spins, and the orthogonal complement of states that break spin-permutation
symmetry.

We now discuss the spin Hamiltonian $\hat{H}_{\text{int}}$ in Eq. (9). The
operator $\bm{\hat{s}}_{p}\cdot\bm{\hat{s}}_{q}$ simply swaps the nuclear spin
states of two atoms pinned to modes $p,q$. The term
$-\bm{\hat{s}}_{p}\cdot\bm{\hat{s}}_{q}$ thereby assigns a definite energy of
$-1$ ($+1$) to a pair of spins that are symmetric (anti-symmetric) under
exchange. In this sense, $\bm{\hat{s}}_{p}\cdot\bm{\hat{s}}_{q}$ is analogous
to the enforcement of SU(2) spin alignment by ferromagnetic interactions,
which similarly assigns different energies to the anti-symmetric spin-0
singlet $\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow}$ and the symmetric
spin-1 triplets
$\left\\{\ket{\uparrow\uparrow},\ket{\downarrow\downarrow},\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow}\right\\}$.
By summing over all pair-wise exchange terms
$\bm{\hat{s}}_{p}\cdot\bm{\hat{s}}_{q}$, the interaction Hamiltonian
$\hat{H}_{\text{int}}$ energetically enforces a permutational symmetry among
all spins, opening an energy gap $u$ between the manifold of all
permutationally symmetric (PS) states and the orthogonal complement of excited
(e.g. spin-wave) states that break permutational symmetry. See Figure 1 for a
summary of this section thus far.

In the case of SU(2), the PS manifold is precisely the Dicke manifold of
collective states $\ket{m_{\text{z}}}$ with total spin $S=\frac{N}{2}$ and
definite spin projection $m_{\text{z}}\in\left\\{S,S-1,\cdots,-S\right\\}$
onto a fixed quantization axis. Equivalently, Dicke states
$\ket{m_{\text{z}}}=\ket{m_{\uparrow},m_{\downarrow}}$ can be labeled by a
definite number of spins $m_{\uparrow}=S+m_{\text{z}}$
($m_{\downarrow}=S-m_{\text{z}}$) pointing up (down) along the spin
quantization axis, with $m_{\uparrow}+m_{\downarrow}=N$. In the general case
of SU($n$), the PS manifold is similarly spanned by states
$\ket{m_{s},m_{s-1},\cdots,m_{-s}}$ with a definite number $m_{\mu}$ of spins
in state $\mu$, and $\sum_{\mu}m_{\mu}=N$. The dimension of the PS manifold is
equal to the number of ways of assigning $N$ identical spins to $n$ distinct
internal states, or ${N+n-1\choose n-1}\sim N^{n-1}$.

External fields or additional interactions that respect permutational symmetry
can induce nontrivial dynamics within the PS manifold. Moreover, additional
terms that explicitly break permutational symmetry can nevertheless lead to
interesting dynamics that can be captured within the PS manifold
perturbatively, as long as the coupling to non-PS states is weak compared $u$
(see Appendix B) [47]. This perturbative regime is thereby efficiently
simulable, as the PS manifold has dimension $\sim N^{n-1}$ (as compared to
$n^{N}$ for the entire spin Hilbert space). Simulating dynamics within the PS
manifold requires calculating matrix elements
$\braket{\ell}{\hat{\mathcal{O}}}{m}$ of spin operators $\hat{\mathcal{O}}$
with respect to PS states $\ket{\ell},\ket{m}$; we discuss this calculation in
Appendix C.

Figure 2:  Whereas the state of a two-level spin (qubit) can be represented by
a point on (or inside) the Bloch sphere, the state of an $n$-level spin is
more generally represented by a probability distribution on the Bloch sphere.
The distribution shown for $n=10$ corresponds to a Haar-random pure state.

Finally, we take a moment to discuss individual $n$-level spins. The state of
a two-level spin, or a qubit, is commonly represented by a point on (or
within) the Bloch sphere. More generally, the state $\ket{\psi}$ of an
$n$-level spin can be represented by a quasi-probability distribution
$Q_{\psi}$ on the Bloch sphere (commonly known as the Husimi-$Q$ function,
e.g. in the spin-squeezing community [48]). The value
$Q_{\psi}\left(\bm{v}\right)$ at a point $\bm{v}$ on the sphere is equal to
the overlap of $\ket{\psi}$ with a pure state $\ket{\bm{v}}$ that is maximally
polarized in the direction of $\bm{v}$:
$Q_{\psi}\left(\bm{v}\right)\equiv\absolutevalue{\braket{\bm{v}}{\psi}}^{2}$
(see Figure 2). In the case of a mixed state $\hat{\rho}$, this distribution
is defined by
$Q_{\hat{\rho}}\left(\bm{v}\right)\equiv\braket{\bm{v}}{\hat{\rho}}{\bm{v}}$.
Closely related spherical representations of multilevel spin states and
operators are discussed in Refs. [49, 50]. In practice, it is conceptually
useful to identify the Hilbert space of a single $n$-level spin with the Dicke
manifold of $n-1$ spin-$\frac{1}{2}$ particles.

##  III External control fields

We now consider the addition of external control fields to address atoms’
internal spin states, which will determine the observables we can access and
initial states we can prepare. Specifically, we consider off-resonantly
addressing an electronic $\ket{\text{g}}\to\ket{\text{e}}$ transition of the
atoms, and then perturbatively eliminating electronic $\ket{\text{e}}$
excitations to arrive at an effective ground-state Hamiltonian addressing
nuclear spins. For simplicity, we will assume that the total spin $s$ of the
ground- and excited-state (hyperfine) manifolds are the same, as e.g. with the
${{}^{1}}\text{S}_{0}\to{{}^{3}}\text{P}_{0}$ transition of alkaline-earth-
like atoms (AEAs). However, the results of this section (namely the general
form of effective nuclear spin Hamiltonians, as well as the corresponding set
of accessible observables and initial states) are the same for transitions
that take $s\to s\pm 1$, so in practice one is free to address the hyperfine
manifolds of the ${{}^{1}}\text{S}_{0}\to{{}^{3}}\text{P}_{1}$ transition of
AEAs.

Figure 3:  Sketch of the three-laser drive used to address nuclear spins on a
one-dimensional lattice. Two counter-propagating lasers with right-circular
polarization and amplitudes $\Omega_{\pm}$ point at an angle $\theta$ to the
lattice axis. A third, linearly polarized laser with amplitude $\Omega_{0}$
points in a direction orthogonal to both the lattice and the other driving
lasers. Absorbing a photon from the laser with amplitude $\Omega_{m}$ induces
a transition $(\text{g},\mu)\to(\text{e},\mu+m)$ for the (electronic, nuclear
spin) state of an atom, where nuclear spin is quantized along the $z$ axis.

We consider a specific three-laser driving scheme with a geometry sketched in
Figure 3. Here the lattice lies in the $y$-$z$ plane at an angle $\theta$ to
the $z$ axis, oriented along $\bm{\ell}=(0,\sin\theta,\cos\theta)$. We set the
spin quantization axis along $z$. The laser setup consists of  1. two counter-
propagating right-circularly polarized lasers with drive amplitudes
$\Omega_{\pm}$ and wavevectors $\kappa\bm{v}_{\pm}$, propagating in opposite
directions along the $z$ axis, $\bm{v}_{\pm}=(0,0,\pm 1)$, and  2. a third
laser linearly polarized along $z$, with drive amplitude $\Omega_{0}$ and
wavevector $\kappa\bm{v}_{0}$, propagating along the $x$ axis
$\bm{v}_{0}=(1,0,0)$.  All driving lasers are detuned by $\Delta$ below an
electronic transition. The full Hamiltonian for this three-laser drive can be
written as

$\displaystyle\hat{H}_{\text{3LD}}^{\text{full}}=\sum_{j,m}\Omega_{m}\left(e^{-\mathrm{i}\mkern
1.0mum\phi
j}\hat{s}_{mj}\otimes\outerproduct{\text{e}}{\text{g}}_{j}+\text{h.c.}\right)+\Delta\hat{N}_{\text{e}},$
(10)

where $m\in\left\\{+1,0,-1\right\\}$ indexes the laser pointing along
$\bm{v}_{m}$; the SOC angle
$\phi\equiv\kappa\bm{v}_{+}\cdot\bm{\ell}=\kappa\cos\theta$ (in units with
lattice spacing $a=1$); $\hat{s}_{\text{z},j},\hat{s}_{+,j}\hat{s}_{-,j}$ are
standard axial, spin-raising, and spin-lowering operators for the spin at
lattice site $j$; $\hat{s}_{0,j}\equiv\hat{s}_{\text{z},j}$ for shorthand;
$\ket{\text{g}}_{j}$ and $\ket{\text{e}}_{j}$ respectively denote the ground
and excited electronic states of atom $j$; and
$\hat{N}_{\text{e}}=\mathds{1}\otimes\sum_{j}\outerproduct{\text{e}}{\text{e}}_{j}$
counts the number of excited atoms (with $\mathds{1}$ the identity operator on
all spin degrees of freedom).

In the far-detuned limit
$\absolutevalue{\Delta}\gg\absolutevalue{\Omega_{m}}$, a second-order
perturbative treatment of electronic excitations ($\ket{\text{e}}$) yields an
effective drive Hamiltonian that only addresses ground-state nuclear spins.
After additionally making the gauge transformation $\hat{s}_{mj}\to
e^{\mathrm{i}\mkern 1.0mum\phi j}\hat{s}_{mj}$ (equivalently
$\hat{c}_{j\mu}^{\dagger}\to e^{\mathrm{i}\mkern 1.0mu\phi\mu
j}\hat{c}_{j\mu}^{\dagger}$), the drive Hamiltonian then becomes

$\displaystyle\hat{H}_{\text{3LD}}=\sum_{j}\hat{H}_{\text{3LD},j}^{\text{single}},$
(11)

where $\hat{H}_{\text{3LD},j}^{\text{single}}$ denotes the action of
$\hat{H}_{\text{3LD}}^{\text{single}}$ on spin $j$:

$\hat{H}_{\text{3LD}}^{\text{single}}=\tilde{\Omega}_{+}\tilde{\Omega}_{-}\hat{s}_{\text{z}}+\tilde{\Omega}_{0}\tilde{\Omega}_{-}\hat{s}_{\text{x}}+\tilde{\Omega}_{0}\tilde{\Omega}_{+}(\hat{s}_{\text{z}}\hat{s}_{\text{x}}+\hat{s}_{\text{x}}\hat{s}_{\text{z}})\\\
-\tilde{\Omega}_{0}^{2}\hat{s}_{\text{z}}^{2}-\tilde{\Omega}_{+}^{2}\hat{s}_{\text{x}}^{2}-\tilde{\Omega}_{-}^{2}\hat{s}_{\text{y}}^{2},$
(12)

with

$\displaystyle\tilde{\Omega}_{0}\equiv-\dfrac{\Omega_{0}}{\sqrt{\Delta}},$
$\displaystyle\tilde{\Omega}_{\pm}\equiv\dfrac{\Omega_{+}\pm\Omega_{-}}{\sqrt{\Delta}},$
(13)

where we have made the simplifying assumption that all drive amplitudes are
real to arrive at the form of $\hat{H}_{\text{3LD}}^{\text{single}}$ in Eq.
(12). We relax the assumption of real drive amplitudes in Appendix D.

Table 1:  Drive Hamiltonians (left column) that can be implemented with different amplitude-matching conditions (right three columns), some of which are specified by an arbitrary sign $\sigma\in\left\\{+1,-1\right\\}$. The drives shown here are equal to that of Eq. (12) up to a possible energy shift of $\hat{s}_{\text{x}}^{2}+\hat{s}_{\text{y}}^{2}+\hat{s}_{\text{z}}^{2}=s(s+1)$, and come in mutually commuting pairs: a drive with $\absolutevalue*{\Omega_{m}}=1$ and $\Omega_{n}=0$ for both $n\neq m$ commutes with the drive in which $\Omega_{m}=0$ and both $\absolutevalue*{\Omega_{n}}=1$.  $\hat{H}_{\text{drive}}^{\text{single}}$ | $\tilde{\Omega}_{0}$ | $\tilde{\Omega}_{+}$ | $\tilde{\Omega}_{-}$
---|---|---|---
$-\hat{s}_{\text{z}}^{2}$ | 1 | 0 | 0
$-\hat{s}_{\text{x}}^{2}$ | 0 | 1 | 0
$-\hat{s}_{\text{y}}^{2}$ | 0 | 0 | 1
$\sigma\hat{s}_{\text{z}}+\hat{s}_{\text{z}}^{2}$ | 0 | 1 | $\sigma$
$\sigma\hat{s}_{\text{x}}+\hat{s}_{\text{x}}^{2}$ | 1 | 0 | $\sigma$
$\sigma\left(\hat{s}_{\text{z}}\hat{s}_{\text{x}}+\hat{s}_{\text{x}}\hat{s}_{\text{z}}\right)+\hat{s}_{\text{y}}^{2}$ | 1 | $\sigma$ | 0
$\pm\hat{s}_{\text{z}}\pm\sigma\hat{s}_{\text{x}}+\sigma(\hat{s}_{\text{z}}\hat{s}_{\text{x}}+\hat{s}_{\text{x}}\hat{s}_{\text{z}})$ | 1 | $\sigma$ | $\pm\sigma$

There are three important observations to make about Eqs. (11) and (12).
First, the fact that $\hat{H}_{\text{3LD}}$ acts identically on all spins
means we can freely replace the site index $j$ with a momentum index $q$ (as
can be verified by substituting
$\hat{c}_{j\mu}=\frac{1}{\sqrt{L}}\sum_{k}e^{-\mathrm{i}\mkern 1.0muq\cdot
j}\hat{c}_{q\mu}$), which is important to ensure that this drive addresses the
same spin degrees of freedom as the spin Hamiltonians previously considered in
Section II. Second, each of
$\tilde{\Omega}_{0},\tilde{\Omega}_{+},\tilde{\Omega}_{-}$ can be tuned
independently by changing the amplitudes of the driving lasers; some
particular Hamiltonians for specific values of these amplitudes are shown in
Table 1. Third, due to the appearance of mutually commuting pairs of
Hamiltonians in Table 1, specifically $-\hat{s}_{\alpha}^{2}$ and
$\pm\hat{s}_{\alpha}+\hat{s}_{\alpha}^{2}$ for
$\alpha\in\left\\{\text{z},\text{x}\right\\}$, the three-laser drive admits
pulse sequences that exactly implement arbitrary SU(2) (spatial) rotations of
the form $e^{-\mathrm{i}\mkern 1.0mu\chi\vec{n}\cdot\vec{s}}$, where $\chi$ is
a rotation angle, $\vec{n}$ is a rotation axis, and
$\vec{s}\equiv(\hat{s}_{\text{x}},\hat{s}_{\text{y}},\hat{s}_{\text{z}})$. The
capability to perform arbitrary spatial rotations, together with the
capability to measure the number of atoms with spin projection $\mu$ onto a
fixed quantization axis, $\braket{\hat{S}_{\mu\mu}}$ (where
$\hat{S}_{\mu\nu}=\sum_{j}\hat{s}_{\mu\nu j}$), implies the capability to
reconstruct all components of the mean collective spin matrix
$\braket{\bm{\hat{S}}}=\sum_{\mu\nu}\braket{\hat{S}_{\mu\nu}}\outerproduct{\mu}{\nu}$
via spin qudit tomography [51, 52]. Moreover, we expect that advanced quantum
control techniques (similar to those of Refs. [53, 54]) can be used to
implement arbitrary SU($n$) rotations by designing suitable time-dependent
drive amplitudes.

If the excited-state manifold $\ket{\text{e}}$ has total spin $s\pm 1$, the
effective ground-state Hamiltonians in Eq. (12) and Table 1 remain almost
identical, but with some additional $n$-dependent factors that do not affect
the general results and discussions above. These results still hold if (for
example) all excited hyperfine manifolds of an electronic
${{}^{1}}\text{S}_{0}\to{{}^{3}}\text{P}_{1}$ transition (with total spins
$s+1,s,s-1$) are addressed simultaneously. See Appendix D for additional
details.

Finally, we comment on the preparation of initial states. Initial states are
nominally prepared in the “lab frame”, and must be transformed according to
the gauge transformation $\hat{c}_{j\mu}^{\dagger}\to e^{\mathrm{i}\mkern
1.0mu\phi\mu j}\hat{c}_{j\mu}^{\dagger}$ prior to evolution under the three-
laser drive $\hat{H}_{\text{3LD}}$ in Eq. (11), which is expressed in the
“gauge frame”. We assume the capability to prepare an initial state in which
all spins are maximally polarized along the $z$ axis, i.e.
$\ket{\text{z}}^{\otimes N}=\ket{s}^{\otimes N}$, which is unaffected by the
gauge transformation (up to a global phase). The three-laser then allows us to
rotate this state into one that is polarized along any spatial axis (in the
gauge frame). In addition, when $n>2$ the three-laser drive allows us to
prepare product states with nontrivial intra-spin correlations. For example,
when $n$ is even we can prepare an $N$-fold product of the “kitten” state

$\displaystyle e^{-\mathrm{i}\mkern
1.0mu\frac{\pi}{2}\left(\hat{s}_{\text{y}}+\hat{s}_{\text{y}}^{2}\right)}\ket{s}\stackrel{{\scriptstyle
n\,\text{even}}}{{\propto}}\ket{s}+\ket{-s}.$ (14)

This state has a vanishing mean spin vector,
$\braket{\hat{s}_{\text{x}}}=\braket{\hat{s}_{\text{y}}}=\braket{\hat{s}_{\text{z}}}=0$,
but variances
$\braket{\hat{s}_{\text{x}}^{2}}=\braket{\hat{s}_{\text{y}}^{2}}=s/2$ and
$\braket{\hat{s}_{\text{z}}^{2}}=s^{2}$.

##  IV Spin-orbit coupling

We now consider the effect of spin-orbit coupling (SOC) induced by the control
fields in Section III. Before discussing SOC for $n$-level fermions, we
briefly review the well-studied case of two-level SOC with a one-dimensional
lattice [41, 43, 44, 38]. In this case, SOC is induced by an external driving
field that imprints a phase $e^{-\mathrm{i}\mkern 1.0mu\phi j}$ on lattice
site $j$, or equivalently imparts a momentum kick $q\to q+\phi$, upon the
absorption of a photonbbbIn order for the drive Hamiltonian
$\hat{H}_{\text{drive}}^{(\phi)}$ to be well-defined, $\phi$ should be
commensurate with the lattice, e.g. $\phi\in\mathbb{Z}\times 2\pi/L$ on a one-
dimensional lattice of $L$ sites.:

$\displaystyle\hat{H}_{\text{drive}}^{(\phi)}=\dfrac{\Omega}{2}\sum_{q}\hat{c}_{q+\phi,\uparrow}^{\dagger}\hat{c}_{q,\downarrow}+\text{h.c.}.$
(15)

Identifying a numerical spin index $\mu=+\frac{1}{2}$ ($-\frac{1}{2}$) with
the state $\uparrow$ ($\downarrow$), this drive Hamiltonian can be
diagonalized in its momentum index $q$ by the gauge transformation
$\hat{c}_{q\mu}^{\dagger}\to\hat{c}_{q-\mu\phi,\mu}^{\dagger}$ (equivalently
$\hat{c}_{j\mu}^{\dagger}\to e^{\mathrm{i}\mkern 1.0mu\phi\mu
j}\hat{c}_{j\mu}^{\dagger}$), which takes

$\displaystyle\hat{H}_{\text{drive}}^{(\phi)}\to\hat{H}_{\text{drive}}\equiv\Omega\hat{S}_{\text{x}},$
$\displaystyle\hat{S}_{\text{x}}\equiv\sum_{q}\hat{s}_{\text{x},q},$ (16)

where
$\hat{s}_{\text{x},q}=\frac{1}{2}\hat{c}_{q,\uparrow}^{\dagger}\hat{c}_{q,\downarrow}+\text{h.c.}$
for two-level spins.

The two-level SOC drive in Eq. (15) has been implemented with an external
laser that couples the two electronic states of nuclear-spin-polarized atoms,
with $\downarrow$ ($\uparrow$) indexing the ground (excited) electronic state
[41, 42, 43, 44, 38]. In contrast, the drive we considered in Section III
addresses electronic excitations off-resonantly, inducing an effective
Hamiltonian in the ground-state hyperfine manifold with spin projections
$\mu\in\left\\{s,s-1,\cdots,-s\right\\}$ (a similar scheme was used to study
SOC in a subspace of the ground-state manifold in Ref. [40]). Nonetheless,
both the two-level drive in Eq. (15) and the $n$-level drive in Eq. (11)
become homogeneous (i.e. independent of the spatial mode index $j$ or $q$) and
independent of the SOC angle $\phi$ after the same spin-symmetric gauge
transformationcccThe “asymmetric” gauge transformation
$(\hat{c}_{j,\uparrow}^{\dagger},\hat{c}_{j,\downarrow}^{\dagger})\to(e^{\mathrm{i}\mkern
1.0mu\phi j}\hat{c}_{j,\uparrow}^{\dagger},\hat{c}_{j,\downarrow}^{\dagger})$,
sometimes performed in the two-state SOC literature, does not generalize as
nicely to $n>2$. $\hat{c}_{j\mu}^{\dagger}\to e^{\mathrm{i}\mkern 1.0mu\phi\mu
j}\hat{c}_{j\mu}^{\dagger}$.

Of course, spin-orbit coupling cannot be “gauged away” entirely. Making a
gauge transformation to simplify the drive comes at the cost of making the
kinetic energy in Eq. (5) spin-dependent, taking

$\displaystyle\hat{H}_{\text{kin}}\to\hat{H}_{\text{kin}}^{(\phi)}\equiv-2J\sum_{q}\cos\left(q+\mu\phi\right)\hat{s}_{\mu\mu
q},$ (17)

as visualized in Figure 4. To better interpret this Hamiltonian, we can write
it in the form

$\displaystyle\hat{H}_{\text{kin}}^{(\phi)}=-2J\sum_{q}\left[\cos\left(q\right)\hat{w}_{+,q}^{(\phi)}-\sin\left(q\right)\hat{w}_{-,q}^{(\phi)}\right],$
(18)

where

$\displaystyle\hat{w}_{+,q}^{(\phi)}$
$\displaystyle\equiv\sum_{\mu}\cos\left(\mu\phi\right)\hat{s}_{\mu\mu q},$
(19) $\displaystyle\hat{w}_{-,q}^{(\phi)}$
$\displaystyle\equiv\sum_{\mu}\sin\left(\mu\phi\right)\hat{s}_{\mu\mu q}.$
(20)

For two-level spins with $\mu=\pm\frac{1}{2}$, $\hat{w}_{+,q}^{(\phi)}$ is
proportional to the identity operator and
$\hat{w}_{-,q}^{(\phi)}=2\sin\left(\phi/2\right)\hat{s}_{\text{z},q}$, so the
kinetic Hamiltonian in the gauge frame describes a (synthetic) inhomogeneous
magnetic field:

$\displaystyle\left.\hat{H}_{\text{kin}}^{(\phi)}\right|_{n=2}=4J\sin\left(\phi/2\right)\sum_{q}\sin\left(q\right)\hat{s}_{\text{z},q}.$
(21)

When $n>2$, an inhomogeneous magnetic field is likewise recovered in the weak
SOC limit $s\phi\ll 1$, in which case

$\displaystyle\left.\hat{H}_{\text{kin}}^{(\phi)}\right|_{s\phi\ll
1}=2J\phi\sum_{q}\sin\left(q\right)\hat{s}_{\text{z},q}+O\left((s\phi)^{2}\right).$
(22)

For larger $\phi$, this Hamiltonian acquires terms with higher powers of
$\hat{s}_{\text{z},q}$, up to $\hat{s}_{\text{z},q}^{n-1}$.

Figure 4:  Spin-orbit coupling for 2-level (a,c) and 4-level (b,d) spins.
Colors indicate different spin projections $\mu$. In the “lab frame” (a,b),
kinetic energy is insensitive to spin, but a spin transition $\mu\to\mu+1$ is
accompanied by a momentum kick $q\to q+\phi$ from the drive. Changing into the
“gauge frame” (c,d), essentially by shifting the momentum label $q$ for each
spin state $\mu$, makes the drive diagonal in the momentum index, but comes at
the cost of making kinetic energy spin-dependent.

Finally, the gauge transformation
$\hat{c}_{q\mu}^{\dagger}\to\hat{c}_{q-\mu\phi,\mu}^{\dagger}$ also transforms
the interaction Hamiltonian. Applying this transformation to Eq. (4) and
keeping only terms that respect coherences that can be imposed on initial
states by the laser drive in Section III (applied to an initially spin-down-
polarized state) again results in an effective spin model. For sufficiently
weak SOC ($s\phi\to 0$) this spin model is still well-approximated by
$\hat{H}_{\text{int}}$ in Eqs. (6) and (9). The validity of this approximation
has been previously benchmarked for SU(2)-symmetric interactions [38, 45], and
we provide additional benchmarking for SU(4) and SU(6) in Appendix A (which
finds that the spin model works well even for large $\phi$). To ensure that
$\hat{H}_{\text{kin}}^{(\phi)}$ does not become trivial as $\phi\to 0$, we can
keep $J\phi/u$ constant, either by increasing $J/U$ or decreasing $N/L$.
Altogether, the interacting spin Hamiltonian in the gauge frame becomes

$\displaystyle\hat{H}_{\text{spin}}=-\dfrac{u}{2N}\bm{\hat{S}}\cdot\bm{\hat{S}}+2J\phi\sum_{q}\sin\left(q\right)\hat{s}_{\text{z},q},$
(23)

consisting of a spin-locking $\bm{\hat{S}}\cdot\bm{\hat{S}}$ term that
energetically favors permutational symmetry, and an inhomogeneous magnetic
field that causes inter-spin dephasing.

##  V Mean-field theory and dynamical phases

We now study the dynamical behavior of the SOC spin Hamiltonian
$H_{\text{spin}}$ in Eq. (23), and henceforth work exclusively in the “gauge
frame” of $\hat{H}_{\text{spin}}$ and the three-laser drive
$\hat{H}_{\text{3LD}}$ in Eq. (11). We use a Ramsey-like setup wherein we
prepare an initial state with the three-laser drive (using fast pule
sequences), then let the state evolve freely for some time under
$\hat{H}_{\text{spin}}$, and finally apply again the three-laser drive to map
observables of interest onto spin projection measurements (e.g. with spin
qudit tomography [51, 52]). At the mean-field (MF) level, the undriven spin
Hamiltonian (neglecting constant energy shifts) becomes

$\displaystyle\hat{H}_{\text{MF}}=u\sum_{q}\left[-\braket{\bar{\bm{s}}}\cdot\bm{\hat{s}}_{q}+h\sin\left(q\right)\hat{s}_{\text{z},q}\right].$
(24)

where $\bar{\bm{s}}\equiv\frac{1}{N}\sum_{q}\bm{\hat{s}}_{q}$ is the average
spin matrix, and $h\equiv 2J\phi/u$ is a dimensionless strength of the
inhomogeneous magnetic field. We assume that all momenta
$q\in\mathbb{Z}_{N}\times 2\pi/N$ are occupied. Fixing the atom number $N$,
the spin Hamiltonian has one free parameter, $h$, which determines the
relative strength of the single-particle and interaction terms. One should
therefore expect distinct dynamical behaviors when $h\ll 1$, in which case
strong spin-locking interactions should give rise to a long-range ordered
phase, as opposed to $h\gg 1$, in which case long-range order should be
destroyed by the strong inhomogeneous magnetic field [45].

To investigate these behaviors quantitatively, we examine time-averaged
observables of the form

$\displaystyle\langle\\!\langle\hat{\mathcal{O}}\rangle\\!\rangle_{\text{MF}}=\lim_{T\to\infty}\dfrac{1}{T}\int_{0}^{T}\text{d}t\braket{\hat{\mathcal{O}}\left(t\right)}_{\text{MF}},$
(25)

where $\braket{\hat{\mathcal{O}}\left(t\right)}_{\text{MF}}$ is the mean-field
value of observable $\hat{\mathcal{O}}$ at time $t$. Specifically, we consider
the time-averaged magnetization

$\displaystyle\sigma_{\text{MF}}\equiv\absolutevalue{\langle\\!\langle\vec{\sigma}\rangle\\!\rangle_{\text{MF}}},$
$\displaystyle\vec{\sigma}\equiv\dfrac{1}{Ns}\times\vec{S},$ (26)

where
$\vec{S}\equiv(\hat{S}_{\text{x}},\hat{S}_{\text{y}},\hat{S}_{\text{z}})$ with
$\hat{S}_{\alpha}\equiv\sum_{q}\hat{s}_{\alpha,q}$, and the time-averaged
(dimensionless) interaction energy

$\displaystyle\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}=\dfrac{1}{N^{2}}\times\langle\\!\langle\bm{\hat{S}}\cdot\bm{\hat{S}}\rangle\\!\rangle_{\text{MF}}.$
(27)

By design, these non-negative quantities are normalized to lie on the interval
$[0,1]$, independent of the system size $N$ or spin dimension $n$. In the
remainder of this section we will assume that $n$ is even, both for the sake
of experimental relevance (most relevant atomic nuclei are fermionic) and to
avoid complications from parity effects.

Our numerical simulations of mean-field dynamics are performed with a
Schwinger boson decomposition of spin operators: $\hat{s}_{\mu\nu
q}=\hat{b}_{\mu q}^{\dagger}\hat{b}_{\nu q}$. This decomposition requires no
approximations, and reduces the number of variables to keep track of by a
factor of $\sim n$. See Appendices E and F for additional details about our
numerical simulations and the Schwinger boson equations of motion.

###  V.1 Initial spin-polarized state

Figure 5:  Time-averaged magnetization $\sigma_{\text{MF}}$ and interaction
energy
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
for different spin dimensions $n$ (indicated in the legend) as determined by
mean-field simulations of $N=100$ spins initially in the x-polarized state
$\ket{\text{X}}$ for a time $T=10^{5}/u$. Insets show same data after
rescaling $h\to h\times(n/2)^{1/3}$, and transforming vertical axes according
to Eq. (29).

Figure 5 shows the time-averages of the magnetization $\sigma_{\text{MF}}$ and
interaction energy
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
as computed by mean-field simulations of $N=100$ spins initially in the
x-polarized state $\ket{\text{X}}\equiv\ket{\text{x}}^{\otimes N}$, where

$\displaystyle\ket{\text{x}}\equiv e^{-\mathrm{i}\mkern
1.0mu\frac{\pi}{2}\hat{s}_{\text{y}}}\ket{s}=\dfrac{1}{2^{s}}\sum_{\mu}{2s\choose
s+\mu}^{1/2}\ket{\mu}.$ (28)

Here ${m\choose k}$ is a binomial coefficient. As expected, the spin model
exhibits a mean-field dynamical phase transition between an ordered phase at
small $h$ and a disordered phase at large $h$. The ordered phase has a non-
zero magnetization $\sigma_{\text{MF}}$ and an interaction energy
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
that asymptotically approach their maximal values as $h\to 0$. The disordered
phase has no (time-averaged) magnetization, $\sigma_{\text{MF}}=0$, but the
interaction energy
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
nonetheless indicates persistent nontrivial inter-spin correlations when
$n>2$. These nontrivial correlations vanish as $h\to\infty$, in which case
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
approaches the minimal value allowed by conservation laws (clarified below).
By minimizing the reduced field $h$ for which $\sigma_{\text{MF}}=0$, we
numerically find that the transition between ordered and disordered phases
occurs at a critical field $h_{\text{crit}}=\left(n/2\right)^{-\alpha}$ with
$\alpha\approx 1/3$ (see Figure 6). When $n=2$, this transition is consistent
with the predictions of a Lax vector analysis [55, 56, 57, 58, 45] that
exploits integrability of $\hat{H}_{\text{spin}}$ to determine long-time
behavior. However, additional theoretical tools are necessary to understand
this transition when $n>2$. We elaborate on this point in Appendix G.

Figure 6:  The critical value of $h_{\text{crit}}$ as determined by mean-field
simulations of $N=100$ spins initially in the x-polarized state
$\ket{\text{X}}$. A single-parameter fit to
$h_{\text{crit}}=\left(n/2\right)^{-\alpha}$ finds $\alpha=0.333(5)$, and
$\alpha=1/3$ is consistent with all mean-field results to within an
uncertainty determined by the resolution of $h$ in mean-field simulations.

As shown in insets of Figure 5, mean-field results for different spin
dimensions $n$ collapse onto each other when normalizing the field $h$ to its
critical value, $h\to h\times\left(n/2\right)^{1/3}$, and rescaling

$\displaystyle\sigma_{\text{MF}}\to\dfrac{\sigma_{\text{MF}}}{\gamma\left(n/2\right)},$
$\displaystyle\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}\to\dfrac{\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}-\gamma\left(n\right)}{1-\gamma\left(n\right)},$
(29)

where

$\displaystyle\gamma\left(k\right)\equiv\dfrac{\Gamma\left(k-\frac{1}{2}\right)}{\sqrt{\pi}\,\Gamma\left(k\right)}\stackrel{{\scriptstyle
k\geq 2}}{{\approx}}\dfrac{1}{\sqrt{\pi(k-1)}}.$ (30)

The rescaling of magnetization and interaction energy can be understood by
considering their limiting behavior as $h\to\infty$ or $h\to 0$.

In the strong-field limit $h\to\infty$, we can ignore interactions and treat
spins as though they simply precess at different rates. The time-averaged
transverse magnetization $\sigma_{\text{MF}}$ then trivially vanishes as
$h\to\infty$. The interaction energy
$\braket{\bar{\bm{s}}\cdot\bar{\bm{s}}}_{\text{MF}}=\braket{\bar{\bm{s}}}_{\text{MF}}\cdot\braket{\bar{\bm{s}}}_{\text{MF}}+O(1/N)$,
meanwhile, has contributions from:  1. the diagonal parts of the mean spin
matrix $\braket{\bar{\bm{s}}}_{\text{MF}}$, which are conserved by
inhomogeneous spin precession, and  2. the off-diagonal parts of
$\braket{\bar{\bm{s}}}_{\text{MF}}$, whose oscillations average to zero when
evaluating the time average in
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$.
Altogether, the interaction energy
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
in the strong-field limit is determined by the time-independent diagonal part
$\operatorname{diag}\braket{\bar{\bm{s}}}_{\text{MF}}=\operatorname{diag}\outerproduct{\text{x}}{\text{x}}$,
namely

$\displaystyle\lim_{h\to\infty}\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}=\Tr\left[\left(\operatorname{diag}\outerproduct{\text{x}}{\text{x}}\right)^{2}\right]=\gamma\left(n\right).$
(31)

The same result can be obtained by computing the time-averaged interaction
energy of two spins precessing at different rates.

In the weak-field limit $h\to 0$, the spin-locking
$\bm{\hat{S}}\cdot\bm{\hat{S}}$ interactions of the Hamiltonian
$\hat{H}_{\text{spin}}$ energetically restrict dynamics to the permutationally
symmetric (PS) manifold. To first order in $h$, the effect of the
inhomogeneous field can be acquired by projecting it onto the PS manifold,
which takes $\hat{s}_{\text{z},q}\to\frac{1}{N}\hat{S}_{\text{z}}$. The first
order effect of the inhomogeneous field thus vanishes, as

$\displaystyle\sum_{q}\sin\left(q\right)\hat{s}_{\text{z},q}\to\sum_{q}\sin\left(q\right)\times\frac{1}{N}\hat{S}_{\text{z}}=0.$
(32)

At second order in $h$, the effective Hamiltonian within the PS manifold is
related to the variance of the inhomogeneous field, rather than its
(vanishing) average. On a high level, the second-order effect of the
inhomogeneous field within the PS manifold thus consists of permutation-
symmetrized products of two spin-z operators,
$\hat{s}_{\text{z},p}\hat{s}_{\text{z},q}$ (with $p,q$ possibly equal).
Altogether, the effective spin Hamiltonian at second order in $h$ is (see
Appendix B)

$\displaystyle\hat{H}_{\text{spin}}^{\text{eff}}=\dfrac{h^{2}u}{2(N-1)}\times\left[\hat{S}_{\text{z}}^{2}-N\sum_{q}\hat{s}_{\text{z},q}^{2}\right],$
(33)

which in the mean-field approximation becomes

$\displaystyle\hat{H}_{\text{MF}}^{\text{eff}}=-\frac{1}{2}h^{2}u\sum_{q}\hat{s}_{\text{z},q}^{2},$
(34)

where we have used the fact that the axial magnetizations
$\braket{\hat{s}_{\text{z},q}}=\frac{1}{N}\braket{\hat{S}_{\text{z}}}$ within
the PS manifold, and the initial value of $\braket{\hat{S}_{\text{z}}}=0$ is
conserved by $\hat{H}_{\text{spin}}$. The weak-field effective Hamiltonian
preserves permutational symmetry, so
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}\to
1$ as $h\to 0$. Moreover, the initial y-magnetization
$\braket{\hat{S}_{\text{y}}}=0$ is conserved by $\hat{H}_{\text{spin}}$, so
the long-time-averaged magnetization $\sigma_{\text{MF}}$ is determined by the
time-average of $\hat{s}_{\text{x}}$ for a single (any) spin:

$\displaystyle\lim_{h\to
0}\sigma_{\text{MF}}=\dfrac{1}{s}\absolutevalue{\lim_{T\to\infty}\dfrac{1}{T}\int_{0}^{T}\text{d}\tau\braket{\text{x}}{\hat{s}_{\text{x}}\left(\tau\right)}{\text{x}}},$
(35)

where

$\displaystyle\hat{s}_{\text{x}}\left(\tau\right)=e^{\mathrm{i}\mkern
1.0mu\tau\hat{s}_{\text{z}}^{2}}\hat{s}_{\text{x}}e^{-\mathrm{i}\mkern
1.0mu\tau\hat{s}_{\text{z}}^{2}}.$ (36)

We can adapt exact analytical results for the dynamics of an infinite-range
Ising model [59]dddSee Appendix K of Ref. [60] for a simpler adaptation of the
analytics in Ref. [59] to the one-axis twisting model
$\hat{H}_{\text{OAT}}=\chi\hat{s}_{\text{z}}^{2}$. to find that

$\displaystyle\braket{\text{x}}{\hat{s}_{\text{x}}\left(\tau\right)}{\text{x}}=s\left(\cos\tau\right)^{n-2},$
(37)

so for even $n$

$\displaystyle\lim_{h\to
0}\sigma_{\text{MF}}=\dfrac{1}{2\pi}\int_{0}^{2\pi}\text{d}\tau\left(\cos\tau\right)^{n-2}=\gamma\left(\dfrac{n}{2}\right).$
(38)

When going beyond mean-field theory, inter-spin correlations generated by
$\hat{S}_{\text{z}}^{2}$ in Eq. (33) will cause $\braket{\hat{S}_{\text{x}}}$
(and thereby the magnetization $\braket{\vec{\sigma}}$) to decay as
$e^{-O(t^{2}/Ns)}$; the timescale of this decay diverges as $N\to\infty$. On a
lattice of linear size $L$ without periodic boundary conditions, additional
corrections to the behavior predicted above will appear on $O(L/J)$
timescales.

###  V.2 Initial kitten states

Figure 7:  A corollary of Figure 5 for the initial state
$\ket{\text{X}\text{X}}$. The inset for interaction energy
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
in Figure 5 subtracts off the minimal value of
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
and rescales to lie on the interval $[0,1]$, as prescribed in Eq. (29). Here
the subtracting and rescaling is identical, but with a minimal value of
$2\gamma\left(n\right)$ rather than $\gamma\left(n\right)$.  Figure 8:  A
corollary of Figure 5 for the initial state
$\ket{\text{X}\text{X}_{\mathrm{i}}\mkern 1.0mu}$. Insets show the same data
shifted and rescaled identically to Figure 5.

We now consider the same setup as above, but with the initial “kitten” states
$\ket{\text{X}\text{X}}\equiv\ket{\text{x}\text{x}}^{\otimes N}$ and
$\ket{\text{X}\text{X}_{\mathrm{i}}\mkern
1.0mu}\equiv\ket{\text{x}\text{x}_{\mathrm{i}}\mkern 1.0mu}^{\otimes N}$,
where

$\displaystyle\ket{\text{x}\text{x}}\equiv\dfrac{\ket{\text{x}}+\ket{-\text{x}}}{\sqrt{2}},$
$\displaystyle\ket{\text{x}\text{x}_{\mathrm{i}}\mkern
1.0mu}\equiv\dfrac{\ket{\text{x}}+\left(-1\right)^{s}\ket{-\text{x}}}{\sqrt{2}},$
(39)

and $\ket{-\text{x}}$ is a state polarized along $-x$, defined similarly to
$\ket{\text{x}}$ in Eq. (28):

$\displaystyle\ket{-\text{x}}\equiv e^{-\mathrm{i}\mkern
1.0mu\frac{\pi}{2}\hat{s}_{\text{y}}}\ket{-s}=\dfrac{1}{2^{s}}\sum_{\mu}\left(-1\right)^{s+\mu}{2s\choose
s+\mu}^{1/2}\ket{\mu}.$ (40)

Similarly to Figure 5, Figures 7 and 8 show the time-averaged magnetization
$\sigma_{\text{MF}}$ and interaction energy
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
throughout mean-field dynamics of the initial states $\ket{\text{X}\text{X}}$
and $\ket{\text{X}\text{X}_{\mathrm{i}}\mkern 1.0mu}$. These figures exclude
the trivial case of spin dimension $n=2$, for which
$\ket{\text{x}\text{x}}=\ket{-s}$ is an eigenstate of $\hat{H}_{\text{spin}}$
and $\ket{\text{x}\text{x}_{\mathrm{i}}\mkern 1.0mu}=e^{-\mathrm{i}\mkern
1.0mu\frac{\pi}{2}\hat{s}_{\text{z}}}\ket{\text{x}}\equiv\ket{\text{y}}$ is
spin-polarized along the $y$ axis. The first and perhaps most interesting
observation to make about Figures 7 and 8 is that they are different,
signifying the importance of intra-spin coherences for the dynamical behavior
of multilevel spin models.

Unlike Figure 5 (for $\ket{\text{X}}$), Figure 7 (for
$\ket{\text{X}\text{X}}$) exhibits no sharp transition between distinct
dynamical phases: the time-averaged magnetization $\sigma_{\text{MF}}=0$ for
all values of the field $h$, and the interaction energy
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
smoothly crosses over from a maximal value of 1 to a minimal value of
$2\gamma\left(n\right)$. The minimal value of
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
approached as $h\to\infty$ can be explained with arguments identical to those
in the paragraph containing Eq. (31), which now imply that

$\displaystyle\lim_{h\to\infty}\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}=\Tr\left[\left(\operatorname{diag}\outerproduct{\text{x}\text{x}}{\text{x}\text{x}}\right)^{2}\right]=2\gamma\left(n\right).$
(41)

The vanishing magnetization $\sigma_{\text{MF}}=0$ in Figure 7 is protected by
symmetries of $\hat{H}_{\text{spin}}$ and $\ket{\text{X}\text{X}}$. For all
initial states that we have considered, the value of
$\braket{\hat{S}_{\text{z}}}=0$ is conserved by the spin Hamiltonian
$\hat{H}_{\text{spin}}$. Moreover, both the spin Hamiltonian
$\hat{H}_{\text{spin}}$ and the state $\ket{\text{X}\text{X}}$ are invariant
(up to global phase) under the action of $\hat{R}_{\text{z}}^{\pi}$, where
$\hat{R}_{\text{z}}^{\theta}\equiv e^{-\mathrm{i}\mkern
1.0mu\theta\hat{S}_{\text{z}}}$, which is to say that

$\displaystyle\hat{R}_{\text{z}}^{\pi}\hat{H}_{\text{spin}}\hat{R}_{\text{z}}^{\pi}{}^{\dagger}=\hat{H}_{\text{spin}}$
$\displaystyle\hat{R}_{\text{z}}^{\pi}\ket{\text{X}\text{X}}\simeq\ket{\text{X}\text{X}},$
(42)

where $\simeq$ denotes equality up to an overall phase. This symmetry implies
that

$\displaystyle\braket{\hat{S}_{\text{x}}}=\braket{\hat{R}_{\text{z}}^{\pi}{}^{\dagger}\hat{S}_{\text{x}}\hat{R}_{\text{z}}^{\pi}}=-\braket{\hat{S}_{\text{x}}}=0,$
(43)
$\displaystyle\braket{\hat{S}_{\text{y}}}=\braket{\hat{R}_{\text{z}}^{\pi}{}^{\dagger}\hat{S}_{\text{y}}\hat{R}_{\text{z}}^{\pi}}=-\braket{\hat{S}_{\text{y}}}=0$
(44)

at all times, so altogether $\sigma_{\text{MF}}=0$.

Turning now to mean-field results for the initial kitten state
$\ket{\text{X}\text{X}_{\mathrm{i}}\mkern 1.0mu}$ in Figure 8, we remark that
the magnetization $\sigma_{\text{MF}}$ and interaction energy
$\langle\\!\langle\bar{\bm{s}}\cdot\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
behave identically to those for the initial spin-polarized state
$\ket{\text{X}}$ in Figure 5. This finding can be understood through the fact
that

$\displaystyle\ket{\text{X}\text{X}_{\mathrm{i}}\mkern
1.0mu}\simeq\hat{R}_{\text{z}}^{\pi/2}\hat{T}_{\text{z}}^{\pi/2}\ket{\text{X}},$
(45)

where $\hat{T}_{\text{z}}^{\theta}\equiv e^{-\mathrm{i}\mkern
1.0mu\theta\hat{S}_{\text{z}}^{2}}$. The operators
$\hat{R}_{\text{z}}^{\theta}$ and $\hat{T}_{\text{z}}^{\theta}$ are generated
by axial fields that respect permutational symmetry, and therefore commute
with the spin Hamiltonian $\hat{H}_{\text{spin}}$, so

$\displaystyle e^{-\mathrm{i}\mkern
1.0mut\hat{H}_{\text{spin}}}\ket{\text{X}\text{X}_{\mathrm{i}}\mkern 1.0mu}$
$\displaystyle\simeq e^{-\mathrm{i}\mkern
1.0mut\hat{H}_{\text{spin}}}\hat{R}_{\text{z}}^{\pi/2}\hat{T}_{\text{z}}^{\pi/2}\ket{\text{X}}$
(46)
$\displaystyle\simeq\hat{R}_{\text{z}}^{\pi/2}\hat{T}_{\text{z}}^{\pi/2}e^{-\mathrm{i}\mkern
1.0mut\hat{H}_{\text{spin}}}\ket{\text{X}}.$ (47)

In turn, expanding $\bar{\bm{s}}\cdot\bar{\bm{s}}$ according to Eq. (8) shows
that

$\displaystyle\hat{T}_{\text{z}}^{\theta}{}^{\dagger}\hat{R}_{\text{z}}^{\theta}{}^{\dagger}\,\bar{\bm{s}}\cdot\bar{\bm{s}}\,\hat{R}_{\text{z}}^{\theta}\hat{T}_{\text{z}}^{\theta}=\bar{\bm{s}}\cdot\bar{\bm{s}},$
(48)

which implies that the interaction energy
$\braket{\bar{\bm{s}}\cdot\bar{\bm{s}}}$ throughout dynamics of the initial
kitten state $\ket{\text{X}\text{X}_{\mathrm{i}}\mkern 1.0mu}$ is the same as
that of the spin-polarized state $\ket{\text{X}}$.

To make sense of why the magnetization $\sigma_{\text{MF}}$ is identical in
Figure 8 for $\ket{\text{X}\text{X}_{\mathrm{i}}\mkern 1.0mu}$ as in Figure 5
for $\ket{\text{X}}$, we follow a four-part argument:

  1. (i)

The time-averaged magnetization vector
$\langle\\!\langle\vec{\sigma}\rangle\\!\rangle_{\text{MF}}$ can be written as
a function of the time-averaged spin matrix
$\langle\\!\langle\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$.

  2. (ii)

The spin matrix $\langle\\!\langle\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$
is only ever nonzero on its diagonal and anti-diagonal, regardless of the
initial state. That is, nonzero components
$\langle\\!\langle\bar{s}_{\mu\nu}\rangle\\!\rangle_{\text{MF}}$ of
$\langle\\!\langle\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$ always have
$\mu=\pm\nu$ (see discussion below).

  3. (iii)

The twist operator $\hat{T}_{\text{z}}^{\theta}$ acts trivially on the
diagonal and anti-diagonal components of $\bar{\bm{s}}$, which together with
point (ii) implies that
$\langle\\!\langle\hat{T}_{\text{z}}^{\theta}{}^{\dagger}\,\bar{\bm{s}}\,\hat{T}_{\text{z}}^{\theta}\rangle\\!\rangle_{\text{MF}}=\langle\\!\langle\bar{\bm{s}}\rangle\\!\rangle_{\text{MF}}$.

  4. (iv)

The rotation operator $\hat{R}_{\text{z}}^{\theta}$ merely rotates the
magnetization vector
$\langle\\!\langle\vec{\sigma}\rangle\\!\rangle_{\text{MF}}$ without changing
its magnitude.

Altogether, points (i)–(iv) imply that the magnetization

$\displaystyle\sigma_{\text{MF}}=\absolutevalue{\langle\\!\langle\vec{\sigma}\rangle\\!\rangle_{\text{MF}}}=\absolutevalue{\langle\\!\langle\hat{T}_{\text{z}}^{\theta}{}^{\dagger}\hat{R}_{\text{z}}^{\theta}{}^{\dagger}\,\vec{\sigma}\,\hat{R}_{\text{z}}^{\theta}\hat{T}_{\text{z}}^{\theta}\rangle\\!\rangle_{\text{MF}}}$
(49)

is the same for the initial state $\ket{\text{X}\text{X}_{\mathrm{i}}\mkern
1.0mu}$ as for $\ket{\text{X}}$.

The only nontrivial step in the above argument is point (ii), which says that
$\langle\\!\langle\bar{s}_{\mu\nu}\rangle\\!\rangle_{\text{MF}}$ is guaranteed
to be zero unless $\mu=\pm\nu$. This observation, nominally a numerical result
of mean-field simulations, can be understood as follows. The eigenstates
$\ket{m,w}$ of $\hat{H}_{\text{spin}}$ are uniquely identified by definite
numbers $m=(m_{s},m_{s-1},\cdots,m_{-s})$ of atoms occupying each internal
spin state $\mu\in\left\\{s,s-1,\cdots,-s\right\\}$, and an auxiliary index
$w$ that encodes how $\ket{m,w}$ transforms under permutations of all spins
(see Appendix B)eeeSeen otherwise, since $\hat{S}_{\mu\mu}$ commutes with
$\hat{H}_{\text{spin}}$, eigenvectors of $\hat{H}_{\text{spin}}$ can be
indexed by eigenvalues of $\hat{S}_{\mu\mu}$. The number $m_{\mu}$ is then the
eigenvalue of $\ket{m,w}$ with respect to $\hat{S}_{\mu\mu}$, i.e.
$\hat{S}_{\mu\mu}\ket{m,w}=m_{\mu}\ket{m,w}$, while $w$ encodes all other
information required to uniquely specify $\ket{m,w}$.. The operator
$\bar{s}_{\mu\nu}=\frac{1}{N}\hat{S}_{\mu\nu}$ with $\mu\neq\nu$ couples the
state $\ket{m,w}$ to states $\ket{m^{\prime},w^{\prime}}$ in which
$(m_{\mu}^{\prime},m^{\prime}_{\nu})=(m_{\mu}+1,m_{\nu}-1)$. Generically,
states $\ket{m,w}$ and $\ket{m^{\prime},w^{\prime}}$ with $m\neq m^{\prime}$
will have different energies, so their coherence oscillates and averages to
zero when evaluating time-averaged expectation values.

However, degeneracies yield stationary (time-independent) coherences that
survive time-averaging. In the weak-field limit $h\to 0$, such a degeneracy
occurs at the mean-field level between PS states differing only in the
populations $m_{\mu},m_{-\mu}$ (with a fixed value of $m_{\mu}+m_{-\mu}$), as
the effective Hamiltonian becomes
$\hat{H}_{\text{MF}}^{\text{eff}}\propto\sum_{\mu}\mu^{2}m_{\mu}$. This
symmetry is preserved at all orders in perturbation theoryfffOnly even powers
of the “perturbation” $\sum_{q}\sin\left(q\right)\hat{s}_{\text{z},q}$ can be
nonzero within the PS manifold, and even powers of this perturbation exhibit
the same mean-field degeneracy between states differing only in the
populations $m_{\mu},m_{-\mu}$., so some coherence between such states is
preserved as $h\to h_{\text{crit}}$, although this coherence decays as
perturbative corrections to degenerate eigenstates cause them to leak out of
the PS manifold (and thereby have a smaller overlap with the initial state
$\ket{\text{X}}$). Note that beyond-mean-field effects break the symmetry
protecting anti-diagonal components of
$\langle\\!\langle\bm{\hat{s}}\rangle\\!\rangle_{\text{MF}}$, causing them to
decay on time scales that should diverge as $N\to\infty$.

##  VI Conclusions and future directions

Starting with an SU($n$) Hubbard model describing ultracold fermionic
alkaline-earth(-like) atoms on an optical lattice, we derived a momentum-space
multilevel spin model with all-to-all SU($n$)-symmetric interactions. We then
introduced external control fields, finding a simple three-laser drive that
homogeneously addresses nuclear spins with a variety of spin Hamiltonians.
Taking a closer look at the effect of the spin-orbit coupling (SOC) induced by
the driving lasers, we found that maintaining the validity of the spin model
requires weak SOC, which in turn gives rise to a (synthetic) inhomogeneous
magnetic field. Finally, we examined dynamical behavior of the SU($n$) spin
model at the mean-field level, finding that long-time observables obey simple
scaling relations with $n$, and that when $n>2$ dynamical behavior can be
highly sensitive to intra-spin coherences.

Our work makes important progress in understanding the SU($n$) Fermi-Hubbard
model in experimentally relevant parameter regimes, and we expect our findings
to be readily testable in experiments with ultracold atoms. Given the
possibility for long-range SU($n$) interactions, we hope our work stimulates
further efforts into simulating SY and SYK-like models [15, 61] in cold atomic
platforms. In follow-up work, it would be interesting to study the
relationship between initial states and dynamical phases of our SU($n$) spin
model more systematically, and to consider the effect of quantum corrections
to mean-field behavior. There is also room to improve on the three-laser drive
introduced in this work, for which it is natural to ask what additional
techniques or ingredients are necessary to implement universal control of
individual nuclear spins. Universal control would allow for an experimental
study of $n$-dependence (including even/odd-$n$ parity effects) in a single
experimental platform, simply by controlling the occupation and coherence of
internal spin states. Finally, one can also study the SU($n$) Hubbard model in
the super-exchange regime that gives rise to a real-space (as opposed to
momentum-space) spin model, where SOC gives rise to chiral multilevel spin
interactions. Unlike our present work, the super-exchange regime does not
require weak SOC, and therefore has a larger parameter space in which to
explore dynamical behavior.

## Acknowledgements

We thank Victor Gurarie, Emil Yuzbashyan, Asier P. Orioli, and Jeremy T. Young
for helpful discussions on this work. This work was supported by the AFOSR
grants FA9550-18-1-0319, FA9550-19-1-027, by the DARPA and ARO grant
W911NF-16-1-0576, ARO W911NF-19-1-0210, DOE (QSA), NSF PHY1820885, NSF JILA-
PFC PHY-1734006, NSF QLCI-2016244, and by NIST.

## References

  * Maldacena [1999] J. Maldacena, The Large-N Limit of Superconformal Field Theories and Supergravity, International Journal of Theoretical Physics 38, 1113 (1999).
  * Ryu and Takayanagi [2006] S. Ryu and T. Takayanagi, Holographic Derivation of Entanglement Entropy from the anti-de Sitter Space/Conformal Field Theory Correspondence, Physical Review Letters 96, 181602 (2006).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a Mott insulator: Physics of high-temperature superconductivity, Reviews of Modern Physics 78, 17 (2006).
  * Read and Sachdev [1989] N. Read and S. Sachdev, Valence-bond and spin-Peierls ground states of low-dimensional quantum antiferromagnets, Physical Review Letters 62, 1694 (1989).
  * Rokhsar [1990] D. S. Rokhsar, Quadratic quantum antiferromagnets in the fermionic large-$n$ limit, Physical Review B 42, 2526 (1990).
  * Kaul and Sandvik [2012] R. K. Kaul and A. W. Sandvik, Lattice model for the $\mathrm{SU}(n)$ néel to valence-bond solid quantum phase transition at large $n$, Physical Review Letters 108, 137201 (2012).
  * Hermele and Gurarie [2011] M. Hermele and V. Gurarie, Topological liquids and valence cluster states in two-dimensional $\mathrm{SU}(n)$ magnets, Physical Review B 84, 174441 (2011).
  * Hermele  _et al._ [2009] M. Hermele, V. Gurarie, and A. M. Rey, Mott Insulators of Ultracold Fermionic Alkaline Earth Atoms: Underconstrained Magnetism and Chiral Spin Liquid, Physical Review Letters 103, 135301 (2009).
  * Chen  _et al._ [2016] G. Chen, K. R. A. Hazzard, A. M. Rey, and M. Hermele, Synthetic-gauge-field stabilization of the chiral-spin-liquid phase, Physical Review A 93, 061601 (2016).
  * Nataf  _et al._ [2016] P. Nataf, M. Lajkó, A. Wietek, K. Penc, F. Mila, and A. M. Läuchli, Chiral spin liquids in triangular-lattice SU($n$) fermionic mott insulators with artificial gauge fields, Physical Review Letters 117, 167202 (2016).
  * Freedman  _et al._ [2004] M. Freedman, C. Nayak, K. Shtengel, K. Walker, and Z. Wang, A class of P,T-invariant topological phases of interacting electrons, Annals of Physics 310, 428 (2004).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-Abelian anyons and topological quantum computation, Reviews of Modern Physics 80, 1083 (2008).
  * Nataf and Mila [2014] P. Nataf and F. Mila, Exact diagonalization of heisenberg SU($N$) models, Physical Review Letters 113, 127204 (2014).
  * Nataf and Mila [2016] P. Nataf and F. Mila, Exact diagonalization of heisenberg SU($n$) chains in the fully symmetric and antisymmetric representations, Physical Review B 93, 155134 (2016).
  * Sachdev and Ye [1993] S. Sachdev and J. Ye, Gapless spin-fluid ground state in a random quantum Heisenberg magnet, Physical Review Letters 70, 3339 (1993).
  * Wu  _et al._ [2003] C. Wu, J.-p. Hu, and S.-c. Zhang, Exact SO(5) symmetry in the spin-$3/2$ fermionic system, Physical Review Letters 91, 186402 (2003).
  * Cazalilla  _et al._ [2009] 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).
  * Gorshkov  _et al._ [2010] A. V. Gorshkov, M. Hermele, V. Gurarie, C. Xu, P. S. Julienne, J. Ye, P. Zoller, E. Demler, M. D. Lukin, and A. M. Rey, Two-orbital SU($n$) magnetism with ultracold alkaline-earth atoms, Nature Physics 6, 289 (2010).
  * Cazalilla and Rey [2014] M. A. Cazalilla and A. M. Rey, Ultracold fermi gases with emergent $\mathrm{SU}(n)$ symmetry, Reports on Progress in Physics 77, 124401 (2014).
  * Hazzard  _et al._ [2012] K. R. A. Hazzard, V. Gurarie, M. Hermele, and A. M. Rey, High-temperature properties of fermionic alkaline-earth-metal atoms in optical lattices, Physical Review A 85, 041604 (2012).
  * Bonnes  _et al._ [2012] L. Bonnes, K. R. A. Hazzard, S. R. Manmana, A. M. Rey, and S. Wessel, Adiabatic loading of one-dimensional $\mathrm{SU}(n)$ alkaline-earth-atom fermions in optical lattices, Physical Review Letters 109, 205305 (2012).
  * Stellmer  _et al._ [2013] S. Stellmer, F. Schreck, and T. C. Killian, Degenerate quantum gases of strontium, in  _Annual Review of Cold Atoms and Molecules_, Annual Review of Cold Atoms and Molecules, Vol. Volume 2 (WORLD SCIENTIFIC, 2013) pp. 1–80.
  * Yip  _et al._ [2014] S.-K. Yip, B.-L. Huang, and J.-S. Kao, Theory of $\mathrm{SU}(n)$ fermi liquids, Physical Review A 89, 043610 (2014).
  * Pagano  _et al._ [2014] G. Pagano, M. Mancini, G. Cappellini, P. Lombardi, F. Schäfer, H. Hu, X.-J. Liu, J. Catani, C. Sias, M. Inguscio, and L. Fallani, A one-dimensional liquid of fermions with tunable spin, Nature Physics 10, 198 (2014).
  * Choudhury  _et al._ [2020] S. Choudhury, K. R. Islam, Y. Hou, J. A. Aman, T. C. Killian, and K. R. A. Hazzard, Collective modes of ultracold fermionic alkaline-earth-metal gases with SU($n$) symmetry, Physical Review A 101, 053612 (2020).
  * Song  _et al._ [2020] B. Song, Y. Yan, C. He, Z. Ren, Q. Zhou, and G.-B. Jo, Evidence for bosonization in a three-dimensional gas of $\mathrm{SU}(n)$ fermions, Physical Review X 10, 041053 (2020).
  * Sonderhouse  _et al._ [2020] L. Sonderhouse, C. Sanner, R. B. Hutson, A. Goban, T. Bilitewski, L. Yan, W. R. Milner, A. M. Rey, and J. Ye, Thermodynamics of a deeply degenerate SU($n$)-symmetric fermi gas, Nature Physics 16, 1216 (2020).
  * Taie  _et al._ [2012] S. Taie, R. Yamazaki, S. Sugawa, and Y. Takahashi, An SU(6) Mott insulator of an atomic Fermi gas realized by large-spin Pomeranchuk cooling, Nature Physics 8, 825 (2012).
  * Hofrichter  _et al._ [2016] C. Hofrichter, L. Riegger, F. Scazza, M. Höfer, D. R. Fernandes, I. Bloch, and S. Fölling, Direct probing of the mott crossover in the $\mathrm{SU}(n)$ fermi-hubbard model, Physical Review X 6, 021030 (2016).
  * Taie  _et al._ [2020] S. Taie, E. Ibarra-Garc\́mathrm{i} a-Padilla, N. Nishizawa, Y. Takasu, Y. Kuno, H.-T. Wei, R. T. Scalettar, K. R. A. Hazzard, and Y. Takahashi, Observation of antiferromagnetic correlations in an ultracold SU($n$) hubbard model, arXiv (2020), arXiv:2010.07730 [cond-mat] .
  * Messio and Mila [2012] L. Messio and F. Mila, Entropy dependence of correlations in one-dimensional $\mathrm{SU}(n)$ antiferromagnets, Physical Review Letters 109, 205306 (2012).
  * Cappellini  _et al._ [2014] 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, Physical Review Letters 113, 120402 (2014).
  * Scazza  _et al._ [2014] F. Scazza, C. Hofrichter, M. Höfer, P. C. D. Groot, I. Bloch, and S. Fölling, Observation of two-orbital spin-exchange interactions with ultracold $\mathrm{SU}(n)$-symmetric fermions, Nature Physics 10, 779 (2014).
  * Zhang  _et al._ [2014] X. Zhang, M. Bishof, S. L. Bromley, C. V. Kraus, M. S. Safronova, P. Zoller, A. M. Rey, and J. Ye, Spectroscopic observation of $\mathrm{SU}(n)$-symmetric interactions in sr orbital magnetism, Science 345, 1467 (2014).
  * Beverland  _et al._ [2016] M. E. Beverland, G. Alagic, M. J. Martin, A. P. Koller, A. M. Rey, and A. V. Gorshkov, Realizing exactly solvable $\mathrm{SU}(n)$ magnets with thermal atoms, Physical Review A 93, 051601 (2016).
  * Goban  _et al._ [2018] 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 (2018).
  * Perlin and Rey [2019] M. A. Perlin and A. M. Rey, Effective multi-body SU($N$)-symmetric interactions of ultracold fermionic atoms on a 3D lattice, New Journal of Physics 21, 043039 (2019).
  * He  _et al._ [2019] P. He, M. A. Perlin, S. R. Muleady, R. J. Lewis-Swan, R. B. Hutson, J. Ye, and A. M. Rey, Engineering spin squeezing in a 3D optical lattice with interacting spin-orbit-coupled fermions, Physical Review Research 1, 033075 (2019).
  * Perlin  _et al._ [2020a] M. A. Perlin, C. Qu, and A. M. Rey, Spin Squeezing with Short-Range Spin-Exchange Interactions, Physical Review Letters 125, 223401 (2020a).
  * Mancini  _et al._ [2015] M. Mancini, G. Pagano, G. Cappellini, L. Livi, M. Rider, J. Catani, C. Sias, P. Zoller, M. Inguscio, M. Dalmonte, and L. Fallani, Observation of chiral edge states with neutral fermions in synthetic Hall ribbons, Science 349, 1510 (2015).
  * Wall  _et al._ [2016] M. L. Wall, A. P. Koller, S. Li, X. Zhang, N. R. Cooper, J. Ye, and A. M. Rey, Synthetic Spin-Orbit Coupling in an Optical Lattice Clock, Physical Review Letters 116, 035301 (2016).
  * Livi  _et al._ [2016] 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, Physical Review Letters 117, 220401 (2016).
  * Kolkowitz  _et al._ [2016] S. Kolkowitz, S. L. Bromley, T. Bothwell, M. L. Wall, G. E. Marti, A. P. Koller, X. Zhang, A. M. Rey, and J. Ye, Spin-orbit-coupled fermions in an optical lattice clock, Nature 542, 66 (2016).
  * Bromley  _et al._ [2018] S. L. Bromley, S. Kolkowitz, T. Bothwell, D. Kedar, A. Safavi-Naini, M. L. Wall, C. Salomon, A. M. Rey, and J. Ye, Dynamics of interacting fermions under spin–orbit coupling in an optical lattice clock, Nature Physics 14, 399 (2018).
  * Smale  _et al._ [2019] S. Smale, P. He, B. A. Olsen, K. G. Jackson, H. Sharum, S. Trotzky, J. Marino, A. M. Rey, and J. H. Thywissen, Observation of a transition between dynamical phases in a quantum degenerate Fermi gas, Science Advances 5, eaax1568 (2019).
  * Lewis-Swan  _et al._ [2021] R. J. Lewis-Swan, D. Barberena, J. R. K. Cline, D. J. Young, J. K. Thompson, and A. M. Rey, Cavity-QED Quantum Simulator of Dynamical Phases of a Bardeen-Cooper-Schrieffer Superconductor, Physical Review Letters 126, 173601 (2021).
  * Bravyi  _et al._ [2011] S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer–Wolff transformation for quantum many-body systems, Annals of Physics 326, 2793 (2011).
  * Ma  _et al._ [2011] J. Ma, X. Wang, C. P. Sun, and F. Nori, Quantum spin squeezing, Physics Reports 509, 89 (2011).
  * Dowling  _et al._ [1994] J. P. Dowling, G. S. Agarwal, and W. P. Schleich, Wigner distribution of a general angular-momentum state: Applications to a collection of two-level atoms, Physical Review A 49, 4101 (1994).
  * Li  _et al._ [2013] F. Li, C. Braun, and A. Garg, The Weyl-Wigner-Moyal formalism for spin, EPL (Europhysics Letters) 102, 60006 (2013).
  * Newton and Young [1968] R. G. Newton and B.-l. Young, Measurability of the spin density matrix, Annals of Physics 49, 393 (1968).
  * Perlin  _et al._ [2020b] M. A. Perlin, D. Barberena, and A. M. Rey, Spin qudit tomography and state reconstruction error, arXiv (2020b), arXiv:2012.06464 [quant-ph] .
  * Anderson  _et al._ [2015] B. E. Anderson, H. Sosa-Martinez, C. A. Riofrio, I. H. Deutsch, and P. S. Jessen, Accurate and Robust Unitary Transformations of a High-Dimensional Quantum System, Physical Review Letters 114, 240401 (2015).
  * Lucarelli [2018] D. Lucarelli, Quantum optimal control via gradient ascent in function space and the time-bandwidth quantum speed limit, Physical Review A 97, 062346 (2018).
  * Yuzbashyan  _et al._ [2005] E. A. Yuzbashyan, B. L. Altshuler, V. B. Kuznetsov, and V. Z. Enolskii, Nonequilibrium cooper pairing in the nonadiabatic regime, Physical Review B 72, 220503 (2005).
  * Yuzbashyan and Dzero [2006] E. A. Yuzbashyan and M. Dzero, Dynamical Vanishing of the Order Parameter in a Fermionic Condensate, Physical Review Letters 96, 230404 (2006).
  * Yuzbashyan  _et al._ [2006] E. A. Yuzbashyan, O. Tsyplyatyev, and B. L. Altshuler, Relaxation and Persistent Oscillations of the Order Parameter in Fermionic Condensates, Physical Review Letters 96, 097005 (2006).
  * Yuzbashyan  _et al._ [2015] E. A. Yuzbashyan, M. Dzero, V. Gurarie, and M. S. Foster, Quantum quench phase diagrams of an s-wave BCS-BEC condensate, Physical Review A 91, 033628 (2015).
  * Foss-Feig  _et al._ [2013] M. Foss-Feig, K. R. A. Hazzard, J. J. Bollinger, and A. M. Rey, Nonequilibrium dynamics of arbitrary-range Ising models with decoherence: An exact analytic solution, Physical Review A 87, 042101 (2013).
  * Perlin and Rey [2020] M. A. Perlin and A. M. Rey, Short-time expansion of Heisenberg operators in open collective quantum spin systems, Physical Review A 101, 023601 (2020).
  * Bentsen  _et al._ [2019] G. Bentsen, I.-D. Potirniche, V. B. Bulchandani, T. Scaffidi, X. Cao, X.-L. Qi, M. Schleier-Smith, and E. Altman, Integrable and Chaotic Dynamics of Spins Coupled to an Optical Cavity, Physical Review X 9, 041011 (2019).

##  Appendix A Numerical benchmarking of the spin model

In this appendix we present numerical evidence to support the validity of the
spin models derived in Sections II and IV. Figures 9 and 10 show a set of
time-averaged observables computed with numerically exact simulations of a
Fermi-Hubbard model and an effective spin model, respectively, with $n=4$
(Figure 9) and $n=6$ (Figure 10) internal levels per spin. Details for these
simulations are provided in the caption of Figure 9. Our main conclusion from
these figures is that the two models show remarkable agreement for the
observables considered in our work. Note that these results are only intended
to benchmark the approximation of a Fermi-Hubbard model by a spin model; these
results are not expected to agree with the mean-field theory in Section V due
to strong finite-size effects.

Figure 9:  Numerical results (analogous to Figures 5, 7, and 8 of the main
text) for the time-averaged interaction energy and magnetization (both
normalized to a maximal value of 1) in a system of $L=5$ lattice sites, for
both a Fermi-Hubbard model (dots) and spin model (lines) with $n=4$ internal
states per spin. The corresponding initial state (defined in Section V of the
main text) is indicated in each panel, and observables are averaged over a
time $tJ=200$. Color indicates the value of $U/J$, and the field $h$
corresponds to $2J\phi/u$ in the case of the Fermi-Hubbard model. Simulations
are performed in real-space, with spin-orbit coupling (SOC) implemented
through a homogeneous drive (with no site or $\phi$ dependence) and nearest-
neighbor tunneling terms that contain factors of $e^{\pm\mathrm{i}\mkern
1.0mu\mu\phi}$. Results for the initial kitten state
$\ket{\text{X}\text{X}_{\mathrm{i}}\mkern 1.0mu}$ are excluded because they
are identical to those of $\ket{\text{X}}$, and magnetization for the initial
state $\ket{\text{X}\text{X}}$ is always $0$. Note that while panels (a) and
(b) are representative of infinite-time behavior, the inset in panel (c) shows
that the Fermi-Hubbard and spin models exhibit different behaviors on very
long time scales, although good agreement is restored by rescaling time in the
spin model, indicating the likelihood of a need to renormalize spin model
parameters. In any case, such time scales are inaccessible in current
experiments and diverge as $N\to\infty$, so these corrections do not affect
the main results of our work.  Figure 10:  Numerical results identical to
Figure 9, but with $L=4$ lattice sites and $n=6$ internal states per spin.

##  Appendix B Perturbation theory for SU($n$) ferromagnets

Here we work out a general perturbation theory for SU($n$) ferromagnets with a
gapped permutationally symmetric (PS) manifold. We begin with an
SU($n$)-symmetric interaction Hamiltonian of the form

$\displaystyle\hat{H}_{0}=\sum_{i<j}g_{ij}\,\hat{\Pi}_{ij},$
$\displaystyle\hat{\Pi}_{ij}\equiv\bm{\hat{s}}_{i}\cdot\bm{\hat{s}}_{j}=\sum_{\mu,\nu}\hat{s}_{\mu\nu
i}\hat{s}_{\nu\mu j},$ (50)

where $g_{ij}$ are (real) scalar coefficients for the permutation operators
$\hat{\Pi}_{ij}$, and $\hat{s}_{\mu\nu i}\equiv\hat{c}_{\mu
i}^{\dagger}\hat{c}_{\nu i}$ is a transition operator for spin $i$. We can
then consider the addition of, for example, an inhomogeneous magnetic field or
Ising couplings,

$\displaystyle\hat{H}_{\text{field}}=\sum_{i}B_{i}\hat{s}_{\text{z},i},$
$\displaystyle\hat{H}_{\text{Ising}}=\sum_{i\neq
j}J_{ij}\hat{s}_{\text{z},i}\hat{s}_{\text{z},j},$ (51)

or more generally an $M$-body operatorgggAt face value, an $M$-body operator
with $M>2$ does not typically appear in experiments. Nonetheless, considering
$M>2$ illuminates the structure of eigenstates (and eigenvalues) of
$\hat{H}_{0}$, and allows us to go to high orders in perturbation theory with
single- and two-body perturbations.

$\displaystyle\hat{\mathcal{O}}(w,\hat{X})=\sum_{k\in\mathcal{D}_{N}\left(M\right)}w_{k}\hat{X}_{k},$
(52)

where $w$ is a dimension-$M$ (i.e. $M$-index) tensor of scalar coefficients
$w_{k}\equiv w_{k_{1}k_{2}\cdots k_{M}}$; $X$ is an $M$-spin operator, e.g.
$\hat{s}_{\text{z}}\otimes\hat{s}_{\text{z}}$ in the case of Ising
interactions with $M=2$; $k\equiv\left(k_{1},k_{2},\cdots,k_{M}\right)$ is a
list of the individual spins
$k_{i}\in\mathbb{Z}_{N}\equiv\left\\{1,2,\cdots,N\right\\}$ that the operator
$\hat{X}_{k}\equiv\hat{X}_{k_{1}k_{2}\cdots k_{M}}$ acts on; and

$\displaystyle\mathcal{D}_{N}\left(M\right)\equiv\left\\{k\in\mathbb{Z}_{N}^{M}:\text{all
entries $k_{i}$ of $k$ are distinct}\right\\},$ (53)

is the strictly “off-diagonal” part of $\mathbb{Z}_{N}^{M}$, which is
necessary to identify for a consistent definition of $\hat{X}_{k}$ as an
$M$-body operator. In this notation, the magnetic field and Ising Hamiltonians
in Eq. (51) respectively become
$\hat{\mathcal{O}}\left(B,\hat{s}_{\text{z}}\right)$ and
$\hat{\mathcal{O}}\left(J,\hat{s}_{\text{z}}\otimes\hat{s}_{\text{z}}\right)$.

If the addition $\hat{\mathcal{O}}(w,\hat{X})$ to the SU($n$)-symmetric
Hamiltonian $\hat{H}_{0}$ in Eq. (50) is sufficiently small, namely with
operator norm $\norm*{\hat{\mathcal{O}}(w,\hat{X})}$ less than half the
spectral gap $\Delta_{\text{gap}}$ of $\hat{H}_{0}$,
$\norm*{\hat{\mathcal{O}}(w,\hat{X})}<\Delta_{\text{gap}}/2$, then we can
treat the effect of $\hat{\mathcal{O}}(w,\hat{X})$ on the ground-state PS
manifold $\mathcal{E}_{0}$ perturbatively. The effective Hamiltonians
$\hat{H}_{\text{eff}}^{(1)}$ and $\hat{H}_{\text{eff}}^{(2)}$ induced by
$\hat{\mathcal{O}}(w,\hat{X})$ on the PS manifold $\mathcal{E}_{0}$ at leading
orders in perturbation theory are [47]

$\displaystyle\hat{H}_{\text{eff}}^{(1)}=\hat{\mathcal{P}}_{0}\hat{\mathcal{O}}(w,\hat{X})\hat{\mathcal{P}}_{0},$
$\displaystyle\hat{H}_{\text{eff}}^{(2)}=-\sum_{\Delta\neq
0}\dfrac{1}{\Delta}\hat{\mathcal{P}}_{0}\hat{\mathcal{O}}(w,\hat{X})\hat{\mathcal{P}}_{\Delta}\hat{\mathcal{O}}(w,\hat{X})\hat{\mathcal{P}}_{0},$
(54)

where $\hat{\mathcal{P}}_{\Delta}$ is a projector onto the eigenspace
$\mathcal{E}_{\Delta}$ of $\hat{H}_{0}$ with interaction energy $\Delta$ above
that of the PS manifold. The first order effective Hamiltonian
$\hat{H}_{\text{eff}}^{(1)}$ simply projects $\hat{\mathcal{O}}(w,\hat{X})$
onto the PS manifold $\mathcal{E}_{0}$, and takes the form

$\displaystyle\hat{H}_{\text{eff}}^{(1)}=\overline{w}\,\underline{X},$ (55)

where the coefficient $\overline{w}$ is the average of all coefficients
$w_{k}$; and $\underline{X}$ is a collective version of $X$:

$\displaystyle\overline{w}\equiv\dfrac{1}{\absolutevalue{\mathcal{D}_{N}\left(M\right)}}\sum_{k\in\mathcal{D}_{N}\left(M\right)}w_{k},$
$\displaystyle\underline{X}\equiv\sum_{k\in\mathcal{D}_{N}\left(M\right)}\hat{X}_{k},$
(56)

with
$\absolutevalue{\mathcal{D}_{N}\left(M\right)}=\prod_{j=0}^{M-1}\left(N-j\right)$.
In the case of a magnetic field $\hat{s}_{\text{z}}$ or Ising interactions
$\hat{s}_{\text{z}}\otimes\hat{s}_{\text{z}}$, for example,

$\displaystyle\underline{\hat{s}_{\text{z}}}=\sum_{i}\hat{s}_{\text{z}}^{(i)}=\hat{S}_{\text{z}},$
$\displaystyle\underline{\hat{s}_{\text{z}}\otimes\hat{s}_{\text{z}}}=\sum_{i\neq
j}\hat{s}_{\text{z}}^{(i)}\hat{s}_{\text{z}}^{(j)}=\hat{S}_{\text{z}}^{2}-N\sum_{i}\hat{s}_{\text{z},i}^{2}.$
(57)

The second order effective Hamiltonian $\hat{H}_{\text{eff}}^{(2)}$ in Eq.
(54) takes more work to simplify due to the presence of a projector
$\hat{\mathcal{P}}_{\Delta}$ onto the manifold $\mathcal{E}_{\Delta}$ of
states with excitation energy $\Delta$. This projector essentially picks off
the part of $\hat{\mathcal{O}}(w,\hat{X})$ that is strictly off-diagonal with
respect to the ground- and excited-state manifolds $\mathcal{E}_{0}$ and
$\mathcal{E}_{\Delta}$. We therefore need to decompose
$\hat{\mathcal{O}}(w,\hat{X})$ into components that generate states of
definite excitation energy when acting on PS states
$\ket{\psi}\in\mathcal{E}_{0}$. The SU($n$) symmetry of $\hat{H}_{0}$ enables
such a decomposition to take the form

$\displaystyle\hat{H}_{0}\hat{\mathcal{O}}(w,\hat{X})\ket{\psi}=\sum_{\Delta}\left(E_{0}+\Delta\right)\hat{\mathcal{O}}(w^{\Delta},\hat{X})\ket{\psi},$
$\displaystyle E_{0}\equiv\sum_{i<j}g_{ij},$ (58)

where $E_{0}$ is the interaction energy of PS states, and thinking of the
tensor $w$ as a $\absolutevalue{\mathcal{D}_{N}\left(M\right)}$-component
vector, the tensor $w^{\Delta}$ can be found by  1. using the coefficients
$g_{ij}$ to construct a matrix $g^{(M)}$ of dimensions
$\absolutevalue{\mathcal{D}_{N}\left(M\right)}\times\absolutevalue{\mathcal{D}_{N}\left(M\right)}\sim
N^{M}\times N^{M}$, and  2. projecting $w$ onto the eigenspace of $g^{(M)}$
with eigenvalue $\Delta$.  We construct $g^{(M)}$ for the single-body ($M=1$)
case below (in Appendix B.1), and provide explicit forms of $g^{(M)}$ with
arbitrary $M$.

Equipped with the decomposition
$\hat{\mathcal{O}}(w,\hat{X})=\sum_{\Delta}\hat{\mathcal{O}}(w^{\Delta},\hat{X})$
with terms $\hat{\mathcal{O}}(w^{\Delta},\hat{X})$ that generate states of
definite excitation energy $\Delta$, we can expand

$\displaystyle\hat{H}_{\text{eff}}^{(2)}=-\sum_{\Delta\neq
0}\dfrac{1}{\Delta}\hat{\mathcal{P}}_{0}\hat{\mathcal{O}}(w^{\Delta},\hat{X})^{2}\hat{\mathcal{P}}_{0}.$
(59)

If $X$ is a single-body operator, then

$\displaystyle\hat{H}_{\text{eff}}^{(2)}=\sum_{\Delta\neq
0}\dfrac{w^{\Delta}\cdot
w^{\Delta}}{N\left(N-1\right)\Delta}\left(\underline{X}^{2}-N\underline{X^{2}}\right),$
(60)

and if furthermore all $g_{ij}=-U/N$, as for $\hat{H}_{\text{int}}$ in Eq.
(9), then the only relevant excitation energy is $\Delta=U$ (see Section B.2),
and

$\displaystyle w^{U}\cdot
w^{U}=\sum_{i}\left(w_{i}-\overline{w}\right)^{2}=N\widetilde{w}^{2}$ (61)

is simply $N$ times the variance $\widetilde{w}^{2}$ of $w$, so

$\displaystyle\hat{H}_{\text{eff}}^{(2)}=\dfrac{\widetilde{w}^{2}}{\left(N-1\right)U}\left(\underline{X}^{2}-N\underline{X^{2}}\right).$
(62)

###  B.1 Generating excitation energy eigenstates

Here we construct the matrix $g^{(M)}$ that enables decomposing $M$-body
operators $\hat{\mathcal{O}}(w,\hat{X})$ into terms
$\hat{\mathcal{O}}(w^{\Delta},\hat{X})$ that generate states of definite
excitation energy $\Delta$ above the PS manifold, as in Eq. (58). We work
through the calculation of $g^{(1)}$ explicitly, and provide the result for
$g^{(M)}$ from a generalized version of the same calculation. To this end, we
consider the action of a single-body operator
$\hat{\mathcal{O}}(w,\hat{X})=\sum_{i}w_{i}\hat{X}_{i}$ on an arbitrary PS
state $\ket{\psi}\in\mathcal{E}_{0}$ and expand

$\displaystyle\hat{H}_{0}\hat{\mathcal{O}}(w,\hat{X})\ket{\psi}=\dfrac{1}{2}\sum_{i\neq
j}\sum_{k}g_{ij}w_{k}\hat{\Pi}_{ij}\hat{X}_{k}\ket{\psi},$ (63)

where strictly speaking $g_{ij}$ has only been defined for $i<j$, so for
completeness we define $g_{ji}=g_{ij}$ and $g_{ii}=0$. The sum in Eq. (63) has
terms with $k\in\left\\{i,j\right\\}$ and terms with
$k\notin\left\\{i,j\right\\}$. In the case of $k\notin\left\\{i,j\right\\}$,
the permutation operator $\hat{\Pi}_{ij}$ commutes with $\hat{X}_{k}$ and
annihilates on $\ket{\psi}$, and we can replace the sum

$\displaystyle\sum_{k\notin\left\\{i,j\right\\}}\to\sum_{k}-\sum_{k\in\left\\{i,j\right\\}},$
(64)

allowing us to simplify

$\displaystyle\dfrac{1}{2}\sum_{i\neq
j}\sum_{k\notin\left\\{i,j\right\\}}g_{ij}w_{k}\hat{\Pi}_{ij}\hat{X}_{k}\ket{\psi}=E_{0}\hat{\mathcal{O}}(w,\hat{X})\ket{\psi}-\dfrac{1}{2}\sum_{i\neq
j}\sum_{k\in\left\\{i,j\right\\}}g_{ij}w_{k}\hat{X}_{k}\ket{\psi},$ (65)

where $E_{0}=\frac{1}{2}\sum_{i\neq j}g_{ij}$ is the interaction energy the PS
state $\ket{\psi}\in\mathcal{E}_{0}$. Switching the order of sums over $i\neq
j$ and $k\in\left\\{i,j\right\\}$ as

$\displaystyle\sum_{i\neq
j}\sum_{k\in\left\\{i,j\right\\}}\to\sum_{k}\sum_{\begin{subarray}{c}i\neq
j\\\ \left\\{i,j\right\\}\ni k\end{subarray}},$ (66)

we can simplify

$\displaystyle\dfrac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
\left\\{i,j\right\\}\ni
k\end{subarray}}g_{ij}=\dfrac{1}{2}\sum_{i}g_{ik}+\dfrac{1}{2}\sum_{j}g_{kj}=g_{k},$
$\displaystyle g_{k}\equiv\sum_{i}g_{ik},$ (67)

which implies that the terms in Eq. (63) with $k\notin\left\\{i,j\right\\}$
are

$\displaystyle\dfrac{1}{2}\sum_{i\neq
j}\sum_{k\notin\left\\{i,j\right\\}}g_{ij}w_{k}\hat{\Pi}_{ij}\hat{X}_{k}\ket{\psi}=E_{0}\hat{\mathcal{O}}(w,\hat{X})\ket{\psi}-\sum_{k}g_{k}w_{k}\hat{X}_{k}\ket{\psi}.$
(68)

The terms in Eq. (63) with $k\in\left\\{i,j\right\\}$, meanwhile, are

$\displaystyle\dfrac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
k\in\left\\{i,j\right\\}\end{subarray}}g_{ij}w_{k}\hat{\Pi}_{ij}\hat{X}_{k}\ket{\psi}=\sum_{i,j}g_{ij}w_{j}\hat{X}_{i}\ket{\psi}.$
(69)

so in total

$\displaystyle\hat{H}_{0}\hat{\mathcal{O}}(w,\hat{X})\ket{\psi}=E_{0}\hat{\mathcal{O}}(w,\hat{X})\ket{\psi}+\sum_{k}\left[\sum_{j}g_{kj}w_{j}-g_{k}w_{k}\right]\hat{X}_{k}\ket{\psi}.$
(70)

The action of the single-body perturbation $\hat{\mathcal{O}}(w,\hat{X})$ on a
permutationally symmetric state therefore generates an eigenstate of
$\hat{H}_{0}$ with interaction energy $E_{0}+\Delta$ if the vector
$w=\sum_{k}w_{k}\ket{k}$ satisfies the eigenvalue equation

$\displaystyle g^{(1)}\cdot w=\Delta w,$ $\displaystyle g^{(1)}\equiv
g-\operatorname{diag}\vec{g},$ (71)

where $g\equiv\sum_{i,j}g_{ij}\outerproduct{i}{j}$ is a matrix of all
couplings $g_{ij}$; the vector
$\vec{g}\equiv\sum_{i,j}g_{ij}\ket{i}=\sum_{i}g_{i}\ket{i}$ is the sum of all
columns of $g$; and the matrix
$\operatorname{diag}\vec{g}\equiv\sum_{i}g_{i}\outerproduct{i}{i}$ has
$\vec{g}$ on the diagonal and zeroes everywhere else.

A similar calculation as above with arbitrary $M$ yields an eigenvalue
equation of the form

$\displaystyle g^{(M)}\cdot w=\Delta w,$ (72)

where we treat $w$ as an
$\absolutevalue{\mathcal{D}_{N}\left(M\right)}$-component vector, and
$g^{(M)}$ is a matrix with dimensions
$\absolutevalue{\mathcal{D}_{N}\left(M\right)}\times\absolutevalue{\mathcal{D}_{N}\left(M\right)}$.
In the case of $M=2$, we have

$\displaystyle
g^{(2)}=\sum_{\left(k,\ell\right)\in\mathcal{D}_{N}\left(2\right)}\ket{k\ell}\left[\sum_{\begin{subarray}{c}i\in\mathbb{Z}_{N}\\\
i\notin\left\\{k,\ell\right\\}\end{subarray}}\left(g_{ik}\bra{i\ell}+g_{i\ell}\bra{ki}\right)+g_{k\ell}\bra{\ell
k}-\left(g_{k}+g_{\ell}-g_{k\ell}\right)\bra{k\ell}\right],$ (73)

and more generally

$\displaystyle
g^{(M)}=\sum_{k\in\mathcal{D}_{N}\left(M\right)}\ket{k}\left[\sum_{a\in\mathbb{Z}_{M}}\sum_{\begin{subarray}{c}i\in\mathbb{Z}_{N}\\\
i\notin
k\end{subarray}}g_{ik_{a}}\bra{k_{a:i}}+\sum_{\left\\{a,b\right\\}\in\mathcal{C}_{M}\left(2\right)}g_{k_{a}k_{b}}\bra{k_{a\leftrightarrow
b}}-\tilde{g}_{k}\bra{k}\right],$ (74)

where $k_{a}\in k=\left(k_{1},k_{2},\cdots,k_{M}\right)$; $k_{a:i}$ a list
that is equal to $k$ except at the $a$-th position, where $k_{a}$ replaced is
by $i$, i.e. $k_{a:i}=\left(\cdots,k_{a-1},i,k_{a+1},\cdots\right)$;
$\mathcal{C}_{L}\left(p\right)$ is the set of all subsets (“choices”) of $p$
elements from $\mathbb{Z}_{L}$; $k_{a\leftrightarrow b}$ is equal to $k$
except at the $a$-th and $b$-th positions, at which $k_{a}$ and $k_{b}$ are
switched; and

$\displaystyle\tilde{g}_{k}\equiv\sum_{\begin{subarray}{c}\left\\{i,j\right\\}\in\mathcal{C}_{N}\left(2\right)\\\
i\in k~{}\text{or}~{}j\in k\end{subarray}}g_{ij}=\sum_{i\in
k}g_{i}-\sum_{\left\\{a,b\right\\}\in\mathcal{C}_{M}\left(2\right)}g_{k_{a}k_{b}}.$
(75)

If the tensor $w$ is permutationally symmetric, meaning that $w_{k}$ is
invariant under arbitrary permutations of $k$, then this symmetry is preserved
by $g^{(M)}$. In this case, we can replace sums over
$k\in\mathcal{D}_{N}\left(M\right)$ in Eqs. (73) and (74) by sums over
$k\in\mathcal{C}_{N}\left(M\right)$, and replace vectors
$\ket{k_{1},k_{2},\cdots,k_{M}}\to\ket{\left\\{k_{1},k_{2},\cdots,k_{M}\right\\}}$,
such that e.g. $\ket{k_{a\leftrightarrow b}}=\ket{k}$. These replacements
reduce the size of $g^{(M)}$ from
$\absolutevalue{\mathcal{D}_{N}\left(M\right)}\times\absolutevalue{\mathcal{D}_{N}\left(M\right)}$
to
$\absolutevalue{\mathcal{C}_{N}\left(M\right)}\times\absolutevalue{\mathcal{C}_{N}\left(M\right)}$,
where
$\absolutevalue{\mathcal{D}_{N}\left(M\right)}=\prod_{j=0}^{M-1}\left(N-j\right)=M!\times{N\choose
M}$ and $\absolutevalue{\mathcal{C}_{N}\left(M\right)}={N\choose M}$.
Additional symmetries of $g$ and $w$, such as translational invariance or
lattice symmetries, can be used to further reduce the computational complexity
of the eigenvalue problem in Eq. (72).

###  B.2 Recovering spin-wave theory

If the interaction Hamiltonian $\hat{H}_{0}$ is translationally invariant,
then the single-body eigenvalue problem in Eq. (71) is solvable analytically.
In this case, the couplings $g_{ij}$ depend only on the separation
$\absolutevalue{i-j}$, so eigenvectors of $g$ are plane waves of the form

$\displaystyle w_{k}\equiv\sum_{d\in\mathbb{Z}_{L}^{D}}e^{\mathrm{i}\mkern
1.0mud\cdot k}\ket{d},$ (76)

where on a $D$-dimensional periodic lattice of $N=L^{D}$ spins, lattice sites
are indexed by vectors $d\in\mathbb{Z}_{L}^{D}$, and wavenumbers take on
values $k\in\mathbb{Z}_{L}^{D}\times 2\pi/L$. The eigenvalues of $g$ can be
determined by expanding

$\displaystyle g\cdot
w_{k}=\sum_{c,d\in\mathbb{Z}_{L}^{D}}g_{cd}e^{\mathrm{i}\mkern 1.0mud\cdot
k}\ket{c}=\sum_{c,d\in\mathbb{Z}_{L}^{D}}g_{c,c+d}e^{\mathrm{i}\mkern
1.0mu\left(c+d\right)\cdot
k}\ket{c}=\sum_{d\in\mathbb{Z}_{L}^{D}}g_{0,d}\cos\left(d\cdot k\right)w_{k},$
(77)

where the imaginary contributions vanish in the sum over $d$ because
$g_{0,d}=g_{0,-d}$. The remainder of Eq. (71) that we need to sort out is
$\operatorname{diag}\vec{g}$, where all
$g_{i}=\sum_{i,j}g_{ij}=\sum_{d}g_{0,d}$ are equal, which implies that
$\operatorname{diag}\vec{g}=\sum_{d}g_{0,d}$ is a scalar. We thus find that

$\displaystyle g^{(1)}\cdot w_{k}=\Delta_{k}w_{k},$
$\displaystyle\Delta_{k}\equiv\sum_{d\in\mathbb{Z}_{L}^{D}}g_{0,d}\left[\cos\left(d\cdot
k\right)-1\right],$ (78)

in agreement with standard spin-wave theory. Excitations generated by the
action of $\hat{\mathcal{O}}\left(w_{k},X\right)$ on PS states
$\ket{\psi}\in\mathcal{E}_{0}$ are known as spin-waves. If $g_{ij}=-U/N$ is
constant, then the spin-wave excitation energies are $\Delta_{k}=U$
independent of the wavenumber $k$.

##  Appendix C Restricting spin operators to the permutationally symmetric
manifold

Here we provide the restriction of a general $M$-body spin operator
$\hat{\mathcal{O}}$ to the permutationally symmetric (PS) manifold of $N$
spins (each with $n$ internal states). Denoting the projector onto the PS
manifold by $\hat{\mathcal{P}}_{0}$, our task is essentially to find the
coefficients of the expansion

$\displaystyle\hat{\mathcal{P}}_{0}\hat{\mathcal{O}}_{M}\hat{\mathcal{P}}_{0}=\sum_{a,b\in\mathcal{A}_{n}\left(N\right)}\braket{a}{\hat{\mathcal{O}}_{M}}{b}\outerproduct{a}{b},$
(79)

where $\mathcal{A}_{n}\left(N\right)$ is the set of all ways to assign $N$
(identical) spins to $n$ (distinct) states, such that for any
$a\in\mathcal{A}_{n}\left(N\right)$ the state
$\ket{a}=\ket{a_{1},a_{2},\cdots,a_{n}}$ is labeled by the occupation number
$a_{\mu}$ of state $\mu$, with $\sum_{\mu}a_{\mu}=N$. Written out explicitly,

$\displaystyle\ket{a}=\dfrac{1}{\sqrt{\mathcal{C}\left(a\right)}}\sum_{\begin{subarray}{c}\text{distinct}\\\
\text{permutations}\\\
\hat{\Pi}\,\text{of}\,\tilde{a}\end{subarray}}\hat{\Pi}\ket{\tilde{a}},$
$\displaystyle\ket{\tilde{a}}\equiv\bigotimes_{\mu}\ket{\mu}^{\otimes
a_{\mu}},$
$\displaystyle\mathcal{C}\left(a\right)\equiv\dfrac{\left(\sum_{\mu}a_{\mu}\right)!}{\prod_{\nu}a_{\nu}!}.$
(80)

Here $\mathcal{C}\left(a\right)$ is a multinomial coefficient that counts the
number of distinct ways to permute the tensor factors of the “standard-
ordered” state $\ket{\tilde{a}}$, enforcing $\braket{a}{a}=1$. Using these
states, with some combinatorics we can expand

$\displaystyle\braket{a}{\hat{\mathcal{O}}_{M}}{b}=\sum_{\begin{subarray}{c}\alpha,\beta\in\mathcal{A}_{n}\left(M\right)\\\
\alpha\leq a,\beta\leq
b\end{subarray}}\delta_{a-\alpha,b-\beta}\sqrt{\dfrac{\mathcal{C}\left(\alpha\right)\mathcal{C}\left(a-\alpha\right)\mathcal{C}\left(\beta\right)\mathcal{C}\left(b-\beta\right)}{\mathcal{C}\left(a\right)\mathcal{C}\left(b\right)}}\,\braket{\alpha}{\hat{\mathcal{O}}_{M}}{\beta},$
(81)

where the restriction $\alpha\leq a$ and the difference $a-\alpha$ are
evaluated element-wise, i.e. $\alpha\leq a\implies\alpha_{\mu}\leq a_{\mu}$
and $\left(a-\alpha\right)_{\mu}=a_{\mu}-\alpha_{\mu}$ for all $\mu$; and
$\delta_{cd}=1$ if $c=d$ and zero otherwise. We sum over both $\alpha$ and
$\beta$ above merely to keep the expression symmetric with respect to
transposition $\left(a,\alpha\right)\leftrightarrow\left(b,\beta\right)$; in
practice, one can simply sum over $\alpha\in\mathcal{A}_{n}\left(M\right)$ and
set $\beta=b-a+\alpha$, throwing out terms with any $\beta_{\mu}<0$. Note
that, by slight abuse of notation, the operator $\hat{\mathcal{O}}_{M}$ on the
left of Eq. (81) acts on an arbitrary choice of $M$ spins (out of $N$),
whereas the operator $\hat{\mathcal{O}}_{M}$ on the right of Eq. (81) is
simply an $M$-spin operator, with matrix elements
$\braket{\alpha}{\hat{\mathcal{O}}_{M}}{\beta}$ evaluated with respect to the
PS $M$-spin states $\ket{\alpha},\ket{\beta}\in\mathcal{A}_{n}\left(M\right)$.

##  Appendix D Relaxing assumptions of the three-laser drive

In order to arrive at the drive Hamiltonian in Eq. (12) of the main text, we
made two simplifying assumptions:  1. that the excited-state hyperfine
manifold had the same total spin $s$ as the ground-state manifold, and  2.
that all drive amplitudes are real (which enforces a phase-locking condition
between the driving lasers).  To derive an effective drive Hamiltonian for the
general case in which the excited-state hyperfine manifold has total spin
$s+r$ with $r\in\left\\{+1,0,-1\right\\}$, we decompose all lasers into their
right- and left-circular polarization components and write the full drive
Hamiltonian in the form

$\displaystyle\hat{H}_{\text{drive}}^{\text{full}}=\sum_{j,\bm{v},\sigma}\Omega_{\bm{v}\sigma}\left(e^{-\mathrm{i}\mkern
1.0mu\kappa\bm{v}\cdot\bm{\ell}j}\hat{s}_{\bm{v}\sigma
j}^{(r)}\otimes\outerproduct{\text{e}}{\text{g}}_{j}+\text{h.c.}\right)+\Delta\hat{N}_{\text{e}},$
(82)

where $\Omega_{\bm{v}\sigma}$ is the amplitude of $\sigma$-polarized light
propagating along axis $\bm{v}$, with $\sigma=+1$ and $-1$ respectively for
right and left circular polarizations; and $\hat{s}_{\bm{v}\sigma j}$ is a
spin-raising/lowering operator for atom $j$ along axis $\bm{v}$, defined by
appropriately rotating the single-atom spin operators

$\displaystyle\hat{s}_{\pm}^{(r)}\equiv-\sqrt{\dfrac{n(n+1)(n-1)}{6}}\times\hat{T}_{\pm}^{(r)},$
$\displaystyle\hat{T}_{\pm}^{(r)}\equiv\mp\sqrt{\dfrac{2(s+r)+1}{2\ell+1}}\sum_{\mu}\braket{s\mu;1,\pm
1}{s+r,\mu\pm 1}\outerproduct{\mu\pm 1}{\mu}.$ (83)

Here $\braket{j_{1}m_{1};j_{2}m_{2}}{j_{3}m_{3}}$ is a Clebsch-Gordan
coefficient, and we have normalized $\hat{T}_{\pm}^{(r)}$ such that
$\tr\left[\hat{T}_{\pm}^{(r)}{}^{\dagger}\hat{T}_{\pm}^{(r)}\right]=1$. Still
assuming real drive amplitudes, the corresponding effective drive Hamiltonian
that replaces Eq. (12) in the far-detuned limit
$\absolutevalue{\Delta}\gg\absolutevalue{\Omega_{\bm{v}\sigma}}$ is then

$\displaystyle\hat{H}_{\text{3LD}}^{\text{single}}=f_{r}^{(1)}\left[\tilde{\Omega}_{+}\tilde{\Omega}_{-}\hat{s}_{\text{z}}+\tilde{\Omega}_{0}\tilde{\Omega}_{-}\hat{s}_{\text{x}}\right]+f_{r}^{(2)}\left[\tilde{\Omega}_{0}\tilde{\Omega}_{+}(\hat{s}_{\text{z}}\hat{s}_{\text{x}}+\hat{s}_{\text{x}}\hat{s}_{\text{z}})-\left(\tilde{\Omega}_{0}^{2}\hat{s}_{\text{z}}^{2}+\tilde{\Omega}_{+}^{2}\hat{s}_{\text{x}}^{2}+\tilde{\Omega}_{-}^{2}\hat{s}_{\text{y}}^{2}\right)\right]-f_{r}^{(3)}\sum_{m}\tilde{\Omega}_{m}^{2},$
(84)

where $f_{r}^{(k)}$ are scalars that depend on the spin dimension $n$:

$\displaystyle f_{0}^{(1)}$ $\displaystyle=1,$ $\displaystyle f_{+1}^{(1)}$
$\displaystyle=-s,$ $\displaystyle f_{-1}^{(1)}$ $\displaystyle=s+1,$ (85)
$\displaystyle f_{0}^{(2)}$ $\displaystyle=1,$ $\displaystyle f_{+1}^{(2)}$
$\displaystyle=-\dfrac{s}{n+2},$ $\displaystyle f_{-1}^{(2)}$
$\displaystyle=-\dfrac{s+1}{n-2},$ (86) $\displaystyle f_{0}^{(3)}$
$\displaystyle=0,$ $\displaystyle f_{+1}^{(3)}$
$\displaystyle=\dfrac{s(s+1)^{2}}{n+2},$ $\displaystyle f_{-1}^{(3)}$
$\displaystyle=\dfrac{s^{2}(s+1)}{n-2}.$ (87)

If additionally the drive amplitudes are complex,
$\Omega_{m}\to\Omega_{m}e^{-\mathrm{i}\mkern 1.0mu\eta_{m}}$ (with real
$\Omega_{m},\eta_{m}$), then

$\hat{H}_{\text{3LD}}^{\text{single}}=f_{r}^{(1)}\tilde{\Omega}_{+}\tilde{\Omega}_{-}\hat{s}_{\text{z}}+\tilde{\Omega}_{0}\sum_{\sigma\in\left\\{\pm
1\right\\}}\dfrac{\tilde{\Omega}_{+}+\sigma\tilde{\Omega}_{-}}{2}\left[f_{r}^{(1)}\sigma\hat{s}_{\tilde{\eta}_{\sigma},\text{x}}+f_{r}^{(2)}\left(\hat{s}_{\text{z}}\hat{s}_{\tilde{\eta}_{\sigma},\text{x}}+\hat{s}_{\tilde{\eta}_{\sigma},\text{x}}\hat{s}_{\text{z}}\right)\right]\\\
-f_{r}^{(2)}\left[\tilde{\Omega}_{0}^{2}\hat{s}_{\text{z}}^{2}+\tilde{\Omega}_{+}\hat{s}_{\tilde{\eta}_{0},\text{x}}^{2}+\tilde{\Omega}_{-}\hat{s}_{\tilde{\eta}_{0},\text{y}}^{2}\right]-f_{r}^{(3)}\sum_{m}\tilde{\Omega}_{m}^{2},$
(88)

where $\hat{s}_{\eta\alpha}\equiv e^{-\mathrm{i}\mkern
1.0mu\eta\hat{s}_{\text{z}}}\hat{s}_{\alpha}e^{\mathrm{i}\mkern
1.0mu\eta\hat{s}_{\text{z}}}$ is a rotated spin-$\alpha$ operator (e.g.
$\hat{s}_{\pi/2,\text{x}}=\hat{s}_{\text{y}}$), and

$\displaystyle\tilde{\eta}_{\pm}\equiv\pm\left(\eta_{\pm}-\eta_{0}\right),$
$\displaystyle\tilde{\eta}_{0}\equiv\dfrac{\eta_{+}-\eta_{-}}{2},$ (89)

are the relative phases of the drive amplitudes.

##  Appendix E Mean-field theory

Here we describe the mean-field theory used to simulate the spin Hamiltonian

$\displaystyle\hat{H}_{\text{spin}}=-\dfrac{u}{2N}\bm{\hat{S}}\cdot\bm{\hat{S}}+2J\phi\sum_{q}\sin\left(q\right)\hat{s}_{\text{z},q}$
(90)

in Eq. (23) of the main text. We begin by decomposing individual spin
operators into Schwinger bosons as $\hat{s}_{\mu\nu q}=\hat{b}_{\mu
q}^{\dagger}\hat{b}_{\nu q}$, such that the spin Hamiltonian becomes

$\displaystyle\hat{H}_{\text{spin}}\to\hat{H}_{\text{boson}}=-\dfrac{u}{2N}\sum_{p,q,\mu,\nu}\hat{b}_{\mu
p}^{\dagger}\hat{b}_{\nu p}\hat{b}_{\nu q}^{\dagger}\hat{b}_{\mu
q}+2J\phi\sum_{q,\mu}\sin\left(q\right)\mu\,\hat{b}_{\mu
q}^{\dagger}\hat{b}_{\mu q}.$ (91)

The Heisenberg equations of motion for the Schwinger boson operators are (see
Appendix F)

$\displaystyle\mathrm{i}\mkern 1.0mu\partial_{t}\hat{b}_{\mu
q}=-\dfrac{u}{N}\sum_{\nu,p}\hat{b}_{\nu p}^{\dagger}\hat{b}_{\mu
p}\hat{b}_{\nu q}+2J\phi\sin\left(q\right)\mu\,\hat{b}_{\mu q}.$ (92)

Our mean-field theory then treats all boson operators in these equations of
motion as complex numbers, $\hat{b}_{\mu q}\to\braket{\hat{b}_{\mu
q}}_{\text{MF}}$, with the initial value $\braket{\hat{b}_{\mu
q}\left(t=0\right)}_{\text{MF}}$ equal to the initial amplitude of spin $q$ in
state $\mu$. Specifically, for an $N$-fold product state of the form
$\ket{\psi}=\bigotimes_{q}\sum_{\mu}\psi_{\mu q}\ket{\mu}$ we set
$\braket{\hat{b}_{\mu q}\left(t=0\right)}_{\text{MF}}=\psi_{\mu q}$. For pure
initial product states, this mean-field treatment of the boson operators
$\hat{b}_{\mu q}$ is mathematically equivalent to a mean-field treatment of
the spin operators $\hat{s}_{\mu\nu q}$, as in Eq. (24), but reduces the
number of variables to keep track of by a factor of $\sim n$.

##  Appendix F Schwinger boson equations of motion for quadratic spin
Hamiltonians

Here we decompose a quadratic spin Hamiltonian into Schwinger bosons, and
derive the equations of motion for the resulting boson operators. We begin
with a general spin Hamiltonian of the form

$\displaystyle\hat{H}=\sum_{\begin{subarray}{c}\mu,\nu,\rho,\sigma\\\
j<k\end{subarray}}g^{\mu\nu j}_{\rho\sigma k}\hat{s}_{\mu\nu
j}\hat{s}_{\rho\sigma k}+\sum_{\mu,\nu,j}\epsilon_{\mu\nu j}\hat{s}_{\mu\nu
j},$ (93)

where $\mu,\nu$ index orthogonal states of an $n$-level spin; $j,k$ index one
of $N$ spins; $g^{\mu\nu j}_{\rho\sigma k}$ and $\epsilon_{\mu\nu j}$ are
scalars; and $\hat{s}_{\mu\nu j}=\outerproduct{\mu}{\nu}_{j}$ is a transition
operator for spin $j$. Strictly speaking, Eq. (93) only defines the couplings
$g^{\mu\nu j}_{\rho\sigma k}$ for $j<k$, so we enforce $g^{\mu\nu
k}_{\rho\sigma j}=g^{\mu\nu j}_{\rho\sigma k}$ and $g^{\mu\nu j}_{\rho\sigma
j}=0$ for completion. Decomposing spin operators into Schwinger bosons as
$\hat{s}_{\mu\nu j}=\hat{b}_{\mu j}^{\dagger}\hat{b}_{\nu j}$, where
$\hat{b}_{\nu j}$ a annihilates a boson of type $\nu$ on site $j$, we can
write this Hamiltonian as

$\displaystyle\hat{H}=\sum_{\begin{subarray}{c}\mu,\nu,\rho,\sigma\\\
j<k\end{subarray}}g^{\mu\nu j}_{\rho\sigma k}\hat{b}_{\mu
j}^{\dagger}\hat{b}_{\nu j}\hat{b}_{\rho k}^{\dagger}\hat{b}_{\sigma
k}+\sum_{\mu,\nu,j}\epsilon_{\mu\nu j}\hat{b}_{\mu j}^{\dagger}\hat{b}_{\nu
j}.$ (94)

The Heisenberg equations of motion for the boson operators are then

$\displaystyle\mathrm{i}\mkern
1.0mu\partial_{t}\hat{b}_{\alpha\ell}=\left[\hat{b}_{\alpha\ell},\hat{H}\right]$
$\displaystyle=\sum_{\begin{subarray}{c}\mu,\nu,\rho,\sigma\\\
j<k\end{subarray}}g^{\mu\nu j}_{\rho\sigma
k}\left[\hat{b}_{\alpha\ell},\hat{b}_{\mu j}^{\dagger}\hat{b}_{\nu
j}\hat{b}_{\rho k}^{\dagger}\hat{b}_{\sigma
k}\right]+\sum_{\mu,\nu,j}\epsilon_{\mu\nu
j}\left[\hat{b}_{\alpha\ell},\hat{b}_{\mu j}^{\dagger}\hat{b}_{\nu j}\right]$
(95) $\displaystyle=\sum_{\mu,\nu,\rho,\sigma,k}g^{\mu\nu\ell}_{\rho\sigma
k}\left[\hat{b}_{\alpha\ell},\hat{b}_{\mu\ell}^{\dagger}\hat{b}_{\nu\ell}\right]\hat{b}_{\rho
k}^{\dagger}\hat{b}_{\sigma
k}+\sum_{\mu,\nu}\epsilon_{\mu\nu\ell}\left[\hat{b}_{\alpha\ell},\hat{b}_{\mu\ell}^{\dagger}\hat{b}_{\nu\ell}\right]$
(96)
$\displaystyle=\sum_{\mu,\nu}\left(\sum_{\rho,\sigma,k}g^{\mu\nu\ell}_{\rho\sigma
k}\hat{b}_{\rho k}^{\dagger}\hat{b}_{\sigma
k}+\epsilon_{\mu\nu\ell}\right)\left[\hat{b}_{\alpha\ell},\hat{b}_{\mu\ell}^{\dagger}\hat{b}_{\nu\ell}\right]$
(97)

where

$\displaystyle\left[\hat{b}_{\alpha\ell},\hat{b}_{\mu\ell}^{\dagger}\hat{b}_{\nu\ell}\right]=\delta_{\alpha\mu}\delta_{\alpha\nu}\hat{b}_{\alpha\ell}+\delta_{\alpha\mu}\left(1-\delta_{\alpha\nu}\right)\hat{b}_{\nu\ell}=\delta_{\alpha\mu}\hat{b}_{\nu\ell},$
(98)

so

$\displaystyle\mathrm{i}\mkern
1.0mu\partial_{t}\hat{b}_{\alpha\ell}=\sum_{\nu}\left(\sum_{\rho,\sigma,k}g^{\alpha\nu\ell}_{\rho\sigma
k}\hat{b}_{\rho k}^{\dagger}\hat{b}_{\sigma
k}+\epsilon_{\alpha\nu\ell}\right)\hat{b}_{\nu\ell}.$ (99)

In the case of uniform SU($n$)-symmetric interactions of the form
$\frac{g}{2}\bm{\hat{S}}\cdot\bm{\hat{S}}$ and a diagonal external field, we
have

$\displaystyle g^{\alpha\nu\ell}_{\rho\sigma
k}=g\times\delta_{\alpha\sigma}\delta_{\nu\rho},$
$\displaystyle\epsilon_{\alpha\nu\ell}=\epsilon_{\alpha\ell}\times\delta_{\alpha\nu}$
(100)

so

$\displaystyle\mathrm{i}\mkern
1.0mu\partial_{t}\hat{b}_{\alpha\ell}=g\sum_{\nu,k}\hat{b}_{\nu
k}^{\dagger}\hat{b}_{\alpha
k}\hat{b}_{\nu\ell}+\epsilon_{\alpha\ell}\hat{b}_{\alpha\ell}.$ (101)

##  Appendix G Lax vector analysis

We start with the spin Hamiltonian

$\displaystyle\hat{H}_{\text{spin}}=-\dfrac{u}{2N}\sum_{\mu,\nu}\hat{S}_{\mu\nu}\hat{S}_{\nu\mu}+2J\phi\sum_{q}\sin\left(q\right)\hat{s}_{\text{z},q},$
(102)

where $\hat{S}_{\mu\nu}=\sum_{q}\hat{s}_{\mu\nu q}$. The single-body operators
that appear in this Hamiltonian have squared norms

$\displaystyle\tr\left(\hat{s}_{\mu\nu q}^{\dagger}\hat{s}_{\mu\nu
q}\right)=1$ and
$\displaystyle\tr\left(\hat{s}_{\text{z},q}^{\dagger}\hat{s}_{\text{z},q}\right)=\sum_{\mu}\mu^{2}=\dfrac{1}{12}(n+1)n(n-1)\equiv\xi^{2}.$
(103)

The Lax formulation (following Refs. [55, 56, 57, 58, 45]) requires all
single-body operators involved to have the same normalization, so we
substitute $\hat{s}_{\tilde{\text{z}},q}\equiv\hat{s}_{\text{z},q}/\xi$ to
expand

$\displaystyle\dfrac{\hat{H}_{\text{spin}}}{u}=-\dfrac{1}{2N}\sum_{\mu,\nu}\hat{S}_{\mu\nu}\hat{S}_{\nu\mu}+\xi
h\sum_{q}\sin\left(q\right)\hat{s}_{\tilde{\text{z}},q},$ where $\displaystyle
h\equiv\dfrac{2J\phi}{u}.$ (104)

The intensive, dimensionless, $(n^{2}-1)$-component Lax vector
$\vec{\ell}\left(z\right)$ associated with $\hat{H}_{\text{spin}}$, which is
defined with an auxiliary complex parameter $z$, has components

$\displaystyle\ell_{\alpha}\left(z\right)=\dfrac{1}{N}\sum_{q}\dfrac{\hat{s}_{\alpha,q}}{z-\sin
q}+\delta_{\alpha,\tilde{\text{z}}}\,\xi h,$ (105)

where $\alpha$ indexes elements of a basis $\left\\{\hat{s}_{\alpha}\right\\}$
of self-adjoint generators of SU($n$), with normalization
$\tr\left(\hat{s}_{\alpha}^{2}\right)=1$. The squared magnitude
$\vec{\ell}\left(z\right)^{2}=\sum_{\alpha}\ell_{\alpha}\left(z\right)^{2}$ is
a constant of motion (for any $z$), and its residues provide $N$ mutually
commuting quantities whose weighted sum recovers $\hat{H}_{\text{spin}}$. When
$n=2$, conservation of these residues provides sufficient dynamical
constraints to make the spin system fully integrable. In this case, dynamical
behavior is governed by the roots of $\vec{\ell}\left(z\right)^{2}$, and the
presence (or absence) of complex roots marks distinct dynamical phases of
$\hat{H}_{\text{spin}}$. However, the size of Hilbert space grows with $n$,
while the number of conserved quantities provided by the Lax analysis (namely,
$N$) does not. When $n>2$, there is therefore no guarantee that the roots of
$\vec{\ell}\left(z\right)^{2}$ will similarly govern dynamical behavior. In
fact, a straightforward generalization of the Lax analysis to $n>2$ makes
predictions that are inconsistent with the mean-field results in Figures 5–8
of the main text. We substantiate this claim with a direct calculation of the
roots of $\vec{\ell}\left(z\right)^{2}$ below.

Within the permutationally symmetric manifold, we can replace
$\hat{s}_{\alpha,q}\to\bar{s}_{\alpha}\equiv\frac{1}{N}\sum_{q}\hat{s}_{\alpha,q}$
at the cost of $O(1/N)$ errors that vanish as $N\to\infty$, so taking this
limit we find

$\displaystyle\ell_{\alpha}\left(z\right)=\mathcal{I}\left(z\right)\bar{s}_{\alpha}+\delta_{\alpha,\tilde{\text{z}}}\,\xi
h,$ (106)

where

$\displaystyle\mathcal{I}\left(z\right)\equiv\lim_{N\to\infty}\dfrac{1}{N}\sum_{q}\dfrac{1}{z-\sin\left(q\right)}=\dfrac{1}{2\pi}\int_{0}^{2\pi}\dfrac{\text{d}q}{z-\sin\left(q\right)}=\dfrac{1}{\sqrt{z^{2}-1}}$
for $\displaystyle z\notin\left[-1,1\right].$ (107)

The squared magnitude of the Lax vector is therefore

$\displaystyle\vec{\ell}\left(z\right)^{2}=\sum_{\alpha}\ell_{\alpha}\left(z\right)^{2}=\mathcal{I}\left(z\right)^{2}\sum_{\alpha\neq\tilde{\text{z}}}\bar{s}_{\alpha}^{2}+\left[\mathcal{I}\left(z\right)\bar{s}_{\tilde{\text{z}}}+\xi
h\right]^{2},$ (108)

where we can define the scalar $Q^{2}\equiv\sum_{\alpha}\bar{s}_{\alpha}^{2}$
to simplify

$\displaystyle\vec{\ell}\left(z\right)^{2}=\mathcal{I}\left(z\right)^{2}\left(Q^{2}-\bar{s}_{\tilde{\text{z}}}^{2}\right)+\left[\mathcal{I}\left(z\right)\bar{s}_{\tilde{\text{z}}}+\xi
h\right]^{2}=\mathcal{I}\left(z\right)^{2}Q^{2}+\xi^{2}h^{2}+2\mathcal{I}\left(z\right)\xi
h\bar{s}_{\tilde{\text{z}}}.$ (109)

For initial states with $\braket{\bar{s}_{\text{z}}}=0$, we thus find that

$\displaystyle\vec{\ell}\left(z\right)^{2}=\dfrac{Q^{2}}{z^{2}-1}+\xi^{2}h^{2},$
(110)

which is zero whenhhhStrictly speaking, the zeros in Eq. (111) occur at values
of $z$ at which $\mathcal{I}\left(z\right)$ is undefined. We avoid this issue
by analytically continuing $\mathcal{I}\left(z\right)^{2}$ to the interval
$z\in(-1,1)$.

$\displaystyle z=\pm\sqrt{1-\left(\dfrac{Q}{\xi h}\right)^{2}}.$ (111)

These roots change character when $z=0$, suggesting that the critical field
$h_{\text{crit}}$ separating dynamical phases satisfies

$\displaystyle
h_{\text{crit}}^{2}\stackrel{{\scriptstyle?}}{{=}}\dfrac{Q^{2}}{\xi^{2}},$
(112)

where we use the relation $\stackrel{{\scriptstyle?}}{{=}}$ to indicate that
this “prediction” of the Lax analysis is not necessarily valid for all $n$.
For a permutationally symmetric state, up to vanishing $O(1/N)$ corrections we
can expand

$\displaystyle
Q^{2}=\sum_{\alpha}\bar{s}_{\alpha}^{2}=\sum_{\mu,\nu}\bar{s}_{\mu\nu}\bar{s}_{\nu\mu}-\dfrac{1}{n}=1-\dfrac{1}{n}=\dfrac{n-1}{n},$
(113)

which implies that

$\displaystyle
h_{\text{crit}}^{2}\stackrel{{\scriptstyle?}}{{=}}\dfrac{n-1}{n}\times\dfrac{12}{n(n+1)(n-1)}=\dfrac{12}{n^{2}\left(n+1\right)}.$
(114)

This Lax analysis correctly predicts that $h_{\text{crit}}=1$ when $n=2$, but
otherwise predicts $h_{\text{crit}}\sim n^{-3/2}$, which is inconsistent with
the finding that $h_{\text{crit}}\sim n^{-1/3}$ in the mean-field results of
the main text (see Figure 6). We emphasize that this inconsistency is not a
failure of the Lax formalism, but rather an indication that new theoretical
tools are necessary to understand multilevel spin models.