\subsection{Experimental system}

\begin{figure*}[!htb]
\centering
\includegraphics[width=140mm]{FigureS1.pdf}
\caption{Spectroscopy of spin exchange in double wells.
(a) Illustration of the spin-dependent potentials. For two spin configurations in a double well, $\ket{\downarrow,\uparrow}$ and $\ket{\uparrow,\downarrow}$, their energies can be well controlled by tuning the offset value $\delta_s$. For instance, at negative side, the state $\ket{\uparrow,\downarrow}$ has lower energy than $\ket{\downarrow,\uparrow}$.
(b) Spin superexchange under different staggered magnetic potentials.
The states are initialized to $\ket{\uparrow,\downarrow}$ in double wells, corresponding to the staggered magnetization of $M_z=-0.5$. Then, we quench the lattice depth to allow spin exchange dynamics at the coupling strength of $J=26.8(5)$ Hz, holding for 20 ms. The spectroscopy shows a symmetric spin-dependent effect. The solid curve is a guide for the eyes. Error bars are standard deviations throughout this supplementary material.}
\label{figS1}
\end{figure*}

Our experiment begins with a nearly pure Bose-Einstein condensate of $^{87}$Rb staying on the hyperfine state of  $5S_{1/2}\left|F=1,m_F=-1\right>$.
Then, we compress the atoms along the $z-$axis by increasing the confinement, followed by which atoms are loaded  into a single node of a 3-$\mu m$ spacing blue-detuned lattice.
Thus, we create a quasi-2D quantum gas.
The confinement in the $x-y$ plane is provided by a red-detuned optical dipole trap propagating along the $z$ direction.
To prepare a low-entropy Mott insulator, we implement a newly-developed cooling technique in optical lattices.
The cooling concept is that the entropy of these environment-isolated atoms is transferred into contacting superfluid reservoirs, which are removed from the system later.
In our region of interest (36 copies of 70-site chains), the probability of unity-filling for the cooled Mott insulator is 0.992(1).



To initialize the N\'eel order along the $x$-direction, we employ a site-selective spin addressing technique to overcome the optical diffraction limit.
The atoms in the Mott insulator state are transferred into an spin-dependent optical superlattice \cite{Yang:2017a}, which breaks the odd-even symmetry of the lattice sites and thereby separates the atoms into two subsystems.
For odd and even subsystems, the lattice potentials are controlled to split their transition frequencies between $\ket{\downarrow}$ and $\ket{\uparrow}$ hyperfine levels.
Applying an edited microwave pulse, we flip the atoms of odd sites into $\ket{\uparrow}$ with an efficiency of 0.995(3) through the rapid adiabatic passage.
For each site, the preparation fidelity of the N\'eel state can be estimated as 0.987(3).
Therefore, in the 70-site chain, a defect-free N\'eel order is achieved with a probability of $\sim$ 40\%.

\subsection{Staggered potentials}



The superlattice naturally forms staggered potentials for our ground-band ultracold atoms.
Here, we make use of two distinct effects of the superlattice, which are the spin-dependent and the  spin-independent potentials.
Such optical potential can be expressed as,

\begin{equation}
V(x) = V_{s,\sigma}\cos^2(k_s x+\theta_\sigma) + V^{\text{offset}}_{s,\sigma} - V_l\cos^2(k_s x/2 + \theta),
\end{equation}
where $k_s=2\pi/\lambda_s$ is the wave vector of the short-lattice laser.
Here, $\theta_\sigma$ is the spin-dependent phase shift of the short-lattice, where $\sigma$ denotes the two spin states.
The original trap depth of short-lattice $V_s$ is also altered to two spin-dependent terms $V_{s,\sigma}$, and $V^{\text{offset}}_{s,\sigma}$.
At $\theta_{\sigma} = 0$, the superlattice is simplified to the spin-independent form, $V^{\text{offset}}_{s,\sigma} =0$.
By tuning the applying voltage of a electro-optical modulator, we can control the spin-dependent energy with a bandwidth larger than 20 kHz \cite{Yang:2017a}.

