Dynamics of an itinerant spin-3 atomic dipolar gas in an optical lattice

TL;DR

This study investigates the spin dynamics of a large ensemble of spin-3 bosonic chromium atoms in an optical lattice, revealing distinct regimes of behavior influenced by lattice depth, interactions, and magnetic field gradients, with implications for quantum simulation.

Contribution

It provides a comprehensive analysis of spin dynamics across the superfluid to Mott insulator transition, highlighting the roles of long-range dipolar interactions and the necessity of beyond mean-field theories at high lattice depths.

Findings

Spin dynamics are governed by particle motion, interactions, and magnetic gradients at low lattice depths.

Long-range dipolar interactions dominate spin behavior at high lattice depths.

Mean-field theory suffices at low depths, but beyond mean-field approaches are needed at high depths.

Abstract

Arrays of ultra-cold dipolar gases loaded in optical lattices are emerging as powerful quantum simulators of the many-body physics associated with the rich interplay between long-range dipolar interactions, contact interactions, motion, and quantum statistics. In this work we report on our investigation of the quantum many-body dynamics of a large ensemble of bosonic magnetic chromium atoms with spin S = 3 in a three-dimensional lattice as a function of lattice depth. Using extensive theory and experimental comparisons we study the dynamics of the population of the different Zeeman levels and the total magnetization of the gas across the superfluid to the Mott insulator transition. We are able to identify two distinct regimes: At low lattice depths, where atoms are in the superfluid regime, we observe that the spin dynamics is strongly determined by the competition between particle motion, onsite interactions and external magnetic field gradients. Contact spin dependent interactions help to stabilize the collective spin length, which sets the total magnetization of the gas. On the contrary, at high lattice depths, transport is largely frozen out. In this regime, while the spin populations are mainly driven by long range dipolar interactions, magnetic field gradients also play a major role in the total spin demagnetization. We find that dynamics at low lattice depth is qualitatively reproduced by mean-field calculations based on the Gutzwiller ansatz; on the contrary, only a beyond mean-field theory can account for the dynamics at large lattice depths. While the cross-over between these two regimes does not correspond to sharp features in the observed dynamical evolution of the spin components, our simulations indicate that it would be better revealed by measurements of the collective spin length.

Tables

Table 1. Table 1: Parameters, for different lattice depths V 0 subscript 𝑉 0 V_{0} measured in units of the recoil energy E r subscript 𝐸 𝑟 E_{r} for λ = 532 𝜆 532 \lambda=532\, nm. ω z subscript 𝜔 𝑧 \omega_{z} is the estimated trapping frequency in the direction of the magnetic field gradient. After the loading, sites are populated at the most with three atoms, which is the upper cutoff of our simulation. Note U 1 = 7.40 × 10 − 2 ​ U 0 subscript 𝑈 1 7.40 superscript 10 2 subscript 𝑈 0 U_{1}=7.40\times 10^{-2}U_{0} , U 2 = 0.795 ​ U 0 subscript 𝑈 2 0.795 subscript 𝑈 0 U_{2}=0.795U_{0} and U 3 = − 4.71 × 10 − 3 ​ U 0 subscript 𝑈 3 4.71 superscript 10 3 subscript 𝑈 0 U_{3}=-4.71\times 10^{-3}U_{0} . The gradient, which does not vary with lattice depth, is given by γ ​ b ​ ( 𝐫 i ) = 29 × 10 6 ​  Hz m − 1 × 𝐫 i ⋅ 𝐮 ^ Z 𝛾 𝑏 subscript 𝐫 𝑖 ⋅ 29 superscript 10 6 superscript  Hz m 1 subscript 𝐫 𝑖 subscript ^ 𝐮 𝑍 \gamma b(\mathbf{r}_{i})=29\times 10^{6}\text{ Hz m}^{-1}\times\mathbf{r}_{i}\cdot\hat{\mathbf{u}}_{Z} where the lattice spacing in the 𝐮 ^ Z subscript ^ 𝐮 𝑍 \hat{\mathbf{u}}_{Z} direction is 269 nm. The values of the quadratic Zeeman shifts used in the simualtions are given in Fig.  3 (o)

	$V_0=3E_r$	$V_0=9E_r$	$V_0=15E_r$
$U_0/h$ (Hz)	$1250$	$2860$	$4190$
${\omega }_Z/2\pi$ (Hz)	$279$	$337$	$387$
$J_X/h$ (Hz)	$917$	$106$	$18.6$
$J_Y/h$ (Hz)	$11.8$	$5.67\times {10}^{-2}$	$1.11\times {10}^{-3}$
$J_Z/h$ (Hz)	$1380$	$297$	$82.3$
$U_{dd}\left(\mathrm{𝟎}\right)/h$ (Hz)	$-6.96$	$-15.9$	$-23.2$
$U_{dd}\left(\pm {𝐮}_X\right)/h$ (Hz)	$2.86$	$3.02$	$3.02$
$U_{dd}\left(\pm {𝐮}_Y\right)/h$ (Hz)	$-0.176$	$-0.175$	$-0.175$
$U_{dd}\left(\pm {𝐮}_Z\right)/h$ (Hz)	$-2.52$	$-2.60$	$-2.62$

Equations

$$\hat{H}=$$

$$+\sum _i{\mathcal{V}}_{\mathrm{tr}}({\mathbf{r}}_i){\hat{N}}_i$$

$$+q\sum _j\sum _m{m}^2{\hat{N}}_{mj}$$

$$-\gamma \sum _{i}b({\mathbf{r}}_i)\sum _{mn}{S}_{mn}^z{\hat{a}}_{m,i}^{†}{\hat{a}}_{n,i}$$

$$+\frac{1}{2}\sum _{m,m^{\prime },n,n^{\prime }}\sum _i{C}^{m{m}^{\prime }n{n}^{\prime }}{\hat{a}}_{m,i}^{†}{\hat{a}}_{m^{\prime },i}^{†}{\hat{a}}_{n^{\prime },i}{\hat{a}}_{n,i}$$

$$+\frac{1}{2}\sum _{m,m^{\prime },n,n^{\prime }}\sum _{i,j}{D}_{ij}^{m{m}^{\prime }n{n}^{\prime }}{\hat{a}}_{m,i}^{†}{\hat{a}}_{m^{\prime },j}^{†}{\hat{a}}_{n^{\prime },j}{\hat{a}}_{n,i}$$

$$C$$

$$D_{ij}^{m{m}^{\prime }n{n}^{\prime }}$$

$$\times U_{dd}({\mathbf{r}}_i-{\mathbf{r}}_j),$$

