\section{Methods}
\textbf{Extended Hubbard model.}
To determine the parameters $J_{xy}$ and $J_{z}$ in Eq.\,\eqref{Heisenberg_eq}
\begin{align}
    J_{xy} = - \frac{4t^2}{U_{\uparrow\downarrow}}, \qquad J_z = \frac{4t^2}{U_{\uparrow\downarrow}}  - \left( \frac{4t^2}{U_{\uparrow\uparrow}} + \frac{4t^2}{U_{\downarrow\downarrow}}\right) \nonumber
\end{align}
we use measurements of the lattice depth $V_0$ and of the three scattering lengths $a_{\uparrow\uparrow}$, $a_{\uparrow\downarrow}$, $a_{\downarrow\downarrow}$. From the calibrated lattice depth, the Hubbard parameters $t^{(0)}$ (in the non-interacting limit) and $U^{(0)}$ (in the single-band approximation) are calculated as \cite{jacksh98}: 
\begin{equation*}
    t^{(0)} = \int dz\, w^*({z}-{a})\left[ -\frac{\hbar^2}{2 m}\frac{d^2}{dz^2} - V_0\sin^2(kz) \right] w({z})
\end{equation*}
\begin{equation*}
    U^{(0)}_{\sigma \sigma^\prime} = g_{\sigma \sigma^\prime} \int d^3r \, |w({\bf r})|^4
\end{equation*}
where $g_{\sigma \sigma^\prime}\,{=}\,{4\pi\hbar^2 a_{\sigma \sigma^\prime}/m}$ and $\sigma,\sigma^\prime\,{=}\uparrow,\downarrow$, using the calculated lowest-band Wannier functions $w(z)$ for each lattice depth \cite{greinerthesis03}. 

For the precision needed to compare experiment to theory, three corrections are applied \cite{mateusz13, munich_spinHelix,will10,Lhmann2012}. (1) Tunneling is modified by the so-called single-band bond-charge \cite{Lhmann2012}, which for single-occupancy is
\begin{equation*}
    t_{\sigma\sigma^\prime} = t^{(0)} - g_{\sigma \sigma^\prime} \int d^3r \ w^*({\bf r}-\delta{\bf r})w^*({\bf r})w({\bf r})w({\bf r}) \nonumber
\end{equation*}
where $\delta{\bf r}\,{=}\,(0,0,a)$ is a displacement by one lattice constant $a$ in the tunneling direction. Through this correction, the tunnelling matrix elements $t_{\uparrow\uparrow}$, $t_{\uparrow\downarrow}$, $t_{\downarrow\downarrow}$ are now slightly spin-dependent. (2) For the on-site interaction, we include admixtures of higher bands, as discussed in \cite{will10,Lhmann2012}. The dominant part is captured by a perturbative correction due to the first and second excited bands
\begin{align}
    U_{\sigma\sigma^\prime} = U^{(0)}_{\sigma\sigma^\prime} -g_{\sigma\sigma^\prime}^2\sum_{{\bf n}_1,{\bf n}_2}\frac{|\braket{{\bf n}_1,{\bf n}_2|{\bf 0},{\bf 0}}|^2}{E_{\text{bg}}} \nonumber
\end{align}
where ${\bf n}_1$ and ${\bf n}_2$ are the three-dimensional band indices of the two atoms and $E_{\text{bg}}$ is the sum of the bandgap energies. Corrections to the tunneling rate $t$ due to population of higher bands is negligible for a Mott insulator with occupation $n=1$, since higher bands are admixed only through virtual doubly-occupied sites. These modifications of $t$ and $U$ both contribute to a modification of superexchange $J_{xy}\,{=}\,t^2/U$ (Extended Data Fig.~\ref{ED-fig:corrections}). The relative correction to $J^{(0)}_{xy}\,{=}\,(t^{(0)})^2/U^{(0)}$, given by $(J_{xy}-J^{(0)}_{xy})/J^{(0)}_{xy}$, is almost independent of lattice depth (for the range of lattice depths studied), is approximately quadratic in scattering length, and in the experiment $J_{xy}$ is typically reduced by $10\,\%$ to $15\,\%$. Note that the correction has the opposite sign as in \cite{munich_spinHelix} since $J_{xy}$ is anti-ferromagnetic in this work. (3) In order to accurately determine $J_z$, one must also consider off-site interactions of the form \cite{Lhmann2012}
\begin{align}
    V_{\sigma\sigma^\prime}=g_{\sigma\sigma^\prime} \int d^3r \ w^*({\bf r}-\delta{\bf r})w^*({\bf r})w({\bf r}-\delta{\bf r})w({\bf r}), \nonumber
\end{align}
where $\delta{\bf r}$ is defined as before. One finds that $J_z$ is modified by the addition of $ 2 ( V_{\uparrow\uparrow} + V_{\downarrow\downarrow} - 2 V_{\uparrow\downarrow} ) $. Depending on the signs and magnitudes of the three interactions, the off-site terms can add to or subtract from the two previously-discussed corrections to $J$ (Extended Data Fig.~\ref{ED-fig:corrections}).

