Tunable spin model generation with spin-orbital coupled fermions in optical lattices

We study the dynamical behaviour of ultracold fermionic atoms loaded into an optical lattice under the presence of an effective magnetic flux, induced by spin-orbit coupled laser driving. At half filling, the resulting system can emulate a variety of iconic spin-1/2 models such as an Ising model, an XY model, a generic XXZ model with arbitrary anisotropy, or a collective one-axis twisting model. The validity of these different spin models is examined across the parameter space of flux and driving strength. In addition, there is a parameter regime where the system exhibits chiral, persistent features in the long-time dynamics. We explore these properties and discuss the role played by the system's symmetries. We also discuss experimentally-viable implementations.


I. INTRODUCTION
Understanding and quantifying the behaviour of interacting quantum particles in lattices is a fundamental goal of modern quantum science. While there is a plethora of research directions, one vital aspect is the response of particles to externallyimposed magnetic fields. Such fields induce an effective flux that threads through the plaquettes of the lattice [1,2], coupling the charge and spin degrees of freedom and modifying the particle dynamics. Interpreting the dynamical response to an applied flux is important for many applications in condensed matter, including ferroelectrics [3], spintronics [4,5] and spin-glass physics [6,7].
One of the best ways to study such phenomena is with ultracold atomic experiments. State-of-the-art ultracold systems provide exceptional levels of cleanliness, isolation and tunability. They allow for pristine implementations of iconic Fermi-or Bose-Hubbard models that describe interacting particles in a lattice, with additional terms to account for the synthetic magnetic fields [8]. A magnetic flux is easy to impose and control with tools such as laser driving, using Raman couplings or direct optical transitions.
There has been a great deal of theoretical work on Fermior Bose-Hubbard models with synthetic gauge fields that ultracold experiments could investigate, exploring ground-state phases [9][10][11] or phenomena such as many-body localization [12]. However, the interplay of magnetic flux together with particle interactions can lead to complex dynamical behaviour that still lacks a good theoretical understanding. In the case of fermions, even non-interacting atoms can exhibit nontrivial behaviour due to Pauli exclusion, while the addition of Hubbard repulsion renders the dynamics even more complex.
In this work, we focus on fermionic atoms loaded into 3D optical lattices. In the Mott insulating limit with one atom per site, each atom acts as an effective spin and the physics can be simplified by mapping to an effective spin model. This spin model emerges via virtual atomic tunneling processes that lead to second-order superexchange interactions between * mikhail.mamaev@colorado.edu these spins. Conventional lattice systems are often captured by isotropic Heisenberg models, while the presence of magnetic flux allows for tunability of the spin interactions. As we show here, the flux leads to more elaborate spin models such as anisotropic XXZ models, collective one-axis twisting Hamiltonians [13], or Dzyaloshinskii-Moriya (DM) interactions [14,15]. The latter is particularly intriguing, as it exhibits chirality, topological features and complex phase diagrams depending on the flux [16,17]. A system that can realize several different spin models with easy tunability can be very useful to the field, as current-generation optical lattice experiments have only recently begun to probe anisotropic interacting spin physics [18,19].
Here we study the case when the internal atomic states are driven by an external laser, which imprints a site-dependent phase that emulates a magnetic flux. We use the drive strength and the magnetic flux as tuning parameters, and show the different types of spin interactions that can be realized. We explore the dynamical properties of these interactions and specify regimes in the parameter space of flux and driving strength where the system's time evolution can be captured by simple models such as Ising, XY, XXZ or collective-spin one-axis twisting. We also study a regime where the corresponding spin model maps to a Heisenberg model in a twisted frame, causing the long-time dynamics to develop non-trivial features such as infinite-time magnetization and chiral spin imbalance. These latter features are not limited to the strongly interacting regime, but also hold in the weakly interacting limit of the Fermi-Hubbard model where atomic motion is relevant. We discuss the role that symmetry plays in preventing the system from relaxing. Our predictions can be readily implemented in many ultracold systems, and are especially relevant for alkaline-earth or earth-like atoms in 3D lattices or tweezer arrays, which provide exceptional coherence times [20][21][22] to overcome the inherent slow interaction rates while avoiding issues from heating given the very low spontaneous emission rates of their low-lying electronic levels.
The outline of the paper is as follows. Section II introduces the underlying Fermi-Hubbard model and derives the corresponding effective spin models that emerge at half filling. A dynamical classification of the spin model behaviour in different parameter regimes is given. Section III focuses on the low drive regime, and discusses the persistent magnetization or chiral features that can be observed due to the additional symmetry of a Heisenberg model in a twisted frame. Section IV provides a detailed discussion on possible experimental implementations.

