Hamiltonian engineering of spin-orbit–coupled fermions in a Wannier-Stark optical lattice clock

Engineering a Hamiltonian system with tunable interactions provides opportunities to optimize performance for quantum sensing and explore emerging phenomena of many-body systems. An optical lattice clock based on partially delocalized Wannier-Stark states in a gravity-tilted shallow lattice supports superior quantum coherence and adjustable interactions via spin-orbit coupling, thus presenting a powerful spin model realization. The relative strength of the on-site and off-site interactions can be tuned to achieve a zero density shift at a “magic” lattice depth. This mechanism, together with a large number of atoms, enables the demonstration of the most stable atomic clock while minimizing a key systematic uncertainty related to atomic density. Interactions can also be maximized by driving off-site Wannier-Stark transitions, realizing a ferromagnetic to paramagnetic dynamical phase transition.

S1 Theoretical Model S1.1 SU(N ) interactions in the Sr optical lattice clock Fermionic 87 Sr atoms have two long-lived electronic orbitals, the 1 S 0 and 3 P 0 clock states, as well as a nuclear spin degree of freedom with I = 9/2. We denote the electronic states as |gi and |ei respectively and the N = 2I + 1 nuclear spin levels as m = I, I + 1, · · · , I. Due to the lack of hyperfine coupling between the nuclear and electronic degrees of freedom, the scattering parameters that describe two-body interactions are independent of the nuclear spin states. This property gives rise to an interaction Hamiltonian invariant under SU(N ) rotations (37,42).
S-wave interactions occur under spatially symmetric collisions. Due to the requirement for fermionic atoms to feature a fully antisymmetric total wave function, to collide under the s-wave channel symmetric nuclear spin states require their electronic orbitals to be antisymmetric. Therefore, the state (|gei |egi)/ p 2 is the only one that can feature s- where M is the mass of a Sr atom, gm (R) and em 0 (R) are fermionic annihilation field operators of nuclear spin m in ground manifold and nuclear spin m 0 in excited manifold respectively.
For spatially antisymmetric p-wave interactions, antisymmetric nuclear spin states require their electronic orbitals to be antisymmetric for an antisymmetric total wavefunction, so the only possible electronic state is (|gei |egi)/ p 2, which interacts via the p-wave scattering volume (b eg ) 3 . Symmetric nuclear spin states require their electronic orbitals to be symmetric, so the three possible electronic states are |ggi, |eei, (|gei + |egi)/ p 2, which interact via the p-wave scattering volumes b 3 gg , b 3 ee and (b + eg ) 3 respectively. The pwave pseudopotential for SU(N ) interaction takes the following form: V gg p (R 12 ) / b 3 gg P + , 3 ]P 12 /2. This leads to the p-wave interaction Hamiltonian in the second quantized form, We focus on the least magnetically sensitive clock transition in 87 Sr, , denoted by |gi and |ẽi respectively. In a large magnetic field, the flip-flop process of nuclear spin states in Eq. (S1) and Eq. (S2) can be ignored, so the interaction Hamiltonian including s-wave and p-wave contributions can be restricted to these two states, so (S3) S1.2 Spin model for the carrier transition As described in the main text, our experimental system is a vertical 1D lattice with magic wavelength ( L = 813 nm), so the external trapping potential V ext (R) is the same for |gi and |ẽi states. To the leading order, we have Here, k L = 2⇡/ L is the wavenumber of the lattice that sets the atomic recoil energy E rec =~2k 2 L /2M , g is the gravitational acceleration, and ! R is the radial trapping frequency. In addition, the clock laser ( c = 698 nm), aligned with the lattice direction, drives the transitions between |gi and |ẽi states with bare Rabi frequency ⌦ and detuning . In the rotating frame of the clock laser, the second quantized Hamiltonian take the following form, where H int is given by Eq. (S3), and Here, k c = 2⇡/ c is the wavenumber of the clock laser, and ↵ (R) is the annihilation field operator for a fermionic atom of internal state ↵.
Our experiment operates in the regime where the collisional rate of relaxation for motional degrees of freedom is slower than internal spin dynamics and trapping frequencies (28,29,(36)(37)(38). This condition ensures only internal levels evolve while atoms remain frozen in their single-particle eigenstates during the dynamics. We first focus on the case of the carrier transition, where the clock laser couples the following two single particle states: n, | " n i ⌘ |ẽ; n X , n Y , W n i and | # n i ⌘ |g; n X , n Y , W n i, where n = {n X , n Y , n}, with n X , n Y denoting the radial harmonic oscillator modes and n the lattice site index of the center of the Wannier-Stark state |W n i. We expand the field operator ↵ (R) in terms of single-particle eigenstates as follows, where c n" and c n# are fermionic annihilation operators for | " n i and | # n i states respectively. Here, the harmonic oscillator wave function is where H n X (X) are Hermite polynomials. The wave function for the Wannier-Stark state where J n (x) are Bessel functions, Under the frozen mode approximation, we treat each atom as a spin-1/2 system spanned by | " n i and | # n i states. Therefore, we define the spin operators, and rewrite the interaction Hamiltonian, where Here, ⌘ |n m| is a dimensionless overlap integral of Wannier-Stark states defined as (S14) U ↵ nm and V ↵ nm are s-wave and p-wave interaction parameters respectively (↵, = {g, e}), where Note that in V ↵ nm we ignore the pwave contributions in theẐ direction, because its leading order terms are overlap matrix elements of gradients of Wannier functions in nearest-neighbor lattice sites based on the expansion in Eq. (S10). These matrix elements are small for parameters used in the experiment.
On the carrier transition, H laser becomes where Here, ' = k c a L = ⇡ L / c is the clock laser phase di↵erence between nearest-neighbor Wannier-Stark states, generating spin-orbit coupling. The Rabi frequency for carrier transition is Fig. S1(A). In the following discussions, it is convenient to remove the phase dependence on lattice sites in the H laser term by a gauge transformatioñ c n" = e in' c n" ,c n# = c n# . Under the gauge transformation the spin operators become: S x n = cos(n')S x n sin(n')S y n ,S y n = sin(n')S x n + cos(n')S y n , S z n = S z n ,Ñ n = N n . (S19) Combining the discussions above, the e↵ective Hamiltonian in the gauged frame be- Now we discuss the dependence of interaction parametersJ ? nm ,˜ nm and D nm on radial harmonic oscillator modes (n X , n Y , m X , m Y ) and the distance along lattice direction (|n m|). As reported in (28,36,37), the overlap integrals are not overly sensitive to the radial modes in consideration, allowing us to simplify the Hamiltonian dynamics in terms of collective spin operators at each lattice site, S x,y,z n = P n X n YS x,y,z n , N n = P n X n YÑ n . Due to the partial delocalization of the Wannier-Stark states along the lattice direction, the dominant terms are on-site and nearest-neighbor interactions. Since the characteristic s-wave interaction strength is much larger than p-wave interaction strength at ultracold temperatures (⇠ 100 nK in our case), we include p-wave interaction only for on-site terms.
All these approximations simplify Eq. (S20) into a large-spin Hamiltonian in a 1D lattice, The interaction parameters for these collective spin operators are calculated by performing a thermal average over radial harmonic oscillator modes, Here, ⌘ 0 and ⌘ 1 are dimensionless overlap integrals for on-site and nearest-neighbor interaction respectively [defined in Eq. (S14)], and the thermal average for s-wave (U ↵ ) and The dependence of interaction parameters 0 , 1 , C 0 on lattice depth V 0 is shown in

S1.3 Density shift of the carrier transition
As described in the main text, we measure the density shift of the carrier transition in Rabi spectroscopy. Note that the clock transition frequency is obtained by the average of two frequencies with the same excitation fraction on the positive and negative detuned side of the ⇡-pulse Rabi spectrum, typically with an excitation fraction near 0.45 (the maximum excitation fraction is near 0.9). The density shift per atom where left and right are the laser detuning from clock transition resonance for the excitation fraction we set on the positive and negative detuned side of the Rabi spectrum, and N loc = 1 2L+1 P L m= L N n+m is the averaged atom number per site in a local region centered around site n. The local region is 2L + 1 ⇠ 15 lattice sites, corresponding to our 6 µm imaging resolution.
To calculate the density shift, we apply a mean-field approximation to Eq. (S21), where Note that we drop the J ? 0 term in Eq. (S21) because this term is a constant for any collective state at each lattice site.
We further simplify Eq. (S25) by assuming all lattice sites share the same atom number N loc in a local region (15 lattice sites). We calculate the spin dynamics in this local region by assuming translationally invariant conditions hS x,y,z In this way, we have a homogeneous field on each site, The mean-field Hamiltonian becomes where B ? ? hSi is the perpendicular component of B, with Dropping the parallel component of B because it does not contribute to the mean-field dynamics, d dt Using the mean-field equations above, we simulate the experimental protocol and obtain theoretical predictions for the density shift. In Rabi spectroscopy, we initialize all the atoms in the ground state (hS z n i = N loc /2) for g ! e case, and all the atoms in the excited state (hS z n i = N loc /2) for e ! g case.
From Eq. (S28), we obtain a simple expression for the density shift by setting it to be the value of at which B ? z = 0: Here, ⌫ s,p ↵! are the s-wave and p-wave contributions to the density shift, & z ↵! is a fitting parameter that accounts for the time evolution of hS z i/N loc during the Rabi dynamics, which depends on the details of the Rabi drive such as the pulse area, excitation fraction, and initial conditions used in the experiment, g ! e or e ! g. Note that ⌫ p ↵! are generated by on-site p-wave interactions, while ⌫ s ↵! are generated by nearest-neighbor s-wave interaction. This dependence allows us to control the density shift by adjusting the spatial extension of the Wannier-Stark states, which is tunable by varying the lattice depth, the key idea to eliminating the density shift presented in the main text.

S1.4 Spin model for o↵-site Wannier-Stark transitions
Apart from the carrier transition, we can also drive transitions to other Wannier-Stark states by using the clock laser to couple two di↵erent internal and motional states of an atom, | " n i ⌘ |ẽ; n X , n Y , W n+l i and | # n i ⌘ |g; n X , n Y , W n i, with l = ±1, ±2, · · · .
Compared to the carrier transition, the subscript n labels di↵erent motional states for | " n i and | # n i states in o↵-site Wannier-Stark transitions. In this case, we expand the field operator ↵ (r) in terms of single-particle eigenstates, Similar to the frozen-mode approximation used for the carrier transition, we treat each atom as a spin-1/2systemspanned by the | " n i and | # n i states defined for the 1 ) compared to the carrier transition parameter ( 0 + 1 ), with the former significantly enhanced compared to the latter. In this case as well the radial temperature at each lattice depth is the reported experimental value S2.1.
Wannier-Stark states coupled by the laser, and rewrite the interaction Hamiltonian in terms of the corresponding spin operators, where Here, ⌘ |n m| , U ↵ nm , V ↵ nm have the same definition as the ones used for the carrier transition The dependence of ⌦ l (l = 1, 2, 3) on lattice depth V 0 is shown in Fig. S1(A).
In the following discussions, we focus on the l = 1 Wannier-Stark transition. Following We can also drop the K l nm term due to the uniform atom population for lattice sites in a local regime. These approximations lead to the following large-spin Hamiltonian in a 1D lattice, where l=1 1 In Fig. S1(D), we compare l=1 1 with its counterpart 0 + 1 in the carrier transition. It is clear that the interaction is significantly enhanced due to site-changing Wannier-Stark transitions.

S1.5 Dynamical phase transition
In the main text we presented theoretical and experimental results on the ferromagnetic to paramagentic dynamical phase transition (DPT) when we address the l = 1 Wannier-Stark transition. Given that Eq. (S37) is a large-spin Hamiltonian, its dynamical phase diagram is well captured by a mean-field approximation. Similar to the carrier transition, we apply the translationally invariant condition hS x,y,z n i = hS x,y,z i in a local regime (15 lattice sites), where hS x,y,z i = 1 2L+1 P L m= L hS x,y,z n+m i. This leads to the following mean-field Hamiltonian, where Writing mean-field equations can be written in terms of normalized expectation value of collective spin operators on a single site s x,y,z = 2hS x,y,z i/N loc , Note that Eq. (S41) takes the same form as the mean-field equations obtained in (35,36), which predicted a DPT between ferromagnetic and paramagnetic phases.
In general terms, a DPT is characterized by the existence of a critical point separating phases with distinct dynamical properties in many-body systems after a sudden quench.
The analog of thermodynamic order parameters is found in long-time average observables, which have a non-analytic dependence on system parameters. We initialize all the atoms in the | #i state, the ground state of our model when 1 ! 1, and then perform a sudden quench of the longitudinal field to its final value 1 . The DPT is signaled by a sharp change in behavior of the long-time average excitation fraction n " = (s z + 1)/2, where In the dynamical ferromagnetic phase, n " ⇡ 0 persists even when the final longitudinal field 1 is varied. In the dynamical paramagnetic phase, n " dynamically adjusts itself following the change of final longitudinal field 1 [See Fig. S2(C)].
In the following, we analyze the critical points for the DPT in our system based on the procedure described in (35,36). Using energy conservation in H MF for an initial state with s z = 1, s x = s y = 0, as well as the identity, In the large-N loc limit, we can eliminate s x and s y , and obtain the following di↵erential equation for s z , 1 2 where V (s z ) = (s z + 1) Our experimental conditions lie in the parameter regime where N loc We interpret Eq. (S44) as the Hamiltonian of a classical particle with position s z moving in the e↵ective potential V (s z ), which is shown in Fig. S2(A). The condition V (s z )=0 determines the physical turnover points of s z . Since V ( 1) = 0, V 0 ( 1)  To count the number of roots in V (s z ), we factor out the known root s z = 1, and then consider the discriminant = 18abcd 4b 3 d + b 2 c 2 4ac 3 27a 2 d 2 of cubic equation /⌦ 1 < 0) dominated by single-particle Rabi flopping where no DPT takes place. Based on our experimental condition (22E rec ,190nK and 2.3s ⇡-pulse), the boundary of these two regimes N loc l 1 =1 /⌦ 1 = 8 p 3/9isequivalent to N loc =63.4, indicated by the black dashed line in Fig. 4(C) of the main text.
These two regimes can also be determined by the asymmetry of the long-time averaged excitation fraction or Rabi lineshape, defined as A LR = (n R n L )/(n R + n L ). Here, where max is the detuning at which the peak value of n " is reached, and we choose f to cover almost the entire frequency range of non-vanishing n " . In Fig. S2(D), we compare the A LR obtained from the long-time averaged excitation fraction and the the one obtained from the Rabi lineshape after a ⇡ pulse. In both cases, the asymmetry A LR becomes more pronounced as the atom number increases in the crossover regime, while A LR saturates near the maximum value in the DPT regime. Note that the many-body decoherence discussed in the next section generally reduces the asymmetry. Nevertheless the saturation behavior observed in the DPT regime is maintained. For convenience, in Fig. 4(C) of the main text we normalize the maximum value of A LR obtained from experimental lineshapes to 1.

