1 Experimental Details↩︎

Each experimental cycle begins by generating an \(F=1\) spinor BEC of up to \(1\times 10^5\) sodium atoms in a crossed optical dipole trap with trapping frequencies \(\omega_{x,y,z} \sim2\pi(130,130,180)~\mathrm{Hz}\). The atoms are prepared in their superfluid ground state, the longitudinally polarized (LP) state where \(\rho_0 = 1\) and \(M = 0\). Here \(\rho_{m_F}\) is the fractional population of the \(m_F\) hyperfine state and \(M = \rho_1 - \rho_{-1}\) is the magnetization. The atoms are then loaded into a cubic optical lattice with a lattice spacing of \(\lambda/2\) and associated recoil energy \(E_R=h^2 /(2M_{\mathrm{Na}}\lambda^2)\). Here \(M_{\mathrm{Na}}\) is the atomic mass of sodium, \(\lambda=1064~\mathrm{nm}\) is the wavelength of the lattice beams, and \(h\) is the Planck constant. Similar to our prior work [1], nonequilibrium spin dynamics are initiated by two distinct quantum quench sequences. For one quench sequence, the lattice depth is quenched across the superfluid to Mott insulator phase transition to the final lattice depth \(u_L\) while the quadratic Zeeman shift is held constant at \(q\). In the other sequence, the lattice depth is adiabatically raised through the superfluid Mott insulator phase transition to \(u_L\), before the quadratic Zeeman shift is quenched from a very high value to its final value \(q\). Both sequences result in few-body nonequilibrium dynamics which are monitored using a spin-resolved two-stage microwave time of flight imaging process [1] after allowing the atoms to evolve in the lattice for a time \(t_{\mathrm{hold}}\).

2 Details on theoretical modeling↩︎

2.1 Determination of \(U_2\), \(V_2\) and \(\chi_n\) from simulations↩︎

This section discusses the optimization procedure used to determine \(U_2\) and \(V_2\) from simulation data with the two body and three body Hamiltonians in Eq. [TwoBody] and [ThreeBody] of the main text. The accuracy of the results is assessed via the distance (cost) \(D = \frac{1}{N_p} \sum_{f=f_{\rm low}}^{f_{\rm high}} |A^{\rm exp}(f) - A^{\rm theory}(f)|\) between the theory and experimental spectral curves, where \(N_p\) is the number of points in the summation. To isolate the effect of two-body interactions, \(U_2\) is computed purely by performing this minimization using \(n=2\) data over a frequency range \(f \in [f_{\rm low}, f_{\rm high}]\) in the vicinity of \(f_2^{\rm harm}\) based on the harmonic theory prediction for \(U_2\). In practice, we take \(f_{\rm low / high } = f^{\rm harm}_2 \pm 2 \Delta f\), where \(\Delta f = (t_{\rm final} - t_{\rm init})/2\pi\) is the frequency resolution. Note that we use zero padding to generate ‘smooth’ Fourier spectra \(A(f)\), discretized on a much finer grid than \(\Delta f\). Fig. 1 (a) and (b) show \(D\) vs \(U_2\) for both types of quench sequence. We find that \(D\) is typically a smooth function of \(U_2\), with a clear minimum at an optimal \(U_2\) value. The curve is well approximated by a Gaussian distribution fit (dotted line) which can be used to estimate the uncertainty on \(U_2\). In order to reduce the impact of the chosen time interval on the Fourier transform, the optimization procedure is carried out at least 10 times, each with a slightly different time interval. The results shown in Fig. 1 correspond to the \(D\) value obtained after averaging over the \(D\) values obtained from the different simulation results. The uncertainties stated in the main text are similarly derived.

