Realization of a Bosonic Antiferromagnet

Our experiment begins with a nearly pure Bose-Einstein condensate of ⁸⁷Rb prepared in the hyperfine state of $5S_{1/2}\left|F=1,m_{F}=-1\right>$. Then, we increase the confinement along the $z-$axis and load the atoms into a single node of a 3-$\mu m$ spacing blue-detuned lattice to 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. In the cooling process, the entropy of the Mott insulator is transferred into the contacting superfluid reservoirs. In our region of interest (36 copies of 70-site chains), the probability of unity-filling for the Mott insulator is 0.992(1).

To initialize the Néel order along the $x$-direction, we employ a site-selective spin addressing technique to overcome the optical diffraction limit. The atomic Mott insulator is transferred into an spin-dependent superlattice 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. By applying an edited microwave pulse, we flip the atoms of odd sites into $\ket{\uparrow}$ with an efficiency of 0.995(3) through a rapid adiabatic passage. For each site, the preparation fidelity of the Néel state can be estimated as 0.987(3). Therefore, in the 70-site chain, a defect-free Néel order is achieved with a probability of $\sim$ 40%.

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 effects. Such an optical potential can be expressed as,

where $k_{s}=2\pi/\lambda_{s}$ is the wave vector of the short-lattice laser; $\theta_{\sigma}$ is the spin-dependent phase shift of the short-lattice, and $\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 an electro-optical modulator, we can control the spin-dependent energy with a bandwidth larger than 20 kHz Yang:2017a .