\textbf{Determination of the Heisenberg parameters.}
We calibrate the lattice depth using amplitude modulation spectroscopy \cite{stofferle04}. We record the excitation spectrum of a Bose-Einstein condensate in a 1D lattice when the depth of the lattice is modulated by $3\,\%$ providing the cloud averaged lattice depth with a statistical uncertainty of of $0.2\,\%$. Due to an asymmetric excitation profile, we estimate a systematic error of $1\,\%$. The trapping potential causes inhomogeneity of the lattice depth (and hence the Heisenberg parameters) across the atom cloud. In our case, the trapping potential is determined solely by the Gaussian curvature of the lattice beams. For a Mott insulator of 44 lattice sites in diameter and lattice beams with $1/e^2$-radius of $125\,\mu\text{m}$, the lattice depth along $z$ varies by $1.3\,\%$. Since in our experiment the lattice depth of the beams along $x$ and $y$ are kept at $V_0\,{=}\,35\,E_R$, this results in a variation at $11\,E_R$ of $t$ by $3.6\,\%$, of $U$ by $0.4\,\%$, and of $J\,{\sim}\,t^2/U$ by $8\,\%$.

The lowest and second-lowest hyperfine states of $^7$Li realize the $\ket{\downarrow}$ and $\ket{\uparrow}$ states. We use our previous measurements \cite{interactionSpectroscopy} of $U_{\uparrow\uparrow}$ and $U_{\downarrow\downarrow}$ (lattice depth modulation) as well as measurements of $U_{\uparrow\uparrow}\,{-}\,U_{\uparrow\downarrow}$ and $U_{\uparrow\downarrow}\,{-}\,U_{\downarrow\downarrow}$ (interaction spectroscopy) to determine the three scattering lengths $a_{\uparrow\uparrow}$, $a_{\uparrow\downarrow}$, $a_{\downarrow\downarrow}$ (under inclusion of higher-band corrections) for several magnetic fields $B$. The determined anisotropies $\Delta$ are shown in Fig.~\ref{fig:setup}b (points). Hyperbolic fits to $a_{\uparrow\uparrow}(B)$, $a_{\uparrow\downarrow}(B)$, $a_{\downarrow\downarrow}(B)$ are used to interpolate the values for the anisotropy (solid line). Extended Data Fig.~\ref{ED-fig:corrections} shows $J_{xy}$, $J_z$ and $\Delta$ with (and without) corrections. A recent detailed theoretical analysis \cite{secker19} of the interaction spectroscopy data also provided precise scattering lengths across several Feshbach resonances.  However, this analysis slightly disagreed with our lattice depth modulation data in the range of magnetic fields studied here, and therefore we relied on the experimental data.

The lattice depth calibration and the experimental determination of the scattering length $a_{\uparrow\downarrow}$ lead to an uncertainty for the spin-exchange times $\hbar/J_{xy}$ of about $\pm10\,\%$. The accuracy of the determined anisotropies $\Delta$ is limited by the experimental determination of all scattering lengths $a_{\uparrow\uparrow}$, $a_{\uparrow\downarrow}$, $a_{\downarrow\downarrow}$. The uncertainty of $\Delta$ is estimated to be about $\pm0.1$.

In the experiment (see Fig.~\ref{fig:setup}b), $\Delta\,{\approx}\,0$ was realized by tuning the magnetic field to $B_0\,{=}\,882.63\,\text{G}$. Here the measured Hubbard parameters result in $\Delta\,{=}\,{-0.12}$ including higher-order correction (and $\Delta\,{=}\,{-0.02}$ without corrections). $B_0^\prime = 842.95\,\text{G}$ is a second magnetic field value, which also realizes $\Delta\,{\approx}\,0$ (actually $\Delta\,{=}\,{-0.13}$ including corrections and $\Delta\,{=}\,0.01$ without corrections). We directly compare these two points $B_0$ and $B_0^\prime$ in Extended Data Fig.~\ref{ED-fig:decay_anisotropies}a~and~b, and observe quantitative agreement. Arbitrary anisotropies were realized by using the magnetic field region in between: $B_0^\prime\,{<}\,B\,{<}\,B_0$. In particular, regions with $\Delta\,{>}\,0$ ($\Delta\,{<}\,0$) were realized by magnetic fields $B\,{<}\,850\,\text{G}$ ($B\,{>}\,850\,\text{G}$). The isotropic point $\Delta\,{\approx}\,1$ was realized at $B_1\,{=}\,847.30\,\text{G}$ (actually $ \Delta\,{=}\,1.01$ and corrections here are negligible). 

\textbf{Experimental Setup.} 
In the experiment, we prepare $4.5\,{\times}\,10^4$ $^7$Li atoms in an optical lattice with spacing $a\,{=}\,532\,\text{nm}$ in the Mott insulating regime with one atom per site \cite{tiltedMottInsulator}. We prepare a far from equilibrium initial spin state and probe the spin dynamics in 1D. The lattice beams in the $x$- and $y$-directions are kept at a large constant value (lattice depth of $35\,E_R$) separating atoms into an array of independent 1D chains, with a typical maximum length of $L_\text{max}\,{=}\,44\,a$ (given by the diameter of the Mott insulator), and with an average length ${\langle L \rangle}\,{=}\,33\,a$ (Extended Data Fig.~\ref{ED-fig:chain_lengths}). Initially, the z-lattice depth is also $35\,E_R$, so that the spin couplings are off. The magnetic field is then ramped to the value required for a desired anisotropy $\Delta$. Using radio frequency pulses and a magnetic field gradient, a helical spin pattern is created where the spin component along the chain winds in the xz-plane of the Bloch sphere with a wavevector $Q\,{=}\,2\pi/\lambda$, where $\lambda$ is the wavelength of the spin helix (see Fig.~\ref{fig:setup}c and Methods section \enquote{Preparation of the spin helix}). For $\lambda$ smaller than the system size, the total magnetization of this state is close to zero.