The spin-dependent energy shift can be precisely controlled in our system.
To calibrate the zero-energy $\delta_s =0$, we perform a spectroscopy measurement of spin-exchange dynamics in isolated double-wells.
We prepare the atoms in $\ket{\uparrow,\downarrow}$ state and then quench the superlattice to $V_s=16.8(1) E_r$, $V_l=10.3(1)$ Er and $\theta=0$.
Here $E_r = h^2/(2 m \lambda_s^2)$ is the recoil energy with $h$ the Planck constant and $m$ the atomic mass.
In this balanced structure, the resonating superexchange frequency within double wells is $26.8(5)$ Hz.
The inter-well coupling is forbidden by the long-lattice barrier.
The degeneracy of the states $\ket{\uparrow,\downarrow}$ and $\ket{\downarrow,\uparrow}$ is lifted by a nonzero spin-dependent effect.
As we scan the $\delta_s$ and monitor the staggered magnetization $M_z$, a clear spectroscopy of the spin exchange dynamics is observed, as shown in Fig. \ref{figS1}.
The position of the central peak can be determined with a precision of $\sim$0.4 Hz.
In addition, the energy offsets of the side peaks agree well with our calculations based on the lattice structure.


\subsection{Establishing the isotropic Heisenberg Hamiltonian}

We design a special 1D superlattice along the $x$-direction for realizing the isotropic antiferromagnetic Heisenberg model.
When $V_l \ll V_s$, the inter-well coupling becomes comparably strong as the intra-well coupling.
In this 1D chain, there are two types of links connecting the neighboring sites, which can be denoted as odd-even (type I) and even-odd (type II) links.
If $\theta = \pi/4$, these two types of links are identical to each other.
While the spin-dependent effect will disappear in this configuration.
To enable the spin-dependent term, we set the superlattice phase to $\theta = \pi/3$ and keep the long-lattice at $V_l = 0.44(1) E_r$.
The tunneling strength on the type II links $t_2$ is roughly 10 percent larger than $t_1$ on the type I links \cite{Walters:2013s}.

In the second-order perturbation theory, the spin exchange coupling along these two types of links are,

\begin{equation}
J_{1}=-\frac{4t_{1}^2U}{U^2-(\delta_0-\Delta)^2}, \ \
J_{2}=-\frac{4t_{2}^2U}{U^2-(\delta_0+\Delta)^2}.
\end{equation}
Here, $\delta_0$ and $\Delta$ represent energy offsets between neighboring sites due to the spin-independent term and a linear potential, respectively.
$t_1$ and $t_2$ are tunneling strengths along the two types of nearest-neighbor links.
Since the lattice depths are almost the same for the two spin states, the tunneling strengths are considered as spin-independent.
$U$ corresponds to the onsite-interaction for all the spin state.
Here we neglect the slight spin-dependence of scattering lengths for $^{87}$Rb and adopt $U_{\uparrow\downarrow}=U_{\uparrow\uparrow}=U_{\downarrow\downarrow}=U$.
When $U < \delta_0 - \Delta \ \text{and} \ \delta_0 + \Delta$, we reverse the sign of $J_1$ and $J_2$ to positive, establishing an antiferromagnetic spin-exchange coupling \cite{Trotzky:2008}.

The difference between the spin exchange couplings $J_1$ and $J_2$ is compensated by the state-independent linear potential.
Since our $x$-axis lies at an angle of 4 degree related to the horizontal plane, the gravity leads to this constant 57-Hz/site gradient along the 1D chain.
Consequently, the discrepancy between these two types of exchange coupling becomes roughly 2 percent [see Figure \ref{figS2}(c)], one order of magnitude lower than the case without this compensation ($J_1$ and $J_2$ would have 20 percent difference in the absence of the linear gradient).