The optimization of \(V_2\) is more complicated. Due to the finite time interval of the experimental signal, the sampling resolution in the frequency space (\(\approx 15\) Hz) is larger than the expected shift in \(f_3\) due to three-body interactions. However, higher atom numbers can have substantially larger shifts, scaling with \(n^3\). We therefore perform independent simulations for each \(n\) value and consider the linear superposition of these signals, i.e. \(A^{\rm theory}(f) = \sum^{n_{\rm cut}}_{n=2} \chi_n A^{\rm theory}_n(f),\) up to a maximum atom number \(n_{\rm cut} = 8\) selected to bound the dominate Mott lobe occupations observed in prior works [1], [2]. The relative population in the \(n=1\) Mott lobe cannot be inferred using this method as it does not contribute to the spin dynamics. Therefore, for the results shown in Fig. [fig:4] of the main text we have normalized \(\chi_n\) such that \(\sum_{n=2}^{n_{\rm cut}} \chi_n = 1\). We select \(f_{\rm low}\) and \(f_{\rm high}\) to encompass all of the band-gap frequencies shown in Fig. [fig:1] of the main text. However, frequencies above \(f_{\rm high}\) are also included, albeit with a weighted contribution to \(D\) via a scale factor of 0.3, up to approximately 1000 Hz. This helps to avoid the over-estimation of the contribution from high Mott lobes and the excited states of lower Mott lobes (\(n \approx n_{\rm cut}\)), i.e. a large \(\chi_3\) may assist in matching an experimental frequency, but this may also introduce weight in high frequencies that are undesirable. While typically the eigenmodes associated with these higher frequencies have less overlap with the initial state, this is not always true.

We minimize \(D\) to determine the optimal values of \(V_2\) and \(\chi_n\) simultaneously. The optimization of \(\chi_n\) is carried out with the Nelder-Mead algorithm using the Optim package [3] in Julia for a given \(V_2\). This allows us to obtain \(D\) as a function of \(V_2\) with optimal \(\chi_n\) choices. Fig. 1 (c) and (d) show \(D\) vs \(V_2\) for the same datasets used for Fig. 1 (a) and (b) respectively. While the cost landscape is more complex than the case for \(U_2\), it still has a clear minimum. Once again, fitting a Gaussian distribution to the curve near the minima allows us to estimate the uncertainty. In cases where the curve near the minima is relatively flat and not entirely convex (though the curve outside this region is typically smooth and convex), we allow the center of the Gaussian to also be a fitting parameter. The off-set of this center from the true minimum is then added to the standard deviation of the Gaussian to give the final uncertainty estimate.