The power in the third lattice beam ($z$-direction) controls the superexchange rate within the chains. Time evolution is  initiated by ramping down the z-lattice depth to a value between $9$ and $13\,E_R$. The ramp time is $0.5\,\text{ms}$, fast compared to superexchange $\hbar/J_{xy}$, but slow compared to tunneling $\hbar/t$. The ensuing coherent dynamics along each chain is governed by a 1D Heisenberg XXZ-model with anisotropy $\Delta$, Eq.~\eqref{Heisenberg_eq}. This is a quantum quench to a far from equilibrium initial state. After a variable evolution time $t$ the dynamics is frozen by rapidly increasing the lattice depth back to $35\,E_R$. The atoms are imaged in the $\ket{\uparrow}$ state via state-selective polarization-rotation imaging (see below).

\textbf{Preparation of the spin helix.} 
A global $\pi/2$-pulse of $75\,\mu\text{s}$ rotates the spin $\ket{\uparrow}_i$ on each site $i$ into the xy-plane of the Bloch sphere ${\ket{\varphi}_i}\,{=}\,{[\ket{\uparrow}_i+\ket{\downarrow}_i]/\sqrt{2}} $. A magnetic field gradient in z-direction causes spin precession at rates which depend linearly on position $z_i$ of the spin thus creating a spin helix ${\ket{\varphi}_i}\,{=}\,{[\ket{\uparrow}_i+e^{-iQz_i} \ket{\downarrow}_i]/\sqrt{2}} $ where the spin winds in the xy-plane. The strength and duration of the gradient determine the wavevector $Q\,{=}\,2\pi/\lambda$ where $\lambda$ is the wavelength of the spin helix. An additional $\pi/2$-pulse rotates the spin helix into a state where the spin winding occurs in the xz-plane ${\ket{\varphi}_i}\,{=}\,{\sin(Q z_i / 2) \ket{\uparrow}_i}+{\cos(Q z_i / 2) \ket{\downarrow}_i}$, so that the full many-body xz-spin helix state is ${|\psi(Q)\rangle}\,{=}\,{\prod_i \ket{\varphi}_i}$. In practice the phase $\theta$ of the winding varies from realization to realization, which amounts to replacing $Qz_i\,{\mapsto}\,{Q z_i}\,{+}\,\theta$. This is caused by small magnetic bias field drifts on the $10^{-5}$ level. The range of $\lambda$ used in the experiment was limited on the short side by optical resolution to $\lambda\,{\geq}\,5.6\,a$ and on the long side by the length of the chains $L_\text{max}\,{=}\,44\,a$. 

\textbf{Imaging.} 
The optical density of the atomic ensemble is too high ($>14$) to allow for in-situ observation of the modulation of $\braket{S^z}$ via absorption imaging. Instead, we employ dispersive imaging, which uses the phase accumulated by the transmitted light in order to form an image of the atomic density distribution. When light at frequency $\omega_L$ is detuned from the atomic resonance $\omega_0$ by many natural linewidths $\Gamma$, it picks up an approximate phase $\theta\,{\approx}\,{-2\delta/\Gamma\,{\times}\,\text{OD}(x,y,\delta)}$, where OD is the optical density at detuning $\delta\,{=}\,\omega_L\,{-}\,\omega_0$, while absorption is suppressed by a sufficiently large detuning $\delta$. In order to form an image, the phase-shifted light must be interfered with a reference beam. In this work, we make use of the fact that the optical transition we use for imaging is driven only by a single polarization component; after passing through the atoms, the shifted and unshifted components are combined on a polarizer. A judicious choice of input and output polarizers yields an interference signal $I$ on the camera which is $I\,{=}\,I_0(1\,{-}\,\sin\theta)/2$ \cite{huletBEC, ketterle_varenna1999}. 

The optical resolution of our imaging system ($\text{NA}\,{\approx}\,0.2$) was determined to have a cut-off at modulation wavelength $\lambda\,{\approx}\,3.0\,\mu\text{m}\,{=}\,5.6\,a$ (330 line pairs per mm). The reduction of the modulation transfer function $\text{MTF}(Q)$ near the cut-off reduces the observed contrast $\mathcal{C}(t)\,{=}\,\text{MTF}(Q)\,{\cdot}\,c(t)$ compared to the real contrast $c(t)$. This does not affect the decay times $\tau$. Assuming that the experimental preparation sequence for the initial spin helix state achieves full contrast $c(0)\,{=}\,1$ for any wavevector $Q$ (based on careful pulse calibration and characterization), we can use $\mathcal{C}(0)$ as a direct measurement of $\text{MTF}(Q)$ and determine the real contrast as $c(t)\,{=}\,\mathcal{C}(t)/\mathcal{C}(0)$.

\textbf{Constant background contrast.} 
For long evolution times, the contrast does not fully decay, but it asymptotically goes to a finite background value. For example, in Fig.~\ref{fig:XX}a this is about $c_0\,{=}\,0.08$. The fitting function $c(t)\,{=}\,{(a_0\,{+}\,b_0\cos\omega t) e^{-t/\tau}}\,{+}\,c_0$  needs to include this offset $ c_0 $, in order to describe the data accurately. The numerical simulations, however, show a decay to zero for long evolution times $t$. Therefore an offset of $c_0\,{=}\,0.08$ had to be manually added in as well, to see agreement between experiment and theory.
In the experiment, the offset is caused by the inhomogeneous density of the atom cloud:  only $90\,\%$ of the atoms are in the Mott insulator state which realizes an  array of 1D spin chains. A small fraction of atoms are in dilute spatial wings, separated by holes which are immobile due the gradient of the trapping potential (which suppresses first-order tunneling as shown in our previous work \cite{tiltedMottInsulator}). These atoms preserve an imprinted spin modulation pattern for long times. We have checked this mechanism by increasing the amount of thermal atoms and clearly observe an increase of the background  $c_0$ (Extended Data Fig.~\ref{ED-fig:holes}). Furthermore, a position sensitive measurement of the contrast confirms, that the main contribution is indeed from atoms in the spatial wings (Extended Data Figs.~\ref{ED-fig:holes}f and \ref{ED-fig:holes}g). In agreement with this model, numerical simulations always show a decay to zero for long spin chains (see Extended Data Figs.~\ref{ED-fig:universal}, \ref{ED-fig:negative_anisotropies}, \ref{ED-fig:initial_phase}, \ref{ED-fig:chain_lengths}).