In our experimental conditions, the Hubbard parameters satisfy such relations $t_1, t_2 \ll \delta_0, U$, as shown in Fig. \ref{figS2}.
This Bose-Hubbard Hamiltonian can be mapped into a spin model as,

\begin{equation}
\hat{H}=\sum_{i}(J_{1}\hat{S}_{2i-1}\cdot\hat{S}_{2i} + J_{2}\hat{S}_{2i}\cdot\hat{S}_{2i+1}) + \delta_s\sum_{j}(S_{2j}^z-S_{2j+1}^z).
\end{equation}
Neglecting the slight difference between $J_1$ and $J_2$, we can adopt an isotropic Heisenberg model with the spin-exchange strength $J$, which is the geometric average of $J_1$ and $J_2$, as $J=\sqrt{J_{1}J_{2}}$.

We can understand the model in two limits.
In the limit of $\delta_s/J \rightarrow \pm \infty$, the two type N\'eel orders are the ground states, which are product states.
We should stress that even though atoms tend to occupy the low-energy site of the staggered potential in the ground state of the one-dimensional Hubbard model, the atoms cannot overcome the staggered potential via tunneling in our state-preparation time and the system can be faithfully described by the spin model.
In another limit, $\delta_s/J =0$, unlike the Ising-type spin model, this Heisenberg model owns a single energetic ground state $\ket{\psi_{\text{AF}}}$.
The SU(2) symmetry of this isotropic Hamiltonian leads to the same rotational symmetry of the ground state.

\begin{figure*}
\centering
\includegraphics[width=150mm]{FigureS2.pdf}
\caption{ From Bose-Hubbard model to Heisenberg spin model.
(a) The blue and black curves represent the tunneling strengths and spin-exchange coupling, respectively.
The lighter colors correspond to the type I links, representing $t_1$ and $J_1$.
While the darker colors correspond to the type II links , representing $t_2$ and $J_2$.
The discrepancy between the exchange coupling strengths has been mostly compensated by the gravity potential.
(b) The blue and black curves shows the on-site interaction $U$ and the spin-independent staggered potential $\delta_0$ under different short lattice depths.
(c) The ramping curves of the spin-exchange strength $J$ and the staggered magnetic potential $\delta_s$ in the adiabatic passage.}
\label{figS2}
\end{figure*}

\subsection{Adiabatic passage}

Starting from the initial N\'eel-ordered state $\left|\psi_0\right>=\left|\uparrow\downarrow\uparrow\downarrow\cdots\right>$, we experimentally optimize the ramping curves to minimize the non-adiabaticity to approach the Heisenberg antiferromagnet $\ket{\psi_{\text{AF}}}$.
The $y-$lattice with 42.0(3) $E_r$ isolates the system into 36 copies of 70-sites 1D chains.
The superlattice phases are switched to $\theta = \pi/3$ and $\theta_{\ket{\uparrow}} = -2 \theta_{\ket{\downarrow}} = 0.005 \pi$ before ramping down the short-lattice potential.

Then, we slowly ramp the lattice potentials in terms of the Hubbard parameters.
At 0 ms, the staggered magnetic gradient is $\delta_s/h=-17.6(3)$ Hz, and the coupling strength is $J/h= 0.24(1)$ Hz.
During the 60 ms ramping, we separate the curve into two stages with 30 ms interval, as shown in Fig. \ref{figS2}(c).
The exchange coupling is turned on gradually, and simultaneously the staggered gradient is also increased.
At 60 ms, $\delta_s$ is slightly larger than 0, but the zero value of the staggered magnetization $M_z = 0.00(2)$ indicates the achievement of a Heisenberg antiferromagnet.
In another 60 ms, we verify the coherence of the system by slowly ramping back the lattice potentials and increasing  $\delta_s$ oppositely.
To detect the property of the state at each time, we freeze the state by quenching the lattice barrier to disable the atom dynamics.

