Observing emergent hydrodynamics in a long-range quantum magnet

The experiments are performed with a string of ⁴⁰Ca⁺ ions confined in a linear radiofrequency trap. For confining long ion strings, the trapping potential is made strongly anisotropic, with oscillation frequencies of 2.93 MHz and 2.898 MHz in the radial plane, and an axial trapping frequency of 126.3 kHz (126.6 kHz) for trapping 25 ions (51 ions). At the above stated confining frequencies, the 25-ion (51-ion) string has a length of 173 $\mu$m (247 $\mu$m). (Pseudo-)spins are encoded in superpositions of 4S${}_{1/2},m=1/2$ and 3D${}_{5/2},m=5/2$ electronic states and coherently manipulated by a frequency-stable laser at 729 nm coupling to the quadrupole transition. A low-noise magnetic field of 4.17 Gauss produced by an assembly of SmCo magnets is used to lift the degeneracy of the Zeeman states.

At the beginning of each experiment, all transverse collective modes of the ion string are cooled close to the ground state by sideband cooling on the S${}_{1/2}\leftrightarrow\text{D}{}_{5/2}$ quadrupole transition. The axial modes are polarization-gradient cooled to sub-Doppler temperatures and the ions are prepared in the S${}_{1/2},m=+1/2$ state by frequency-resolved optical pumping [36].

At the end of each experimental run, a quantum state measurement projecting onto the (pseudo-)spin basis states is carried out by spatially resolved measurements of the ions’ fluorescence emitted on the S_1/2 to P_1/2 transition. The state assignment error probability is below $10^{-3}$ per ion.

Preparation of the (pseudo-)spins in the desired initial product states is achieved by using a spin-echo sequence that combines resonant laser pulses coupling to all spins with nearly the same strength and off-resonant pulses addressed to individual spins that are sandwiched between the global pulses. For a 25-ion crystal, the state preparation starts with a set of composite $\pi/2$ pulses that prepares all calcium ions in an equal superposition of the 4S${}_{1/2},m=1/2$ and 3D${}_{5/2},m=5/2$ electronic states. A tightly focused laser beam, with beam radius of about 2 $\mu$m, is steered over the ion string (ion spacing $>$ 4 $\mu$m) with an acousto-optic deflector. Pre-calibrated addressing pulses are placed in a specific order within the two halves of the spin-echo sequence such that errors resulting from cross-talk while addressing neighbouring ions are minimized. The pulse sequence is completed by rotating the unaddressed (addressed) ions to the S (D) state by the second global $\pi/2$ pulse. For a 51-ion crystal, we perform the above described sequence but on 4S${}_{1/2},m=1/2$ and 3D${}_{5/2},m=-3/2$ electronic states before restoring the population to the desired qubit states. Due to the restriction of the polarization angle and beam geometry, we only achieve optimum off-resonant coupling, thus short interaction times, with the addressing beam while working on 4S${}_{1/2},m=1/2$ and 3D${}_{5/2},m=-3/2$ states.

Composite laser pulses: Ions are collectively excited with an elliptically shaped laser beam with a beam diameter of about 470 $\mu$m, which is shone over the ion string of length 247 $\mu$m from the transverse direction of the ion string. This leads to a non-uniform illumination, resulting in spatially dependent rotation angles of the ions when performing the global single qubit rotations. In our experiment, we reduce the impact of these unequal rotations of the spins by employing Composite pulse sequences [37]. The Composite pulse sequences for $\theta=\pi/2$ and $\theta=\pi$ angle rotations are given as

	$\displaystyle U^{x}_{\theta=\pi/2}=e^{-i\frac{\pi}{4}\sigma_{y}}e^{-i\frac{\pi}{4}\sigma_{x}},$		(5)
	$\displaystyle U^{x}_{\theta=\pi}=e^{-i\frac{\pi}{4}\sigma_{x}}e^{-i\frac{\pi}{2}\sigma_{y}}e^{-i\frac{\pi}{4}\sigma_{x}}.$		(6)

