Adiabatic state preparation of stripe phases with strongly magnetic atoms

TL;DR

This paper presents a protocol for preparing stripe phases in two-dimensional spin models using strongly magnetic atoms in optical lattices, validated by simulations that estimate realistic preparation times.

Contribution

It introduces a new adiabatic state preparation protocol for stripe phases with validation through advanced simulation methods.

Findings

Successful simulation of dynamics for up to 64 sites.

Estimated high-fidelity preparation times for stripe phases.

Validation of the protocol with exact small-system results.

Abstract

We propose a protocol for realizing the stripe phase in two spin models on a two-dimensional square lattice, which can be implemented with strongly magnetic atoms (Cr, Dy, Er, etc.) in optical lattices by encoding spin states into Zeeman sublevels of the ground state manifold. The protocol is tested with cluster-mean-field time-dependent variational ans\"atze, validated by comparison with exact results for small systems, which enable us to simulate the dynamics of systems with up to 64 sites during the state-preparation protocol. This allows, in particular, to estimate the time required for preparation of the stripe phase with high fidelity under real experimental conditions.

Equations

$$H=\sum _{i<j}^N{C}_{ij}{\sigma }_z^{(i)}{\sigma }_z^{(j)}+\Omega \sum _i^N{\sigma }_x^{(i)}-\Delta \sum _i^N{\sigma }_z^{(i)},$$

$$H$$

$$+$$

$$M_z=\frac{1}{N}\sum _i⟨S_z^{(i)}{S}_z^{(i+1_z)}⟩,M_y=\frac{1}{N}\sum _i⟨S_z^{(i)}{S}_z^{(i+1_y)}⟩,$$

$$|\Phi ⟩=\prod _{i=1}^N({\alpha }_i(t){|\uparrow ⟩}_i+{\beta }_i(t){|\downarrow ⟩}_i),$$

$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{\alpha }}_i^{\ast }}\right)=\frac{\partial L}{\partial {\alpha }_i^{\ast }},$$

$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{\beta }}_i^{\ast }}\right)=\frac{\partial L}{\partial {\beta }_i^{\ast }},$$

$$L=\frac{i}{2}\left(⟨\Phi |\dot{\Phi }⟩-⟨\dot{\Phi }|\Phi ⟩\right)-⟨\Phi |H|\Phi ⟩.$$

$$-i{\dot{\alpha }}_i=\left[-\Delta (t)+\sum _j{C}_{ij}(|{\beta }_j{|}^2-|{\alpha }_j{|}^2)\right]{\alpha }_i-\Omega (t){\beta }_i,$$

$$-i{\dot{\beta }}_i=\left[\Delta (t)+\sum _j{C}_{ij}(|{\alpha }_j{|}^2-|{\beta }_j{|}^2)\right]{\beta }_i-\Omega (t){\alpha }_i.$$

$${|\psi ⟩}_{\mu }=\sum _{l,l^{\prime }}{\alpha }_{l{l}^{\prime }}^{(\mu )}{|l⟩}_{1\mu }{|l^{\prime }⟩}_{2\mu },$$

$$|\Phi ⟩=\prod _{\mu =1}^{N/2}{|\psi ⟩}_{\mu }.$$

$$H_{\mathrm{SB}}=\sum _{\mu }ϵ(t)\left({|1⟩}_{1\mu }{⟨1|}_{1\mu }-{|1⟩}_{2\mu }{⟨1|}_{2\mu }\right)$$

$$V_{ij}=\frac{c_d}{r_{ij}^3}\left({\mathbf{J}}^{(i)}.{\mathbf{J}}^{(j)}-3({\mathbf{J}}^{(i)}.{\hat{\mathbf{r}}}_{ij})({\mathbf{J}}^{(j)}.{\hat{\mathbf{r}}}_{ij})\right)$$

$$H_{\mathrm{int}}=H_{q=0}+H_{q=1}+H_{q=2},$$

$$H_{q=0}=\sum _{i<j}{J}_{ij}\left(J_z^{(i)}{J}_z^{(j)}-\frac{1}{4}(J_{+}^{(i)}{J}_{-}^{(j)}+h.c.)\right)$$

$$H_{q=1}=-\frac{3}{2}\sum _{i\ne j}\frac{c_d}{r_{ij}^3}\sin {\theta }_{ij}\cos {\theta }_{ij}(J_{+}^{(i)}{J}_z^{(j)}{e}^{-i{\varphi }_{ij}}+h.c.),$$

$$H_{q=2}=-\frac{3}{4}\sum _{i<j}\frac{c_d}{r_{ij}^3}{\sin }^2{\theta }_{ij}(J_{+}^{(i)}{J}_{+}^{(j)}{e}^{-2i{\varphi }_{ij}}+h.c.),$$

$$h$$

$$h_J=\sum _{m=-J}^J(\alpha m+{ϵ}_m)|m⟩⟨m|,$$

$$|\Phi ⟩=\prod _{i=1}^N\sum _l{\alpha }_{l,i}{|l⟩}_i,$$

$$-i{\dot{\alpha }}_{-1}^i$$

$$-i{\dot{\alpha }}_0^i$$

$$-i{\dot{\alpha }}_1^i$$