The spin-dependent energy shift can be precisely controlled in our system. To calibrate the zero-energy point $\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)E_{r}$ and $\theta=0$. Here, $E_{r}=h^{2}/(2m\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 Extended Data Fig. 1. 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.

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. At $\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 depth 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 Walters:2013s .

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

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 independent of spin components. $U$ corresponds to the on-site interaction for all the spin states. Here we neglect the slight spin-dependence of scattering lengths for ⁸⁷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 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 a constant 57-Hz/site gradient along the 1D chain. Consequently, the discrepancy between $J_{1}$ and $J_{2}$ is less than 2 percent [see Extended Data Fig. 2(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 such a linear potential).

Under our experimental conditions, the Hubbard parameters satisfy such relations $t_{1},t_{2}\ll\delta_{0},U$, as shown in Extended Data Fig. 2. We should stress that even though atoms tend to occupy the low-energy sites owing to the non-uniform potential along the 1D chain, the atoms cannot overcome the staggered potential via tunneling in our experimental relevant time scale. Therefore, the Bose-Hubbard Hamiltonian can be mapped into a spin model as,

By neglecting the slight difference between $J_{1}$ and $J_{2}$, we can adopt the geometric mean of them as $J=\sqrt{J_{1}J_{2}}$, and write the Hamiltonian as an isotropic Heisenberg model.

We can understand the spin model in two extreme cases. In the limit of $\delta_{s}/J\rightarrow\pm\infty$, the two types of Néel orders are the ground states of the system. In another limit, $\delta_{s}/J=0$, unlike the Ising-type spin model, this Heisenberg model has a single energetic ground state $\ket{\psi_{\text{AF}}}$. The SU(2) symmetry of this isotropic Hamiltonian leads to the same rotational symmetry of this ground state.

Starting from the initial Néel-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-site 1D chains. The superlattice phases are switched to $\theta=\pi/3$ 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 intervals, as shown in Extended Data Fig. 2(c). The coupling strength and the staggered gradient are increased gradually during the sweep. 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 ramp back the lattice potentials slowly and increase $\delta_{s}$ oppositely. The system is driven towards the opposite Néel order with $M_{z}=0.27(1)$, implying the low-entropy property of the antiferromagnet. To detect the property of the state at each time, we freeze the state by quenching the lattice barrier to disable the atom dynamics.

Starting from the prepared antiferromagnet, we apply another adiabatic protocol to the system by ramping reversely $\delta_{s}$ after 60 ms to a large negative value (see Extended Data Fig. 3). Finally, $M_{z}$ can return to $-0.24(2)$, which means the decoherence and nonadiabaticity of the system can be well distinguished.

We quantify the translational symmetry of the antiferromagnet. Extended Data Fig. 4(a) shows that the measured $M_{z}$ does not gain any spatial periodicity. Such phenomena is consistent with the translational-invariant property of the antiferromagnet. Besides, we hold the atoms in the strong coupling regime after the 60-ms adiabatic passage, where $J/h=26.0(5)$ Hz and $\delta_{s}=0$. As shown in Extended Data Fig. 4(b), $M_{z}$ maintains almost constant and does not exhibit any oscillations, which suggests that the created state is close to thermal equilibrium.

The $in\ situ$ atomic densities are detected via an absorption imaging with $\sim$1 $\mu$m spatial resolution. Even though the imaging process is destructive, we can record the density distributions of both spin states ($\ket{\uparrow}/\ket{\downarrow}$) in two successive images. We choose a region of interest (ROI) containing 36 copies of 70-site 1D chains in the center of the atomic cloud for the statistical analysis.

The nearest-neighbor correlation function is defined as,

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}$. The first term in the sum is obtained by averaging the correlated measurements on two types of links, as $\frac{1}{2}\left(\langle\hat{S}_{2j-1}^{z}\hat{S}_{2j}^{z}\rangle+\langle\hat{S}_{2j}^{z}\hat{S}_{2j+1}^{z}\rangle\right)$. The correlation functions are given as following,

Here, $P^{\text{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. These two types of links are measured at two kinds of superlattice configurations, one has $\theta=0$ and the other has $\theta=\pi/2$. Depending on the atom number on the combined site, we can derive the desired probabilities.

We probe the parity of site occupations by applying a photoassociation (PA) laser to the atomic sample. The laser is red-detuned to the D2 line of ⁸⁷Rb by 13.6 $\rm{cm}^{-1}$, which has an intensity of 0.56 $\rm{W/cm}^{2}$ and lasts for 20 ms. 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 Fioretti:2001s . This unstable molecule would decay back to free atoms and gain 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 Extended Data Fig. 5(a), where two sets of measurements are carried out for these two types of links. In this way, we obtain the probabilities of $P_{\ket{\downarrow,\downarrow}}^{\text{I,II}}$, as shown in Extended Data Fig. 5(b).

For the other 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 Yang:2017a . For instance, we flip all of the atoms on odd sites for probing the probabilities $P_{\ket{\uparrow,\downarrow}}^{\text{I}}$ and $P_{\ket{\downarrow,\uparrow}}^{\text{II}}$. Finally, the probabilities of these four states allow us to calculate the spin correlation $C(d=1)$. However, all of the measurement errors would propagate into the final results. 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 Extended Data Fig. 5(a) into calculations.

In Extended Data Fig. 5(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 is mainly caused by the motional excitations in Hubbard model and the lattice heating during the adiabatic passage. Besides, the inefficiency of the state manipulations accumulates and finally reduces the detected spin correlations.

Besides spin excitations, the particle-hole excitations emerge when the atoms tunnel to their neighboring sites. Such motional excitations are detected by counting the ratio of doublons $\ket{\downarrow\downarrow}$ making use of the parity-projection technique. Using a global spin-flip, we can acquire the doublons on the state $\ket{\uparrow\uparrow}$ as well. However, atoms have to tunnel to their next-nearest-neighbor sites in order to create such identical-spin doublons. Extended Data Fig. 6 shows that the excitations to the state $\ket{\uparrow\uparrow}$ and $\ket{\downarrow\downarrow}$ are generally rare [1(1) percent] in the adiabatic passage.

The more probable motional excitations are original from the atoms that tunnel to their neighboring sites and therefore result in the $\ket{\uparrow\downarrow}$ doublons. From the normalization condition, the ratio of these excitations can be deduced as $1-\sum_{\mu,\nu}P_{\ket{\mu,\nu}}$, with $\mu,\nu=\ \uparrow,\downarrow$. As shown in Extended Data Fig. 6, the ratio of this type of motional excitations grows gradually during the state evolution, which is 0.04(3) for our target state. The propagation of such excitations is strongly suppressed by the spin-independent staggered offset $\delta_{0}$. Besides, the linear potential also suppresses the long-range motional transport Yang:2020s .

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

However, without a site-resolved imaging, the superlattice structure normally only provides the access to the nearest-neighbor spin correlation because of its $Z_{2}$ translational symmetry 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.

Extended Data Figure 7(a) shows the procedure of detecting the probability of spin state $\ket{\downarrow;\downarrow}$ in next-nearest-neighboring links. 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,

where $P_{\ket{\uparrow;\uparrow}}$ and $P_{\ket{\downarrow;\downarrow}}$ represent the probabilities of spin state $\left|\uparrow;\uparrow\right>$ and $\left|\downarrow;\downarrow\right>$ on next-nearest-neighbor links. These probabilities equal to four times of the corresponding atom losses in the PA collision. For determining the probabilities accurately, we subtract the contribution of the motional excitations in the PA collision.

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

Based on the Heisenberg spin model, the numerical results for both the post-quench dynamics and the adiabatic passage are calculated using the $t$-DMRG algorithm in the framework of matrix product states. In the numerics, we neglect the slight difference between the spin-exchange couplings of the two types of links and adopt their geometric mean $J=\sqrt{J_{1}J_{2}}$ in the spin model. We implement the methods using the OpenMPS library Jaschke:2018s and fix the maximal bond dimension to $D=1500$. Convergence are achieved at a time step of $10^{-2}\ h/J$ and $3\times 10^{-4}$ s for the quench dynamics and the adiabatic passage, respectively. The truncation threshold is $10^{-6}$ per time step. We investigate the finite-size effects by calculating the spin dynamics in 1D systems at several sizes $N$. As shown in Extended Data Fig. 8, the discrepancy between the curves for $N=40$ and $N=70$ is on the order of $10^{-2}$, indicating the volume convergence of our calculations.

Temperature $T$ and entropy $\mathcal{S}/N$ of the state are estimated by comparing the experimental results with finite-temperature quantum Monte Carlo methods (QMC) Bauer:2011s . All the simulations are carried out in the grand-canonical ensemble on a 70-site open chain. In such a large system, the finite-size effect can be eliminated. After $10^{5}$ times of thermalization, we continue to perform the QMC sweeps for at least $5\times 10^{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 Extended Data Fig. 9(a). In Extended Data Fig. 9(b), we present the next-nearest-neighbor correlation as a function of the temperature. At a temperature of 0.53(7)$J$, $C^{\prime}(d=2)$ reaches 0.57(1) in this thermal state. The entropy per particle is obtained using another QMC method, which is called quantum Wang-Landau algorithm Wang:2001s . The cutoff of series expansion is set to be $10^{3}$. As shown in Extended Data Fig. 9(c), the entropy per particle deduced from the temperature of the spin chain is 0.42(5)$k_{B}$.

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 Extended Data Fig. 9(d).

From the $t$-DMRG simulations, we plot the spin correlations of the antiferromagnet up to 15 sites. As shown in Extended Data Fig. 9(e), the correlation function exhibits a power-law decay as the distance $d$ increases, which is consistent with the property of the Tomanaga-Luttinger Liquid Haldane:1981 ; Giamarchi:2004 ; Yang:2017s .

The spin fluctuations in terms of $\delta M_{z}$ can reflect non-local spin correlations of many-body states. However, the detection of $\delta M_{z}$ is obstructed by the detection noise, for instance the shot noise of our absorption imaging. Here, we develop a method to resolve the spin fluctuations on the basis of the two-image measurements. The total density $n_{A}+n_{B}$ corresponds to a unit-filling Mott insulator with vanishing density fluctuations. Since the staggered magnetization is $M_{z}=n_{A}/(n_{A}+n_{B})-1/2$, where the spin fluctuations only exist on $n_{A}$. In this respect, we can extract $\delta M_{z}$ directly from the fluctuations of $n_{A}$. For reference, we consider the background fluctuations due to the photon shot-noise by analyzing the total density $n_{A}+n_{B}$. We resolve the spin fluctuations $\delta M_{z}$ by subtracting the shot-noise background, as $\delta M_{z}\propto\sqrt{\sigma(n_{A})^{2}-\sigma_{bg}^{2}}$.

As Extended Data Fig. 10 (a) shows, the density histogram of $n_{A}$ throughout the adiabatic sweep are obtained by counting the densities on 2.37 copies 70-site chains (one pixel equals to 2.37 sites). Based on 100 experimental repetitions, we acquire 1500 groups of measurement for each histogram. Here, $\sigma(n_{A})$ represents the Gaussian waist of the histogram of $n_{A}$ (the Gaussian function is used for fitting), as shown in Extended Data Fig. 10 (b). Taking into account the imaging resolution and the pixel size, the fluctuation of staggered magnetization $\delta M_{z}$ can be extracted.