\subsection{Nearest-neighbor spin correlation}

The $in\ situ$ atomic densities are detected via an absorption imaging with the optical resolution of $\sim$1 $\mu$m.
Even though it is one kind of destructive imaging techniques, we can take two successive images for recording both the spin states, the first image for $\ket{\uparrow}$ and the second image $\ket{\downarrow}$.
We choose a region of interest (ROI) containing 36 copies of 70-site 1D chains in the center of the atom cloud for the statistical analysis.

The nearest-neighbor correlation function is defined as
\begin{equation}
C(d=1)= \frac{4}{N}\sum_{j}\left(\left<\hat{S}_j^z\hat{S}_{j+1}^z\right>-\left<\hat{S}_j^z\right>\left<\hat{S}_{j+1}^z\right>\right),
\end{equation}
where $N$ represents the number of sites.
In our spin-balanced system, the second term in the sum is equal to the square of the staggered magnetization, as $-1/N\sum_{j}\langle\hat{S}_j^z\rangle\langle\hat{S}_{j+1}^
z\rangle= M_z^2$.
For obtaining the first term in the sum, we perform correlated measurements on these two type of links respectively, as

\begin{equation}
\begin{split}
&\text{Type I :} \ \ \\
&-\frac{4}{N}\sum_{j}\langle\hat{S}_{2j-1}^z\hat{S}_{2j}^z\rangle =\frac{1}{2} \left(P_{\ket{\uparrow,\downarrow}}^{I} + P_{\ket{\downarrow,\uparrow}}^{I} -P_{\ket{\uparrow,\uparrow}}^{I}- P_{\ket{\downarrow,\downarrow}}^{I}\right),\\
\\
&\text{Type II :}  \ \ \\
&-\frac{4}{N}\sum_{j}\langle\hat{S}_{2j}^z\hat{S}_{2j+1}^z\rangle =\frac{1}{2} \left(P_{\ket{\uparrow,\downarrow}}^{II} + P_{\ket{\downarrow,\uparrow}}^{II} -P_{\ket{\uparrow,\uparrow}}^{II}- P_{\ket{\downarrow,\downarrow}}^{II}\right).
\end{split}
\end{equation}
Here, $P^{I(II)}$ represent the probabilities of spin states on the corresponding links.
Making use of the superlattice, we can combine the atoms on neighboring sites into a single lattice well.
Depending on the atom number on the combined site, we can derive the desired probabilities.
These two types of links are measured at two kinds of superlattice configurations, one has $\theta=0$ and the other has $\theta=\pi/2$.

\begin{figure*}
\centering
\includegraphics[width=140mm]{FigureS3.pdf}
\caption{Detection of the nearest-neighbor spin correlation.
(a) The detection scheme for the spin components on two types of nearest-neighbor links. In step (i), all $\ket{\uparrow}$ atoms are removed by a resonating imaging light pulse.
Then in (ii), the atoms on type I or type II links are respectively loaded into superlattice double-well units with $\theta=0$  and $\theta=\pi/2$.
In step (iii), the atoms within each double-well are combined and doublons in the long-lattice site get lost after the photoassociation collision.
Finally, in step (iv), we count the remaining atoms and deduce the probabilities of the state $\ket{\downarrow,\downarrow}$.
(b) The upper graph shows the state probabilities on type I links corresponding to the tunneling strength $t_1$, and for type II links the results are depicted in the lower graph.
In these two graphs, the red circles, the blue squares, the green triangles and the yellow inverted triangles represents the probabilities of $\ket{\uparrow,\downarrow},\ket{\downarrow,\uparrow},\ket{\uparrow,\uparrow}$ and $\ket{\downarrow,\downarrow}$, respectively.
The solid curves are numerical caculations of a 70-site spin chain.}
\label{figS3}
\end{figure*}