$$U_{dd}(\mathbf{R})$$

$$|\psi (t){⟩}_G$$

$$|\psi (t){⟩}_i$$

$$i\hslash \frac{d}{dt}|\psi (t){⟩}_i={\hat{H}}_G^i|\psi (t){⟩}_i,$$

$${\hat{H}}_G^i$$

$$-\sum _{\alpha }{J}_{\alpha }\sum _{j_{\alpha }^{(i)}}\sum _m\left(⟨{\hat{a}}_{m,j}⟩\left[{\hat{a}}_{m,i}^{†}-\frac{1}{2}⟨{\hat{a}}_{m,i}^{†}⟩\right]+h.c.\right)$$

$$+\sum _{m,m^{\prime },n,n^{\prime }}\sum _{j\ne i}{D}_{ij}^{m{m}^{\prime }n{n}^{\prime }}⟨{\hat{a}}_{m^{\prime },j}^{†}{\hat{a}}_{n^{\prime },j}⟩$$

$$\times \left({\hat{a}}_{m,i}^{†}{\hat{a}}_{n,i}-\frac{1}{2}⟨{\hat{a}}_{m,i}^{†}{\hat{a}}_{n,i}⟩\right),$$

$$\hat{H}$$

$$+\frac{1}{2}\sum _{i,j\ne i}{U}_{dd}({\mathbf{r}}_i-{\mathbf{r}}_j)[\frac{1}{2}({\hat{S}}_i^x{\hat{S}}_j^x+{\hat{S}}_i^y{\hat{S}}_j^y)-{\hat{S}}_i^z{\hat{S}}_j^z],$$

$$N_m(t)$$

$$+(1-\eta (t))N_m^{\mathrm{S}=3}(t),$$

$$|\psi (0)⟩$$

$$N_m^{(2)}(t)=⟨$$

$$\displaystyle-it\biggl{\langle}\left[\hat{N}_{m},\hat{H}_{2}\right]\biggr{\rangle}-\frac{t^{2}}{2}\biggl{\langle}\left[\left[\hat{N}_{m},\hat{H}_{2}\right],\hat{H}_{2}\right]\biggr{\rangle},$$

$$N_0^{(2)}(t)$$

$$⟨\bar{S}⟩=\sqrt{⟨{\hat{S}}^x{⟩}^2+⟨{\hat{S}}^y{⟩}^2+⟨{\hat{S}}^z{⟩}^2},$$

$$\displaystyle\langle\bar{S}\rangle\sim\sum_{i}N^{i}\bigg{(}3-\frac{3t^{2}}{2}\bigg{[}$$

$$\displaystyle+\frac{27}{2}\sum_{j^{(i)}}\left(U_{dd}(\mathbf{r}_{i}-\mathbf{r}_{j})\right)^{2}\bigg{]}\bigg{)},$$

$$i\hslash \frac{\partial f_{\mathbf{N}}^j}{\partial t}=F_{\mathbf{N}}^j(\{f_{\mathbf{N}}^i\}),$$

$${\tilde{c}}_0=c_0-\frac{c_3}{7},{\tilde{c}}_1=c_1-\frac{5c_3}{84}{\tilde{c}}_2=c_2-\frac{5c_3}{3},{\tilde{c}}_3=\frac{c_3}{126}$$

$$C^{m{m}^{\prime }n{n}^{\prime }}=$$

$$+$$

$$|\psi (0)⟩=$$

$$+$$

$$N_m^{(2)}(t)$$

$$N_{0,\{1,1\}}^{(2)}(t)$$

$$N_{0,\{2,0\}}^{(2)}(t)=\frac{5}{8}-$$

$$\times \left(U_1+\frac{1}{2}{U}_{dd}^0\right)\left(q+\frac{3}{4}{U}_{dd}^0\right),$$

Full text

I Introduction

Ultra-cold gases provide an excellent platform to study strongly correlated out-of-equilibrium quantum matter. So far, a broad range of atomic, molecular, and optical systems Gross and Bloch (2017); Bloch (2008); Bohn et al. (2017) including trapped-ions Blatt and Roos (2012); Gärttner et al. (2017); Bohnet et al. (2016); Neyenhuis et al. (2017) polar molecules Yan et al. (2013); Hazzard et al. (2014), Rydberg atoms Labuhn et al. (2016); Zeiher et al. (2017); Bernien et al. (2017); Barredo et al. (2018); Guardado-Sanchez et al. (2018), magnetic atoms Lepoutre et al. (2019); Baier et al. (2016); de Paz et al. (2016a, 2013, b); Tang et al. (2018); Böttcher et al. (2019); Chomaz et al. (2018), and cavity QED arrays Norcia et al. (2018); Davis et al. (2019) have been used to realize quantum many-body systems with long-range interactions and to probe equilibrium properties and out-of-equilibrium dynamics both in pinned and itinerant systems.

Magnetic atoms trapped in optical lattices  Lahaye et al. (2009) naturally form a quantum simulator for complex $S>1/2$ models due to the exponential growth in Hilbert space Aharonov et al. (2009); Hallgren et al. (2013). In these systems, the large number of spin degrees of freedom, as well as the ability to create lattices with itinerant particles, where both motional and interaction effects can not be neglected, quickly limit the capability of current state-of-the-art numerical methods to tackle the complex quantum dynamics. The high level of control and tunability in these simulation platforms has already resulted in numerous pioneering experiments in three dimensional optical lattices. These include studies of extended Bose-Hubbard models Baier et al. (2016), and spin lattice models Patscheider et al. (2019) with erbium (Er) atoms using $S=6$ bosonic and $F=19/2$ fermionic isotopes respectively, as well as the spin dynamics of $S=3$ bosonic chromium (52Cr) atoms in both the superfluid and the Mott insulator regimes de Paz et al. (2013, 2016b, 2016a); Lepoutre et al. (2019).

In this work, we present experimental results together with an extensive numerical study of the spin dynamics seen in an array of bosonic 52Cr atoms in a 3D lattice. By tuning the lattice depth, we explore the itinerant regime, where spin dynamics and tunneling occur over similar timescales. This regime is expected to be the most relevant for quantum simulations, since the complexity is then such that exact calculations are intractable by classical computers. We focus on lattice depths that span the Mott-to-superfluid transition. The dynamics is initialized by rapidly rotating a fully polarized equilibrium state, initially pointing along the magnetic field direction, to orient the atomic spins in a direction perpendicular to the magnetic field, corresponding to a coherent superposition of all the magnetic sublevels. The system is then left to freely evolve at a fixed lattice depth. A residual magnetic field gradient and both dipolar and contact interactions generate spin dynamics beyond the expected simple spin precession. Notably, the populations of the seven different Zeeman sublevels are observed to evolve following the spin rotation in a non-trivial way. We also theoretically investigate the dynamics of the total magnetization (i.e. collective spin length) of the gas which tends to decay as the system evolves.