Off-resonant rotations with a tightly focused beam: The tightly focused off-resonant laser beam is used in conjunction with the global laser beam. Pulse length and laser intensity of the addressing beam are chosen such that the specific ions are prepared in the desired quantum states with a single ion fidelity of about $99\%$. Addressing of 51(25) ions, required about 0.5(0.25) ms of interaction time, which is much shorter than the coherence time ($\sim$ 50 ms) of our system, hence the initial state is well protected against any decoherence effects. Slow spatial drift of the ion string along the weak confining axis, which impacts the quality of addressing operations, is corrected on-the-fly by routinely measuring the ion positions and consecutively applying compensating potentials to the trap electrodes. Laser beam intensities of both global and addressing laser beams are stabilized via a sample and hold intensity stabilization technique. Furthermore, the addressing pulse lengths are routinely calibrated and corrected against any unwanted slow variations.

Product state fidelities: Fig. 5 shows initial state preparation fidelities of the 25-ion and 51-ion strings. Each spin is prepared in the desired up/down basis states with about 98.8$\%(99.3\%)$ fidelity, which corresponds to a fidelity of about $73$%$(70\%)$ of the desired product state for a 25(51) ion string.

Spin-spin interactions are engineered with a two-tone laser beam inducing Mølmer-Sørensen type interactions; i.e. the laser off-resonantly couples to both the lower and upper vibrational sidebands of all $2N$ collective motional modes describing the ion motion in the plane that is normal to the direction of an $N$-ion string. The tunable long-range coupling is provided by a specific choice of laser parameters that generates a long-range Ising-type spin-spin interaction with a coupling strength approximately described by a power-law [38, 39]. The laser beam is elliptically shaped and elongated along the axis of the ion string to provide near-uniform illumination over the ion string. Unwanted light shifts, which occur due to the presence of other electronic transitions, are compensated with a third laser beam that is mixed with the two-tone laser beam [30]. The laser beam propagates from the transverse direction of the axis of the ion string and imparts a near flat waveform over the ion string, which otherwise can have adverse effects in the laser-ion interaction.

The engineered spin-spin interaction, for $B\gg J$ gives rise to a flip-flop type interaction described by

$$\hat{H}=\sum_{i<j}\frac{J}{|i-j|^{\alpha}}(\hat{\sigma}^{+}_{i}\hat{\sigma}^{-}_{j}+\hat{\sigma}^{-}_{i}\hat{\sigma}^{+}_{j}),$$

(7)

where exponent $\alpha$ and strength of the spin-spin interaction $J$ can be tuned by changing the detuning and intensity of the entangling laser beam. The experimentally realized Hamiltonian slightly differs from the power law expression given in equation 7. The exact spin-spin matrix connecting all spins is given by

$$J_{ij}=\frac{\Omega_{i}\Omega_{j}}{2}\sum_{m=1}^{2N}\frac{\eta_{im}\eta_{jm}}{\Delta_{m}},$$

(8)

where $\eta_{jm}$ is the Lamb-Dicke parameter of the $j^{th}$ ion and the $m^{th}$ motional mode. $\Delta_{m}$ is the detuning of the two-tone laser field from the motional mode of interest and $\Omega_{i}$ is the Rabi frequency of the $i^{th}$ ion. The detuning of the laser beams from the vibrational sideband ranges from a few tens of kHz for the center-of-mass (COM) mode to about 1.5 MHz for the lowest-frequency mode, resulting in effective spin-spin interactions that decay with distance between the spins. The structure of the coupling matrix can be experimentally controlled by varying the laser detuning from the COM mode and the anisotropy of the trapping potential. A comparison between the approximated power law and experimentally realized interaction Hamiltonian is shown in Fig. 6 for all three $\alpha$ values that are used in the main manuscript. The figure displays a small mismatch between the power-law interaction of eq. 7 and the effective interaction realized in the trapped ion systems [39], which is inherently present in any the laser-ion interaction. Additionally, due to technical reasons, a non-uniform illumination of the laser beam over the ion chain alters the spin-spin interaction and slightly weakens the interactions of the spins at the ends of the chain. The numerical simulations presented in the main manuscript are performed for the experimentally realized spin-spin interaction which includes the non-uniform laser intensity over the ion string.