$$\frac{d}{dt}{\alpha }_{\mu }=C{\alpha }_{\mu },$$

$$C=\left[\begin{array}[]{ccccccccc}C_{0}&J_{1}A^{*}&0&J_{1}B^{*}&0&0&0&0&0\\ J_{1}A&C_{1}&J_{2}A^{*}&\frac{J_{1}^{2}}{4}&J_{1}B^{*}&0&0&0&0\\ 0&J_{2}A&C_{2}&0&\frac{J_{1}J_{2}}{4}&J_{1}B^{*}&0&0&0\\ J_{1}B&\frac{J_{1}^{2}}{4}&0&C_{3}&J_{1}A^{*}&0&J_{2}B^{*}&0&0\\ 0&J_{1}B&\frac{J_{1}J_{2}}{4}&J_{1}A&C_{4}&J_{2}A^{*}&\frac{J_{1}J_{2}}{4}&J_{2}B^{*}&0\\ 0&0&J_{1}B&0&J_{2}A&C_{5}&0&\frac{J_{2}^{2}}{4}&J_{2}B^{*}\\ 0&0&0&J_{2}B&\frac{J_{1}J_{2}}{4}&0&C_{6}&J_{1}A^{*}&0\\ 0&0&0&0&J_{2}B&\frac{J_{2}^{2}}{4}&J_{1}A&C_{7}&J_{2}A^{*}\\ 0&0&0&0&0&J_{2}B&0&J_{2}A&C_{8}\end{array}\right],$$

$$A$$

$$B$$

$$C_0$$

$$C_1$$

$$C_2$$

$$C_3$$

$$C_4$$

$$C_5$$

$$C_6$$

$$C_7$$

Full text

1 Introduction

Experiments with ultra-cold quantum gases of atoms and molecules have recently witnessed tremendous progresses in realising and probing many-body physics with long-range interacting spin systems [1]. Several works have reported the observation of coherent dynamics with both magnetic and Rydberg atoms [2, 3] and polar molecules [1], ranging from the time evolution of spin models after a quantum quench [4, 5, 6, 7], to the realisation of Hubbard and extended Bose-Hubbard models featuring both local and non-local interactions [8, 9, 10, 11]. These experimental studies are paving the way to the investigation of quantum magnetism in atomic and molecular quantum systems, in particular, about the stability of ordered magnetic phases in frustrated systems [1, 2].

One example of such phases is the stripe phase - supporting ferromagnetic order along one direction, and antiferromagnetic along the perpendicular one. Stripe order has been widely discussed in quantum magnetism, both in the context of spin (Heisenberg and Ising type) [12, 13] and fermionic (Hubbard and t-J) models (with possible implications on the finite-doping phases of high-temperature superconductors) [14]. Moreover, stripe ordering represents a natural counterpart to anti-ferromagnets in the presence of angular dependent interactions [12]: these interactions do naturally occur in dipolar gases of magnetic atoms such as Cr, Dy, and Er, all of which are by now being observed in quantum degenerate regimes [15, 16, 17, 18, 19]. The main virtues of such atomic species is that they can combine dipolar interactions with close-to-ideal initial state preparation (such as Mott insulators with unit filling) and long coherence times (of the order of 1 $s$) which compare well with the interaction timescales (typically of the order of 10 to 100 Hz in units of $\hbar$) [1, 2]. These two features make magnetic atoms an ideal platform for the realisation and investigation of stripe phases in atomic many-body systems, an approach very much complementary to recent experimental findings [20] in the context of spin orbit coupled Bose gases [21, 22].

One of the key questions toward the realisation of such magnetic phases is the identification of clear-cut state preparation protocols which can warrant experimental observability of specific signatures (i.e. correlation functions) under realistic experimental conditions. In this work, we address the preparation of the stripe phase in spin-1 Heisenberg-type and spin-$1/2$ Ising-type models on the square lattice, with antiferromagnetic interactions between neighbours in one direction, and ferromagnetic along the other, which can be realised with magnetic atoms trapped in optical lattices. For both models, the proposed adiabatic protocols for preparing the stripe phase are tested using a time-dependent variational principle, directly inspired by (cluster) mean field theory. This method agrees well with the exact treatment of small systems, and allows us to tackle the full-time dependent dynamics for lattices containing of the order of 100 sites - which is comparable to the largest, defect-free Mott insulator states available in current experiments [23]. Our results strongly indicate that the realisation of the stripe phases is within the reach of current experiments with strongly magnetic atoms such as Dy or Er, with strong signatures of the stripe formation already visible for moderate speeds of changes in the system’s parameters during the preparation protocols. Beyond the magnetic atom implementation discussed here, our findings are also applicable to other physical systems, such as polar molecules and Rydberg atoms in optical lattices [2], and arrays of superconducting circuits [24, 25], which are described by very similar lattice spin Hamiltonians.

The structure of the paper goes as follows. In section 2, we discuss the two dipolar spin models of interest. In section 3, we introduce the time-dependent variational principle, and discuss in detail the state preparation protocols for stripe phases in both spin-$1/2$ and spin-1 models. Section 4 contains two implementation schemes using magnetic atoms in optical lattices, and provides a glimpse of the corresponding experimental timescales for both coherent and incoherent (dissipation) dynamics. Finally, we draw our conclusions in section 5.