Through our theoretical analysis we find that the system dynamics falls into two general regimes of behavior. (i) At low lattice depths, where the system is in the superfluid regime, spin transport is important as the gradient field drives the magnetic sublevels to spatially separate. The spin dynamics happens in a way that is strongly affected by the onsite interactions. We develop a mean-field Gutzwiller model Kimura et al. (2005); Yamashita and Jack (2007) to study this regime. The Gutzwiller model is able to qualitatively describe the experimentally observed population dynamics as well as the inhibition of the demagnetization process due to the spin dependent onsite interactions. This is similar to the protection against demagnetization and the persistence of ferromagnetic textures induced by the spin dependent contact interactions. These behavior was recently observed in a Bose Einstein condensate in a dipole trap Lepoutre et al. (2018a).(ii) In deeper lattices transport is inhibited as the system enters the Mott-Insulator regime. Here the dynamics of spin populations are driven primarily by long range dipolar interactions. We use a generalized discrete Truncated Wigner Approximation (GDTWA) Lepoutre et al. (2019); Schachenmayer et al. (2015a, b); Polkovnikov (2010) to describe the dynamics in this frozen atom regime. Similar observations were reported before at deep lattices Lepoutre et al. (2019), however what is surprisingly observed in this study is that the spin population dynamics remains almost independent of the lattice depth as the system approaches the Mott insulating regime and that we are able to reproduce well the observed dynamics using the GDTWA method over a broad range of moderate depth lattices.

The Gutzwiller method predicts a strong reduction of spin dynamics as the lattice is raised and crosses the transition to the Mott-insulating state. However, this is at odds with the experimental data, where a pronounced signature of the underlying quantum critical point in the spin population dynamics is absent. We attribute this to the quantum fluctuations which are not taken into consideration in the Gutzwiller treatment. Instead of the abrupt change obtained within the Gutzwiller approach, what is observed is a gradual change between the two previously studied regimes: a classical ferrofluid without lattices Lepoutre et al. (2018b), and a correlated spin model at large lattice depths Lepoutre et al. (2019). We note that a clearer signature of the transition appears in the theoretically calculated spin length, which we expect could be measured in future experiments.

The paper is structured as follows. We discuss the experimental system and introduce the corresponding generalized Bose-Hubbard model that describes the dynamics. We then introduce the Gutzwiller and the GDTWA methods and present the results of a second order perturbative treatment of the spin dynamics. Next, we study the dynamics starting from a simplified two-site model, which can be simulated exactly. We use the two-site model to identify key processes and to also benchmark the validity of the Gutzwiller approximation and the perturbative treatments. We then use the different numerical approaches to model the experimentally observed population dynamics at various lattice depths. We also study the dynamics of the collective spin length (i.e. the dynamical evolution of the the total magnetization). We use both the populations and magnetization dynamics to look for signatures of the superfluid to Mott insulator transition and to understand the underlying physics in the two regimes.

II System and Hamiltonian

Chromium atoms are loaded into a three dimensional (3D) optical lattice as illustrated schematically in Fig. 1, and described by the lattice vectors $\{\mathbf{u}_{X},\mathbf{u}_{Y},\mathbf{u}_{Z}\}$. The lattice depth along each of the three lattice vectors is proportional to $V_{0}$. The spatial extent of the system is determined by the harmonic confinement $\mathcal{V}_{\mathrm{tr}}=\frac{1}{2}m\sum_{\alpha=x,y,z}\omega_{\alpha}^{2}x_{\alpha}^{2}$, and the total number of atoms $N\approx 3\times 10^{4}$. The spin degree of freedom is encoded in the Zeeman sublevels of the $S=3$ ground state of the 52 Cr atoms (see inset to Fig. 1). An external magnetic field in the direction $\hat{\mathbf{B}}$ sets the quantization axis, with the gradient in this field lying approximately along the $\mathbf{u}_{Z}$ direction (see inset to Fig. 1).

In the tight-binding regime Chromium atoms in the lowest band of the optical lattice are modeled as occupying Wannier state spatial orbitals $\mathrm{w}(\mathbf{r}-\mathbf{r}_{i})$, centered on each lattice site $\mathbf{r}_{i}$. Here we consider the lattice site index $i=(i_{X},i_{Y},i_{Z})$ a triad of integers, such that $\mathbf{r}_{i}=\sum_{\alpha=X,Y,Z}i_{\alpha}\mathbf{u}_{\alpha}$. The dynamics of our system is described by the generalized Bose-Hubbard model

[TABLE]

where the indices $m,m^{\prime},n,n^{\prime}$ label the Zeeman sublevels of the atoms, and $\langle..\rangle_{\alpha}$ is used to indicate nearest neighbours along direction $\mathbf{u}_{\alpha}$, with $J_{\alpha}$ being the tunneling amplitude in this direction. The operator $\hat{a}_{m,j}$ $\left(\hat{a}^{\dagger}_{m,j}\right)$ destroys (creates) a spin-3 bosonic particle in the $m$-Zeeman sublevel at site $j$, $\hat{N}_{mj}=\hat{a}_{m,j}^{\dagger}\hat{a}_{m,j}$, and the total number of atoms on each site is given by $\hat{N}_{j}=\sum_{m}\hat{N}_{mj}$. The linear magnetic field gradient is described by $b(\mathbf{r})\approx b\,\mathbf{r}\cdot\mathbf{u}_{Z}$, with $S_{mn}^{x,y,z}$ being the $x,y,z$-components of the spin-3 matrices respectively and $\gamma=g\mu_{B}/\hbar$ where $g\simeq 2$ is the Landé-g factor and $\mu_{B}$ is the Bohr Magneton. The atoms also experience an effective quadratic Zeeman energy shift $q$ that arises from tensorial light shifts of the atomic levels. Both $q$ and $b$ vary with the lattice depth (e.g. see Ref.Lepoutre et al. (2019)).

The last two terms in the Hamiltonian (II) describe the interactions. The contact interaction is of the form Kawaguchi and Ueda (2012)

[TABLE]

where $U_{n}=\tilde{c}_{n}\int d\mathbf{r}|\text{w}(\mathbf{r})|^{4}$. The contact interaction can also include higher order spin terms (for full form see Appendix B), but for the states accessible during the dynamics in consideration these have negligible effect compared to the contributions from the $U_{0}$ and $U_{1}$ terms. Finally, the dipolar interactions are described by (see Kawaguchi et al. (2010))