\textbf{Time-reversal invariance of spin dynamics.}
For an xz-spin helix initial state and time evolution via the XXZ Hamiltonian, the contrast is time-reversal symmetric: $c(t)\,{=}\,c(-t)$, which follows because the state, Hamiltonian and observable (the local magnetization $S^z_i$) can be all expressed real. This also implies invariance against the overall sign of the Hamiltonian $H\,{\mapsto}\,{-H}$. The same argument holds for the system with holes evolving under the bosonic $t$-$J$ model. The initial dynamics of the contrast (in the ideal scenario) in both cases is therefore quadratic, $c(t)\,{=}\,1\,{-}\frac{1}{2}\Gamma^2 t^2\,{+}\,\cdots$, with 
\begin{equation*}
    \Gamma^2 = -\frac{2}{L} \sum_i  \cos(Qz_i + \theta) \langle \psi(Q) | [H,[H,S^z_i] | \psi(Q)\rangle
\end{equation*}
Therefore the timescale of the initial quadratic decay $|\Gamma|^{-1}$ is the superexchange timescale $\hbar/J_{xy}$ (XXZ model) or $\hbar/t$ ($t$-$J$ model). The fact that experimentally we do not observe an initial quadratic behavior indicates either (i) the presence of holes, but we are not resolving the fast timescale $\hbar/t$, or (ii) that the initial state is not time-reversal invariant (i.e.~cannot be expressed real, in the same basis that  the Hamiltonian is written in). The latter could arise from pulse imperfections, or the fact that the ramp-down of the optical lattice takes place over a finite duration of time, leading to deviations from the ideal initial state. Nevertheless, we expect that the overall behavior of the  decay of the visibility, e.g.~scaling behavior of dynamics with wavevector $Q$, is not strongly affected by (i) or (ii).
 
\textbf{Numerical simulations.} 
In the numerical simulations we consider: (i) a spin helix quench under XXZ Hamiltonian (Eq.~\eqref{Heisenberg_eq}) dynamics, and (ii) a spin helix with $2.5\,\%$ to $5\,\%$ hole probability evolving under the bosonic $t$-$J$ model (i.e. assuming no doubly occupied sites), given by
\begin{align*}
    H & = \sum_{\langle ij \rangle } J_{xy}(S^x_i S^x_j + S^y_i S^y_j) + J_z S^z_i S^z_j  + H_d \nonumber \\
\end{align*}
with
\begin{align*}
    H_d & = -\sum_{\sigma,\langle ij \rangle } t a_{\sigma i}^\dagger a_{\sigma j} - \sum_{\sigma, \langle i jk \rangle} \Bigg[\frac{t^2}{U_{\uparrow \downarrow}} a_{\sigma i}^\dagger n_{\bar\sigma j} a_{\sigma k}    \nonumber \\ 
    &  + \frac{t^2}{U_{\uparrow \downarrow}} a_{\bar\sigma i}^\dagger S_j^\sigma a_{\sigma k} + \frac{2t^2}{U_{\sigma\sigma}} a_{\sigma i}^\dagger n_{\sigma j} a_{\sigma k}+ \text{h.c.} \Bigg]
    \label{eqn:tJ}
\end{align*}
where spin ${\sigma}\,{=}\,{\uparrow,\downarrow}$. Here $a_{\sigma i}$, $a^\dagger_{\sigma i}$ are bosonic lowering and raising operators at site $i$,  $S_i^\uparrow\,{\equiv}\,S_i^+$ is defined as $a_{\uparrow i}^\dagger a_{\downarrow i}$, while $S_i^\downarrow\,{\equiv}\,S_i^-$ is defined as $a_{\downarrow i}^\dagger a_{\uparrow i}$. 
We use parameters from experiments and focus on a lattice depth of $11\,E_R$, in which case we have $U_{\uparrow\downarrow}/U_{\uparrow\uparrow}\,{=}\,1.206$, $1.406$, $1.401$, $1.398$, $1.397$, $1.392$, $U_{\uparrow\downarrow}/U_{\downarrow\downarrow}\,{=}\,{-0.188}$, $ 0.264$, $0.459$, $0.575$, $0.659$, $0.862$ and $U_{\uparrow\downarrow}/t\,{=}\,{-17.94}$, $-24.32$, $ -24.17$, $-24.08$, $-24.05$, $-23.94$ for anisotropies $\Delta\,{=}\,0.020$, $0.670$, $0.860$, $0.973$, $1.055$, $1.256$ respectively. In the absence of holes, the action of the term $H_d$ vanishes, and the Hamiltonian reduces to the XXZ Hamiltonian.

In both cases we employ the real-time time-evolving block decimation (TEBD) method with matrix product states (MPS) on an open spin chain. In case (i): local Hilbert-space dimension$\,{=}\,2$, length $L\,{=}\,40$ sites and bond dimensions up to $800$; case (ii):  local Hilbert-space dimension$\,{=}\,3$, $L\,{=}\,40$ for $\Delta\,{=}\,0$ or $L\,{=}\,20$ for all other $\Delta$. Simulations are cut off in simulated time due to the rapid increase in entanglement of the state, requiring ever increasing computational times. For the special case $\Delta\,{=}\,0$, without holes, we alternatively employ free fermionic methods, significantly speeding up the calculations. For technical reasons, instead of simulating the full distribution of holes we simply average over the situations in which there are either exactly one or two holes in the chain. For $L\,{=}\,20$ ($40$), a single hole corresponds to an average of $5\,\%$ ($2.5\,\%$) holes.

To obtain the contrast, we determine the local magnetization $\langle S^z_i(t) \rangle\,{=}\,\langle n_{i,\uparrow}(t) \rangle\,{-}\,\langle n_{i,\downarrow}(t) \rangle$ and determine the Fourier component at wavevector $Q$ via $c(t)\,{:=}\,\frac{2}{L} \sum_i \langle S^z_i(t) \rangle \cos(Q z_i + \theta)$, assuming that the initial spin helix has wavevector $Q$ with a given  phase $\theta$ (see Methods section \enquote{Creation of the spin helix}); we then compare this to the experimentally extracted contrast. However as noted in the Methods section \enquote{Time-reversal invariance of spin dynamics}, the numerically simulated contrast always has an initial quadratic decay, unlike that seen in the experiments. We therefore determine the decay timescales and its power law scaling with $Q$ with one of the following methods: for positive anisotropies $\Delta$, we define the decay time constant as the time it takes for the contrast $c(t)$ to decay from 1 to 0.6, multiplied by ${-\ln(0.6)}$ to convert to a \enquote{$1/e$} time (see Extended Data Fig.~\ref{ED-fig:universal}). For negative anisotropies, for regime I, we take the time to decay from 1 to 0.9, while for regime II we fitted a simple exponential decay profile to obtain the decay timescales  once the curves start \enquote{peeling off} when plotted in time rescaled by $\lambda$ (see Fig.~\ref{fig:negative}c and Extended Data Fig.~\ref{ED-fig:negative_anisotropies}). The resulting power law scalings (Extended Data Fig.~\ref{ED-fig:powerlaw_theory}) and resulting exponents (Fig.~\ref{fig:powerlaw}c) are in reasonable agreement with the experimental results.

\textbf{Finite size effects.}
Numerical simulations can study the effect of different chain lengths and initial phases of the spin helix on the dynamics. This is relevant because in the experiment, the atoms are distributed uniformly over a sphere with a typical diameter of $L_\text{max}\,{=}\,44\,a$ in three-dimensions, leading to an ensemble of 1D chains of varying lengths. The experimentally measured contrast $c(t)$ is an average over all chain lengths with a probability distribution shown in Extended Data Fig.~\ref{ED-fig:chain_lengths}. Furthermore, due to drifts of the applied magnetic field, the initial phase $\theta$ of the spin helix state varies from shot to shot. Here we study numerically both effects.

We concentrate on the XX model ($\Delta\,{=}\,0$) without holes.  Extended Data Fig.~\ref{ED-fig:initial_phase} shows that the strong dependence of the  visibility $c(t)$ on the initial phase is due to reflection of magnetization off the boundaries of the chain. This suggests that averaging over various chain lengths or averaging over initial phases should give similar results which is confirmed in Extended Data Fig.~\ref{ED-fig:chain_lengths}.  It is even sufficient to average over only two phases, $\theta\,{=}\,0$ and $\pi/2$, to achieve insensitivity to initial and boundary conditions. The distribution of chains, and the simulation of a magnetization profile averaged over the distribution of chains, are illustrated in Extended Data Fig.~\ref{ED-fig:chain_lengths}b-d.

The conclusion is that the experiment is naturally performing an average over different phases and different chain lengths, washing out the sensitive dependence of the spin dynamics on initial conditions.  For comparison with simulations, it is sufficient to use a system with a  fixed chain length $L\,{=}\,40\,a$, and average over only the two phases $\theta\,{=}\,0$ and $\pi/2$.

\setcounter{figure}{0}
\renewcommand{\figurename}{Extended Data Fig.} 

\onecolumngrid

\newpage

\begin{figure}[h]
    \includegraphics[width=0.8\linewidth,keepaspectratio]{EDFig1.pdf}
    \caption{
        \textbf{Corrections to Hubbard parameters}. \textbf{a}, Corrections for tunnelling $ (t-t^{(0)})/t^{(0)} $ (green), on-site interactions $ (U-U^{(0)})/U^{(0)} $ (blue), superexchange $ (J-J^{(0)})/J^{(0)} $ (purple) and off-site interactions $ -2 V / J^{(0)} $ (red), where $ t^{(0)} $, $ U^{(0)} $, $ J^{(0)} = 4(t^{(0)})^2/U^{(0)}$ are the uncorrected values and $ t $, $ U $, $ J = 4t^2/U $ include corrections, at a lattice depth of $13\,E_R$ (solid line) and $9\,E_R$ (dashed line). \textbf{b},\textbf{c}, As a function of magnetic field $ B $ we show the transverse coupling constant $J_{xy}$ (yellow), the longitudinal coupling constant $J_{z}$ (purple) and the anisotropy $\Delta\,{=}\,J_z/J_{xy} $ (black), without corrections (dashed) and including corrections (solid) for a lattice depth of $11\,E_R$. The excluded region ($|a_{\uparrow\uparrow}|\,{>}\,700\,a_0$) is around a Feshbach resonance in the $\ket{\uparrow}$ state near $845.4\,\text{G}$.}
	\label{ED-fig:corrections}
\end{figure}

\newpage

\begin{figure}[h]
    \includegraphics[width=1\linewidth,keepaspectratio]{EDFig2.pdf}
    \caption{
        \textbf{Contrast measurement}. The images \textbf{a}-\textbf{d} show the distribution of atoms in the $\ket{\uparrow}$ state. Every pixel is a local measurement of the column density (number of atoms per unit area). The y- and z-axes are displayed in units of lattice spacings $a\,{=}\,0.532\,\mu\text{m}$. The images are projected (integrated) along the y-direction to obtain the linear density (number of atoms per unit length). The resulting 1D distributions are fitted with $f(z)\,{=}\,{g(z)\,{\cdot}\,[{1}\,{+}\,{\mathcal{C} \cos(Qz+\theta)}]/2}$ (solid line) where $ g(z) $ is a Gaussian envelope (dashed line). Examples \textbf{a}-\textbf{d} were measured at different evolution times (\textbf{a}) $t\,{=}\,0\,\hbar/J_{xy}$, (\textbf{b}) $t\,{=}\,2.3\,\hbar/J_{xy}$, (\textbf{c}) $t\,{=}\,6.3\,\hbar/J_{xy}$, (\textbf{d}) $t\,{=}\,12.0\,\hbar/J_{xy}$, for anisotropy $\Delta\,{\approx}\,0$ and wavelength $\lambda\,{=}\,10.4\,a$. The obtained contrast $\mathcal{C}(t)$ is shown in Fig.~2a. In general, we also normalize by the initial measured contrast $\mathcal{C}(0)$ to correct for finite optical imaging resolution. This is important for shorter wavelengths $\lambda$ close to the optical resolution of $3\,\mu\text{m}$, where the measured contrast $\mathcal{C}(t)$ is reduced compared to the real contrast $c(t)\,{=}\,\mathcal{C}(t)/\mathcal{C}(0)$. }
	\label{ED-fig:imageFit}
\end{figure}

\newpage

\begin{figure}[h]
    \includegraphics[width=1\linewidth,keepaspectratio]{EDFig3.pdf}
    \caption{
        \textbf{Dispersion relations}. For all positive anisotropies $\Delta\,{\geq}\,0$, the time evolution of the contrast $c(t)$ shows a damped oscillatory component, in addition to the overall exponential decay. For larger $ \Delta $, the oscillations become smaller. \textbf{a}, Decay and weak oscillation of the contrast  at the isotropic point $\Delta\,{\approx}\,1$. The observed period is $T\,{=}\,2\pi/\omega\,{=}\,37(2)\,\hbar/J_{xy}$ for $ \lambda = 10.4a $. \textbf{b}, Decay curves for different wavelengths $\lambda\,{=}\,15.7\,a$, $13.4\,a$, $11.7\,a$, $9.4\,a$, $7.8\,a$, $6.3\,a$ show how the measured oscillation frequencies $\omega$ vary with $\lambda$. \textbf{c}, They follow linear dispersion relations $\omega(Q)\,{=}\,v Q$ shown for $\Delta\,{=}\,{-0.12}$ (red), $0.35$ (orange), $0.78$ (yellow), $1.01$ (blue) and $1.27$ (light blue). \textbf{d}, The obtained velocities $v$ decrease with increasing anisotropy $\Delta$. For $\Delta\,{=}\,1.58$ (open symbol) oscillations were very small and the measurement was limited to large values of $Q$, which precluded recording a full dispersion relation.}
	\label{ED-fig:dispersion}
\end{figure}

\newpage

\begin{figure}[h]
    \includegraphics[width=1.0\linewidth,keepaspectratio]{EDFig4.pdf}
    \caption{
        \textbf{Effect of finite hole concentration} shown for $\Delta\,{\approx}\,0$ and $\lambda\,{=}\,10.4\,a$.  By varying the thermal fraction $N_\text{th}/N$ of the Bose-Einstein condensate before it is loaded into the optical lattice, we vary the energy and entropy of the atoms in the spin chain, and therefore the concentration of holes. (For our conditions, doubly occupied sites have higher energies than holes).  \textbf{a}, Decay curves $c(t)$ for varying hole concentrations ranging from low (blue) to high (red) thermal fraction. Solid lines are fits $c(t)\,{=}\,\left(a_0\,{+}\,b_0\cos\omega t\right) e^{-t/\tau}\,{+}\,c_0$. \textbf{b}, The background contrast $c_0$ increases monotonously with thermal fraction $N_\text{th}/N$. A linear fit (solid line) extrapolates to $c_0\,{=}\,0.01(2)$, consistent with zero, for $N_\text{th}/N\,{=}\,0$. This suggests that all of the background contrast is due to hole excitations.  Higher hole concentrations suppress (\textbf{c}) the oscillating fraction $b_0/(a_0\,{+}\,b_0)$, (\textbf{d}) don't affect the oscillation period $T\,{=}\,2\pi/\omega$ and (\textbf{e}) decrease the decay time $\tau$ only slightly.  The behavior shown in \textbf{c} and \textbf{e} is most likely caused by mobile holes in the central part of the Mott insulator -- indeed, numerical simulations of the $t$-$J$-model reproduce such effects (Fig.~2a). On the other hand, a finite background contrast (\textbf{b}) is likely caused by immobile holes located in the outer parts of the atom distribution where first-order tunneling is suppressed by the gradient of the (harmonic) trapping potential \cite{tiltedMottInsulator}. Since they disrupt spin transport, we expect that the imprinted spin modulation in these regions will not (or only very slowly) decay. This region is visible as a shell of low atomic density surrounding the Mott insulator in the in-situ images for (\textbf{f}) large hole concentration and absent for (\textbf{g}) low hole concentration. For different evolution times $t$, the plots in \textbf{f} and \textbf{g} show (on the left) the local contrast as function of distance $r$ from the center of the atom cloud, obtained from in-situ images (on the right). From top to bottom, the evolution times are $t\,{=}\,0\,\hbar/J_{xy}$ (top image), $2.7\,\hbar/J_{xy}$ and $21.7\,\hbar/J_{xy}$ (bottom image). The dashed line indicates a contour of constant radius (\textbf{f}) $r\,{=}\,30\,a$ and (\textbf{g}) $r\,{=}\,20\,a$. }
	\label{ED-fig:holes}
\end{figure}

\newpage

\begin{figure}[h]
    \includegraphics[width=0.7\linewidth,keepaspectratio]{EDFig5.pdf}
    \caption{
        \textbf{Decay behavior as a function of anisotropy} ranging from (\textbf{a}) negative to (\textbf{b}) positive, for one fixed wavelength $\lambda\,{=}\,10.4\,a$. Using $\Delta\,{\approx}\,0$ as a reference point, we show how the temporal profile of the decay curve $c(t)$ changes when we introduce positive or negative interactions. Every data point is an average of two measurement at lattice depths $9\,E_R$ and $11\,E_R$. \textbf{a}, from bottom to top: $\Delta\,{=}\,{-0.12}$ (red), ${-0.59}$ (pink), ${-0.81}$ (yellow), ${-1.02}$ (blue), ${-1.43}$ (green), ${-1.79}$ (purple). \textbf{b}, from bottom to top: $\Delta\,{=}\,{-0.13}$ (red), $0.08$ (purple), $0.35$ (pink), $0.55$ (orange), $0.78$ (yellow), $1.01$ (blue), $1.27$ (light blue), $1.58$ (green). Regardless of the sign, for increasing $|\Delta|$ the decay always slows down and the revivals damp out more quickly. However, there is a big difference in how this slowdown happens: \textbf{b}, For increasing positive interactions $\Delta\,{>}\,0$, the initial rate of decay increases continuously. \textbf{a}, In contrast, for all negative interactions $\Delta\,{<}\,0$, the initial rate of decay stays constant (and is ballistic), coinciding with the $\Delta\,{\approx}\,0$ case. It is only after a critical time $t_0$ that the decay suddenly starts slowing down (and becomes diffusive) for times $t\,{>}\,t_0$. This critical time $t_0$ decreases with increasing negative interaction strength $|\Delta|$.}
		\label{ED-fig:decay_anisotropies}
\end{figure}

\newpage

\begin{figure}[h]
    \includegraphics[width=1\linewidth,keepaspectratio]{EDFig6.pdf}
    \caption{
        \textbf{Collapse of decay curves} for positive anisotropies $\Delta\,{>}\,0$. All decay curves $c(t)$ for different wavelengths $\lambda\,{=}\,15.7\,a$, $13.4\,a$, $11.7\,a$, $10.4\,a$, $9.4\,a$, $8.5\,a$, $7.8\,a$, $7.2\,a$, $6.7\,a$  collapse very well into a single curve for all evolution times $t$, when time units are rescaled by $\lambda^\alpha$, where the exponent $\alpha$ is a function of anisotropy $\Delta$, both for experiment (points) and theory (solid lines). Experimental points were measured for lattice depths: $9\,E_R$ (red), $11\,E_R$ (blue), $13\,E_R$ (yellow). {\textbf{a}, \textbf{b},} Ballistic regime ($\alpha\,{=}\,1$), \textbf{c}, super-diffusion ($\alpha\,{=}\,1.5$), \textbf{d}, diffusion ($\alpha\,{=}\,2$), \textbf{e}, \textbf{f}, sub-diffusion ($\alpha\,{=}\,2.5,\,3$ for experiment and $\alpha\,{=}\,3.5,\,4.5$ for numerical simulations. In \textbf{f}, experiments covered a reduced range $\lambda\,{\leq}\,10.4\,a$). However, the experimentally measured oscillation frequencies $\omega$ still follow a linear dispersion relation for all anisotropies $\Delta\,{\geq}\,0$ (Extended Data Fig.~\ref{ED-fig:dispersion}). Therefore the collapse is not perfect. Nevertheless, away from the ballistic regime $\alpha\,{\approx}\,1$, such oscillations are small and therefore only lead to a small deviation from the collapse behavior. Note also the rather different timescales in experiments and simulations for $\Delta\,{>}\,1$.}
	\label{ED-fig:universal}
\end{figure}

\newpage

\begin{figure}[h]
    \includegraphics[width=1\linewidth,keepaspectratio]{EDFig7.pdf}
    \caption{
        \textbf{Collapse at short times} for negative anisotropies $\Delta\,{<}\,0$. All decay curves $c(t)$ for different wavelengths $\lambda$ collapse into a single curve at early times, when time units are rescaled by $\lambda$ (indicating ballistic behavior). For later times the decay is diffusive with different scaling. \textbf{a}-\textbf{c}, Theory (from top to bottom: $\lambda\,{=}\,31.3\,a$, $23.5\,a$, $18.8\,a$, $15.7\,a$, $13.4\,a$, $11.7\,a$, $10.4\,a$, $9.4\,a$, $8.5\,a$, $7.8\,a$, $7.2\,a$, $6.7\,a$, $6.3\,a$). The dotted lines are exponential fits $e^{-t/\tau_\text{II}}$ to the diffusive regime and the time constants $\tau_\text{II}$ are shown in Extended Data Fig.~\ref{ED-fig:powerlaw_theory}a. \textbf{d}-\textbf{f}, Experiment (from top to bottom: $\lambda\,{=}\,18.8\,a$, $13.4\,a$, $10.4\,a$, $8.5a\,$, $7.2\,a$, $6.3\,a$) shown for a lattice depth of $11\,E_R$. The dashed line indicates the ballistic case $\Delta\,{\approx}\,0$ (see Fig.~2c).}
	\label{ED-fig:negative_anisotropies}
\end{figure}

\newpage

\begin{figure}[h]
    \includegraphics[width=0.9\linewidth,keepaspectratio]{EDFig8.pdf}
    \caption{
        \textbf{Power law scalings (theory) and diffusion coefficients}. \textbf{a},\textbf{b}, Decay time constants $\tau$ for different anisotropies $\Delta$ ranging from (\textbf{a}) negative to (\textbf{b}) positive. Numerical results are shown in \textbf{a} for $\Delta\,{=}\,{-1}$ (blue), ${-1.5}$ (green), ${-2}$ (purple) and in \textbf{b} for $\Delta\,=\,0$ (red), ${0.5}$ (orange), ${0.85}$ (yellow), ${1}$ (blue), ${1.5}$ (green). Solid lines are power law fits (to the filled symbols). Open symbols are excluded from the fit due to finite size effects. Crossed symbols are results from $t$-$J$-model simulations including $5\,\%$ hole fraction. Fitted power law exponents are shown in Fig.~4c of the main text. For positive anisotropies $\Delta\,{\geq}\,0$ the decay time $\tau$ is defined as $\tau\,{=}\,\tau'/\ln(1/0.60)$ with $c(\tau')\,{=}\,0.60$. For negative anisotropies $\Delta\,{\leq}\,0$, the decay time $\tau_\text{I}$ for (I) short times is defined as $\tau_\text{I}\,{=}\,10\tau_\text{I}'$ with $c(\tau_\text{I}')\,{=}\,0.90$. For (II) longer times, the decay time $\tau_\text{II}$ is obtained from exponential fits $e^{-t/\tau_\text{II}}$ (see dotted curves in Extended Data Fig.~\ref{ED-fig:negative_anisotropies}). \textbf{c}, Diffusion coefficients obtained from theory (open symbols) and experiment (filled symbols). For negative anisotropies $\Delta\,{<}\,0 $, values were determined from quadratic power law fits $1/\tau\,{=}\,DQ^2$ to the data points in \textbf{a} (theory) and Fig.~3a (experiment) for the diffusive regime (II). Note that for $\Delta\,{\geq}\,0$ the system is only diffusive for $\Delta\,{=}\,{+1}$, as shown in \textbf{b} (theory) and Fig.~3b (experiment). }
	\label{ED-fig:powerlaw_theory}