2 Model Hamiltonians

In this section we introduce two spin models of interacting dipoles whose phase diagrams include stripe phases. The setup we have in mind is shown in figure 1(a). It contains magnetic atoms at fixed positions on a square lattice, interacting via dipole-dipole interactions (DDI) and with an external magnetic field along the $z$-axis. This is the common starting point for the realisation of both spin-$1/2$ and spin-1 lattice models. In the following we discuss the corresponding Hamiltonians and the phase diagrams. For the implementation of the models, we refer the reader to section 4.

2.1 Spin-$1/2$ model

The first model we consider is obtained by assuming the two involved atomic states as states $\left|{\uparrow}\right\rangle,\left|{\downarrow}\right\rangle$ of pseudo-spin $1/2$ particles. The corresponding Hamiltonian reads

[TABLE]

where $N$ is the number of atoms, $\sigma_{z}^{(i)}=\left|{\uparrow}\right\rangle_{i}\left\langle{\uparrow}\right|_{i}-\left|{\downarrow}\right\rangle_{i}\left\langle{\downarrow}\right|_{i}$ and $\sigma_{x}^{(i)}=\left|{\uparrow}\right\rangle_{i}\left\langle{\downarrow}\right|_{i}+\left|{\downarrow}\right\rangle_{i}\left\langle{\uparrow}\right|_{i}$ are the Pauli matrices acting at site $i$. The first term in this Hamiltonian describes the Ising-type DDI between two spins on sites $i$ and $j$ with $C_{ij}=c_{d}(1-3\cos^{2}\theta_{ij})/r_{ij}^{3}$, where the interaction constant $c_{d}$ depends on the specific implementation of the spin models, $r_{ij}=|\textbf{r}_{i}-\textbf{r}_{j}|$ is the distance (in units of the lattice spacing $a$) between the sites, and $\theta_{ij}$ the corresponding angle with respect to the quantisation axis $z$ [c.f. figure 1(a)]. In order to have the stripe phase as the ground state, we place the quantisation axis in the plane of the atoms such that DDI is attractive in one direction ($\theta_{ij}=0$) and repulsive in the other direction ($\theta_{ij}=\pi/2$) 111If the quantisation axis is perpendicular to the lattice, the angle $\theta_{ij}$ is fixed to $\pi/2$ for all atomic pairs.. Finally, the last two terms in equation (1) correspond to transverse ($\Omega$) and longitudinal ($\Delta$) magnetic fields, which will be used for the state-preparation.

The zero-temperature phase diagram of the Hamiltonian (1) is schematically shown in figure 1(c). It hosts a stripe phase for \left|\Omega\right|,\,\left|\Delta\right|$$<\left|c_{d}\right| with spins aligned along the $z$-axis and anti-aligned along the $y$-axis [c.f. figure 1(b)]. This phase is destroyed by increasing magnetic fields when \left|\Omega\right|,\,\left|\Delta\right|$$\sim\left|c_{d}\right| (see [26, 27, 28] in the context of Rydberg atoms), and for \left|\Omega\right|,\,\left|\Delta\right|$$\gg\left|c_{d}\right| the system is in the ferromagnetic phase with spins oriented along the total magnetic field [29]. In particular, for $\Omega\gg\Delta,|C_{ij}|$, all spins align preferentially along the $x$-axis resulting in the ground state $\left|{G}\right\rangle=\left|{\rightarrow}\right\rangle_{1}...\left|{\rightarrow}\right\rangle_{N}$, with $\left|{\rightarrow}\right\rangle=(\left|{\uparrow}\right\rangle-\left|{\downarrow}\right\rangle)/\sqrt{(}2)$.

2.2 Spin-1 model

The second model is for spin-1 particles, with three internal atomic states representing the spin-one states $\left|{-1}\right\rangle$, $\left|{0}\right\rangle$ and $\left|{1}\right\rangle$ (see section 4 for details), described by the following Hamiltonian

[TABLE]

where $N$ is the number of spins, $S_{z}^{(i)}=\left|{1}\right\rangle_{i}\left\langle{1}\right|_{i}-\left|{-1}\right\rangle_{i}\left\langle{-1}\right|_{i}$, $S_{+}^{(i)}=\sqrt{2}(\left|{1}\right\rangle\left\langle{0}\right|+\left|{0}\right\rangle\left\langle{-1}\right|)$ and $S_{-}^{(i)}=(S_{+}^{(i)})^{\dagger}$ are spin one operators, $\beta$, $\gamma$ and $\delta$ are constants depending on chosen atomic states in the implementation scheme (see section 4.3). The first term in this Hamiltonian corresponds to an external “magnetic” field which acts only on the [math]-th component of the spin and serves as a controlling parameter during state preparation, while the last two terms represent the DDI including a XXZ-type Hamiltonian (the second term) and higher order spin operators (the last term). Obviously [c.f. figure 1], the state $\left|{G^{\prime}}\right\rangle=\left|{0}\right\rangle_{1}...\left|{0}\right\rangle_{N}$ minimises the energy of the Hamiltonian for large negative $\Delta E$, $\Delta E\ll-\left|C_{ij}\right|$. For small $\Delta E$ ($\left|\Delta E\right|\ll\left|C_{ij}\right|$), the spin exchange term of the Hamiltonian leads to a disordered phase. Finally, for large positive $\Delta E$, $\Delta E\gg\left|C_{ij}\right|$, the DDI arranges spins into the stripe phase.