The engineered entangling interaction is experimentally tested by time evolving an initial state where the middle ion (remaining ions) is initialized in a spin-up (spin-down) state. In Fig. 7(a) we show dynamics of spin excitation for $\alpha=1.1$. We compare the experimental results with numerical simulations that we carry out in a subspace where the magnetization is conserved, see Fig. 7(b) and find a good agreement at early times whereas some discrepancy is visible at late times. This discrepancy is attributed to decoherence effects that arise from non-perfect initial state preparation and non-conservation of magnetisation, discussed in the following section.

While experimentally realizing the spin-spin Hamiltonian, the system undergoes some unwanted decoherence effects that can alter the dynamics. In this section, we describe the two major decoherence channels that are important for the initial state preparation and dynamics with long ion strings. The first one is dephasing noise, which arises from laser frequency/phase noise and fluctuating magnetic fields. Jointly, they reduce the coherence time of our spins. Most importantly, the superposition, which is an essential step for individual control of the spins for preparing arbitrary product states, decays due to the dephasing noise. The laser frequency noise is reduced by locking the laser to an ultra-stable reference cavity. The trap setup is placed inside a mu-metal shield to suppress the magnetic field noise. The overall coherence time of the system is about 50 ms. Furthermore, we reduce the impact of the dephasing noise by applying spin-echo sequences. While these effects play an important role for the state preparation, the dynamics is not affected, as we work in the decoherence-free-subspace in which the spins are protected against the dephasing noise.

The second type of decoherence that spins experience can be described by amplitude damping and depolarization channels. Due to the finite lifetime of the metastable state D_5/2 ($\tau\sim$1 s), in which our spin-up state is encoded, decays to the S_1/2 electronic state. In addition to this decay channel, the spins experience unwanted excitation/decay due to the bichromatic laser field coupling to the motional sidebands. These processes violate the conservation of magnetization that would be expected from the flip-flop interaction of eq. (7). We experimentally measure the spin-flip rates for the experimental conditions that are used in the main manuscript. In Fig. 8 we show magnetization variation as a function of interaction time for a 51-ion chain. Here, we model two spin flip rates: spontaneous decay of the excited state $\Gamma$ and the laser induced spin flips rate $\gamma_{\text{flip}}$, similar to what has been described in reference [10]. The model assumes that the spin-flip rate due to incoherent laser-ion coupling does not depend on the spin orientations. The flip rates that we extract by fitting the experimental data are $\gamma_{\text{flip}}=0.78$/s and $\Gamma=0.91$/s. These two competing processes in the end give rise to a dominating decay of the spins that are encoded as spin up over those are encoded as spin down. For the measurements that are presented in this manuscript, we remove the biased decay of the excitation to leading order by preparing a second set of initial states, where the central spin is prepared in a $\ket{\downarrow}$-state and average over the two distinct sets of initial states.

The magnetization decay model, which we assumed for fitting of the experimental data, considers that all the ions have the same spin flip rate. However, we experimentally observed that this rate depends slightly on the position of the ions in the string. The exact mechanism underpinning this effect is still under investigation. We speculate that the spatially dependent laser-ion interactions, coupling to the low-frequency motional modes of the ion chain, could induce such adverse effects. Most importantly, cross-frequency components that arise while producing the entangling laser beam could excite motional modes which in an ideal scenario should be suppressed. Fig. 9 displays one of these artifacts which is observed while time evolving a spin state with all spins prepared in spin-down orientation. It may be noted that the magnetization varies the most for the middle ions in comparison to the outer ions. Ideally, the Hamiltonian does not couple this state to any other states hence we expect the state to be unaltered. However, due to the incoherent spin flips the total magnetization of the time evolved state changes with the interaction time. We further investigate the impact of spatially dependent incoherent processes in our current study of emergent hydrodynamics. In Fig. 10, we show the correlation function that we experimentally measure at $t=40$ ms. The measurement sequence is identical to that has been described in Fig. 1 of the main manuscript, except that we exchange the encoding of ion number 26${}^{\text{th}}$ with the first ion in the 51-ion string. Note that the current measurements reveal no signature of unwanted variation of the magnetization which signifies the robustness of the measurement protocol against magnetization errors.