[TABLE]

where

[TABLE]

is the kernel of the time averaged dipole-dipole interaction between atoms in Wannier states separated by distance $\mathbf{R}$, with $c_{dd}=\mu_{0}(g\mu_{B})^{2}/{4\pi}$. Here we have used $\mathcal{F}$ and $\mathcal{F}^{-1}$ to denote the Fourier transform and its inverse, respectively. Because the dipole-dipole interactions decay rapidly with spatial separation, we restrict the summation to pairs of particles separated by up to one lattice site in each direction. This restriction reduces the computational difficulty, and we find that it does not affect the main results of our paper, as long as only a qualitative agreement is sought for.

To understand the system dynamics it is important to quantify the many microscopic parameters of the system which appear in the Hamiltonian. We summarize these values for three cases of the lattice depth $V_{0}$ in Table 1. We note that the 3D lattice considered here is non-separable and we approximate the Wannier functions of the ground-band as Gaussians in order to obtain estimates of the interaction terms $\{U_{0},U_{dd}(\mathbf{R})\}$. The initial distribution of the atoms in the lattice depends on lattice depth because it varies the strength of the onsite interactions ($U_{n}$) but also because the focused lasers used to make the lattice contribute to the harmonic confinement which increases with lattice depth.

In our experiments an initial state is produced by loading a Bose-Einstein condensate of Cr atoms spin polarized in the $m=-3$ sublevel into the optical lattice at the desired depth $V_{0}$. We assume that the atoms are initially in equilibrium. For our parameters the superfluid to Mott insulator transition is predicted to occur according to the Gutzwiller model at approximately $V_{c}\approx 8E_{r}$, where $E_{r}$ is the recoil energy for $\lambda=532\,$nm. A fast $\pi/2$ microwave pulse ($\sim 5\mu$s long) is then applied to rotate the atomic spins to be along the spin-$x$ axis, and thus in a superposition of all $7$ Zeeman sublevels. The corresponding populations $N_{m}(t)$ are then measured as a function of hold time after the pulse, using Stern-Gerlach imaging. The observed evolution in these populations (see data on Fig. 3) motivates the theoretical analysis that we develop in the following section.

III Theoretical methods

III.1 Gutzwiller method

The Gutzwiller method is a meanfield technique suited to describing bosons in an optical lattice. This approach has been applied to spin-1 bosons (e.g. see Kimura et al. (2005); Yamashita and Jack (2007)), but as far as we are aware our work is the first to consider the extension of this method to describe spin-3 bosons with dipole-dipole interactions. The Gutzwiller method treats onsite terms exactly and inter-site couplings (due to tunneling and interactions) at the meanfield level. The system state in this method is written as a product of states at each lattice site, i.e.

[TABLE]

where $|\mathbf{N}\rangle_{i}=|N_{3},N_{2},N_{1},N_{0},N_{-1},N_{-2},N_{-3}\rangle_{i}$ is the spin-3 Fock state basis at site $i$ and $f_{\mathbf{N}}^{i}$ are the respective amplitudes of the onsite expansion. It can then be shown that the equation of motion for the state at site $i$ is given by

[TABLE]

where $H_{G}^{i}$ is the Hamiltonian at lattice site $i$. This contains the exact onsite Hamiltonian $\hat{H}_{i}$ with all the onsite terms of Eq.(II), and the mean field contribution from nearest neighbours, giving Jreissaty et al. (2011)

[TABLE]

where $j^{(i)}_{\alpha}$ are the nearest neighbours to site $i$ in the $\alpha$ direction, and we also choose to only include the sum over nearest neighbour $j$ in the DDI term.

The Gutzwiller ansatz reduces the Hilbert space dimension significantly and couples sites only through meanfield terms. As such, if every site is identical, the 3D problem only requires the solution at one site, i.e.  $f_{\mathbf{N}}^{i}(t)\to f_{\mathbf{N}}(t)$ for all $i$. However, in our problem the magnetic field gradient and the harmonic confinement break the translational symmetry. Thus, the system can no longer be reduced to a single site within the Gutzwiller approximation. To produce a computationally tractable case we use a model with spatially varying coefficients only along the $\mathbf{u}_{Z}$ direction and assume the other two perpendicular directions to be spatially invariant, i.e. we use the ansatz $f_{\mathbf{N}}^{i_{Z}}(t)$ to describe the dynamics which assumes the sites $i=(i_{X},i_{Y},i_{Z})$ behave identically to $(0,0,i_{Z})$. This choice leaves us explicitly modelling an effective one-dimensional problem along the direction of the gradient field, as schematically shown in the inset to Fig. 1, yet allows us to retain the 3D character of the long range DDIs and tunnelling.

Due to the spin-degrees of freedom the local Hilbert space grows rapidly with the number of atoms per site. We restrict $f_{\mathbf{N}}^{i_{Z}}$ to $\{\mathbf{N}:\sum_{m}N_{m}\leq 3\}$, i.e. up to 3 atoms per site, which should be adequate as the experimental initial state typically has some occupation of doublons (i.e. sites with 2 atoms), but the number of triply occupied sites is negligible both initially and during the dynamics. To match the experimental system we consider a system of 30 lattice $i_{Z}$-sites (see Appendix A for additional details). The initial state is computed by finding the Gutzwiller ground state for $N_{\rm tot}=18$ atoms (occupying the 1D line of sites along the gradient direction), all restricted to the $m=-3$ sublevel. We note that this choice for $N_{\rm tot}$ does not reproduce the initial number of doublons observed experimentally. However, in the experiment the total number of atoms is observed to decrease rapidly with time as $N(t)\approx N(0)\exp(-\gamma t)$, where $\gamma\approx 100s^{-1}$ (see Appendix B) due to dipolar relaxation. Our Gutzwiller theoretical model does not include losses, but our choice of $N_{\rm tot}$ gives a doublon fraction closer to the average value observed experimentally over the $\lesssim 10\,$ms time period. This accounts for the rapidly decaying population of doubly occupied sites and our effective trapping conditions(the choice of $N_{\rm{tot}}$ is discussed further in Sec. IV.2). The initial state density distribution for various lattice depths is shown as the black line in Fig. 3(p), with the result for $V_{0}=15\,E_{r}$ demonstrating a uniformly filled Mott insulator state. To initiate the dynamics we apply a $\pi/2$ spin rotation about the $y$ axis to this initial state.

III.2 GDTWA