3 Dynamical preparation of the stripe phase

We now proceed with the description of the stripe-phase preparation protocols for both the above models, which is based on the adiabatic evolution of the ground state [28, 30, 31, 32]. This is achieved by changing adiabatically the magnetic field from the values at which the ground state is simple and easily realisable in experiments, to the values at which the stripe phase is the ground state. To characterise the state of the system during the protocol, we introduce the stripe order parameters along the $z$ and $y$ lattice directions [c.f. figure 1(b)]

[TABLE]

where $S_{z}^{(i)}$ are the spin operators of the corresponding model ($S_{z}^{(i)}\rightarrow\sigma_{z}^{(i)}$ for the spin-$1/2$ case), and $1_{\alpha}$ represents a shift by one site along the $\alpha$-direction. The system is in the stripe phase when $M_{z}\rightarrow 1$ and $M_{y}\rightarrow-1$.

3.1 Dynamical state preparation of stripe phases for the spin-$1/2$

model

For the spin-$1/2$ model, we consider initially (at time $t=t_{0}$) the system with $\Delta(t_{0})=0$, $\Omega(t_{0})\gg\left|c_{d}\right|$ in its ground state $\left|{G}\right\rangle$, which can be experimentally prepared with a high fidelity, and then adiabatically decrease $\Omega$ to reach the final point $\Delta(t_{f})=0$, $\Omega(t_{f})=0$ at time $t_{f}$, where the stripe phase is the ground state. Note, however, that during the parameter ramp, we cross the phase transition point. Therefore, strictly speaking the criterion of adiabaticity can only be fulfilled in a finite system where the finite-size effects ensure the existence of a finite gap even at the transition point. Nevertheless, keeping in mind that experimental realisations with cold atoms are always dealing with finite systems, this does not cause any real problems, but requires the entire time $t_{f}-t_{0}$ of the preparation protocol be sufficiently large in order to minimise the population of excited states along the ramp (see below).

We test the protocol with an approach based on a time dependent variational principle (TDVP) which allows us to simulate the dynamics of the system with large number of atoms [28, 33]. This makes possible, in particular, to estimate the minimal coherence time needed to obtain the stripe phase experimentally with a good fidelity. The justification of our approach is based on the comparison with the results of exact diagonalisation (ED) for small systems.

Within the TDVP approach, the time-evolution of a many-body quantum state of the spin-$1/2$ model is described via a product state wave-function

[TABLE]

where $\alpha_{i}(t)$ and $\beta_{i}(t)$ are the time-dependent coefficients satisfying the normalisation condition $|\alpha_{i}|^{2}+|\beta_{i}|^{2}=1$. Note that the ground state of the Hamiltonian (1) with $\Delta=0$ in the classical limit $\Omega\rightarrow 0$ (the stripe phase), is described by this ansatz with the red and blue circles in figure 1(b) corresponding to $|\alpha_{i}|^{2}=1,\,\beta_{i}=0$ and $\alpha_{i}=0,\,|\beta_{i}|^{2}=1$, respectively. Our mean field (MF) ansatz, however, neglects the quantum pair correlation in the non-classical limit, $\Omega\neq 0$.

The time evolution of variational coefficients $\alpha_{i}(t)$ and $\beta_{i}(t)$ with changing of $\Omega$ and $\Delta$ can be calculated by the Euler-Lagrange equations [28]

[TABLE]

[TABLE]

where $L$ is the Lagrangian of the system

[TABLE]

After substituting the MF ansatz (4) and the Hamiltonian (1) into (7), the Euler-Lagrange equations (5) and (6) give

[TABLE]

[TABLE]

For $C_{ij}=0$, equations (8) and (9) describe the single-spin Rabi oscillations with the frequency $\sqrt{\Delta^{2}+\Omega^{2}}$. We see that the interactions $C_{ij}$ in the MF description generate nonlinear corrections to the longitudinal magnetic field $\Delta$ so that the rate of local spin flips depends on the spin configuration on the neighbouring sites.

The ansatz in equation (4) can also be used to investigate the ground state of the system by considering the parameters $\alpha_{i}$ and $\beta_{i}$ as independent variational parameters (subjected to the normalisation constraint only) and minimising the energy of the system. In figure 2(a) we present the results of such calculations for a $6\times 6$ lattice by showing the quantity $M_{z}M_{y}$ as a function of $\Delta$ and $\Omega$. As expected, the stripe phase ($M_{z}M_{y}\rightarrow-1$) is the variational ground state of the system when the DDI is larger than the transverse and longitudinal magnetic fields $\Delta$ and $\Omega$ [c.f. figure 2(b)]. With increasing magnetic fields, the stripe phase disappears [c.f. figure 2(c)], and for large $\Omega$ and $\Delta$ the ground state is the ferromagnetic one (for the parameters in figure 2 (c), the spins are mostly oriented in the $x$-direction, $S_{x}=N^{-1}\sum_{i}\left\langle\sigma_{x}^{(i)}\right\rangle=-0.96$).