Not limited by our optical resolution, we probe the parity of atom number on lattice sites with a photoassociation (PA) laser.
The laser is red-detuned to the D2 line of $^{87}$Rb by 13.6 $\rm{cm}^{-1}$.
Assisted by this laser, atom pairs in hyperfine state $\ket{\downarrow}$ can be excited to the $\nu=17$ vibrational state of the $0_g^-$ long-range molecular channel \cite{Fioretti:2001s}.
This unstable molecule would decay back to free atoms and gains kinetic energy to escape from the trap.
Therefore, the ratio of atom loss is equivalent to the probability of the state $\ket{\downarrow,\downarrow}$ before the combining operation.

The detection scheme is illustrated in Fig. \ref{figS3}(a), where two sets of measurements are carried out for these two types of links.
There are several steps in the detection procedure.
The main concern is removing the atoms on $\ket{\uparrow}$ state before applying the atom combination.
The PA light has an intensity of 0.56 $\rm{W/cm}^2$ and lasts for 20 ms.
In this way, we obtain the probabilities of $P_{\ket{\downarrow,\downarrow}}^{I}$ and $P_{\ket{\downarrow,\downarrow}}^{II}$, as shown in Fig. \ref{figS3}(b).

For the other three spin-states $\ket{\uparrow,\uparrow},\ \ket{\uparrow,\downarrow},\ \ket{\downarrow,\uparrow}$, we employ the site-selective spin addressing technique to transfer them into the  $\ket{\downarrow,\downarrow}$ state before implementing the detection procedures \cite{Yang:2017a}.
For instance, we flip all of the atoms on odd sites for probing the probabilities $P_{\ket{\uparrow,\downarrow}}^{I}$ and $P_{\ket{\downarrow,\uparrow}}^{II}$.
The probabilities of these four states allow us to calculate the spin correlations $C(d=1)$.
However, all of the measurement errors would propagate into the final spin correlations.
We use the basic normalization condition $P_{\ket{\uparrow}}+P_{\ket{\downarrow}} =1$ to reduce the errors by only taking the residual atoms in the last step of Fig. \ref{figS3}(a) into calculations.

In Fig. \ref{figS3}(b), the probabilities of the four spin state for type I and type II nearest-neighbor links are depicted.
The deviation from the numerical calculations can be caused by the motional excitations in Hubbard model and the lattice heating during the adiabatic passage.
Besides, after steps of superlattice manipulation in the detection method, the accumulative inefficiency of the detection  also contributes to the reduction of the measured spin-state probabilities,


\subsection{Motional excitations}

\begin{figure*}[!htb]
\centering
\includegraphics[width=70mm]{FigureS4.pdf}
\caption{Motional excitations during the adiabatic passage. The blue (red) circles represent the ratios of the doublon excitations with identical (opposite) spin components. The solid curves are linear fits to the experimental data.}
\label{figS4}
\end{figure*}

Besides spin excitations in the state evolution, we also measure the ratios of motional excitations in these lattice sites, which are absent in the spin model.
In the Hubbard model, one particle-hole excitation emerges when an atom tunnels to its neighboring site.

The particle-hole excitations with doublons $\ket{\downarrow\downarrow}$ are detected by the same parity-projection technique making use of the PA laser.
After a global spin-flip, we can also acquire the doublons on the state $\ket{\uparrow\uparrow}$.
The particle excitations on the state $\ket{\uparrow\uparrow}$ and $\ket{\downarrow\downarrow}$ generated in the adiabatic passage are rare, since the atoms need to tunnel to its next-nearest neighboring sites.
Fig. \ref{figS4} shows that the total ratio of these excitations keeps constant at 0.01(1) in the state evolution.