S2.1 Sample preparation
Characteristic axial scans at 300 E rec are shown in Fig. S3. Without axial cooling, the red trace, the sample thermally populates many axial modes. Due to the anharmonicity of the trapping potential, transitions between di↵erent oscillator levels are resolvable as distinct peaks in both the positively detuned blue and negatively detuned red sidebands.
We verify that sideband cooling populates the ground oscillator level by observing the elimination of the red sideband.
The radial temperature T r over a range of lattice depths is shown in Fig. S4. The piecewise function Eq. 7 is plotted alongside the red experimental data points.  Figure S3: Sideband Spectroscopy. Axial sideband spectroscopy with a radially cooled sample at 300 E rec . We scan the detuning of the clock laser from the carrier transition using high intensity clock light to probe the axial mode filling. The red trace was measured without axial sideband cooling. The blue trace was measured after implementing axial sideband cooling, illustrating near perfect sample preparation in the lowest motional band. The di↵erent axial transitions are apparent in the axial sidebands enlarged in the lower plot.

S2.2 Coherence time of o↵-site Wannier-Stark transitions
We utilize a similar method as in Ref. (6) to extract a coherence time on the |g ; W n i ! |e ; W n±1 i transition. Using a Ramsey sequence with a randomly sampled phase for the second pulse, we fit the contrast decay between two spatially resolved regions of the sample as a function of dark time, see Fig. S5.
Due to very strong s-wave interactions on this transition, we observe a significant dependence on local density. While this dependence merits further study, here we operate the system in a relatively low density regime compared to our density shift measurements to determine a near optimal coherence. Fitting a single exponential time decay to the contrast, we measure an atomic coherence time of 20 (1) Figure S5: Coherence Time. To determine the coherence time of a site-changing Wannier-Stark transition, we use a Ramsey sequence with a randomly sampled phase for the second rotation. As in Ref. (6), randomly sampling a phase for a given dark time traces out an ellipse. We fit the ellipse to extract a contrast measurement, reported here as purple dots. The contrast decay is fit with a single exponential with decay time ⌧ = 20(1) s.