We now analyse the dynamical preparation of the stripe phase along the trajectory shown in figure 1(b) with $\Delta=0$ for all times and $\Omega$ decreasing from its initial value $\Omega=\Omega_{0}$ to 0, as represented in figure 2(d). For the system initialised in the state $\left|{G}\right\rangle=\left|{\rightarrow}\right\rangle_{1}...\left|{\rightarrow}\right\rangle_{N}$, with $\left|{\rightarrow}\right\rangle=(\left|{\uparrow}\right\rangle-\left|{\downarrow}\right\rangle)/\sqrt{2}$, figures  2(e) and (f) present the final values of $M_{z}$ and $M_{y}$ as a function of $\Omega_{0}$ and $t_{f}$. We see that the larger initial value $\Omega_{0}$ the longer it takes to reach the stripe phase with high fidelity. For $\Omega_{0}=3c_{d}$ the stripe phase is reached already for $t_{f}\sim 10^{2}t_{d}$, where $t_{d}=\hbar/c_{d}$.

To validate our MF approach we compare it with the exact numerical solutions of the time-dependent Schrödinger equation for the same initial state $\left|{G}\right\rangle=\left|{\rightarrow}\right\rangle_{1}...\left|{\rightarrow}\right\rangle_{N}$ on the $4\times 4$ lattice ($N=16$). The comparison is presented in figure 3, and shows that the two methods are in a very good agreement when comparing the energy. The order parameters (spin-spin correlators) $M_{z,y}$ are in agreement within 10% for all times except those close to the transition point, where the TDVP simulations show strong oscillations. They correspond to a large-amplitude motion of the system over a practically flat potential landscape in the low-energy manifold of states near the transition. Away from the transition point, when the system chooses one of the two possible stripe configurations (see also the discussion below), the oscillations disappear. In contrast, the numerical solutions of the Schrödinger equation do not show any noticeable oscillations. We attribute this difference to the fact that the evolution within the ansatz (4) allows to reach only a very limited (as compared to the true quantum-mechanical evolution) part of the Hilbert space, making impossible the smearing out of the total contribution of various oscillating terms in the spin-spin correlation functions.

In figure 4, we have plotted the evolution of the occupation probabilities of the states $\left|{\uparrow}\right\rangle$ (red curve) and $\left|{\downarrow}\right\rangle$ (blue curve) to illustrate the dynamics of the formation of the stripe pattern within the TDVP. It is important to mention that within the TDVP the system chooses spontaneously one of the two possible stripe configurations. We believe that this spontaneous symmetry breaking is a result of finite numerical precision during the calculation. In contrast, our numerical integration of the time-dependent Schrödinger equation always results in a superposition of the two possible stripe configurations, showing therefore no sign of spontaneous symmetry breaking. Such “robustness” of the numerical solutions of the Schrödinger equation is due the quantum correlations (entanglement) over a large part of the Hilbert space, which is built in the system during the purely quantum-mechanical time evolution. In real physical systems, however, decoherence processes (due to finite temperature, noise, etc.) destroy most of these correlations and cause spontaneous symmetry breaking.

As a final comment we note that, as we have already seen in figure 2, the larger $\Omega_{0}$ is, the longer ramps are required to prepare the stripe phase.

3.2 Dynamical state preparation of stripe phases for the spin-1 model

For the preparation of the stripe phase within the spin-$1$ model (LABEL:eq:hamil_one), we start with $\Delta E_{0}=\Delta E(t_{i})\ll-\left|C_{ij}\right|$ in the Hamiltonian and initialise the system in the state $\left|{G^{\prime}}\right\rangle=\left|{0}\right\rangle_{1}...\left|{0}\right\rangle_{N}$. We then increase $\Delta E$ adiabatically to a large positive value $\Delta E_{f}=\Delta E(t_{f})\gg|C_{i,j}|$ [c.f. figure 1(d)] to end up in the stripe phase formed on the manifold of spin states $\left|{1}\right\rangle_{i}$ and $\left|{-1}\right\rangle_{i}$ (the spin states $\left|{0}\right\rangle_{i}$ are effectively eliminated by the large positive $\Delta E_{f}$).

Before proceeding to the analysis of the time-evolution of the system under the ramp we note that the MF variational ansatz used above for the spin-$1/2$ model, can not be applied here because it is unable to describe the flip-flop processes of the type $\left|{00}\right\rangle\longleftrightarrow\left|{1,-1}\right\rangle$ involving two spins, see A. We thus simulate the dynamics of the spin-1 model via a time-dependent cluster mean field (CMF) approach with clusters containing two spins along the $y$-axis [c.f. figure 1(a)]. The wave function of the $\mu$-th cluster is written as

[TABLE]

where $l,l^{\prime}=-1,\,0,\,1$, the indices $1$ and $2$ label the particles in the cluster, and $\alpha_{ll^{\prime}}^{(\mu)}$ are the variational parameters subjected to the normalisation constraint $\sum_{l,l^{\prime}}\left|\alpha_{ll^{\prime}}^{(\mu)}\right|^{2}=1$. The corresponding ansatz for the wave function of the system is a product state of all clusters

[TABLE]