For the present measurements, the experimental error bars for an individual ion are estimated from the quantum projection noise model [40]. In our experiment, we repeat the measurements $N_{m}$ times for each initial state. Each initial state is chosen from a uniform random distribution function that samples the individual spin state with an equal probability of spin-ups and spin-downs, while restricting the total magnetization of the whole system to be $+1$. The measurements are performed $N_{u}$ times. The mean standard deviation of measurements for $i^{\text{th}}$ ion, while preparing the initial state $\ket{\Phi_{u}}=\bigotimes_{i=1}^{N}\ket{\phi_{u}^{i}}$, is given by

$$\sigma_{u}^{i}=\sqrt{p_{u}^{i}(1-p_{u}^{i})/N_{m}},$$

(9)

where $p_{u}^{i}$ is the probability of finding the $i^{th}$ ion, prepared in $\ket{\phi_{u}^{i}}$ state, in the measurement basis. We further propagate this formula for the auto-correlation function described in Eq. (3) of the main text and estimate the mean standard deviation to be

$$\sigma_{corr}=2\frac{\sqrt{\sum_{u=1}^{N_{u}}\sigma_{u}^{2}}}{N_{u}}.$$

(10)

.

To extract transport coefficients from the experimental data, we perform a least-square fitting with the spatio-temporal function

$$C_{j}(t)=(D_{\alpha}t)^{-\beta}F_{\alpha}\left(\frac{|j|}{(D_{\alpha}t)^{\beta}}\right),$$

(11)

where $D_{\alpha}$ represents the transport coefficient and the exponent $\beta$ determines the rate at which auto-correlations decay with time. $F_{\alpha}$ are the stable symmetric distributions defined in Eq. 19 below.

Fitting of scaling functions: In Fig. 2 of the main text we fitted a normalized Lorentzian and Gaussian to the data points which are finite within error bars. Hence, the full profile was fitted for $\alpha=1.5$, while the fit was restricted to the central 27 sites for $\alpha=1.1$.

Extraction of diffusion coefficients: In Fig. 4 of the main text, we used the expected values for the exponent $\beta=1/(2\alpha-1)$ to rescale the data. We then fitted the Ansatz in Eq. 11 to the rescaled data to extract the diffusion coefficient $D_{\alpha}$. The error bars are deduced by taking into account the uncertainty in extracting the $\alpha$ parameters.

Spatio-temporal evolution of correlations: The unscaled spatio-temporal evolution of correlations for all three $\alpha$ values are shown in Fig. 11.

In Fig. 12 we compare the experimentally measured auto-correlation function with an exact simulation for the 25 ions chain, finding excellent agreement.

Transport in the long-range XY model (7) is described by Lévy flights [29]. The Hamiltonian connects states with an up spin at site $i$ and a down spin at site $j$, $\ket{\uparrow_{i}\downarrow_{j}}$, with the opposite configuration $\ket{\downarrow_{i}\uparrow_{j}}$. As these states have equal energy, we can use Fermi’s golden rule to estimate the transition probability for a spin exchange as

$$W_{i\rightarrow j}\equiv|\braket{\uparrow_{i}\downarrow_{j}}{\hat{H}}{\downarrow_{i}\uparrow_{j}}|^{2}=\frac{\lambda}{|i-j|^{2\alpha}},$$

(12)

with $\lambda\propto J^{2}$. Moreover, $W_{i\rightarrow j}=W_{j\rightarrow i}\equiv W_{ij}$. Hence, we can construct a classical Master equation for the spin excitation probability density $f_{i}\equiv(\braket{\hat{\sigma}^{z}_{i}}+1)/2$,

$$\partial_{t}f_{i}=\sum_{j}W_{ij}(f_{j}-f_{i}),$$

(13)

where the first and second term on the right hand side are the “gain” and “loss” terms, respectively.

We can solve this equation by Fourier transforming to momentum $k$ and taking both a continuum and infinite system size limit while retaining a short-distance cutoff of unity, leading to

$$\partial_{t}f_{k}=W_{k}f_{k},$$

(14)

with the transition rates

$$W_{k}=\left(-c_{\alpha}|k|^{2\alpha-1}+\frac{k^{2}}{3-2\alpha}\right)\lambda,$$

(15)