While the aforementioned procedure determines the optimal \(\chi_n\), the uncertainty of these parameters in principle requires determining \(D\) in an \(n_{\rm cut}-1\) dimensional space; information that is not directly obtained by the Nelder-Mead algorithm. As such, we obtain estimates for the uncertainty by looking at how \(D\) varies as a function of a single \(\chi_{n'}\), with all other \(\chi_n\) (\(n \neq n'\)) fixed to their optimal values \(\chi^{\rm opt}_n\), i.e. a one dimensional cut of the full landscape. Fig. 2 shows \(D\) vs \(\chi_n\) for all \(n\) in the case of a quench at \(q/U_2\approx 0.85\). It can be seen that the cost functions are smooth and well described by a Gaussian fit with standard deviation \(\sigma_n\) (red line), which provides an estimate for the uncertainty. For simplicity, this procedure is carried out with un-normalized values of \(\chi_n\), as this allows the distribution to be symmetric and therefore better described by the fit. In cases where \(\chi_n \approx 0\) the Gaussian fit is not a good approximation, and we refrain from estimating an error. This tends to occur only when a given \(A_n(f)\) is incapable of approximating the data. An example is given in the lower-right plot for \(\chi_6\). Once \(\sigma_{n'}\) is determined we calculate a lower and upper bound \(\chi^{\rm \pm}_{n'}= \chi^{\rm opt}_{n'} \pm \sigma_{n'}\). Finally, we re-calculate the normalization at each of these bounds, \(\mathcal{N}^{\rm \pm}_{n'} = \sum^{n_{\rm cut}}_{n\neq n'} \chi^{\rm opt}_n + \chi_{n'}^{\rm \pm},\) which defines the uncertainty for the normalized parameters as \(\Delta \chi^{\pm}_{n'}= \sigma_{n'}/\mathcal{N}^{\rm \pm}_{n'}.\) This is the uncertainty shown with the error bars in Fig. [fig:4] of the main text. Note that this approximates the naive estimate of \(\sigma_{n'} / \sum_{n=2}^{\rm n_{cut}} \chi_n\) when the uncertainty is small (\(\sigma_n \ll \chi_n\)).

2.2 Example spectral densities \(A_n(f)\)↩︎

Examples of the individual signal contributions from each \(n\) simulation are shown in Fig. 3 (\(q/U_2 \approx 0.85\)) and in Fig. 4 (\(q/U_2 \approx 0.60\)) for the quenches shown in Fig. [fig:4] of the main text. Data for both the real (left column) and imaginary (right column) parts of the spectral density are displayed. For a given \(n\), it can be seen that the three-body and two-body interactions have peaks which are slightly offset, in conformity with Fig. [fig:1](b) of the main text. While we restrict the plots to the most relevant frequency range, in some instances higher frequency oscillations are also visible. The aforementioned fitting procedure does not neglect these contributions. For some frequencies, the experimental signal can only be explained by a single \(n\) value. However, where different \(n\) have comparable frequencies, interference may occur between the different signals in the superposition. Signals that do not contribute towards a viable approximation of the experimental result will be weighted by a very small \(\chi_n\) value.

2.3 Oscillation frequencies from quantum theory↩︎

The quantum state at time \(t\) is expressed in terms of the eigenmodes \(\ket{\nu}\) and eigenvalues \(\epsilon_v\) as \[\begin{align} \ket{\psi(t)} = \sum_{v} c_v e^{-i \epsilon_v t} \ket{v}. \end{align}\] where \(c_v = \braket{\nu}{\psi(0)}\). An arbitrary observable \(\mathcal{O}(t) = \bra{\psi(t)}\hat{O}\ket{\psi(t)}\) can be easily expressed in terms of the eigenmodes as \[\begin{align} \mathcal{O}(t) = \sum_{\nu } |c_{\nu}|^2 \bra{\nu} \hat{O} \ket{\nu} + \sum_{\mu > \nu} 2 \text{Re} \Big( c_{\nu} c^*_{\mu} \bra{\mu} \hat{O} \ket{\nu} e^{-i (\epsilon_{\nu} - \epsilon_{\mu} )t} \Big) , \end{align}\] where all the time-dependence is contained in the second term on the right, which contains the off-diagonal matrix elements. Therefore, the deviation \(\delta \mathcal{O}(t) =\mathcal{O}(t) - \bar{\mathcal{O}}\) from the time-averaged value \(\bar{\mathcal{O}}\) is given by \[\begin{align} \delta \mathcal{O}(t) = \sum_{\mu > \nu} \Big( 2 \text{Re} \big( A_{\mu \nu} \big) \cos \{(\epsilon_{\nu} - \epsilon_{\mu} )t\} + 2 \text{Im}\big(A_{\mu \nu}\big) \sin \{(\epsilon_{\nu} - \epsilon_{\mu} )t\} \Big) \end{align}\] where \(A_{\mu \nu} = c_{\nu} c^*_{\mu} \bra{\mu} \hat{O} \ket{\nu}\). Taking the Fourier transform \(\delta \mathcal{O}(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty} dt e^{i \omega t} \delta \mathcal{O}(t)\) gives \[\begin{align} \delta \mathcal{O}(\omega) = \sum_{\mu > \nu} \Big\{ \text{Re} \big( A_{\mu \nu} \big) \Big( \delta(\omega - \epsilon_{\nu} + \epsilon_{\mu} ) + \delta(\omega + \epsilon_{\nu} - \epsilon_{\mu} ) \Big) + i \text{Im}\big(A_{\mu \nu}\big) \Big( \delta(\omega - \epsilon_{\nu} + \epsilon_{\mu} ) + \delta(\omega + \epsilon_{\nu} - \epsilon_{\mu} ) \Big) \Big\} \end{align}\] Considering only the solution for positive frequencies, the real and imaginary parts of the Fourier transform are \[\begin{align} {\rm Re}[ \delta \mathcal{O}(\omega)] & = \sum_{\mu > \nu} \text{Re} \big( A_{\mu \nu} \big) \delta(\omega - |\epsilon_{\nu} + \epsilon_{\mu}|) , \\ {\rm Im}[ \delta \mathcal{O}(\omega)] & = \sum_{\mu > \nu} \text{Im} \big( A_{\mu \nu} \big) \delta(\omega - |\epsilon_{\nu} + \epsilon_{\mu}|) \end{align}\] This illustrates that the spectra in Fig. [fig:2] of the main text, which show the real and imaginary parts of the Fourier transform of \(\rho_0\), contain information about the off-diagonal matrix elements and initial state overlap factors.

References↩︎