Based on this ansatz and the Hamiltonian (LABEL:eq:hamil_one), we construct the Lagrangian (7) and the corresponding Euler-Lagrange equations describing the time evolution of the variational parameters $\alpha_{ll^{\prime}}^{(\mu)}$, which are rather lengthy and presented in B.

We now apply this scheme to study the dynamical preparation of the stripe phase with the linear ramp of $\Delta E(t)$. We consider in the following the parameters of the Hamiltonian LABEL:eq:hamil_one as $\beta=-2.75$, $\gamma=-1.85$ and $\delta=-1.25$ corresponding to the case of our implementation with Erbium atoms. In doing the calculations, we also add the term

[TABLE]

to the Hamiltonian (LABEL:eq:hamil_one), where $\epsilon(t)$ is the small energy shift [$\epsilon(t_{f})=0$] between the two sites (and, therefore, between the states $\left|{-1,1}\right\rangle$ and $\left|{1,-1}\right\rangle$) in each clusters, to break the up-down symmetry and obtain the stripe phase [c.f. B]. Further details and the results of the calculations on a $8\times 8$ lattice are presented in figure 5. Panel (a) shows an example of the function $\Delta E(t)$ we used for the state preparation. In panel (b) we present the occupation of the state $\left|{-1}\right\rangle$ at the end of evolution starting initially from the state $\left|{G^{\prime}}\right\rangle$, which clearly shows the formation of the stripe phase. For the final value $\Delta E_{f}=50c_{d}$, the spin-spin correlations $M_{z}$ and $M_{y}$ in the final state are shown in figures 5(c) and 5(d), respectively, as a function of $\Delta E_{0}$ and final time of evolution $t_{f}$. One clearly sees the effect of the speed of the ramp on the fidelity of the state preparation.

Finally in figure 6 we compare the results of the time-dependent CMF (solid lines) with those obtained by direct numerical integration (NI) of the time-dependent Schrödinger equations for the system on a $3\times 4$ lattice. In this figure, the dashed red (dotted blue) lines are the results of the NI with (without) the symmetry breaking term $H_{\mathrm{SB}}$ used in the CMF. For all the calculations we used the linear ramp with $\Delta E$ linearly changing from $\Delta E_{0}=-50c_{d}$ to $\Delta E_{f}=50c_{d}$ between $t_{i}=0$ and $t_{f}=5$.

Concerning the final states of the evolution, the results show good agreement between the NI and CMF methods: the energies match within less than one percent, whilst the spin correlations are overestimated (around twenty percent) by the CMF method. During the evolution, however, the deviations between the two methods are much more pronounced. This is not a surprise, keeping in mind the simplicity of the wave function in the CMF ansatz, equations (10) and (11), and have the same origin (very limited part of the Hilbert space used in the CMF ansatz) as in the case of spin-$1/2$ model discussed above.

4 Microscopic realisation of spin models with magnetic atoms

In this section we present an implementation of aforementioned spin models on the basis of magnetic atoms, with spin states being encoded in the Zeeman sublevels of the atomic ground state. We note that, being based on generic properties of the DDIs, the models can also be realised in other experimental platforms such as polar molecules [5, 34], Rydberg atoms [35, 36], and trapped ions [37, 38, 39].

The setup we have in mind is represented schematically in figure 7 with magnetic atoms placed on a square lattice and prepared in their electronic ground state manifold. Considering the case of a fine structure manifold with angular momentum $J$ 222We consider for simplicity the absence of hyperfine interactions., our two models will involve a restricted number of levels ($2$ and $3$), compared to the total number of states $2J+1$ [c.f. figure 7 (b)].

We start with reviewing basic properties of the DDI between magnetic atoms, then show how to select energetically a particular subset of Zeeman levels to implement our two spin models. Finally, we discuss the resulting relevant time scales for the stripe phase preparation.

4.1 Dipole-dipole interactions between magnetic atoms

For two magnetic atoms, which are located at sites $i$ and $j$ of a two-dimensional square lattice, the operator of the DDI reads [2]

[TABLE]

with $c_{d}=\mu_{0}(g_{J}\mu_{B})^{2}/(4\pi a^{3})$, where $\mu_{0}$ is the magnetic permeability of the vacuum, $\mu_{B}$ is the Bohr magneton, $g_{J}$ is the Landé factor, $J^{i(j)}$ denotes the operator of the total angular momentum of the particle $i\,(j)$, and $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$ is the relative position of two particles ($\hat{\bm{r}}_{ij}=\bm{r}_{ij}/r_{ij}$) in units of the lattice spacing $a$. Assuming the presence of a strong external magnetic field $\mathbf{B}=B\mathbf{\hat{z}}$ and choosing its direction as the quantisation axis, the DDI Hamiltonian can be written as [40]

[TABLE]

where the first term

[TABLE]

conserves the total angular momentum of two atoms. Here $J_{ij}=(c_{d}/r_{ij}^{3})(1-3\cos^{2}\theta_{ij})$ with $\theta_{ij}$ being the angle between $\bm{r}_{ij}$ and the $z-$axis. The second and the last terms of the Hamiltonian (14) are

[TABLE]

[TABLE]