where $c_{\alpha}=-2\Gamma(1-2\alpha)\sin(\alpha\pi)$ and $\Gamma$ is the Gamma function. Solving Eq. (14) and Fourier transforming back for an initially localized excitation $f_{j}=\delta_{j,0}$, we get

$$f_{j}(t)=\int\frac{dk}{2\pi}e^{ikj+W_{k}t}.$$

(16)

The exact form of the solution crucially depends on the value of $\alpha$. In the following, we discuss the three cases shown in the dynamical phase diagram in Fig. 1, main text.

Diffusion for $\alpha\geq 1.5$: In this case, the second term in Eq. (15) dominates. Inserting $W_{k}$ into Eq. (16), we find

$$f_{j}^{\alpha>1.5}(t)=(D_{\alpha}t)^{-1/2}G\left(\frac{|j|}{(D_{\alpha}t)^{1/2}}\right),$$

(17)

with a diffusive Gaussian profile $G(y)=\exp(-y^{2}/4)/\sqrt{4\pi}$ and diffusion coefficient $D_{\alpha}=\lambda/(2\alpha-3)$. The long-range character of the transition rates only shows up in subleading corrections, leading to algebraic tails of the distribution function [41, 29].

Superdiffusion for $0.5<\alpha<1.5$: In this case, we get

$$f_{j}^{0.5<\alpha<1.5}(t)=(D_{\alpha}t)^{-1/(2\alpha-1)}F_{\alpha}\left(\frac{|j|}{(D_{\alpha}t)^{1/(2\alpha-1)}}\right),$$

(18)

with $D_{\alpha}=\lambda c_{\alpha}$. The distribution function $F_{\alpha}$ is given by the stable symmetric distributions defined by

$$F_{\alpha}(y)=\int\frac{dk}{2\pi}e^{iyk-|k|^{2\alpha-1}},$$

(19)

which only has a known elementary solution for $\alpha=3/2$ (Gaussian) and $\alpha=1$ (Lorentzian).

Mean-field regime for $\alpha<0.5$: For $\alpha<0.5$, the system enters the mean-field (or unstable transport) regime. Here, a consistent hydrodynamic description is no longer possible in the thermodynamic limit, as revealed by the diverging transition rates in Eq. (15) for $k\rightarrow 0$. Instead, the system forms a zero-dimensional “quantum dot” rather than a one-dimensional chain, and relaxes to thermal equilibrium instantaneously [29, 42], due to the non-local nature of the interactions. This behavior can be best understood from the limit of all-to-all interactions $\alpha=0$, allowing to rewrite the Hamiltonian in terms of a single global spin operator $\vec{\hat{\sigma}}_{\rm global}=\sum_{i=1}^{L}\vec{\hat{\sigma}}_{i}$. In the thermodynamic limit, the diverging spin eigenvalue associated with this operator gives rise to a vanishing relaxation time.

Here, we calculate the short time dynamics of

$$\operatorname{Tr}\left[\hat{\sigma}^{z}_{i}(t)\hat{\sigma}^{z}_{j}(0)\right]$$

(20)

for time evolution under the Hamiltonian

$$\hat{H}=\sum_{i<j}J_{ij}(\hat{\sigma}^{+}_{i}\hat{\sigma}^{-}_{j}+h.c.).$$

(21)

We expand $\hat{\sigma}^{z}_{i}(t)$ in the Heisenberg picture,

	$\displaystyle\hat{\sigma}^{z}_{i}(t)=e^{i\hat{H}t}\hat{\sigma}^{z}_{i}e^{-i\hat{H}t}$
	$\displaystyle=\hat{\sigma}_{i}^{z}+iJt\left[\frac{\hat{H}}{J},\hat{\sigma}_{i}^{z}\right]+\frac{(iJt)^{2}}{2!}\left[\frac{\hat{H}}{J},\left[\frac{\hat{H}}{J},\hat{\sigma}_{i}^{z}\right]\right]+\mathcal{O}((Jt)^{3}).$