II. TUNABLE SPIN DYNAMICS
A. Fermi-Hubbard and spin model The system we describe is a three-dimensional (3D) optical lattice loaded with ultracold fermionic atoms, each possessing two internal states σ ∈ {g, e} corresponding to a spin-1/2 degree of freedom. The effective system dimensionality is freely tuned by changing the lattice confinement strengths. We focus on the case where the system is effectively 1D by considering a strong confinement in two directions that suppresses tunneling along them. Along the direction atoms can tunnel, we assume a lattice of L sites populated by N atoms. The system is depicted in Fig. 1(a). The Hamiltonian of the system is, Hereĉ j,σ annihilates an atom with spin σ on site j. Atoms tunnel at rate J and exhibit onsite repulsion of strength U . In addition to the standard Fermi-Hubbard termĤ FH , we include a laser-driving termĤ Ω that induces the desired flux φ through a spatially-dependent phase e ijφ , providing model tunability. The differential phase imprinted by the laser implements a net spin-orbit coupling (SOC) by generating a momentum kick to the atoms while flipping their spin [1,2,[23][24][25][26]. The realization of this effective flux in 1D is depicted in Fig. 1(b). The drive can be implemented with a direct interrogating laser (if the internal states g, e are split by optical frequency), or a Raman coupling (see Section IV for details). We assume φ ∈ [0, π] without loss of generality. For sufficiently strong interactions U/J 1 at half filling N/L = 1 with one atom per site, double occupancies of lattice sites are strongly suppressed. The system dynamics in this regime may be approximated with a spin model, as depicted in Fig. 1(c). We define a dressed spin basis at different lattice sites by making a gauge transformation, defining new fermionic operators {â j,↑ ,â j,↓ } to remove the SOC phase, and define conventional spin operators for these dressed  1. (a) Schematic of Fermi-Hubbard dynamics in an optical lattice, with nearest-neighbour tunneling rate J, onsite repulsion U , and on-site driving Ω. When the laser drive has a wavelength that is incommensurate with the underlying optical lattice wavelength, the laser induces spin-orbit coupling (SOC) through a spatiallydependent phase e ijφ in the drive term, with φ controlled by the drive implementation (e.g. laser alignment or lattice spacing, see Section IV). (b) The SOC phase creates effective magnetic flux on plaquettes of the ladder state structure in 1D, with lattice index j along the length of the ladder and internal states e, g corresponding to the individual rungs. An atom tunneling around one plaquette picks up a total phase of φ. (c) At half filling and strong interactions U/J 1, the atomic spin σ at different lattice site can be dressed by the laser drive, which modifies the resulting spin dynamics dominated by second-order virtual superexchange processes with rate ∼ J 2 /U , up to possible normalization from the drive Ω.
This is just the Abrikosov pseudo-fermion representation [27,28]. Standard second-order perturbation theory then leads to the following general superexchange (SE) spin model (see Appendix A for derivation), The first two terms correspond to an XXZ model. The third term is a DM interaction with a plane axis ofẑ [thus also taking the form ofẑ · ( σ j × σ j+1 )]. The interaction coefficients are, These share the conventional J 2 /U superexchange energy scale, with different normalization factors coming from the interplay between the drive and flux. This model is valid for all flux values φ, and all drive frequencies Ω far from the resonance point |Ω| = |U |, requiring a spacing of ||U | − |Ω|| J to prevent higher-order effects. Note that there can also be other resonances such as |Ω| = |U |/2, as considered in e.g. Ref. [29], though the relevant width is far smaller as will be seen in the next section. These spin interactions are anisotropic even if the drive is turned off, Ω = 0, because we are using a dressed basis to absorb the SOC phase (see Appendix A). In the context of the underlying Fermi-Hubbard model, one can use a fast pulse of a SOC drive to prepare a desired initial state such as a product state in the dressed basis, after which the drive may be turned off if desired (see Section IV for details). There is also a single-particle term J Ω corresponding to the drive. This term contains an additional single-particle superexchange contribution, but this is typically negligible compared to the bare drive, and so we can approximate J Ω ≈ Ω/2.
If the drive is off, Ω = 0, this spin model commutes with and thus conserves total Ŝ z , whereŜ γ = 1 2 jσ γ j . When the drive is instead very strong, Ω J , J ⊥ , J DM , total Ŝ x is approximately conserved instead because the drive imposes an energy penalty to flipping spins along the ±x Bloch sphere direction. Note that while the superexchange coefficients can be modified by the drive, as a first rough estimate the highdrive condition Ω J , J ⊥ , J DM may be interpreted as a drive faster than the bare superexchange rate, Ω J 2 /U .

B. Spin model regimes
While the general spin model of Eq. (4) is complex, there are parameter regimes where its form simplifies, allowing the dynamical emulation of other more conventional spin models. For a strong drive Ω J DM , the DM interaction is averaged out and can be neglected. Furthermore, assuming the drive also satisfies Ω J , J ⊥ , theσ y (6) Note that this model is not the same as simply taking the XXZlike piece from the first line of Eq. (4), as here the singleparticle drive term commutes with the XXZ term, making it easy to account for in unitary evolution. One could also write an XXZ model for the no-drive limit Ω = 0 (see Appendix B), for which we would just keep the XXZ portion of Eq. (4); this can be valid in the no-drive limit if the J DM term vanishes parametrically due to its sin(φ) factor for φ ≈ 0, π.
As one limiting regime of the XXZ model, in the strongdrive regime Ω J , J ⊥ ,J DM if we also have φ ≈ π, the coefficients J and J ⊥ are approximately equal and opposite (J ≈ −J ⊥ ), causing them to cancel each other out and leave an Ising model, As another special regime of the XXZ model, there is a line in parameter space where the J coefficient vanishes, requiring Ω = ±U cos(φ) sec(φ/2), which causes the (σ x jσ x j+1 + σ y jσ y j+1 ) terms of Eq. (4) to vanish. Assuming that we also still have a strong drive Ω J , J ⊥ ,J DM , the DM interaction remains averaged out as well, leaving only theσ z jσ z j+1 term along with the drive (which is now transverse to the interaction). If we again make the rotating-wave approximation σ z jσ z j+1 ≈ 1 2 (σ y jσ y j+1 +σ z jσ z j+1 ), we arrive at an XY model, We emphasize that this XY model with a strong (commuting) field, and the underlying Ising model with a strong transverse field are only equivalent under unitary time-evolution [30].
For the strong drive regime Ω J , J ⊥ ,J DM with small flux, the coefficients J , J ⊥ are almost equal to J 2 /U . We can thus collect the XXZ model of Eq. (6) into an isotropic Heisenberg term and a perturbativeσ x jσ x j+1 component. In this regime the system can be approximated with a collectivespin one-axis twisting (OAT) model, whose dynamical properties can be explored using restrictions to the fully-symmetric Dicke manifold [31][32][33][34]. The model is written as, The flux must be small compared to 1/L rather than to 1, because the collective regime validity depends on the Dicke manifold gap, which shrinks with system size for nearest-neighbour interactions. Note that while the above model requires a strong drive, like the XXZ model we can also write a one-axis twisting model in the Ω = 0 regime as well (see Appendix B). Similar low-drive collective physics were explored in Ref. [35] for the weakly-interacting regime; in contrast, here we have a strongly-interacting model that nonetheless allows us to map the nearest-neighbour superexchange interactions to a collective-spin model through gap protection. Moreover, as discussed in Ref. [34], the collective behaviour can be more robust when mapping from an XXZ model with anisotropy slightly below unity (on the easyplane side), which this model can realize as discussed in the next section.
Finally, for a small drive Ω J , J ⊥ ,J DM , we must use the full model of Eq. (4). If the drive is turned off completely (Ω = 0), the resulting interaction can be written as, We label the spin interaction in this regime as Heisen-berg+twist (Heisen+T). While a DM interaction is present, its effect for the above parameters can be simplified to Heisenberg model physics in a twisted frame of reference, as will be discussed in Section III. Having a small non-zero Ω J DM will not significantly change this picture aside from adding the corresponding single-particle term J Ω jσ x j , since the coefficients of the model are only weakly dependent on Ω in this regime. For a larger Ω J DM , the DM term will be rotated out.
In addition to all of the above, we can also add an extra field through the use of the laser drive's detuning δ, which would take the form of δ 2 j (n j,e −n j,g ) in the basis of the bare Fermi-Hubbard model, thus adding a spin term of the form ∼ δ 2 jσ z j . While such a term does not commute with the rest of the Hamiltonian, if |δ| |U 2 − Ω 2 |, the superexchange model will remain the same to good approximation. This permits the addition of an extra single-particle term without changing the spin interactions. We do not explicitly do so in this work, but such a detuning nonetheless provides yet another tuning parameter that can be implemented without the need for additional experimental ingredients.

C. Dynamical model comparisons
It is useful to know where the various simplified models discussed in the previous section are applicable. The most rigorous metric of dynamical model agreement is state fidelity, but such a comparison tends to be unnecessarily harsh because the fidelity can drop with increasing system size while experimentally-relevant observables remain in agreement. We instead evaluate the validity of the spin models through comparison of simple collective observables.
We examine two typical time evolutions of the system, starting from product initial states. The first evolution is, The observable C (Z) is the spin contrast, chosen such that there are no fast single-particle oscillations coming from the drive. The second evolution is, where |→ j = (|↑ j + |↓ j )/ √ 2. Here we label our contrast C (X) as the magnetization. For a large drive Ω J , J ⊥ , J DM only the first C (Z) evolution will see non-trivial dynamics, as C (X) is approximately conserved by the drive. On the other hand, for a weak drive Ω J , J ⊥ , J DM we can have non-trivial dynamics for both evolutions depending on the model in question. In principle a full rigorous comparison should consider all possible high-energy state properties rather than the two selected above. Since we mainly seek a qualitative understanding of model regimes, we have chosen a pair of experimentally-simple evolutions for which at least one will have non-trivial dynamics at every point in parameter space.
To determine the validity of a particular spin model, we compute the time-dependence of C (α) for both evolutions α ∈ {Z, X} using a spin model C which is four times the timescale of the fastest superexchange interaction strength at any given point in parameter space. This timescale is used for every model except the OAT, for which we instead use t f = 8/|J − J ⊥ | (16 time units of the twisting term (Ŝ x ) 2 without the L-dependence, which is sufficient to observe entanglement properties such as spinsqueezing). A dynamical error metric for a given spin model is defined as, This metric is a root-mean-square error giving the average difference between a contrast measurement of the Fermi and spin model being considered over the time interval [0, t f ]. For every choice of φ, Ω in the parameter space there will be two error metrics for the two evolutions α ∈ {Z, X}, and we take the worse of the two, Of course, the disagreement between the models will tend to grow at longer times, and to an extent this analysis is qualitative. We choose four times the superexhange rate because this should be sufficient to see non-trivial contrast decay, and be useful for applications such as spin-squeezing or entangled state generation. Fig. 2(a) shows a color plot of the resulting error metrics for the different spin models. The color scheme is chosen such that a given model's color completely vanishes when its error (e)C (X) ϕ=π/2 Ω/J=0 1.

2.
4. where Ω = U , and the spin model description is invalid due to resonance. Purple points indicate parameters for which the contrast does not decay below 1/e at any time. This occurs trivially for φ = 0 (where the spin model is a pure Heisenberg model and no dynamics occur for product states), and for the special regime of φ 1, Ω = 0 (which is discussed in Section III). metric reaches ∆ Spin ≥ 0.25, corresponding to an average error of 0.25 in measurements of C (X) or C (Z) . The Hubbard interaction strength is chosen to be U/J = 50, well into the Mott insulating regime to ensure that it does not contribute additional error. We see the parameter regimes as described in Section II B. Large flux φ ≈ π and strong drive corresponds to the Ising model (red). The line in parameter space satisfying J = 0 corresponds to the XY model (green). There is a narrow red line adjacent to this regime where the Ising model also looks to be valid, but this is a spurious effect caused by the Ising evolution happening to align with the Fermi-Hubbard out to the specific timescale t f we use. For small flux φL 1 and large drive, we have the OAT model (yellow). Underlying the other models is the XXZ model in white, which is valid throughout most of the parameter regime. The XXZ anisotropy parameter ∆ = 2J /(J + J ⊥ ) is plotted in Fig. 2(f); we can attain any ferromagnetic or easy-plane value for Ω/U < 1, and any antiferromagnetic value for Ω/U > 1.
Finally, for low drive Ω = 0 we have the Heisenberg+twist model (blue). Note that the system in this regime can also be described by a low-drive version of the XXZ or OAT models (see Appendix B), which is why the corresponding colors are also present for Ω = 0. The regions in blue are points in parameter space where the chiral DM interaction exclusive to the Heisenberg+twist is necessary to capture the correct timeevolution.
Figs. 2(b-e) show snapshots of the relevant models' timeevolution from different points of the diagram. The dark region near Ω = U is the resonance point where the superexchange denominators vanish and second-order perturbation theory breaks down, causing no spin model to be valid. There is an additional resonance point near Ω = U/2 corresponding to a second-order resonant process not captured by the spin model [29], although the width of this resonance is smaller. We also note that the associated timescales speed up as the resonance point Ω = U is approached. Fig. 2(g) plots a characteristic timescale t d defined as the time needed for the contrast to decay down to 1/e, using the C (Z) contrast for Ω/J J 2 /U (the high drive limit, all data points except the Ω = 0 line) and C (X) for Ω = 0. This allows for an evaluation of experimental tradeoff, where one can move closer to the resonance for faster timescales at the cost of weaker model agreement. There is also a region of φ 1, Ω = 0 where the contrast does not fully decay at any time despite no obvious conservation law protecting it. This persistent magnetization effect will be discussed in Section III.
While these simulations are for relatively small system size L = 8, in general the regimes of validity do not undergo significant change as the size increases because the interactions are nearest-neighbour. The only exception is the OAT model, for which the regime will shrink with increasing L because it requires φL 1 (the yellow region looks relatively large here because we use a small L). We also note that open boundary conditions are used for the above simulations with all models except the OAT; this is to ensure that the chiral properties of the DM interaction are captured, as will be explored in Section III. Comparisons with the OAT use a periodic Fermi-Hubbard model because open boundaries can lead to a minor but non-zero offset to the contrast even in the thermodynamic limit.
For the simpler regimes such as the Ising or one-axis twisting model, the dynamics are well understood and can have analytic solutions [36]. More general XXZ-type dynamics can be complex to treat, as even exact 1D Bethe ansatz techniques are difficult for full dynamical evolution. However, the parameter regime of low drive Ω = 0 where the Heisen-berg+twist model is valid offers a special case. The dynamics there are non-trivial, but exhibit special long-time features that can be understood from even the non-interacting limit, offering analytic tractability while still simulating the dynamical behaviour of a strongly interacting model. In the next section, we will focus on this regime in more detail.

A. Long-time magnetization profiles
Having shown the different regimes of spin models that can be realized with laser-driven SOC optical lattice systems, we focus on the regime of Ω = 0 where the DM interaction plays a role. Conventionally, the ground-state properties of systems including DM interactions can already be quite complex [17]. In our case, we can work directly with the Heisenberg+twist model of Eq. (10). For clarity, we write it again, Here ∆ Ω=0 is an XXZ anisotropy (note that it differs from the anisotropy of the model in Eq. (6) because we are in the low-drive limit here), and D describes the relative strength of the DM term. In our system the coefficients are (again, maintaining Ω = 0), The non-zero flux causes this model to deviate away from a conventional Heisenberg model by a set of position-dependent local rotations of the spin variables [37] because we are in a dressed basis, see Eq.
(2). However, we can map back to a Heisenberg model at the price of changing the initial state. More concretely, the dynamics of H Heisen+T evolving |ψ which correspond to the low-drive evolution from the prior section, are equivalent to: with the initial state becoming a spiral in the plane of the DM interaction, rotating by an angle φ per lattice site under a unitary transformationÛ = e i 2 j jφσ z j . The mapping is exact for open boundary conditions and Ω = 0. Periodic boundaries can lead to incommensurate mismatch if the flux is not a multiple of 2π/L (see discussion in Appendix D). Hereafter, we refer to the dynamics underĤ Heisen+T as the gauged frame (working in a dressed basis), and the equivalent dynamics un-derĤ Heisen as the un-gauged frame (for which the quantization axis of the Bloch sphere is set by the bare atomic states g, e with no site-dependent phases).
The spiral structure from non-zero flux causes the system's dynamics to exhibit non-trivial features due to the additional underlying symmetry of the Heisenberg model. As the simplest example of such features, in Fig. 3(a-d) we plot the timeevolution of the collective magnetization Ŝ x in the gauged frame (starting from |ψ . We find that in all cases, after an initial decay there is a non-zero mean magnetization that persists to infinite time, while the other in-plane component Ŝ y averages to zero. This dependence of this mean magnetization on the flux is shown in Fig. 3(e). If we had started with a product state along theŷ-direction of the Bloch sphere we would have a non-zero mean Ŝ y and zero mean Ŝ x instead. Interestingly enough, we find that while the amplitude of fluctuations about the mean depends strongly on the Hubbard parameter U/J (recall that the spin model is only valid for U/J 1), As an even more intriguing feature, in Fig. 3(f) we find that at small flux φL 1 there is a scaling behaviour as a function of φL, or flux times total system size, which eventually peels off once φ gets large enough. The first minimum of this scaling long-time magnetization occurs at φL = 2π, which corresponds to a full-period twisting of a spiral state in the ungauged frame. We see a non-zero magnetization even at a full period twist because we use open boundary conditions. Periodic boundaries see a similar profile, except with Ŝ x t→∞ falling to zero at φ = 2πn/L for any n ∈ N (see Appendix D for further discussion on boundary effects).
An infinite-time non-zero magnetization independent of U/J for this model is surprising. Naively, one would expect a relaxation to zero magnetization, as the system has rotational symmetry in thex-ŷ plane of the Bloch sphere, and there is no obvious conservation law that discriminates Ŝ x from Ŝ y . One could expect this to be a consequence of 1D integrability, since the model maps to a Heisenberg model in the ungauged frame, but we also find similar effects in equivalent 2D systems (see Sec. III B). One may also consider this to be a fine-tuned regime, but in addition to its persistence for all U/J Hubbard repulsion strengths, we find that this non-zero magnetization is robust to perturbative effects such as harmonic trapping or imperfect filling fraction N/L < 1, which only slightly change the outcome (See Appendix C). Furthermore, the scaling behaviour maintaining a non-zero average at φL = 2π implies that boundary conditions play a non-trivial role even in the thermodynamic limit (this is unique to the interacting spin model, see Appendix D for details).
Aside from simple observables like magnetization, the system can also develop long-time lattice-wide chiral spin imbalances. Fig. 3(g) shows the other site-resolved in-plane magnetization σ y j across the lattice for different snapshots of the time-evolution (with open boundaries). While the mean Ŝ y is zero, we find that the system establishes a tilt in the spinprojection along theŷ direction, indicating that the DM interaction maintains a non-zero spin current at all times, opposed by the relaxation dynamics of the XXZ model. This tilt can be quantified by, which is plotted in Fig. 3(h) showing the same characteristic scaling behaviour as the long-time magnetization. The system generates an extensive spin imbalance depending on the scaled flux φL. We can again approximate it using the noninteracting limit U/J = 0, This imbalance is again surprising, especially because it manifests as an extensive chiral feature resulting from an initial excitation with non-extensive energy in the thermodynamic limit. Furthermore, for φL ≤ 2π the spin is imbalanced in one direction, whereas for φL slightly higher than that the direction is reversed, even though the first 2π only makes a fullperiod revolution of the spins and should not set a preferential spin pumping direction.

B. Symmetry-restricted features
The reason that non-trivial long-time behaviour occurs is because of the underlying exact SU(2) symmetry of the Heisenberg+twist model, together with its associated reduction of available phase space at small values of the flux. While the model we study is in a twisted (gauged) frame, the local basis rotations still preserve the associated conserved quantities, just in a twisted form. With no flux φ = 0 we have a Heisenberg model, which is a critical point of the XXZ model (∆ Ω=0 = 1, D = 0). The addition of the DM term causes this critical point to extend into a line D = ± ∆ 2 Ω=0 − 1 (with a corresponding branch for the ∆ Ω=0 = −1 critical point). For the parameters in Eq. (17), the system remains on the critical line at some position determined by the flux. The additional symmetries of the Heisenberg model, while twisted, still cause the Hilbert space to break into symmetry sectors and restrict the number of states the system can relax into. By comparison, a model sitting off the critical line of ∆ Ω=0 = sec(φ), D = − tan(φ) cannot be mapped to the Heisenberg and will have dynamics that are chaotic even in 1D [38], causing the magnetization to decay to zero for all φ = 0.
Recall that in general, the Heisenberg model conserves both total angular momentum S 2 [with eigenvalues S(S + 1)] and angular momentum projectionŜ x (with eigenvalues M ), splitting the Hilbert space into S 2 shells andŜ x sectors within each shell: With φ = 0, the spiral initial state |ψ Spiral 0 in the ungauged frame will be distributed among these symmetry sectors. For sufficiently small flux, only the highest angular momentum shells are populated. These include the Dicke manifold S = L 2 and the spin-wave manifold S = L 2 − 1, followed by S = L 2 − 2, etc. The higher angular momentum shells have few states per symmetry sector [1 in Dicke, L−1 in spin-wave, then O(L 2 ), O(L 3 ) and so on]. When enough of the initial state population sits in these highest shells, there are insufficient states for the system to relax and an infinite-time magnetization is generated. The equivalence between the Fermi-and spin models can also be understood from this argument; the undriven Fermi-Hubbard model in the un-gauged frame also has SU(2) symmetry regardless of the value of U/J, meaning that the lack of relaxation should persist. Spin-insensitive perturbations such as external harmonic trapping or imperfect filling fraction likewise maintain SU(2) symmetry and preserve the magnetization or imbalance. Non-negligible boundary effects in the thermodynamic limit are also sensible, as the structure of the highest angular momentum shells depends strongly on the boundaries (sinusoidal vs plane-wave), causing the relevant populations as a scaling function of φL to be different. Note that the non-negligible boundary effects are only maintained in the thermodynamic limit for the interacting spin model, however, as discussed in Appendix D. In a sense, the system thermalizes within a restricted set of Hilbert space manifolds that prevent full relaxation in the conventional manner.
To help quantify the above arguments, in Fig 4(a) we plot the overlap of the initial state wavefunction |ψ Spiral 0 in the un-gauged frame with the different symmetry sectors of the Heisenberg modelĤ Heisen , where |φ S,M,n is the n-th eigenstate ofĤ Heisen within the symmetry sector of angular momentum S and projection M . At φ = 0, all population sits in the maximally-polarized Dicke state j |→ j . Increasing φ causes the deeper shells to become populated, increasing the number of states that the system can explore. We give a metric of this property by defining, which is a weighted average of the population in each shell times the number of states per symmetry sector in that shell N S , normalized by the number of states in the largest sector N max = N S=1 (for even L). The metric ζ can be understood as a measure of how big a Hilbert space the system can explore. When ζ 1, the wavefunction has most of its weight in symmetry sectors much smaller in dimension N S than the largest-size ones (with size N max ), and the system will have trouble relaxing. For ζ approaching of order one, most the wavefunction weight is in the biggest possible sectors and we can expect a more conventional decay of magnetization/imbalance. Note that N S is not simply the number of states in the S-shell, but the number of states per sector of fixed M as well. For example, the Dicke manifold has L + 1 states in total, but each symmetry sector of M = − L 2 , . . . , L 2 within it only has 1 state, thus N S=L/2 = 1. In Fig. 4(b) we plot this metric as a function of φL, finding a characteristic scaling crossover in behaviour. The regime where mean infinite time magnetization falls near zero corresponds to the regime where ζ saturates to a value near one.
To connect with more conventional metrics, in the inset of Fig 4(b) we also plot the dynamics of the bipartite entanglement entropy (partitioning the lattice into left/right halves) for different values of φL. As ζ ≈ 1 is approached, the entanglement entropy saturates at its maximum permitted value based on the system size, while for small flux φL 1 it never reaches that value.
The unusual dynamics described above appear reminiscent to other kinds of unusual dynamics associated with integrability, and indeed the 1D nearest neighbor Heisenberg model is integrable. However here we find that integrability is not a necessary (and is in general not a sufficient) ingredient for the observed long-time dynamics. Beyond the mechanism we propose above, which is unrelated to integrability, additional evidence for the unimportance of integrability here is our surprising observation of analogous long-time behaviour in an equivalent 2D Heisenberg model, which is not thought to be integrable. Fig. 5(a) plots the infinite-time magnetization in 2D, using a similar spiral initial state with a 2D structure [i.e. a unitary transformation of the form e iφ(i+j) 2σ z i,j for lattice coordinates (i,j)] . While the qualitative profile is changed, we still find non-zero persistent averages, which actually remain higher out to longer values of φL (with L = L x × L y for lattice length L x and width L y ). This occurs because 2D systems retain more population in the highest shells for the same flux (since the energy gaps between shells scale with coordination number), and those shells have the same symmetry sector sizes N S independent of dimension. Fig. 5(b) confirms this prediction by plotting the same metric ζ, showing that it saturates at ζ ≈ 1 at the same flux that we see the infinitetime magnetization drop to zero. There have also been studies of similar physics in 3D using approximate numerical methods [39], although there persistent infinite-time magnetization was found for φ = 0, π.

IV. EXPERIMENTAL IMPLEMENTATIONS
The general system described in Eq. (1) can be realized in several ways. The most straightforward implementation is with a 3D optical lattice. Dynamics can be restricted to 1D as explored in this work by increasing the transverse lattice depths, e.g. a deep lattice along directionsŷ,ẑ and a shallower depth along a tunneling axisx (though still deep enough to maintain U/J 1). Spin-orbit coupled driving may be realized through a direct optical transition between long-lived internal states (g, e), such as clock states used in conventional atomic clock protocols, which will ensure the coherence times needed for observing spin dynamics. Flux is generated whenever the drive laser wavevector k L has some projection along the tunneling axisx, and will take the value φ = cos(θ)(2πa)/λ L , with a the lattice spacing, λ L the driving laser wavelength and θ the angle between k L and the tunneling axisx.
The most promising platform candidates are alkaline earth or earth-like atom experiments, as they provide long-lived optically separated internal clock states and magnetic field insensitivity [24][25][26]. Alternatively, one may emulate these types of spin physics with nuclear-spin states [40] such as hyperfine states within a given manifold where the spin-orbit coupling is implemented via Raman transitions. The relevant flux in this case will come from the difference in the overall projection of the two Raman beams, e.g. φ = [cos(θ 1 ) − cos(θ 2 )](2πa)/λ L with θ 1 , θ 2 the angles of the beams to the tunneling axisx. Depending on the duration of spin dynamics one wishes to emulate, other atomic platforms such as alkali atoms may also prove useful, especially in regimes near resonance |U | ≈ |Ω| where the spin model still holds, but the timescales are faster. There have also been discussions on the use of Lanthanide atoms [41], which can avoid some of the heating issues typically found in alkali atoms.
Preparing the desired product initial states is straightforward. A product state |ψ (Z) 0 = j |↑ j in the dressed basis is trivial to prepare, as it is equal to a product state of all atoms in the bare atomic basis up to an overall phase, and can thus be initialized with standard optical pumping techniques. Creating a product state |ψ (X) 0 = j |→ j requires a little more effort, because it is an eigenstate of the drive and cannot be generated with the same drive alone. However, such a state can be prepared by skipping the laser phase. Fig. 6 shows such a protocol. One initializes all atoms in the bare atomic ground-state j |g j , implements a π/2 pulse with the same laser used for driving, then skips its phase ahead by π/2. The pulse will create a site-dependent rotation that transforms the state into the desired dressed product state j |→ j , e − iπ 4 j e − iπ 2 e ijφĉ † j,eĉ j,g +h.c.
The same techniques may be used for measuring collective observables. Total Ŝ z is simply measured from the bare atomic excitation fraction, as Ŝ z = 1 2 j ( n j,e − n j,g ). Total Ŝ x can be measured by reversing the above state preparation protocol; one skips the phase by −π/2, makes a π/2 pulse, then measures excitation fraction. Initialize π/2 pulse Skip phase Ωⅇ ⅈjϕ ⅇ -ⅈπ/2 Ωⅇ ⅈjϕ FIG. 6. Schematic for preparing product states in the dressed eigenbasis of the drive. The system is initialized in a bare atomic product state j |g j . The spin-orbit coupled drive is turned on, and a π/2 pulse is made, preparing a spiral state on the equator of the Bloch sphere (the red lines indicate the drive axes of rotation on each lattice site). The phase of the drive is then skipped by π/2, shifting the axes to match the current spin direction on each site (blue lines). This results in a product state j |→ j . Measuring Ŝ x can also be done by reversing this protocol, then measuring j ( nj,e − nj,g ).
In addition to optical lattices, much of the physics can also be done with optical tweezer arrays placed close enough to allow tunneling. These have seen significant recent development due to their tunability and control [42][43][44][45][46][47], and offer an interesting alternative platform for spin dynamics experiments. If using atoms with long-lived internal clock states, a single interrogating laser can be applied in exactly the same way as for the lattice; one simply needs to control either the angle or tweezer spacing to realize the desired φ. Spinorbit coupled Raman schemes are likewise straightforward to adapt. As further possibility, higher-dimensional systems are also simple to generate; one makes the lattice shallower along more than one tunneling direction, or uses a 2D tweezer array. Each direction's corresponding flux will be controlled by the projection of the driving laser(s) onto the relevant axis. This permits studies of the same types of spin models in 2D, as well as interplay between interactions of various types along different dimensions (such as for example an isotropic interaction along one dimension and an anisotropic one along another), which can be relevant to emulating real condensed matter materials.

V. CONCLUSIONS
We have shown that using a single laser drive to induce magnetic flux in a fermionic optical lattice system can realize a wide variety of different spin models across the parameter regime of flux and driving strength. This system is readily implementable in modern optical lattice or tweezer experiments using highly coherent atomic states. It opens a path to greatly improve quantum simulation capabilities using tools already in reach in current experiments. In addition to studying well-understood models such as the Ising or one-axis twisting models, more exotic physics such as lack of relaxation imposed by symmetry constraints can also be explored with this setup.
There are many possible future directions to explore on both experimental and theoretical fronts. One may inquire further as to what happens to similar long-time magnetization behaviour in higher-dimensional systems; we find that 2D also exhibits non-zero steady state averages, but there are qualitative differences in the resulting profiles. Boundary effects can be probed, as we find discrepancies between interacting and free models (see Appendix D), which could become more complex in higher dimensions. A more in-depth study of the persistent magnetization behaviour's breakdown from external perturbations can also be done. Finally, one can study the steady state spin imbalances in the context of spin transport. This can be valid if the DM's prefactor ∼ sin(φ) vanishes for flux close to 0 or π. For the OAT, we project the above model into the Dicke manifold,Ĥ For the dynamical regime diagram in main text Fig. 2, the color scheme is determined as follows. The error metric ∆ Spin is obtained for each of the spin models Spin ∈ {Ising, XY, OAT, XXZ, SE}, and truncated to a chosen maximum error η = 0.25 via ∆ Spin → min(∆ Spin , η), which is the threshold at which the model's color will vanish. Note that the last model "SE" is not the Heisen+T model, but the full spin model, which we use to determine the relevance of the Heisen+T by comparing the SE error to an XXZ error, i.e. checking whether the DM term is relevant. Each point in parameter space is assigned an RGB color coordinate (r, g, b) with color values r, g, b ∈ [0, 1]. The different spin models are likewise assigned colors, with red (1,0,0) for Ising, green (0,1,0) for XY, yellow (1,1,0) for OAT, white (1,1,1) for XXZ and blue (0,0,1) for Heisen+T (which, again, will be determined by the difference in error of the SE and XXZ models). The color coordinates are then determined via, We essentially subtract (1 − ∆ Spin /η) from every color that does not use the corresponding spin model. For example, since the Ising is red, we subtract one minus its (scaled) error metric from green and blue. The only exception is the Heisen+T model, for which we instead subtract the difference in scaled error between the XXZ and SE models (∆ XXZ − ∆ SE )/η, as error in the XXZ but not in the SE means that the DM term is relevant. Note that for the row of points with Ω = 0 (which are the only points where the DM plays a non-trivial role due to the scale of the plot), we use the undriven versions of the XXZ and OAT models,Ĥ XXZ | Ω=0 andĤ OAT | Ω=0 . After computing the colors, we scale them by the overall quality of the full spin model, ν → ν(1 − ∆ SE /η) for ν ∈ {r, g, b}, to capture the resonances where no spin model description works.

C. FILLING FRACTION AND HARMONIC TRAPPING ROBUSTNESS
The persistent magnetization/imbalance properties we see are robust to typical optical lattice imperfections like reduced filling fraction, or harmonic trapping coming from lattice beam curvature. The filling fraction can be modeled by directly replacing some of the lattice site atoms with holes in the product initial-state. Fig. C1(a) shows the magnetization for different filling fractions with a small system of L = 8 using the Fermi-Hubbard model. We find that while there is a reduction of Ŝ x t→∞ , the non-zero infinite-time average persists, in agreement with our symmetry-based arguments.
A harmonic trap can be modeled by adding an additional Hamiltonian term, H Trap = Ω ext j (j − j 0 ) 2 (n j,e +n j,g ) , with Ω ext the trap energy. Fig. C1(b) shows the magnetization with the trap present, showing that the long-time averages likewise persist.

D. BOUNDARY CONDITIONS AND INTERACTION EFFECTS FOR LONG-TIME MAGNETIZATION
As discussed in the main text, the dynamics of the Heisenberg+twist modelĤ Heisen+T evolving from a product state |ψ (X) 0 = j |→ j lead to a non-zero infinite-time magnetization Ŝ x t→∞ . The analogous dynamics of the Fermi-Hubbard model in the same basis lead to similar infinite-time magnetization values. However, there are some deviations depending on the boundary conditions and system size.
We first consider the infinite-time magnetization Ŝ x t→∞ for the Fermi-Hubbard model with open boundaries in the noninteracting limit U/J = 0. The exact result in this regime is given in Eq. (21), which we re-write here for clarity: In Fig. D1(a) we plot this value for a few increasing finite values of L, as well as the asymptotic form in the thermodynamic limit. It is evident that while there is a non-zero mean at the commensurate flux values φL = 2π, 4π, etc. due to the boundaries, this mean vanishes in the thermodynamic limit, which should be expected from a non-interacting model. The situation becomes more complex when considering the U/J → ∞ limit (i.e. the spin model). Fig. D1(b) plots the infinite-time magnetization for the spin model for a few system sizes as a scaling function of φL for both open and periodic boundary conditions. We find that there is a clear difference between the two that persists in the thermodynamic limit due to the scaling behaviour. The periodic boundary case sees the magnetization drop to zero at φL = 2π, 4π, matching the noninteracting model's thermodynamic limit behaviour (as shown by the blue line in comparison). However, the spin model with open boundaries maintains a non-zero magnetization difference even in the thermodynamic limit, which can be viewed as the difference between the green and red points in the figure. We thus see a purely spin-interaction contribution on top of the persistent magnetization behaviour. The SU(2) symmetry creates a lack of conventional relaxation, leading to non-zero values. However, open boundaries cause the interacting system to have non-zero values even at commensurate φL = 2π, 4π, where the thermodynamic limit should expect to see Ŝ x t→∞ = 0, since we have a full-period spiral that should not favor onex-direction over another.
As a final note, we point out that the spin models used in the above plot are in the un-gauged frame, meaning that we evolve under a Heisenberg modelĤ Heisen from a spiral initial state. For open boundary conditions there is no difference between this