where $J_{\pm}=J_{x}\pm iJ_{y}$ and $\phi_{ij}$ is the azimuthal angle of $\bm{r}_{ij}$ in the $x-y$ plane. These terms do not conserve energy and transfer one ($H_{q=1}$) or two ($H_{q=2}$) units of the internal angular momentum of the atoms to or from their relative orbital motion. In this work, we only consider the term $H_{q=0}$ assuming the Zeeman splitting induced by the magnetic field to be much larger than the DDI, which makes the contribution of non-angular momentum conserving terms $H_{q\neq 0}$ negligible.

4.2 State-dependent energy shifts

To exclude energetically all unwanted states from the dynamics, we use AC Stark shifts in combination with the linear Zeeman shifts. To this end we consider two polarised laser beams with polarisation $\sigma_{\pm}$ and same frequency $\omega$, which couple off-resonantly the states $\left|{m}\right\rangle\equiv|J,m_{J}\rangle$ of the ground state manifold to a manifold of excited states $\left|{m^{\prime}}\right\rangle\equiv|J^{\prime},m_{J}^{\prime}\rangle$, [c.f. figure 7]. The corresponding Hamiltonian within the rotating wave approximation and in the frame rotating with the laser frequency $\omega$ is given by

[TABLE]

where different Landé factors in the two manifolds give rise to different Zeeman shifts $\alpha^{(\prime)}=\mu_{B}g_{J^{(\prime)}}B$. We note that the Rabi frequencies $\Omega_{\pm}^{(m)}=\Omega_{\pm}C_{\pm}^{(m)}$ depend on m due to the different Clebsch-Gordan coefficients $C_{\pm}^{(m)}$ involved in the electric dipole transitions. Finally, $\Delta=\omega-\omega_{0}$ is the detuning where $\omega_{0}$ denotes the resonant frequency between the ground and the first excited states in the absence of the magnetic field. Taking the laser coupling as a perturbation ($|\Delta-\alpha^{\prime}(m\pm 1)+\alpha m|\gg\Omega_{\pm}^{(m)}$), the effective Hamiltonian for the ground state manifold in the second order perturbation theory can be written as

[TABLE]

where $\epsilon_{m}=\frac{\Omega_{+}^{(m)^{2}}}{\Delta-\alpha^{\prime}(m+1)+\alpha m}+\frac{\Omega_{-}^{(m)^{2}}}{\Delta-\alpha^{\prime}(m-1)+\alpha m}$. (Note that we neglected the two photon Raman transitions $|m\rangle\to|m+2\rangle$ assuming $\alpha,\Delta\gg\Omega_{\pm}^{(m)}$.) The AC stark-shifts $\epsilon_{m}$ depend non-linearly on $m$, which will allow energy-conserving dynamics only on a restricted set of two (three) Zeeman states in the ground state manifold for the spin $1/2$ (spin $1$) model. We also note that an alternative way to realise the nonlinear energy shifts is by using the quadratic Zeeman effect.

4.3 Implementation of the spin-1 model

Following the discussion above, the effective Hamiltonian governing the dynamics on the ground state manifold has the form $H=H_{q=0}+h_{J}$. To implement the spin-1 model from section 2.2, we consider, for example, the three adjacent magnetic levels $\left|{m=m_{0}-1,m_{0},m_{0}+1}\right\rangle$ with $m_{0}=-J+1$ [c.f. figure 7] (another option would be $m_{0}=J-1$), and impose the condition $|2\epsilon_{m_{0}+1}-\epsilon_{m_{0}}-\epsilon_{m_{0}+2}|\gg c_{d}$, which suppresses the process $\left|{m_{0}+1,m_{0}+1}\right\rangle\to\left|{m_{0},m_{0}+2}\right\rangle$ (due to $H_{q=0}$) energetically. This allows us to project the Hamiltonian $H$ onto the submanifold with $m=m_{0}-1,m_{0},m_{0}+1$ defining the spin-one states $\left|{k}\right\rangle$ ($k=-1,0,1$) as $\left|{k}\right\rangle=\left|{m=m_{0}+k}\right\rangle e^{-i(\alpha(m_{0}+k)+\epsilon_{m_{0}+k}-\Delta E/2\delta_{k,0})t}$ with $\Delta E=2\epsilon_{m_{0}}-\epsilon_{m_{0}-1}-\epsilon_{m_{0}+1}$. One then gets the parameters $\beta=-J_{1}^{2}/8$, $\gamma=-J_{1}(J_{1}-J_{2})/8$ and $\delta=-(J_{1}-J_{2})^{2}/8$ with $J_{1}=\sqrt{J(J+1)-m_{0}(m_{0}+1)}$ and $J_{2}=\sqrt{J(J+1)-m_{0}(m_{0}-1)}$. Note that the energy shift $\Delta E$ can be varied in time by using the time-dependent Rabi frequencies $\Omega_{\pm}$ and the detuning $\Delta$. As an example, we consider ground state bosonic Erbium atoms with $J=6$ 333The hyperfine-structure is absent in this case. ($g_{J}=1.1638$) [41], which are trapped in a square lattice with lattice spacing $a=266$ nm [11]. This results in the dipole-dipole coupling $c_{d}=2\pi\times 0.9~{}\hbar\mathrm{Hz}$. In order to satisfy the requirements above for preparing the stripe phase for the case of Erbium atoms with $m_{0}=-5$, one can apply two laser beams with the Rabi frequencies $\Omega_{\pm}=2\pi\times 0.5~{}\mathrm{MHz}\pm x$ and the detuning $\Delta=0.8\alpha$, which couple the ground state manifold to the manifold with $J^{\prime}=7$ ($g_{J^{\prime}}=1.070$). In the presence of a magnetic field $B=5~{}\mathrm{G}$ ($\alpha=51.186~{}\hbar\mathrm{MHz}$ and $\alpha^{\prime}=47.06~{}\hbar\mathrm{MHz}$), we get figure 8, which shows all the requirements are fulfilled.