In the deep lattice regime we describe the system with a spin Hamiltonian of the form:

[TABLE]

where the first two terms account for the linear and quadratic Zeeman fields and the last term describes the long-range dipolar interactions that couple Cr atoms in different sites. A fraction of sites are doubly occupied in the experiment and we model each of these sites as a pseudo-atom of spin-6, i.e. the maximum total spin of two spin-3 Cr atoms, and their dynamics are still governed by Eq. (9), but with $\hat{S}^{x,y,z}$ replaced by operators corresponding to $S=6$. This treatment is motivated by the fact that the initial spin state for each doubly occupied site is in the symmetric $S=6$ manifold, and the dominant on-site contact interactions do not change this total spin (see Eq.(2)). We have also checked numerically that this is a good approximation. We solve the spin dynamics under the Hamiltonian Eq. (9) by applying the generalized discrete Truncated Wigner Approximation (GDTWA), a numerical approach first introduced in Ref. Lepoutre et al. (2019). In the GDTWA approach, quantum dynamics is obtained by proper sampling the initial quantum state in phase space and averaging over the ensuing trajectories Zhu et al. . The GDTWA method has been shown to be capable of capturing quantum correlations developed during spin dynamics Lepoutre et al. (2019); Patscheider et al. (2019). With this method, we are able to perform calculations including the effect of quantum fluctuations in a relatively large system ( lattice size $13\times 6\times 13$ ) and check that a convergence in system size is reached within experimental uncertainties. The population of a spin state $m$ is obtained via combining the contributions from both singly and doubly occupied sites

[TABLE]

where $N_{M}^{S=3,6}(t)$ is the averaged population on a spin state $M$ calculated for the singly and doubly occupied sites, respectively. The term $\langle 3,m,3,m^{\prime}|6,M\rangle$ denotes the Clebsch-Gordan coefficients, and $\eta$ denotes the fraction of atoms in doubly occupied sites. Due to atom loss from doubly occupied sites, $\eta$ varies with time. The time dependence that we use for $\eta$ is extracted from experimental measurements of the total atom number $N(t)$ (see Appendix B).

III.3 Perturbative treatment

Due to the complexity of the Hamiltonian (II) many factors play a role in the system evolution following the spin rotation. In order to have a better grasp of the contributions of the various terms in the Hamiltonian we first consider the dynamics of the system in the short-time limit using perturbation theory.

We consider the dynamics of an initial state with $N=2$ atoms occupying two neighbouring lattice sites, $i=1,2$, described by the initial state

[TABLE]

where $|0\rangle_{i}$ is the vacuum state in site $i$. The $f_{m}$ coefficients specify the different spin amplitudes after the spin rotation, $f_{0}=-\frac{\sqrt{5}}{4}$ and $f_{\pm 1}=\frac{\sqrt{15}}{8}$ , $f_{\pm 2}=-\sqrt{\frac{3}{32}}$, and $f_{\pm 3}=\frac{1}{8}$. To second order in $t$ the evolution of the population of the different Zeeman levels is given by:

[TABLE]

where $\hat{H}_{2}$ is the Hamiltonian (II) specialized to two sites separated by $\mathbf{u}_{Z}$.

We find that the term proportional to $t$ vanishes, and to second-order in $t$, the spin population dynamics for $m=0$ is given by:

[TABLE]

The evolution of the population of the other sublevels is of a similar form and presented in Appendix C. We have also ignored in the perturbative treatment the dynamics of doubly occupied sites. As explained in the Appendix C.3, the perturbative treatment tends to break down very quickly when those terms are included and a better approximation is obtained when they are excluded. We also note that up to this order, our perturbative result immediately generalizes to the many-body system by simply summing over all pairs of particles.

Equation (12) provides valuable insight into the short-time dynamics of the system. First, it shows that the dipole-dipole interactions between neighboring sites and the quadratic Zeeman term are what drive the dynamics at short time for all lattice depths. Tunneling or linear magnetic field gradients, when present, only cause redistribution of the spin populations at a higher order (See Appendix C.2).

IV Analysis of Dynamics

IV.1 Double-well Dynamics

The computational complexity of simulating Hamiltonian (II) makes any exact treatment beyond a few lattice sites challenging. In order to benchmark our methods we compare the predictions of the Gutzwiller dynamics and the perturbative expressions obtained above to the exact dynamics of a double-well system with two atoms. In Fig. 2(a) we plot the proportion of doubly occupied sites (doublons) in the initial states. This shows that the Gutzwiller approximation predicts that the sites become single occupied for $V_{0}\gtrsim 8E_{r}$, indicating the existence of a superfluid to Mott-insulator transition at $V_{0}\approx 8\,E_{r}$. In contrast, the exact double-well result shows that the doublon fraction varies smoothly with lattice depth. Nevertheless, the range of lattice depths where doubly occupied sites are dominant is well captured by the Gutzwiller model. The panels (b)-(f) compare the perturbative, exact and the Gutzwiller dynamics for the $m=0$ population. As a meanfield method the Gutzwiller approach is not expected to provide an accurate description of a two-particle system, however the comparison reveals that Gutzwiller results are qualitatively correct for low lattice depths $V_{0}\lesssim 8E_{r}$. In this regime the populations tend to oscillate. Furthermore, we note that we have not included the role of doubly occupied sites in the above expressions (see Appendix C for the full results). We find that by including the doublon sites the solution matches the initial curvature of the exact solution (see the inset of Fig. 2 (b)). However this leads to a divergence between the two approaches for longer times (see Appendix C.3). The slower dynamics of the second order perturbation theory for singly occupied sites follow the exact solution for a longer time as shown throughout the panels. For deeper lattices $V_{0}>9E_{r}$, where the system is in the insulating regime and doublons are negligible, the $N_{0}$ population tends to decay, and is well described by the perturbative result. In this regime on the contrary the Gutzwiller model shows suppressed population dynamics. This reveals the failure of the method to account for the necessary quantum fluctuations that mainly drive the dynamics in the frozen atom limit.

To quantify the validity of the different methods, in panel (g) we plot the accumulated error during the time evolution vs lattice depth. It explicitly shows that while the Gutzwiller gives a better description at low lattice depths, the perturbative formula describes well the short time dynamics in the Mott regime. Note that the GDTWA results were not shown in this two site comparison as GDTWA is a method based on semi-classical trajectories, which is more valid with larger $N$. Therefore it does not work for the extreme case of only two atoms.

IV.2 Full system dynamics

In this section we use the Gutzwiller ansatz and the GDTWA to study the dynamics of larger systems in the parameter regimes used for the experimental measurements of the spin population dynamics.