S2.3 Density shift measurement
As described in the main text, we use extended 'clock locks' to measure the spatially dependent average density shift. Fig. S6 illustrates how we measure the density shift coe cients. Using our camera imaging, we construct a density and frequency map of the sample, shown in Fig. S6A and B respectively. Finally, as illustrated in Fig. S6C, we fit this frequency as a function of atom number with a linear function, with the slope the density shift coe cient.

S3.1 Many-body decoherence in o↵-site Wannier-Stark transitions
Our theoretical model is based on the frozen-mode approximation, which restricts the accessible Hilbert space of each atom into a spin-1/2 degree of freedom spanned by the | " n i and | # n i states. Our spin model is valid in the collisionless regime, breaking down at long times or at high enough densities where mode relaxation is not negligible. Since the interaction strength is significantly enhanced when interrogating site-changing WS transitions, as discussed in previous sections, the mode relaxation rate is expected to be more significant. We take into account the mode-changing collisions phenomenologically by adding a density-dependent dephasing term ( z ) into our mean-field equations for the We use z as a fitting parameter and find it has a linear dependence on N loc as expected from mode changing decoherence. In Fig. S7

S3.2 Scattering parameters
In Ref. (28,37), the relation between the p-wave interaction and p-wave scattering length was missing a factor of 1/2, with the correct coe cient being 3⇡~2b 3 ↵ /2M . Using the past notation, in Ref.  Table 1 that includes the updated s-wave and p-wave scattering lengths that are used in this work.