Finally, we emphasise that the spin-$1$ model only involves magnetisation conserving terms and is thus robust (in first order in the Zeeman shifts) against time-dependent fluctuations of the magnetic field. This is not the case for the spin 1/2 model discussed below.

4.4 Implementation of the spin-$1/2$ model

The spin-$1/2$ model can be implemented by choosing two states $\left|{m=-J}\right\rangle$, $\left|{m=-J+2}\right\rangle$ in the ground state manifold as spin states, $\left|{-J}\right\rangle=\left|{\uparrow}\right\rangle$ and $\left|{-J+2}\right\rangle=\left|{\downarrow}\right\rangle$ [c.f. figure 7]. To exclude the other Zeeman states from the dynamics, we require the conditions $|\epsilon_{-J}+\epsilon_{-J+2}-2\epsilon_{-J+1}|\gg c_{d}$ and $|\epsilon_{-J+3}+\epsilon_{-J+1}-2\epsilon_{-J+2}|\gg c_{d}$. The resulting projection of the Hamiltonian $H_{q=0}$ onto the submanifold $\left|{-J+2}\right\rangle,\left|{-J}\right\rangle$ leads to the first term in (1). To obtain the other terms we can use an additional Raman coupling of the two states $\left|{-J}\right\rangle$ and $\left|{-J+2}\right\rangle$. In the rotating wave-approximation, this coupling takes the form of the two magnetic fields $\Omega$ and $\Delta$ which appear in (1) 444Note that we can neglect the transition between other ground state levels $|m\rangle$ assuming that the Landé factors between ground ($J$) and excited manifolds allow to make state-selective Raman transition between $\left|{-J+2}\right\rangle$ and $\left|{-J}\right\rangle$. For the purpose of the adiabatic state preparation, these fields can be controlled by changing the Rabi-frequency and the detuning of the additional Raman coupling (see, for example, [32]).

4.5 Time scales

With the above implementation schemes in mind, we can now estimate the time scales required for our state-preparation protocol under realistic experimental conditions. Coming back to the previously discussed example, i.e. Erbium atoms, one gets $2\pi\times 3.66~{}\hbar\mathrm{Hz}$ for the exchange interaction. This gives $t_{f}<1\mathrm{s}$ for the time required to adiabatically prepare the stripe phase in the spin-1 case ($t_{d}\sim 0.177\mathrm{s}$), which is well compatible with experimental time scales. (For the spin 1/2 case, the required time is order of magnitude larger, which does not look realistic.) We also note that the Raman couplings involved in the implementation of the models give rise to decoherence with the rate $\Gamma\sim(\frac{\Omega}{\Delta})^{2}\Gamma_{r}$, where $\Gamma_{r}$ is the line width of the excited states involved. This rate should be much smaller than the exchange interaction, which requires $\Delta\gg\Omega,\Gamma_{r}$ [16, 42, 43].

5 Conclusion

We have shown how systems of magnetic atoms in optical lattices can be used to realise stripe phases in current experiments. On the methodological side, we have employed time-dependent variational ansätze which, as also reported in [28], seem particularly well suited to address adiabatic state preparation for magnetically ordered states. Most importantly, they enable the investigation of systems sizes comparable with current cold atom experiments, that is well beyond what can be reached with exact methods. We remark that our many-body results can find immediate applications to other physical settings, such as polar molecules and Rydberg atoms in optical lattices [2], superconducting qubits [24, 25], where both Heisenberg-type and Ising Hamiltonians with dipolar interactions can be realised.

Our results show how the combination of long coherence times and flexibility in initialising and manipulating magnetic atoms represent an ideal tool to explore quantum magnetism, despite the interaction strengths being weaker (in absolute value) with respect to polar molecules and Rydberg gases. Moreover, the tunability of the magnetic field orientation/strength and the properties of the DDIs also allow in principle to realise other interesting quantum phases, such as integer or fractional Chern insulators [44, 45, 46, 47, 48], for which it would be intriguing to formulate a proper time-dependent variational principle to correctly reproduce the exact dynamics.

Acknowlegments

We acknowledge inspiring and fruitful discussions with P. Zoller. We also thank S. Baier, J. H. Becher, L. Chomaz, F. Ferlaino, M. J. Mark, and G. Natale for helpful discussions on several aspects of the implementation schemes. The exact diagonalisation simulations were performed using the QuTiP toolbox [49]. Research in Innsbruck is supported by the European Research Council (ERC) Synergy Grant UQUAM, EU H2020 FET Proactive projects RySQ, and by the Austrian Science Fund through SFB FOQUS (FWF Project No. F4016-N23).