As discussed above, in the Gutzwiller model we do not include the particle-loss (due to dipolar relaxation) dynamically. Instead, we account for it by sweeping over a range of particle numbers $N_{\rm tot}$ and find that the experimental dynamics are best captured by $N_{\rm tot}=18$. We note that the Gutzwiller treatment is not adequate in the deep lattice limit, where we use the GDTWA formalism, including the doublons and losses, to account for the observed dynamics.

The second fitting parameter of the theory is the quadratic Zeeman shift $q$. Since it arises from tensorial light shifts, it is expected to vary with lattice depth $V_{0}$, and while its precise value is not known ab initio, experimental evidence Lepoutre et al. (2019) suggests that its magnitude is bounded by $|q/h|\leq 6$ Hz in the range of lattice depths considered. The observed dynamics, most notably at short times, depends strongly on the sign and magnitude of $q$ (e.g. see Eq. (12)). We have determined the values of $q$ that give the best fit to the initial dynamics. The optimized values of $q$ are shown in Fig. 3(o) for the Gutzwiller and GDTWA methods.

For $V_{0}<10E_{r}$ the optimal value of $q$ determined from the Gutzwiller method is within the expected range of values. However for deeper lattices a much larger $q$ value is required to achieve dynamics of a similar magnitude of that observed experimentally. This suggests that the Gutzwiller dynamics is failing in the deep lattice regime, consistent with our observations of the double-well dynamics. In fact, we know that quantum fluctuations, which play a key role to drive the dynamics in the Mott-insulator regime, are not accounted for in the Gutzwiller model. Hence, the artificially large $q$ value obtained when fitting Gutzwiller results to the experimental data is compensating for the absence of these fluctuations in the simulation. For the GDTWA calculations a value of $q/h=2$ Hz is found to provide good agreement for all lattice depths $>8E_{r}$.

In Fig. 3(p) we show the spatial distributions of the different levels at $t=5$ ms and various lattice depths in the Gutzwiller model, which is a feature that has not been resolved in our experiment. At low lattice depths (see $V_{0}=2\,E_{r}$ result) the Gutzwiller method shows that the spin states spatially separate in wavepackets that are appreciably narrower than the initial density distribution (black). The spin-dependent transport is driven by the magnetic field gradient and causes the center of mass of the various wavepackets to undergo spatial oscillations. The packets tend to be well separated at the time when there is a large dip in the $N_{0}$ population in Fig. 3(a) (and corresponding peak in the higher $|m|$-populations). This observation emphasizes the relevant role played by the interplay between the magnetic field gradient, which drives the spin transport, and the spin dependent contact interactions which lead to a redistribution of the spin populations as observed in the Gutzwiller predictions (see Fig. 3(a)-(d)). As the lattice depth increases spin transport is inhibited, and the wavepacket oscillations gradually become less visible. The Gutzwiller results for the density distributions at $V_{0}=8\,E_{r}$ reveals a small center of mass separation of the spin components (here the wavepackets are broad and mostly overlapping). The Gutzwiller solution predicts the system to become fully insulating, with one atom per site and all transport dynamics freezing out, as the lattice depth exceeds a critical value of $V_{c}\approx 8E_{r}$. For $V_{0}>V_{c}$ the Gutzwiller model fails to capture the experimentally observed spin population dynamics. On the contrary in this regime the GDTWA is able to reproduce very well the observed dynamics, even for $V_{0}$ as shallow as $V_{c}$, as shown in panels Fig. 3(h)-(n).

The sharp change of behavior of the Gutzwiller dynamics as the system crosses the Mott insulator transition is further illustrated in Fig 4(a)-(c), which shows the $m=0$ populations with lattice depth at three times. There we can observe that the Gutzwiller calculations show almost no dependence as a function of lattice depth as soon as the Mott transition is reached. This can also be seen in Fig 3, which shows that dynamics almost disappears in the Gutwiller approach above $V_{c}$. Such a change in behavior around $V_{c}=8E_{r}$ is also observable in the experimental result but in a much less pronounced way. This is because, above the Mott transition, spin dynamics mainly occurs due to inter-site dipole-dipole interactions, and involves the growth of quantum correlations, a feature which can be reproduced by the GDTWA approach Lepoutre et al. (2019) but not by the Gutzwiller simulations. As a consequence, there still exists significant dynamics at large lattice depths, which the Gutwiller ansatz fails to capture.

Although we only observe a gradual and smooth evolution of spin dynamics as the lattice spans the Mott to superfluid transition, it is worth pointing out that the observations are qualitatively different above and below the transition. Below the transition an oscillation is observed, and above the transition the population follows a monotonous evolution, with weak lattice depth dependence. To better quantify this claim, we plot in Fig. 4 (d) the time at which the $m=0$ fractional populations decay from their initial value to 0.25. A clear difference is observed below and above the Mott transition. While the Gutzwiller model overestimates the amplitude of the oscillation for $V_{0}<V_{c}$ it correctly captures its shape. In particular, it is able to capture the slowing down of the dynamics which is experimentally observed as the lattice depth is raised, as highlighted in Fig. 4 (d). To physically understand this slowing down, it is worth referring to our recent results Lepoutre et al. (2018a). Indeed, when the dipoles are tilted compared to the magnetic field axis by $\pi/2$, no spin dynamics can arise (either due to contact interactions or to dipole-dipole interactions) in absence of magnetic field gradients. To leave the initially polarized stationary state, a coupling between spin and motional degrees of freedom is needed which is here provided by the gradient of the magnetic field. As shown in Fig. 3(p), an important effect of the lattice potential is to reduce the relative motion of the different spin states, therefore reducing the effect of the magnetic field gradient, and damping spin dynamics.

On the other hand, only the GDTWA successfully captures the magnitude of the population decay for $V_{0}>V_{c}$. In fact, an unexpected outcome of the comparison between the experimental data and the GDTWA is that GDTWA reproduces well the dynamics even for lattice depths barely above the Mott transition. At large lattice depths, the disagreement of the experimental data with Gutzwiller results together with the good agreement with the GDTWA approach confirm that a new regime is reached, where spin dynamics is governed by different physical processes. As mentioned above, the large lattice depth regime is impacted by strong quantum correlations. We also observe in the simulations that, while spin dynamics is mostly driven by contact interactions at shallow lattice depths, it is almost entirely driven by dipole-dipole interactions at large lattice depths ($V_{0}>V_{c}$) .

V Total Magnetization and Gap Protection