\end{figure}

\newpage

\begin{figure}[h]
    \includegraphics[width=\linewidth,keepaspectratio]{EDFig9.pdf}
    \caption{
        \textbf{Finite size effects from the initial phase} of the spin helix state, illustrated here for $\Delta\,{=}\,0$ and $\lambda\,{=}\,10.4\,a$. \textbf{a}, The time evolution of the contrast $c(t)$ depends strongly on the initial phase $\theta$. The dynamics of the local magnetization $\langle S^z_i(t) \rangle$ for phases (\textbf{b}) $\theta\,{=}\,0$ and (\textbf{c}) $\theta\,{=}\,\pi/2$ reveals that this arises due to the reflection of ballistically propagating magnetization off the ends of the chain. Depending on the initial phase of the spin helix, the reflected magnetization interferes constructively or destructively with the  pattern of the bulk magnetization. }
	\label{ED-fig:initial_phase}
\end{figure}

\newpage

\begin{figure}[h]
    \includegraphics[width=\linewidth,keepaspectratio]{EDFig10.pdf}
    \caption{
        \textbf{Finite size effects from the chain length}. \textbf{a}, Contrast $c(t)$ obtained after a weighted average over all different chain lengths between $L = 0$ and $44$, as illustrated in (\textbf{b}), for $\Delta\,{=}\,0$ and $\lambda\,{=}\,10.4\,a$. The dynamics shows almost no dependence on the phase $\theta$ (red, yellow, blue) in contrast to a single chain length $L\,{=}\,40$ (Extended Data Fig.~\ref{ED-fig:initial_phase}a). Also overlaid are the contrast for a fixed chain length ($L\,{=}\,40$) averaged over all initial phases $0\,{\leq}\,\theta\,{<}\,2\pi$ (black solid line), and averaged over only the two phases $\theta\,{=}\,0$ and $\pi/2$ (black dashed line). The close agreement implies that averaging over either chain lengths or phases suppresses the dependence on initial or boundary conditions. \textbf{b}, A cut through the spherical Mott insulator with diameter $L_\text{max}\,{=}\,44\,a$ (as in the experiment) illustrates the distribution of different chain lengths (oriented along the z-direction). Averaging the local magnetization $S^z$ over the x- and y-direction provides a 1D magnetization profile (bottom panel), which is an average over all chains. \textbf{c}, The number of chains with length $L$ is given by $(\pi/2)(L/a)$.  The total number of chains is $\pi(L_\text{max}/2a)^2\,{\approx}\,1500$. \textbf{d}, The number of atoms in chains with length $L$ is given by $(\pi/2)(L/a)^2$. Since the contribution of each chain to the imaging signal is proportional to the atom number in the chain, the relevant average over chain lengths is weighted by the atom number and is $\langle L \rangle\,{=}\,3/4\,{\cdot}\, L_\text{max}\,{=}\,33\,a$. }
	\label{ED-fig:chain_lengths} 
\end{figure}


\end{document}