When evaluating the commutators and taking the trace, the first order contribution vanishes and we get

	$\displaystyle\operatorname{Tr}[\hat{\sigma}^{z}_{i}(t)\hat{\sigma}^{z}_{j}(0)]=$	$\displaystyle t^{2}\sum_{k\neq i}J_{ik}^{2}\left(\operatorname{Tr}[\hat{\sigma}_{k}^{z}\hat{\sigma}_{j}^{z}]-\operatorname{Tr}[\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}]\right)$
	$\displaystyle+$	$\displaystyle\operatorname{Tr}[\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}].$

Furthermore, $\operatorname{Tr}[\hat{\sigma}_{k}^{z}\hat{\sigma}_{m}^{z}]=\delta_{km}$ and hence

$$\operatorname{Tr}[\hat{\sigma}^{z}_{i}(t)\hat{\sigma}^{z}_{j}(0)]=\bigg{\{}\begin{array}[]{lr}1-t^{2}\sum_{i\neq k}J_{ik}^{2}&\text{for }i=j\\ t^{2}J_{ij}^{2}&\text{for }i\neq j\\ \end{array}.$$

(22)

In Fig. 3 of the main text we used the explicit $J_{ij}$ matrix determined from the laser parameters in the experiment to obtain the short time expansion.

The infinite temperature correlation function in Eq. (3) of the main text can be written as

	$\displaystyle C_{j}(t)$	$\displaystyle=\braket{\hat{\sigma}^{z}_{j}(t)\hat{\sigma}^{z}_{0}(0)}_{T=\infty}$		(23)
		$\displaystyle=\frac{1}{Z}\mathrm{Tr}\left[e^{-\hat{H}/T}\hat{\sigma}^{z}_{j}(t)\hat{\sigma}^{z}_{0}(0)\right]\bigg{|}_{T=\infty}$		(24)
		$\displaystyle=\frac{1}{Z}\mathrm{Tr}\left[\hat{\sigma}^{z}_{j}(t)\hat{\sigma}^{z}_{0}(0)\right],$		(25)

with the partition sum $Z=\mathrm{Tr}[e^{-\hat{H}/T}]|_{T=\infty}=2^{L}$ for systems of size $L$. We expand the trace in the $\hat{\sigma}^{z}$ basis to obtain

	$\displaystyle C_{j}(t)$	$\displaystyle=\frac{1}{Z}\sum_{\sigma_{-\frac{L}{2}},\dots,\sigma_{\frac{L}{2}}}\braket{\sigma_{-\frac{L}{2}}\dots\sigma_{\frac{L}{2}}}{\hat{\sigma}^{z}_{j}(t)\hat{\sigma}^{z}_{0}(0)}{\sigma_{-\frac{L}{2}}\dots\sigma_{\frac{L}{2}}}$		(26)
		$\displaystyle=\frac{1}{Z}\sum_{\sigma_{-\frac{L}{2}},\cdots,\sigma_{\frac{L}{2}}}\sigma_{0}\braket{\sigma_{-\frac{L}{2}}\cdots\sigma_{\frac{L}{2}}}{\hat{\sigma}^{z}_{j}(t)}{\sigma_{-\frac{L}{2}}\dots\sigma_{\frac{L}{2}}},$		(27)

where $\sigma_{j}=\pm 1$. We have therefore reduced the evaluation of a two-time correlation function to the weighted sum of single-time functions using product initial states, which can be realized in the experiment. In order to evaluate Eq. (27) exactly, all $2^{L}$ product states in the z-basis would have to be prepared, time evolved and measured to then obtain $C_{j}(t)$ from the weighted sum. In practice, however, it suffices to sample a small number $M\approx 10-100$ [29] of product states to evaluate Eq. (27) while replacing $Z\rightarrow M$. In the main text, we considered the largest sector with magnetization $\sim 0$. In Fig. 13 we compare the sampling of the largest sector to the exact results of the thermal expectation value. The full trace is evaluated using the typicality approach [43], which involves evolving Haar-random states $\ket{\psi_{r}}$ in the entire Hilbert space and then averaging over $5-10$ such states.

To further improve the convergence of the sampling of Eq. (27) we prepare the conjugate configuration for each randomly drawn product state. In the conjugate state, all spins are flipped except for the central one. This choice of initial states reduces statistical noise by reproducing the exact property $C_{j}(0)=\delta_{j,0}$ of the full trace for each pair of initial states.