The dynamical evolution of the populations of the Zeeman levels, as presented in Fig. 3(a)-(n), varies smoothly when crossing the underlying superfluid to Mott transition. In order to better understand the spin dynamics both in the Superfluid and Mott regimes, we now study the dynamics of the collective spin length, corresponding to the total magnetization, which we find better reveals the underlying abrupt transition. It is given by:

[TABLE]

Here $\hat{S}^{x,y,z}=\sum_{i}\sum_{m,n}S_{mn}^{\alpha}\hat{a}^{\dagger}_{m,i}\hat{a}_{n,i}$ are collective spin observables. For the case of an initial state following a $\pi/2$ rotation, the quantity $\langle{\hat{S}^{z}}\rangle$ which is conserved during the dynamics is equal to zero, and therefore $\langle\bar{S}\rangle=\sqrt{\langle{\hat{S}^{x}}\rangle^{2}+\langle{\hat{S}^{y}}\rangle^{2}}$.

In Fig. 5 we plot the dynamics of the collective spin length obtained by the Gutzwiller predictions ($V_{0}<8E_{r}$) and GDTWA ($V_{0}\geq 8E_{r}$). The dashed green line also shows the short time dynamics obtained from perturbation theory of the double-well system up to terms quadratic in time, $O(t^{2})$, summed over pairs of atoms at neighbouring lattice sites:

[TABLE]

where $N^{i}$ is the population at site $i$ and the sum over $j^{(i)}$ includes all of site $i$’s nearest neighbours. This is dominated by the large magnetic gradient term, and is therefore largely independent of lattice depth.

The collective spin length provides direct information of the spin coherence and it is experimentally accessible in a Ramsey sequence Yan et al. (2013); Hazzard et al. (2014); Lepoutre et al. (2018a) performed by applying a $\pi/2$ pulse after the free evolution before measuring the population. Two leading processes are expected to generate magnetization decay: the field gradient and the interactions. While the magnetic gradient generates single particle dephasing since it causes different lattice sites to precess at slightly different rates, interactions entangle the spins leading to a loss of information when one traces out over a part of the total system as one does when computing local observables such $\langle\hat{s}_{i}^{x,y}\rangle$. The quadratic Zeeman term can also lead to magnetization decay in this case due to the development of intra-spin correlations, i.e. correlations between the individual electrons inside each atom Fernholz et al. (2008).

While generically, interactions and inhomogeneities both can lead to magnetization decay, the spin dependent interactions proportional to $U_{1}$ in Eq. (2) counter-intuitively can favor spin alignment for weakly interacting atoms. This is because these interactions open a gap in the energy spectra that can suppress dephasing processes as experimentally demonstrated in recent work Lepoutre et al. (2018a); Norcia et al. (2018); Smale et al. (2019). Signatures of the gap protection can be observed in Gutzwiller predictions for shallow lattices Fig. 5 (a-f). The protection is present when there are more than one atom per lattice site and when the $U_{1}$ term can lock the spins favoring alignment. This manifests as an oscillatory behavior in the collective spin length instead of rapid decay. The protection however enters as a higher order process and it is not observable in the perturbative analysis that neglects third and higher order terms in time. The perturbative analysis however does provide a relatively good description of the short time dynamics.

In the Mott regime on the other hand interactions and magnetic field gradient cooperate and both lead to a fast decay of the contrast as can be observed in the GDTWA simulations Fig. 5 (h-n). Similar behavior is observed in the exact solution of a double well system Fig. 5 (o), validating the behavior observed in the many-body system.

The relatively sharp transition between oscillatory and overdamped behavior around the critical point ($V_{0}=7E_{r}$ to $V_{0}=8E_{r}$) might be overestimated in the Gutzwiller approximation but might survive in the full quantum system. Experimental measurements of the contrast will be needed to test if this is the case.

VI Conclusions

In this work we presented extensive theoretical and experimental comparisons of the dynamics of itinerant spin-3 Cr atoms in a 3D optical lattice and subject to harmonic trapping along all three directions. The microscopic Hamiltonian governing the dynamics of the system is complex and features single particle motion, as well as contact, and long-range dipolar interactions. Exact modeling of the experimental dynamics exceeds the capabilities of classical computation, so to compare to our experimental results, we have developed a variety of approximate models. First, we studied the exact population dynamics of the Zeeman levels of a two-site system and compared our results to a Gutzwiller treatment, as well as a short-time perturbative treatment. We thus demonstrated that for shallow lattice depths below the superfluid-insulator phase transition at $V_{c}$, where a significant portion of the sites are doubly occupied, the Gutzwiller description provides a good qualitative description of the exact dynamics. However, we found that the Gutzwiller description fails reproducing the exact dynamics for large lattice depths, above the Mott transition. For large lattice depths, on the other hand, we have used an effective spin model description (GDTWA) whose short-time dynamics was shown in Lepoutre et al. (2019) to match the exact solution for large enough lattice depths, where tunneling and double occupancy are suppressed.

Armed with the above intuition, we applied the Gutzwiller approximation to lattice depths $V_{0}<V_{c}$ and the GDTWA method, which incorporates the effect of quantum fluctuations, to $V_{0}>V_{c}$ to benchmark the experimental observations. We observed qualitative agreement between the experimental results and the theoretical studies, which confirms that these approximate methods can be trusted in their respective domain of validity. In turn, the comparison with these two models provides unique physical insights regarding the physics at play for different lattice depths. While that for $V_{0}>V_{c}$ dipolar interactions and quantum effects play an essential role in the observed dynamics, we found that for $V_{0}<V_{c}$ transport and contact interactions both play an essential role (while intersite quantum correlations can then be neglected). Our analysis thus shows that the two different regimes of low and high lattice depths are qualitatively different. This can indeed be seen by contrasting the behavior of spin dynamics in these two regimes: while spin dynamics is oscillatory and lattice depth dependent for $V_{0}<V_{c}$, the behavior is mostly monotonous and lattice depth independent for $V_{0}>V_{c}$. The cross-over between these two behaviors is however smooth, and does not reveal the sharp underlying Mott to superfluid transition. Therefore, we propose an experimental measurement of the spin-length which can be readily implemented experimentally, as it should display a more pronounced change in behaviour as the system crosses from the superfluid to Mott insulator regime.

Acknowledgements.