S3.3 Corrections in the tunneling rate from Gaussian beam geometry
In previous sections, we assume a separable confinement potential and tunneling only along the direction of gravity. However in the experimental system, the Gaussian geometry of the laser beams inevitably couple the axial and radial wave functions. This coupling leads to corrections in the nearest-neighbor tunneling rate which now depends on the thermal distribution of the radial modes. Notice that the Gaussian beam profile of a 1D lattice leads to the following trapping potential, where k L = 2⇡/ L is the lattice wave number, w L is the beam waist, and V 0 > 0 is the lattice depth. Expanding the trapping potential to second order of X, Y , we have where the radial trapping frequency is given by ! R = p 4V 0 /M w 2 L . Based on Eq. (S48), an atomic gas with radial temperature T r feels an e↵ective lattice depth given by V 0 k B T r . Although k B T r ⌧ V 0 , it may still lead to non-negligible corrections to the nearestneighbor tunneling rate, which shows exponential dependence on lattice depth. Note that in the large-spin Hamiltonian discussed in previous sections, the interaction parameters are determined by thermal average over radial modes. To take into account the leading order e↵ects of the thermal distribution, we replace the ground band tunnel coupling by This correction leads to ⇠ 40% increase of nearest-neighbor s-wave interaction strength near the zero-crossing point.

S4 Data File 1
The data presented in Figs. 1-4 is available for download as a separate supplementary file.