However, the more probable particle-hole excitations result in the $\ket{\uparrow\downarrow}$ state, owing to the tunnelling of atoms to their neighboring sites.
From the normalization condition, the ratio of these excitations can be inferred as $1-\sum_{\mu,\nu}P_{\ket{\mu,\nu}}$, with $\mu,\nu=\ \uparrow,\downarrow$.
As depicted in Fig. \ref{figS4}, the ratio of these doublons excitations grows gradually during the 120 ms state evolution, which is 0.04(3) for our Heisenberg antiferromagnet.
The propagation of excitations is strongly suppressed by the spin-independent staggered offset $\delta_0$.
Besides, the linear potential also suppress the long-range motional transport \cite{Yang:2020s}.

\subsection{Next-nearest-neighboring spin correlation}

\begin{figure*}[!htb]
\centering
\includegraphics[width=150mm]{FigureS5.pdf}
\caption{Detection of next-nearest-neighbor spin correlation.
(a) The detection scheme for resolving the probabilities of state $\ket{\downarrow;\downarrow}$ on two types of links.
For probing the correlations on odd (even) sites, in step (i), we remove all the $\ket{\uparrow}$ atoms as well as the atoms residing on even (odd) sites.
In step (ii), the remaining atoms are transferred into the long-lattice potential.
Then in step (iii), the atoms are split into adjacent short-lattice sites with equal possibilities.
Finally, in step (iv), we combine the atoms and count the loss rate of doublons by removing them in the long-lattice potential.
(b) Probabilities of the states $\ket{\uparrow;\uparrow}$ and $\ket{\downarrow;\downarrow}$.
The upper and lower graph show the evolution of spin components for neighboring odd and even sites, respectively.
The blue circles represent the probabilities of state $\ket{\uparrow;\uparrow}$, and the red circles correspond to the state $\ket{\downarrow;\downarrow}$.
The solid curves are the numerical calculations of a 70-site spin chain.}
\label{figS5}
\end{figure*}

According to the Mermin-Wagner-Hohenberg theorem \cite{Mermin:1966s}, the one-dimensional isotropic Heisenberg Hamiltonian cannot support a long-range spin order \cite{Yang:2017s}.
As the distance grows, the magnetic correlation decays rapidly in the 1D Heisenberg antiferromagnet.
So the occurrence of nonzero next-nearest-neighbor spin correlation would be a distinctive signature of approaching a low-entropy antiferromagnet.

However, without the single-site imaging, the superlattice structure normally only provides the access to the nearest-neighbor spin correlation because of its $Z_2$ symmetry \cite{Trotzky:2010s,Dai:2016s}.
Here, owing to our precise control of the superlattice phase, we can extend the capability of the correlated spin measurement to next-nearest-neighbor links.

Figure \ref{figS5}(a) shows the procedure to detect the probability of spin state $\ket{\downarrow;\downarrow}$ along next-nearest-neighboring links, implemented on odd sites and even sites separately.
After a selective atom removal, the superlattice serves as an array of atomic splitters to redistribute the remaining atoms into their adjacent sites with equal probabilities.
Subsequently, changing the superlattice phase $\theta$ by $\pi/2$, we can transfer the $\ket{\downarrow}$ atoms originally lying two sites apart to the same long-lattice unit in a probabilistic way.
The atom loss in a PA collision can be used to deduce the probability of spin state $\ket{\downarrow;\downarrow}$.
Similarly, the $\ket{\uparrow; \uparrow}$ state on next-nearest-neighbor sites can be detected with an additional global spin-flip pulse before the step (i).

We define a next-nearest-neighbor correlator as,

\begin{equation}
C'(d=2)=P_{\ket{\uparrow;\uparrow}}+P_{\ket{\downarrow;\downarrow}},
\end{equation}
Here $P_{\ket{\uparrow;\uparrow}}$ and $P_{\ket{\downarrow;\downarrow}}$ represent the probabilities of spin state $\left|\uparrow;\uparrow\right>$ and $\left|\downarrow;\downarrow\right>$ along next-nearest-neighbor links, which equal to four times the corresponding atom density loss in the photoassociation collision.
For obtaining accurate probabilities, we subtract the contribution of the motional excitations in the PA collision.