We thank Johannes Schachenmayer for useful discussions and Asier Pineiro and Andrew Wilson for reviewing the manuscript. Funding: A.M.R is supported by the AFOSR grant FA9550-18-1-0319 and its MURI Initiative, by the DARPA and ARO grant W911NF-16-1-0576, the DARPA DRINQs grant, the ARO single investigator award W911NF-19-1-0210, the NSF PHY1820885, NSF JILA-PFC PHY-1734006 grants, and by NIST. The Villetaneuse group acknowledges financial support from Conseil Régional d’Ile-de-France under DIM Nano-K / IFRAF, CNRS, Ministère de l’Enseignement Supérieur et de la Recherche within CPER Contract, Université Sorbonne Paris Cité (USPC), and the Indo-French Centre for the Promotion of Advanced Research - CEFIPRA under the LORIC5404-1 contract. A.S.N was supported by the Australian Research Council Centre of Excellence for Engineered Quantum Systems (project number CE170100009), and funded by the Australian Government. B.Z. is supported by the NSF through a grant to ITAMP.

Appendix A Numerical Methods

The Gutzwiller dynamical equations can be obtained by variationally minimizing $\langle\psi_{G}|i\hbar\frac{\partial}{\partial t}-\hat{H}|\psi_{G}\rangle$. This yields a set of nonlinear differential equations for the evolution of the Gutzwiller coefficients $\{f_{\mathbf{N}}^{i}\}$, i.e.

[TABLE]

where $F_{\mathbf{N}}^{j}(\{f_{\mathbf{N}}^{i}\})=_{j}\langle\mathbf{N}|\hat{H}_{G}^{j}|\psi_{G}(t)\rangle.$ The form of the $F_{\mathbf{N}}^{j}(\{f_{\mathbf{N}}^{i}\})$ is analytically cumbersome (e.g. see Asaoka et al. (2016) for the spin-1 case), but is easily evaluated numerically by taking expectations of the various operators terms appearing in $\hat{H}_{G}^{j}$ in terms of the Gutzwiller coefficients. From the initial condition the system of equations (15) is evolved using an adaptive step Runge-Kutta method, with the tolerance set sufficiently low that the solution converges.

Appendix B Parameters

The coefficients for the contact terms are given by Ref. Kawaguchi and Ueda (2012) as

[TABLE]

where $c_{0}=71g_{c}a_{B}$, $c_{1}=3g_{c}a_{B}$, $c_{2}=-15g_{c}a_{B}$, and $c_{3}=-46g_{c}a_{B}$ with $g_{c}=4\pi\hbar^{2}/M$ and $a_{B}$ is the Bohr radius. This gives $U_{n}=\tilde{c}_{n}\int d\mathbf{r}|\text{w}(\mathbf{r})|^{4}$ . The contact interaction also contains a $U_{2}$ and $U_{3}$ term

[TABLE]

which are included in our model, but have negligible effect for the states accessible during the dynamics in consideration compared to the larger $U_{0}$ and $U_{1}$ contributions.

In this paper, with the exception of the GDTWA treatment, we use the nearest-neighbour as the cut-off for the dipolar interactions. This is justified as in shallower lattices where doubly occupied sites exist, the on-site dipolar interactions are dominant. Generally, we find that for the range of lattice depths considered in this paper, the effect of second nearest neighbour interactions is small compared to that of nearest neighbour interactions. We provide a comparison between these terms in Fig. 6.

Atom loss is significant in the experiment, and in Fig. 7 we show the experimental result for the population over time, along with the $N(t)=N(0)\exp(-\gamma t)$ approximation for two values of $\gamma$. At low lattice depths $\gamma=100\,s^{-1}$ fits the experimental data well, while at larger lattice depths there is an initial fast decay of $\gamma=170\,s^{-1}$, which then slows and returns to the $\gamma=100\,s^{-1}$ value by $10$ ms.

Appendix C Perturbative Expressions for Substate Population Dynamics

C.1 Short-time dynamics at $O(t^{2})$

In the main text, Eq. (12) provides a representative expression for the perturbative (short-time) dynamics of the substate populations for an initial state of singlons. Here we provide the rest of the expressions which take a similar form to those presented in the main text. More generally we can write

[TABLE]

where $p_{1}$ is the amplitude for the state with one atom per site and $p_{2}=\sqrt{1-p_{1}^{2}}$ is the amplitude of the state with two atoms on one site and zero on the other, both of which are chosen to be real. To second order in $t$ we find the substate population dynamics is given by,

[TABLE]

where we have used $N^{(2)}_{m,\{1,1\}}$ and $N^{(2)}_{m,\{2,0\}}$ to indicate the expectation value of $\hat{N}_{m}$ for pairs of sites with singlons ($n_{s}$) and the number of atoms in a doublon ($n_{d}$) configuration, respectively. These are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $U_{dd}^{0}=U_{dd}(\mathbf{0})$ is the onsite DDI and $U_{dd}^{1}$ is the DDI to the neighbouring site. In this derivation we neglected the $U_{2}$ and $U_{3}$ terms as their contribution is negligible.

C.2 Short-time dynamics at $O(t^{4})$

In a double-well system, the $t^{4}$ contribution to the short-time dynamics is given by

[TABLE]

where $A^{m}_{1},A^{m}_{2}$, and $A^{m}_{12}$ are the terms due to state with only singlons, only doublons, and the superposition of the two, respectively. Setting $m=0$, we find

[TABLE]

where we have used $B\equiv\gamma\left[b(\mathbf{r}_{1})-b(\mathbf{r}_{2})\right]$ for two nearest neighbouring sites at $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$. If the double well starts in a superposition of $\{1,1\}$ (singlons) and $\{2,0\}$ (doublon-hole) states then we get an additional contribution to the fourth order dynamics given by,

[TABLE]

Similarly for the $m\pm 1$ we find

[TABLE]

[TABLE]

and

[TABLE]

The remaining terms are given by $A_{n}^{\pm 2}=-\frac{1}{2}A_{n}^{0}$ and $A_{n}^{\pm 3}=-A_{n}^{\pm 1}$.

C.3 Divergence Time Scales

In the inset of Fig. 2(b) we see that the second order results including doublons rapidly diverge from the exact simulation. To predict the timescales over which this divergence occurs, we consider the largest term in the doublon expansion

[TABLE]

The $y$ axis in the inset covers $\Delta N_{0}=6\times 10^{-5}$, and we find $\Delta N_{0}^{(2)}$ changes by this amount in $t=0.03$ ms, which agrees well with the figure. Divergence on a similar order of magnitude to the main figure (i.e. $\Delta N_{0}=0.15$), occurs when the higher order terms become important. Considering the largest $4^{th}$ order term, which depends on tunnelling, we get

[TABLE]

Note that this predicts $\Delta N_{0}^{(4)}=0.15$ at $t=0.29$ ms, which indicates that the $2^{nd}$ order results will breakdown for time scales on the order of $t=0.29$ ms.