Figure \ref{figS5}(b) shows the evolution of the probabilities $P_{\ket{\uparrow;\uparrow}}$ and $P_{\ket{\downarrow;\downarrow}}$ for odd sites and even sites separately.
Since only a quarter of atoms contribute to the signal, the probabilities measured here have significantly larger errors compared to their counterparts on nearest-neighbor links.

\subsection{Numerical simulation}

\begin{figure*}[!htb]
\centering
\includegraphics[width=140mm]{FigureS6.pdf}
\caption{Numerical calculations of the Heisenberg spin chains. (a) spin relaxation dynamics after quenching the N\'eel order to the Heisenberg spin model.
(b)(c) and (d) are respectively the numerical results of the evolution of $M_z$, $C(d=1)$ and $C'(d=2)$ in the adiabatic passage.
Different blue colors represent the chain length of 10, 12, 16, 20, 30, 40, 50 and 70, respectively.
The inset graphs show the discrepancies between numerics for different system lengths and the curve with the largest length (N = 70). The finite-size convergence is well below our detection errors.}
\label{figS6}
\end{figure*}

\begin{figure*}[!htb]
\centering
\includegraphics[width=154mm]{FigureS7.pdf}
\caption{A 70-site Heisenberg spin chain is simulated with the finite-temperature Monte Carlo method.
The nearest-neighbor spin correlation $C(d=1)$ and the entropy per particle versus temperatures are depicted in (a) and (b), respectively.
The red circles represent the case of our prepared antiferromagnet.
(c) Probabilities of the spin states in the nearest-neighbor and next-nearest-neighbor links.
The circles represent the measured probabilities of spin components, while the blue bars are the theoretical predictions of a thermalized Heisenberg antiferromagnet at $C(d=1)=-0.44$.}
\label{figS7}
\end{figure*}

The numerical data for the post-quench dynamics and the adiabatic passage of the Heisenberg models are obtained using the the $t$-DMRG algorithm in the framework of matrix product states.
In the numerical simulations, we neglect the slight difference between the two types of nearest-neighbor links and adopt the geometric average $J=\sqrt{J_{1}J_{2}}$ as the spin-exchange strength throughout the chain.
For the implementation, we use the OpenMPS library \cite{Jaschke:2018s} and fix the maximal bond dimension to $D=1500$.
Convergence are achieved at a time step of $3\times10^{-4}$ s and $10^{-2}\ h/J$ for the adiabatic passage and the quench dynamics, respectively, with truncation threshold of $10^{-6}$ per time step.
Finite-size effects are also investigated for several chain lengths.
As show in Fig. \ref{figS6}, the results for system size from $L=40$ to $70$ reveal similar post-quench dynamics, and converge with our experimental results in statistical errors.

For estimating the temperature and the entropy of our state, we numerical calculate the Heisenberg model using a finite-temperature quantum Monte Carlo method (QMC) \cite{Bauer:2011s}.
All the simulations are carried out in the grand-canonical ensemble on open chains of size with $N=70$, which is checked to be large enough to eliminate the finite-size effect.
After $10^5$ times of thermalization, we continue to perform the QMC sweeps for at least $5\times10^6$ times to get the results under each setting.
We choose the temperature $k_B T/J$ from 0.2 to 0.8 with a 0.03 interval.
From the nearest-neighbor correlation of interest [$C(d=1)=-0.44(3)$], we derive the temperature of our system to be 0.53(7)$J$, see Fig. \ref{figS7}(a).
The entropy per particle is obtained using another QMC method, which is called quantum Wang-Landau algorithm \cite{Wang:2001s}.
The cutoff of series expansion is set to be 1000.
As shown in Fig. \ref{figS7}(b), the entropy per particle is 0.42(5)$k_B$, which is deduced from the temperature of the spin chains.

Based on the comparison between the measured probabilities of the spin state and the prediction of the thermalized state, our system is inferred to be a thermalized state, see Fig. \ref{figS7}(c).


