Self-generated quantum gauge fields in arrays of Rydberg atoms

Due to their strong and non-local interaction, their experimental accessibility and the high degree of tunability and control that can be exerted, Rydberg atoms have become a powerful tool for simulating strongly correlated quantum many-body systems [1]. The strong interaction exhibited by Rydberg atoms extending over μm distances, while at the same time featuring high stability, allows to use optical tweezers and atom-by-atom assembly to prepare arbitrary one- and two-dimensional arrays of atoms [2–4]. The repulsive van-der-Waals interaction between pairs of Rydberg atoms then provides a natural way to implement spin models with Ising-type zz interactions [5] and to study them with full coherent control. Experiments include the observation of quantum phase transitions to ordered phases in one-dimensional quantum models with Ising interactions extending beyond the nearest neighbor (NN) [6–8] and two-dimensional square lattices [9–12], as well as the probing of topological spin liquids in dimer models based on Rydberg atoms in Kagome lattices [13]. Here the spin degree of freedom is formed by a low-lying level and a Rydberg state of the atoms. If the spin is made up out of two Rydberg states with a dipole-allowed transition, excitation exchange between pairs of atoms corresponds to an xy spin interaction, which has been used to experimentally investigate symmetry protected topological lattice models of interacting bosons [14].

The realization of external magnetic fields poses a challenge to all platforms based on neutral atoms. To circumvent this problem, different physical effects that mimic the behavior of a magnetic field on charged particles have been proposed including laser assisted tunneling [15], periodic driving [16] or the transfer of angular momentum from light [17]. Recently, in [18] it was demonstrated that spin–orbit coupling in systems of Rydberg atoms can be used to generate Peierls phases in the hopping matrix elements of Rydberg excitations without external light fields, laser assisted tunneling or lattice shaking. The authors showed that the Peierls phases, which are the manifestation of external magnetic fields in lattice models, lead to chiral currents in a triangle-configuration of three atoms. Moreover the complex amplitudes of the hopping process between two Rydberg atoms depended on the occupation of Rydberg states of a third atom in between the two atoms. Such density-controlled hopping processes form one elementary building block for a direct implementation of even more complex many-body models such as dynamical gauge theories [19, 20].

In the present paper we theoretically investigate a simple one-dimensional extension of the triangle system studied in [18] with density-dependent complex hopping of Rydberg excitations competing with repulsive density–density interactions. Interacting bosons on a lattice in an external effective magnetic field can give rise to a transition between a Meissner and a vortex phase [21–23] upon tuning the magnetic flux. If the magnetic field becomes density-dependent this transition can be driven by density alone [24]. Additionally, we show that in the Rydberg system a different type of transition between a Meissner phase and vortex phase occurs. We will show that the Rydberg model features besides a static, classical gauge potential a dynamical quantum component linked to density modulations. Although being a dynamical quantity, the quantum component of the gauge field is gauge invariant as it is a transverse field. Nevertheless, it substantially modifies the ground-state phase diagram, which we numerically investigate using exact diagonalization (ED) methods. It leads to two new liquid phases in addition to a trivial superfluid, which are characterized by alternating vortex currents of Rydberg excitations. In an ordered phase where the repulsive interactions dominate a regular array of current vortices emerges. Surprisingly the quantum gauge field is spontaneously generated even without additional density–density interaction terms, reminiscent of the Higgs-mechanism.

The paper is organized as follows: in section 2 we shortly illustrate the microscopic origin of the model. In the subsequent section 3 we discuss the Rydberg spin model and show its phase diagram as well as the emergence of chiral currents. In order to gain some understanding of the phase diagram we transform the spin model, equivalent to hard core bosons, to a fermion model and discuss a mean-field approximation of the fermionic model in section 4. In section 5 we reformulate the model in terms of operator valued transverse gauge fields which provides an explanation for the vortex-currents found in the numerical simulations. In section 6 we make a few remarks about possible experimental signatures and in section 7 we summarize our results and give an outlook to further studies.

To understand how Rydberg interactions can be tuned to exhibit density-dependent Peierls phases, we first consider the physical behavior of a minimal triangle setup of atoms, shown in figure 1. In this setup, three identical Rydberg atoms are placed on an equilateral triangle in the xy-plane. Furthermore, two Rydberg states are chosen to represent the two spin-states $\left\vert 1\right\rangle$ and $\left\vert 0\right\rangle$. The third atomic state, $\left\vert +\right\rangle$, is detuned and will be adiabatically eliminated, yielding non-resonant processes that scale as the inverse detuning 1/Δ. The level scheme shown in figure 1 is created by external electric and magnetic fields which provide the required energy splittings. This minimal scheme has been experimentally realized in [18], where the authors explain under what conditions the manifold of Rydberg states in the 60S_1/2, 60P_1/2, and 60P_3/2 levels can be reduced to these three states and provide concrete values of detuning and magnetic field strength. For a short summary regarding the restriction to three atomic states, see appendix A. When the magnetic field is orthogonal to the plane of the atoms, the microscopic Hamiltonian coupling atoms i and j can be written as

$$\begin{equation}{\hat{V}}_{ij}=\frac{1}{4\pi {{\epsilon}}_{0}{R}_{ij}^{3}}\left[{\hat{d}}_{i}^{z}{\hat{d}}_{j}^{z}+\frac{1}{2}\left({\hat{d}}_{i}^{+}{\hat{d}}_{j}^{-}+{\hat{d}}_{i}^{-}{\hat{d}}_{j}^{+}\right)-\frac{3}{2}\left({\hat{d}}_{i}^{+}{\hat{d}}_{j}^{+}{e}^{-2i{\phi }_{ij}}+{\hat{d}}_{i}^{-}{\hat{d}}_{j}^{-}{e}^{2i{\phi }_{ij}}\right)\right].\end{equation} \tag{ 1 }$$

Here, ${R}_{i,j}=\left\vert {\boldsymbol{r}}_{i,j}\right\vert$ and ϕ_ij refer to the interatomic distance and angle in polar coordinates, respectively. Additionally, the operators ${\hat{d}}_{i}^{\pm }=\mp \frac{1}{\sqrt{2}}({\hat{d}}_{i}^{x}\pm i{\hat{d}}_{i}^{y})$ represent ladder operators raising and lowering the angular momentum of the ith Rydberg atom. In (1) we recognize the flip–flop terms ${\hat{d}}_{i}^{+}{\hat{d}}_{j}^{-}$ as well as terms ${\hat{d}}_{i}^{+}{\hat{d}}_{j}^{+}{e}^{-2i{\phi }_{ij}}$ that increase the angular momentum of both involved atoms and are multiplied by a complex phase.

Zoom In Zoom Out Reset image size

As a result of the interaction Hamiltonian (1) it is useful to consider two different types of processes. In this section, we will only provide a qualitative explanation, the microscopic derivation is given in appendix A.

We consider the states $\left\vert 1\right\rangle$ and $\left\vert 0\right\rangle$ as the two states of a spin-1/2 system and transform to hard-core bosons for a more intuitive picture. The first and simplest process is the resonant dipole–dipole interaction, where the term ${\hat{d}}_{i}^{+}{\hat{d}}_{j}^{-}$ in (1) causes an exchange of population in the $\left\vert 0\right\rangle$ and $\left\vert 1\right\rangle$ states of the atoms i and j. This mechanism yields for the hard-core bosonic Hamiltonian ${\hat{d}}_{i}^{+}{\hat{d}}_{j}^{-}\to {\hat{b}}_{j}^{{\dagger}}{\hat{b}}_{i}$, where the operators ${\hat{b}}_{i}$ and ${\hat{b}}_{i}^{{\dagger}}$ are annihilation and creation operators for hard-core bosons (spin-1/2) at site i, respectively, which satisfy the hard-core constraint ${\hat{b}}_{i}^{2}={\left({\hat{b}}_{i}^{{\dagger}}\right)}^{2}=0$ at each site i.

The second, more interesting process is the non-resonant excitation transfer involving three sites of the lattice. For this process, illustrated in figure 2, the terms quadratic in raising or lowering operators in (1) are responsible. An important feature is here that de-excitation of the blue atom drives the non-resonant process in the red atom only if the latter is in the ground state $\left\vert 0\right\rangle$. In this manner, the hopping process shown in figure 2 connects two sites (blue and green atoms) depending on the density (occupation of state $\left\vert 1\right\rangle$) of another atom (red). Since the corresponding terms of (1) also include a complex phase, we obtain the following term in hard-core boson language ${\hat{b}}_{j}^{{\dagger}}{\hat{b}}_{i}\left(1-{\hat{n}}_{k}\right){e}^{-4i{{\epsilon}}_{i,j,k}\alpha }$. Here, _i,j,k is the Levi-Civita symbol which distinguishes between a clockwise and counter-clockwise motion of the excitation, and α = ϕ_ij is the interatomic angle for all i ≠ j.

Zoom In Zoom Out Reset image size

Lastly, we need to take into account a special case of the second process, where the initial and the final atoms are identical. In that situation, we no longer have a hopping process, but in fact a NN density–density interaction ${\hat{b}}_{i}^{{\dagger}}{\hat{b}}_{i}\left(1-{\hat{n}}_{k}\right)$, where the complex phase drops out.

Summarising all the terms discussed so far, we can write down the total Hamiltonian for the triangle setup discussed in figure 1

$$\begin{equation}\hat{H}=-J\sum\limits _{i\ne j=1}^{3}{\hat{b}}_{i}^{{\dagger}}{\hat{b}}_{j}+\text{h.c.}-2gJ\sum\limits _{i\ne j\ne k}^{3}{\hat{b}}_{i}^{{\dagger}}{\hat{b}}_{j}\left(1-{\hat{n}}_{k}\right){e}^{-4i{{\epsilon}}_{i,j,k}\alpha }+\text{h.c.}-2gJ\sum\limits _{i\ne k}^{3}{\hat{n}}_{i}\left(1-{\hat{n}}_{k}\right).\end{equation} \tag{ 2 }$$

Here, J is the natural energy scale of the system and $g=\frac{27J}{2{\Delta}}$ is the rescaled detuning of the atomic state $\left\vert +\right\rangle$. The numerical factors in g are explained in the microscopic derivation given in appendix A. This model has been experimentally investigated by the authors of [18], where they confirmed the density-dependence of the chiral motion of excitations on the triangular setup.

In the following we consider an extension of this model that extends in 1D in the form of a zig-zag chain, as shown in figure 1. As a result, more density-dependent terms and longer-range processes enter the Hamiltonian given in (2), but the underlying principle remains unchanged.

3.1. Model

We now consider a one-dimensional zig-zag chain formed by the elementary triangles of Rydberg atoms from the previous section, see figure 1(b). We set the physical angle α = π/3 in which case the triangles are equilateral, and each site j now has four equivalent NNs. As a result, the indirect (density-dependent) hopping terms as well as the NN interaction outlined in the previous section now include all possible NN-combinations as shown in figure 3. In order to account for the fact that more internal Rydberg states may be relevant in the interaction we allow for additional NN density–density interactions resulting from off-resonant van-der Waals interactions. Their relevance can in fact be controlled by proper choice of Rydberg states. For example, the total C₆ interaction coefficients for the Rubidium $\left\vert 60s,60s\right\rangle ,\left\vert 60p,60s\right\rangle$ and $\left\vert 60p,60p\right\rangle$ states are different [25], which leads to a modification of the density–density interaction from the one in (2). We here quantify this effect by a dimensionless parameter η, where η = 1 corresponds to the case when no additional van-der Waals interactions are present (i.e. the case of a three-level atom), and η ≫ 1 means that van-der Waals interactions dominate over the second order processes discussed in the previous section. Note that η depends on the principal quantum number of the Rydberg level and thus is not a continuously tuneable parameter. The complete many-body Hamiltonian for the zig-zag chain then reads up to irrelevant homogeneous single-particle energy shifts

$$\begin{align}\hfill \hat{H}& =-J\sum\limits _{j}{\hat{b}}_{j+1}^{{\dagger}}{\hat{b}}_{j}\left[1+2g\left({e}^{\mp \frac{2\pi i}{3}}(1-{\hat{n}}_{j-1})+{e}^{\pm \frac{2\pi i}{3}}(1-{\hat{n}}_{j+2})\right)\right]+\text{h.c.}\hfill \\ \hfill & \quad -J\sum\limits _{j}{\hat{b}}_{j+2}^{{\dagger}}{\hat{b}}_{j}\left[1+2g\enspace {e}^{\pm \frac{4\pi i}{3}}(1-{\hat{n}}_{j+1})\right]+\text{h.c.}\hfill \\ \hfill & \quad -J\sum\limits _{j}{\hat{b}}_{j+3}^{{\dagger}}{\hat{b}}_{j}\enspace 2g\left[{e}^{\mp \frac{2\pi i}{3}}(1-{\hat{n}}_{j+1})+{e}^{\pm \frac{2\pi i}{3}}(1-{\hat{n}}_{j+2})\right]+\text{h.c.}-J\sum\limits _{j}{\hat{b}}_{j+4}^{{\dagger}}{\hat{b}}_{j}\enspace 2g(1-{\hat{n}}_{j+2})+\text{h.c.}\hfill \\ \hfill & \quad -J\sum\limits _{j}{\hat{b}}_{j}^{{\dagger}}{\hat{b}}_{j}\enspace 2g\eta \left[(1-{\hat{n}}_{j-1})+(1-{\hat{n}}_{j+1})+(1-{\hat{n}}_{j-2})+(1-{\hat{n}}_{j+2})\right].\hfill \end{align} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

The upper (lower) sign in the exponents applies to even (odd) sites j corresponding to the upper (lower) sub-chain. ${\hat{b}}_{j},{\hat{b}}_{j}^{{\dagger}}$ are annihilation and creation operators of hard-core bosons and ${\hat{n}}_{j}={\hat{b}}_{j}^{{\dagger}}{\hat{b}}_{j}$ is the number operator of particles acting at site j. The individual hopping processes are illustrated in figure 3. $g=27J/\left(2{\Delta}\right)$ is the coupling constant of the second order processes relative to the direct hopping. Its magnitude determines the relevance of both the complex, density-dependent hopping processes and the density–density interaction as compared to the direct hopping processes occurring at rate J.

The Hamiltonian (3) conserves the number of hard-core bosons (Rydberg excitations in state $\left\vert 1\right\rangle$), and of particular interest is the case of an excitation number commensurate with the lattice size L. Consequently, in the following we consider half filling, i.e. the total number of atoms in Rydberg states $\left\vert 0\right\rangle$ and $\left\vert 1\right\rangle$ is each one half the number of lattice sites. In [14] the authors explain how to prepare such a state experimentally for the case of a one-dimensional chain with alternating hopping. In the present case of a zig-zag chain half filling could alternatively be achieved e.g. by selectively applying external laser fields only to the lower sub-chain.

3.2. Phase diagram

In the present paper we are interested in the properties of the many-body ground state of the zig-zag chain (3) for a half-filled system of hard-core bosons with repulsive interactions, i.e. g > 0 and η ⩾ 0. We start by considering the two trivial limiting cases.

If all second-order processes are weak, i.e. for g → 0, which is the case for large detuning |Δ| ≫ J, only direct hopping processes between NNs j → j + 1 and j → j + 2 are relevant and one expects a simple bosonic (quasi) superfluid (phase I in figure 4). Controlled by the parameter η, there is a competition between the density–density interaction, preferring an ordered state for ηg ≫ 1 (phase IV in figure 4), and hopping processes which drive the system into a liquid state (phases I–III). For small values of η there is furthermore a competition between direct hopping processes and second-order hopping processes with density-dependent Peierls phases, controlled by g. Of particular interest are effects caused by the complex, density-dependent hopping terms.

Zoom In Zoom Out Reset image size

To attain a qualitative overview of the transition points between these phases we have calculated the dimensionless ground-state fidelity [26]

$$\begin{equation}f(\lambda )=\frac{2}{N}\frac{1-\left\vert \left\langle {{\Phi}}_{0}(\lambda )\right\vert \left\vert {{\Phi}}_{0}(\lambda +\delta \lambda )\right\rangle \right\vert }{\delta {\lambda }^{2}},\quad \delta \lambda \to 0,\end{equation} \tag{ 4 }$$

where |Φ₀(λ)⟩ is the many-body ground state of the system with N particles, and λ denotes a system parameter (g or η) that is varied, while keeping the other constant. In the thermodynamic limit L → ∞ while keeping L/N = 2 and within the same phase f is a smooth function of the system parameter(s) λ but develops a singularity at a phase transition. In a finite system the ground-state fidelity remains finite but shows a peak, which can still be used as a good indicator of a phase transition.

Figure 4 shows a qualitative phase diagram as a function of g and η. Shown are different first-order correlations in a false color plot. The crosses/plusses correspond to peaks of the ground-state fidelity f obtained from ED of Hamiltonian (3) in a finite system of length L = 20 with periodic boundary conditions (PBCs), whereas the inset shows the result of a finite size expansion. The latter demonstrates that the intermediate, liquid regime II survives in the thermodynamic limit. Since we are studying an experimentally motivated system, some comments about the achievable range of values of g are in order. For large g, the adiabatic elimination of state |+⟩ becomes incorrect and our model breaks down. For an estimate we use the values of detuning and hopping strength given in [18], where the minimal triangle system in figure 1(c) was studied. These values give g ≈ 0.46, which is very close to the transition point to regime III for η = 1 and an increase by a factor of 2 seems feasible. Thus all liquid phases at η = 1 are within reach of the experiment.

In order to characterize the different phases we discuss in the following density correlations and currents.

3.3. Density correlations

As can be seen from figure 4 the critical value of η needed for a transition into a phase with density order is larger than unity and on the order of 3. Thus further density–density interactions are needed in addition to the one associated with the second order hopping processes discussed in section 2.

In figure 5 we plotted the density–density correlation ${g}_{j}^{(2)}=\langle {\hat{n}}_{0}{\hat{n}}_{j}\rangle -\langle {\hat{n}}_{0}\rangle \langle {\hat{n}}_{j}\rangle$ for different values of η at g = 2. One recognizes the emergence of long-range density order in both sub-chains when η > 3. In an infinite system or for PBCs and a number of lattice sites that is a multiple of 4, the ordered ground state is four-fold degenerate:

$$\begin{align}\hfill \vert {\psi }_{0}^{(1)}\rangle & =\vert 1\rangle \vert 1\rangle \vert 0\rangle \vert 0\rangle \vert 1\rangle \vert 1\rangle \dots ,\hfill \\ \hfill \vert {\psi }_{0}^{(2)}\rangle & =\vert 0\rangle \vert 1\rangle \vert 1\rangle \vert 0\rangle \vert 0\rangle \vert 1\rangle \dots ,\hfill \\ \hfill \vert {\psi }_{0}^{(3)}\rangle & =\vert 0\rangle \vert 0\rangle \vert 1\rangle \vert 1\rangle \vert 0\rangle \vert 0\rangle \dots ,\hfill \\ \hfill \vert {\psi }_{0}^{(4)}\rangle & =\vert 1\rangle \vert 0\rangle \vert 0\rangle \vert 1\rangle \vert 1\rangle \vert 0\rangle \dots .\hfill \end{align} \tag{ 5 }$$

In these configurations every excited site has exactly one excited neighbor, minimizing the energy from density–density interactions. In any other configuration the number of excited neighbors is increased, leading to a larger interaction energy per particle. Thus ${g}_{j}^{(2)}$ oscillates with a period of 4.

Zoom In Zoom Out Reset image size

3.4. Liquid phases and chiral average currents

More interesting than the liquid–insulator transition from phase III into phase IV seen in figure 5 are those within the liquid phases. Figure 6 shows a plot of f along a vertical cut of figure 4 for η = 1. One recognizes two further peaks, one at g ≈ 0.5 and another one at g ≈ 0.38. For small g (phase I) the system is a simple superfluid. This phase is well characterized by an approximate Hamiltonian (3) where a mean-field decoupling is used in all hopping terms

$$\begin{align}\hfill {\hat{H}}_{\text{MF}}^{\text{B}}& =-J\sum\limits _{j}\left[{\hat{b}}_{j+1}^{{\dagger}}{\hat{b}}_{j}(1-g)+{\hat{b}}_{j+2}^{{\dagger}}{\hat{b}}_{j}r\enspace {e}^{\pm i\phi (g)}-{\hat{b}}_{j+3}^{{\dagger}}{\hat{b}}_{j}g\right]\enspace +\text{h.c.}\hfill \\ \hfill & \quad -2Jg\eta \sum\limits _{\langle j,l\rangle }{\hat{n}}_{j}(1-{\hat{n}}_{l}),\hfill \end{align} \tag{ 6 }$$

where $r{e}^{\pm i\phi (g)}=1+g{e}^{\pm \frac{4\pi i}{3}}$ describes an effective magnetic field and ⟨j, l⟩ denotes all NNs. For such a system one expects a Meissner-like phase for small magnetic fields with anti-parallel screening currents in the upper and lower sub-chain as predicted in [21] and first observed in [22].

Zoom In Zoom Out Reset image size

Since the system is particle-number conserving, currents can be derived from the continuity equation ${\partial }_{t}{\hat{n}}_{j}={\sum }_{i}\enspace {I}_{i\to j}-{I}_{j\to i}$ which contains a sum of incoming and outgoing terms. In the full model the current along the chain thus reads

$$\begin{align}\hfill {I}_{j\to j+2}& =2J\enspace \mathrm{I}\mathrm{m}\langle {\hat{b}}_{j}^{{\dagger}}{\hat{b}}_{j+2}\rangle -Jg\enspace \mathrm{I}\mathrm{m}\langle {\hat{b}}_{j}^{{\dagger}}(1-{\hat{n}}_{j+1}){\hat{b}}_{j+2}\rangle \hfill \\ \hfill & \quad \pm \sqrt{3}Jg\enspace \mathrm{R}\mathrm{e}\langle {\hat{b}}_{j}^{{\dagger}}(1-{\hat{n}}_{j+1}){\hat{b}}_{j+2}\rangle ,\hfill \end{align} \tag{ 7 }$$

where the upper (lower) sign corresponds to even (odd) j for the upper (lower) chain. Note that the appearance of both imaginary and real parts of the correlations is due to the presence of the Peierls phase. At the same time the average current between the two chains must vanish for PBCs.

In figure 6(b) we have plotted I_j→j+2 for the upper sub-chain as a function of g and for η = 1. The current along the lower sub-chain is exactly opposite. For small g there is a linear growth which is fully consistent with the approximate model, equation (6), in which the magnetic flux per elementary triangle in units of the flux quantum, Φ/Φ₀ = ϕ(g)/2π, increases approximately linear with g.

For larger values of g two new phases (II and III) with distinctly different current patterns emerge. One notices two interesting features:

• (a)  
  The direction of the chiral currents reverses when g becomes larger than unity.
• (b)  
  The magnitude of the current is substantially amplified in phase II as compared to phases I and III, with a discontinuous transition.

Furthermore, when η is increased the system crosses from the liquid regime into a regime with density order. Consequently, the currents are suppressed with a power law $\sim {\eta }^{-1}$ and disappear in the limit of large η (see appendix B).

In order to get a qualitative understanding of the emergence of different phases in the liquid regime we first discuss a mean-field approximation accounting for the hard-core character of the bosons. To this end we express the bosonic Hamiltonian in terms of fermions using a Jordan–Wigner transformation [27] along a zig-zag string (see figure 3) via

$$\begin{equation}{\hat{c}}_{i}={\text{e}}^{\text{i}\pi {\sum }_{l< i}\enspace {\hat{n}}_{l}}{\hat{b}}_{i}.\end{equation} \tag{ 8 }$$

Here the ${\hat{b}}_{i}$ are the hard-core boson operators from before and the ${\hat{c}}_{i}$ are fermionic ones. The transformed Hamiltonian then reads

$$\begin{align}\hfill \hat{H}& =\left\{-J\sum\limits _{j}\;{\hat{c}}_{j+1}^{{\dagger}}{\hat{c}}_{j}+{\hat{c}}_{j+2}^{{\dagger}}{\hat{c}}_{j}(1-2{\hat{n}}_{j+1})\right.-2gJ\sum\limits _{j}{\hat{c}}_{j+2}^{{\dagger}}{\hat{c}}_{j}\left({e}^{\mp \frac{2\pi i}{3}}(1-{\hat{n}}_{j-1})+{e}^{\pm \frac{2\pi i}{3}}(1-{\hat{n}}_{j+2})\right)\hfill \\ \hfill & \quad -2gJ\sum\limits _{j}\;{\hat{c}}_{j+2}^{{\dagger}}{\hat{c}}_{j}\left[(1-2{\hat{n}}_{j+1})+2g\enspace {e}^{\pm \frac{4\pi i}{3}}(1-{\hat{n}}_{j+1})\right]\hfill \\ \hfill & \quad -2gJ\sum\limits _{j}\;{\hat{c}}_{j+3}^{{\dagger}}{\hat{c}}_{j}\left[{e}^{\mp \frac{2\pi i}{3}}(1-2{\hat{n}}_{j+2})(1-{\hat{n}}_{j+1})+{e}^{\pm \frac{2\pi i}{3}}(1-2{\hat{n}}_{j+1})(1-{\hat{n}}_{j+2})\right]\hfill \\ \hfill & \left.\quad -2gJ\sum\limits _{j}{\hat{c}}_{j+4}^{{\dagger}}{\hat{c}}_{j}(1-2{\hat{n}}_{j+1})(1-{\hat{n}}_{j+2})(1-2{\hat{n}}_{j+3})+\text{h.c.}\right\}-2\eta gJ\sum\limits _{\langle \;j,l\rangle }{\hat{n}}_{j}(1-{\hat{n}}_{l}),\hfill \end{align} \tag{ 9 }$$

where we have used the identity ${e}^{\pm i\pi {\hat{n}}_{j}}=\left(1-2{\hat{n}}_{j}\right)$. Now we apply a mean-field approximation, replacing the particle number operators in the higher-order terms by expectation values, assuming a uniform half filling $\left\langle {\hat{n}}_{j}\right\rangle =0.5$. The corresponding mean-field Hamiltonian reads

$$\begin{equation}{\hat{H}}_{\text{MF}}=-J\sum\limits _{j}{\hat{c}}_{j+1}^{{\dagger}}{\hat{c}}_{j}(1-g)-gJ\sum\limits _{j}{\hat{c}}_{j+2}^{{\dagger}}{\hat{c}}_{j}\enspace {e}^{\pm \frac{4\pi }{3}i}+\text{h.c.}-2\eta gJ\sum\limits _{j}{\hat{n}}_{j}.\end{equation} \tag{ 10 }$$

One recognizes that the hoppings along the upper and lower sub-chain are complex and contain fixed Peierls phases ±4π/3 of opposite sign, which cannot be gauged away. These phases correspond to a homogeneous, constant effective gauge field for the fermions, which leads to chiral currents along the upper and lower chains. Furthermore, the direct NN hopping j → j + 1 changes its sign when g becomes larger than unity. We may write ${\hat{H}}_{\text{MF}}$ in momentum space

$$\begin{equation}{\hat{H}}_{\text{MF}}=-J\sum\limits _{k}{\left(\begin{matrix}\hfill {\hat{c}}_{Ak}^{{\dagger}} \hfill \\ \hfill {\hat{c}}_{Bk}^{{\dagger}}\enspace \hfill \end{matrix}\right)}^{ \top } \left(\begin{matrix}\hfill 2g\enspace \mathrm{sin}\left(k-\frac{\pi }{6}\right)\hfill & \hfill (1-g)(1+{e}^{-ik})\hfill \\ \hfill (1-g)(1+{e}^{ik})\hfill & \hfill -2g\enspace \mathrm{sin}\left(k-\frac{\pi }{6}\right)\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\hat{c}}_{Ak} \hfill \\ \hfill {\hat{c}}_{Bk}\enspace \hfill \end{matrix}\right),\end{equation} \tag{ 11 }$$

where $k\in \left[-\pi ,\pi \right]$ is the lattice momentum and the indices A and B denote the upper and lower chain, respectively. The corresponding spectrum is shown in figure 7 where one observes that there is no gap opening between the two bands and thus in mean-field approximation the system is not insulating away from unit filling. For g < 0.5 only momentum states at the edge of the Brillouin zone touch the Fermi energy. For g ⩾ 0.5, where the system enters phase III, low-momentum states cross the Fermi energy, changing the properties of the many-body ground state.

Zoom In Zoom Out Reset image size

The nature of the intermediate phase II, in between the superfluid phase I, g → 0, and phase III, where second-order hopping processes are dominant, g ⩾ 0.5, cannot be explained within the simple mean-field approximation and requires to consider the interactions included in the full Hamiltonian (3).

We have seen in section 3.3 that the NN density–density interaction drives the system into a density wave of period two in both sub-chains. In the presence of a weak density wave the mean-field Hamiltonian gets perturbed. The density-dependent hopping terms are then modulated corresponding to a unit cell of four lattice sites. The impact of such a modulation can be understood when folding the mean-field spectrum of figure 7 into a reduced Brillouin zone from −π/2 to π/2. This is illustrated in figure 8 for values of g between 0.3 and 0.45. One recognizes that for g ⩾ 0.36 there is a band crossing of the lower sub-bands, which modifies the nature of the many-body ground state when the density-modulation is taken into account. The emergence of phase II can be attributed to this.

Zoom In Zoom Out Reset image size

From this mean-field description one naively expects that phase II disappears in the absence of density–density interactions, i.e. for η = 0. Quite surprisingly this is not the case. When calculating the ground-state fidelity of the full model (3) for η = 0 we still find two peaks at g ≈ 0.5 and g ≈ 0.375. This must then be attributed to the presence of the density-dependent complex hopping terms alone. We will return to this point at a later point section 5.2.

To further understand the impact of the density-dependent hopping terms included in the Rydberg-spin Hamiltonian (3) we once again consider the currents between NN sites between and along the chains. To simplify the methodology we use open boundary conditions (OBCs) to select a fixed realization of density-order. In figure 9 we show the corresponding current for OBC spatially resolved, where we have used equation (7) and

$$\begin{align}\hfill {I}_{j\to j+1}& =2J\enspace \mathrm{I}\mathrm{m}\langle {\hat{b}}_{j}^{{\dagger}}{\hat{b}}_{j+1}\rangle -Jg\enspace \mathrm{I}\mathrm{m}\langle {\hat{b}}_{j}^{{\dagger}}(2-{\hat{n}}_{j+2}-{\hat{n}}_{j-1}){\hat{b}}_{j+1}\rangle \hfill \\ \hfill & \quad \pm \sqrt{3}Jg\enspace \mathrm{R}\mathrm{e}\langle {\hat{b}}_{j}^{{\dagger}}({\hat{n}}_{j+2}-{\hat{n}}_{j-1}){\hat{b}}_{j+1}\rangle .\hfill \end{align} \tag{ 12 }$$

In the limit g → 0 all currents naturally vanish, as does the effect of the gauge field (left-most column of figure 9). In the middle column we see that for 2g ≈ 1 the currents connecting the sub-chains remain small, whereas the current along the sub-chains shows the behavior, characteristic for a Meissner phase, inasmuch as we find a homogeneous current moving right in sub-chain A (j even) and moving left in sub-chain B (j odd). However, as soon as the density-dependent hopping become dominant over the direct processes (g ≫ 1) we observe that both currents feature a periodicity of four sites. We can draw these currents as arrows between the sites (bottom row of figure 9) to see that the current flows around plaquettes of four sites in alternating manner forming a regular array of vortices. For hardcore bosons hopping on a square ladder in an external effective magnetic field, a transition from a Meissner phase into a phase with current vortices when increasing the magnetic flux has been predicted in [21] and was observed in [22]. Thus the emergence of vortices in figure 9 seems at first glance consistent with these predictions based on the approximate boson Hamiltonian, equation (6). However, one recognises that for larger values of g the current vortices have the opposite direction as expected from equation (6), which is an effect of the density-dependence of the Peierls phases.

Zoom In Zoom Out Reset image size

To understand this effect, we rewrite the fermion Hamiltonian (9) in the following form

$$\begin{align}\hfill \hat{H}& =-J\sum\limits _{j}\left[{\hat{c}}_{j+1}^{{\dagger}}{\hat{c}}_{j}+{\hat{c}}_{j+2}^{{\dagger}}{\hat{c}}_{j}(1-2{\hat{n}}_{j+1})\right]+\text{h.c.}\hfill \\ \hfill & \quad -2gJ\enspace \sum\limits _{j}\left[{\hat{c}}_{j+1}^{{\dagger}}{\hat{c}}_{j}\enspace {\hat{U}}_{j+1,j}\left(1-{\hat{n}}_{j-1}{\hat{n}}_{j+2}\right)+{\hat{c}}_{j+2}^{{\dagger}}{\hat{c}}_{j}\enspace {U}_{j+2,j}(1-{\hat{n}}_{j+1})\right.\hfill \\ \hfill & \left.\quad +{\hat{c}}_{j+3}^{{\dagger}}{\hat{c}}_{j}\enspace {\hat{U}}_{j+3,j}\left(1-{\hat{n}}_{j+1}{\hat{n}}_{j+2}\right)-{\hat{c}}_{j+4}^{{\dagger}}{\hat{c}}_{j}\enspace (1-2{\hat{n}}_{j+1})(1-{\hat{n}}_{j+2})(1-2{\hat{n}}_{j+3})\right]+\text{h.c.}\hfill \\ \hfill & \quad -2\eta gJ\sum\limits _{j}{\hat{c}}_{j}^{{\dagger}}{\hat{c}}_{j}\enspace \left[(1-{\hat{n}}_{j-1})+(1-{\hat{n}}_{j+1})+(1-{\hat{n}}_{j-2})+(1-{\hat{n}}_{j+2})\right].\hfill \end{align} \tag{ 13 }$$

Here we have introduced unitary link operators

$$\begin{align}\hfill {\hat{U}}_{j+1,j}& =\mathrm{exp}\left(i{\hat{\phi }}_{j+1,j}\right)=-\mathrm{exp}\left(\pm \frac{\pi i}{3}({\hat{n}}_{j+2}-{\hat{n}}_{j-1})\right),\hfill \\ \hfill {\hat{U}}_{j+2,j}& =\mathrm{exp}\left(i{\hat{\phi }}_{j+2,j}\right)=\mathrm{exp}\left(\mp \frac{\pi i}{3}\right),\hfill \\ \hfill {\hat{U}}_{j+3,j}& =\mathrm{exp}\left(i{\hat{\phi }}_{j+3,j}\right)=-\mathrm{exp}\left(\pm \frac{2\pi i}{3}({\hat{n}}_{j+1}-{\hat{n}}_{j+2})\right),\hfill \end{align} \tag{ 14 }$$

where the upper (lower) sign corresponds to even (odd) indices j. Thus besides the direct hopping amplitudes at rate J there are complex hopping terms proportional to Jg corresponding to an effective gauge field. In contrast to lattice models describing particles coupled to a fixed background field, such as the Harper–Hofstadter model, the effective gauge fields are here operator valued. Only ${\hat{U}}_{j+2,j}$ corresponds to a fixed, classical gauge field.

Despite the fact that the gauge fields causing the operator-valued Peierls phases are dynamical quantities, there is no additional conserved charge associated with them and the gauge freedom is the same as for a classical background field. In the case of a classical background field, the unitary link operators are just exponential phase factors

$$\begin{equation}{U}_{lj}=\mathrm{exp}\left\{i{\int }_{j}^{l} \mathrm{d}\boldsymbol{r}\cdot \boldsymbol{A}\right\},\end{equation} \tag{ 15 }$$

where A is the external, classical vector potential. Under a $\mathsf{U}(1)$ gauge transformation of the classical field, the link operators transform as

$$\begin{equation}{U}_{lj}\enspace \to \enspace {U}_{lj}^{\prime }=\mathrm{exp}\left\{i{\int }_{j}^{l} \mathrm{d}\boldsymbol{r}\cdot (\boldsymbol{A}+\nabla \alpha )\right\}=\mathrm{exp}\left\{i{\int }_{j}^{l} \mathrm{d}\boldsymbol{r}\cdot \boldsymbol{A}+i\left({\alpha }_{l}-{\alpha }_{j}\right)\right\},\end{equation} \tag{ 16 }$$

which can be compensated by a unitary transformation of the fermion operators, ${\text{e}}^{-\text{i}{\alpha }_{j}{\hat{n}}_{j}}\enspace {\hat{c}}_{j}\enspace {\text{e}}^{\text{i}{\alpha }_{j}{\hat{n}}_{j}}={\hat{c}}_{j}\enspace {\text{e}}^{\text{i}{\alpha }_{j}}$, such that the Hamiltonian is invariant under the total transformation. We should remember that only the longitudinal component of the vector potential $\boldsymbol{A}={\boldsymbol{A}}_{\left.{\Vert}\right.}+{\boldsymbol{A}}_{\perp }$ changes under a gauge transformation, while the transversal part A_⊥ remains the same. Upon inspecting equation (14) one recognizes that the operator-valued contributions $\hat{\boldsymbol{A}}$ to the gauge field are transverse:

$$\begin{equation}{\hat{\phi }}_{l,j}={\int }_{j}^{l} \mathrm{d}\boldsymbol{r}\cdot (\boldsymbol{A}+\hat{\boldsymbol{A}})={\int }_{j}^{l} \mathrm{d}\boldsymbol{r}\cdot ({\boldsymbol{A}}_{\left.{\Vert}\right.}+{\hat{\boldsymbol{A}}}_{\perp }).\end{equation} \tag{ 17 }$$

To see this we note that the operator terms in equation (14) are given by differences of occupation numbers at lattice sites left and right of the hopping path, which is the lattice equivalent of

$$\begin{equation}\hat{\boldsymbol{A}}\sim \nabla \times \hat{n}(x,y){\boldsymbol{e}}_{z}=\left(\frac{\partial \hat{n}}{\partial y},-\frac{\partial \hat{n}}{\partial x},0\right),\quad \mathrm{a}\mathrm{n}\mathrm{d}\enspace \mathrm{t}\mathrm{h}\mathrm{u}\mathrm{s}\quad \nabla \cdot \hat{\boldsymbol{A}}=0.\end{equation} \tag{ 18 }$$

The Peierls phases ${\hat{\phi }}_{l,j}$ themselves have no direct physical relevance as a classical gauge transformation of the longitudinal component, equation (16), leads to ${\hat{\phi }}_{l,j}\to {\hat{\phi }}_{l,j}+{\alpha }_{l}-{\alpha }_{j}$. Quantities that are invariant under these gauge transformations are plaquette fluxes, e.g. in a rhombus of lattice sites j, j + 1, j + 2, j + 4 consisting of two adjacent triangles, ${\hat{{\Theta}}}_{j}^{\square }$ (mod 2π). ${\hat{{\Theta}}}_{j}^{\square }$ contains a fixed classical component ${{\Phi}}_{c}^{\square }=2\pi /3$ and a quantum component ${\hat{{\Phi}}}_{j}^{\square }$

$$\begin{align}\hfill {\hat{{\Theta}}}_{j}^{\square }& =\pm \left({\hat{\phi }}_{j+2,j}+{\hat{\phi }}_{j+3,j+2}-{\hat{\phi }}_{j+3,j+1}-{\hat{\phi }}_{j+1,j}\right),\hfill \\ \hfill & =-\frac{\pi }{3}\left({\hat{n}}_{j+2}+{\hat{n}}_{j+1}-{\hat{n}}_{j-1}-{\hat{n}}_{j+4}\right)+\frac{2\pi }{3}={\hat{{\Phi}}}_{j}^{\square }+{{\Phi}}_{c}^{\square },\hfill \end{align} \tag{ 19 }$$

where the upper (lower) sign correspond to even (odd) sites j. Note that we have defined the sign of the fluxes by a clockwise sequence of sites in the +z direction in figure 3, i.e. for even j: j → j + 2 → j + 3 → j + 1 → j and for odd j: j → j + 1 → j + 3 → j + 2 → j. In the absence of density modulations the quantum gauge fields ${\hat{{\Phi}}}_{j}^{\square }$ have all a vanishing magnitude and one recovers the mean-field behavior. This changes however if density correlations are taken into account. Then the plaquette fluxes are modulated and they induce currents between the two chains I_j→j+1 and modulate the currents I_j→j+2 along the chains, as seen in figure 9.

5.1. Dynamical gauge fields driven by density–density interactions

As shown in figure 5, density–density interactions drive the many-body ground state toward an ordered state with a density wave for large values of η. For the ground states (5) one would find e.g. in state |ψ₀⟩ = ...|0⟩₀|1⟩₁|1⟩₂|0⟩₃|0⟩₄|1⟩₅...

$$\begin{equation*}\langle {\hat{{\Phi}}}_{0}^{\square }\rangle =\frac{2\pi }{3},\hspace{25.0pt}\quad \langle {\hat{{\Phi}}}_{1}^{\square }\rangle =0,\hspace{25.0pt}\quad \langle {\hat{{\Phi}}}_{2}^{\square }\rangle =-\frac{2\pi }{3},\hspace{25.0pt}\quad \langle {\hat{{\Phi}}}_{3}^{\square }\rangle =0,\quad \mathrm{e}\mathrm{t}\mathrm{c}.\end{equation*}$$

Thus the total flux of the gauge field oscillates around the classical value ${{\Phi}}_{c}^{\square }=2\pi /3$.

To describe the transition to the flux-ordered phase in a system with OBCs, we introduce the order parameter

$$\begin{equation}\chi =\left\vert \frac{1}{L}\sum\limits _{j=2n}{\left(-1\right)}^{n}\left(\left\langle {\hat{{\Phi}}}_{j}^{\square }\right\rangle +\left\langle {\hat{{\Phi}}}_{j+1}^{\square }\right\rangle \right)\right\vert .\end{equation} \tag{ 20 }$$

In this expression, we exclude the first and last plaquette at the edges of the system to avoid boundary effects. In figure 10(a) we show χ as function of g in the regime of large η. We recognize a transition into a phase with long-range order of the flux lattice. The origin of the regular pattern of staggered effective fluxes is simply the density order generated by the repulsive interaction proportional to η. Since in the ordered phase (IV) first-order correlations are not strictly zero but are suppressed proportional to η⁻¹ (see appendix), the flux pattern induces a lattice of current vortices similar to those in figure 9. In figure 10(b) we have plotted the local current between the two sub-chains inside the ordered phase IV for OBCs.

Zoom In Zoom Out Reset image size

5.2. Spontaneous gauge field generation

When η is reduced and below the critical value of the transition to the ordered phase, density correlations show a power-law decay as function of distance. In this case the order parameter χ vanishes in the thermodynamic limit, but the local flux ${\hat{{\Phi}}}_{j}^{\square }$ shows non-vanishing fluctuations, i.e. $\left\langle {\left({\hat{{\Phi}}}_{j}^{\square }\right)}^{2}\right\rangle \ne 0$ and liquid-like spatial correlations.

However, when the density–density interaction is switched off entirely, i.e. for η = 0, one naively expects no density modulations in the absence of density-dependent complex hoppings. In fact due to the triangular lattice structure Friedel oscillations do lead to density modulations in a finite system with OBCs even in the mean-field limit. Only at g = 1, where the mean-field Hamiltonian (10) separates into two independent chains, the corresponding density distribution is flat, see figure 11(a), and we expect the quantum gauge field to be in a vacuum state. We note that the classical Peierls phases can be gauged away in this case. The full model does however show density modulations, see figure 11(a), which are entirely caused by the density-dependent complex hopping terms. Thus there is a spontaneous generation of the gauge field with non-vanishing amplitude of $\langle \left({\hat{{\Phi}}}_{j}^{\square }\right)\rangle$ and fluxes that form an oscillating pattern for OBCs. The oscillations decay away from the boundary, however, showing that the flux correlations are liquid-like and there is no long-range order, in contrast to figure 10(a) for η = 10.

Zoom In Zoom Out Reset image size

5.3. Gutzwiller ansatz

The spontaneous generation of a non-zero expectation value of the staggered gauge field results from the competition between an increase of kinetic energy of direct hopping processes on the one hand and a lowering of energy by increase of the imaginary second-order hopping amplitudes on the other. To see this, we performed a variational approach to determine the ground state using a Gutzwiller ansatz. Guided by the numerical results we assume a unit cell of four sites and write the following ansatz wavefunction, illustrated in figure 12

$$\begin{equation}\left\vert \psi \right\rangle =\frac{1}{\mathcal{N}}\prod\limits _{j=4n}\left\vert {\psi }_{j}\right\rangle \left\vert {\psi }_{j+1}\right\rangle \left\vert {\psi }_{j+2}\right\rangle \left\vert {\psi }_{j+3}\right\rangle ,\end{equation} \tag{ 21 }$$

$$\begin{equation}\left\vert {\psi }_{j}\right\rangle =(1-\varepsilon )\left\vert 0\right\rangle +(1+\varepsilon ){\mathrm{e}}^{\text{i}(\theta -\phi )}\left\vert 1\right\rangle ,\end{equation} \tag{ 22 }$$

$$\begin{equation}\left\vert {\psi }_{j+1}\right\rangle =(1+\varepsilon )\left\vert 0\right\rangle +(1-\varepsilon ){\mathrm{e}}^{-\mathrm{i}(\theta -\phi )}\left\vert 1\right\rangle ,\end{equation} \tag{ 23 }$$

$$\begin{equation}\left\vert {\psi }_{j+2}\right\rangle =(1+\varepsilon )\left\vert 0\right\rangle +(1-\varepsilon ){\mathrm{e}}^{\text{i}(\theta +\phi )}\left\vert 1\right\rangle ,\end{equation} \tag{ 24 }$$

$$\begin{equation}\left\vert {\psi }_{j+3}\right\rangle =(1-\varepsilon )\left\vert 0\right\rangle +(1-\varepsilon ){\mathrm{e}}^{-\mathrm{i}(\theta +\phi )}\left\vert 1\right\rangle .\end{equation} \tag{ 25 }$$

Here we use three parameters to accommodate the physics of the model. The parameter ɛ accounts for the observed density modulation, θ distinguishes the two sub-chains and ϕ accounts for the complex phase that is picked up when a particle hops along the same sub-chain (to the right for the upper chain, to the left for the lower one).

Zoom In Zoom Out Reset image size

Using a standard optimization algorithm to minimize the variational energy, we obtain the following results for the three parameters (up to numerical precision) for L = 16, η = 0 and g = 1:

$$\begin{equation}{\varepsilon }_{\text{opt}}\approx 0.134\hspace{25.0pt}{\theta }_{\text{opt}}\approx \frac{\pi }{2}\hspace{25.0pt}{\phi }_{\text{opt}}\approx 0.317.\end{equation} \tag{ 26 }$$

For the case of very small density-dependent hopping, g → 0, all three optimal parameter values are found to be zero as in this case only direct hopping processes survive. In figure 13 we show the expectation value of the ground-state energy in the Gutzwiller state, where we have fixed the third parameter θ to its optimal values, θ = 0 (g = 0.001) and θ = π/2 (g = 1). For small density-dependent hopping the alternating currents quantified by χ vanish. For g = 1 however there exist two degenerate minimum energy configurations point-symmetric to zero which represent the two different chain configurations and feature non-vanishing gauge fields, as indicated by the value of χ for the horizontal cut in figure 13.

Zoom In Zoom Out Reset image size

We note that the Gutzwiller ansatz cannot capture a liquid-type phase with short-range correlations of the gauge field and thus a finite value of the order parameter χ is obtained. It is instructive to consider the optimum value of the density imbalance, characterized by the variational parameter ɛ as function of the strength g of second-order hopping processes. This is shown in figure 14. One recognizes that within the Gutzwiller approach there is a phase transition at g = 0.5 from a phase with vanishing gauge field to a phase with a finite amplitude. Below the critical value of g the energy increase due to a density wave outweighs the energy decrease due to the gauge-field mediated hopping. Above the critical value the situation is reversed.

Zoom In Zoom Out Reset image size

While vortex currents of Rydberg excitations are difficult to measure in experiments, first order correlations can be accessed. Applying a global microwave pulse to all Rydberg atoms allows to rotate all spins formed by states |0⟩ and |1⟩, see figure 1(c), from the x- or y-direction on the Bloch-sphere onto the z axis. A subsequent spatially resolved measurement of correlations of Rydberg excitations thus yields xx or yy spin correlations, from which one can obtain the real part of first-order boson correlations

$$\begin{equation}\mathrm{R}\mathrm{e}\left[\langle {\hat{b}}_{k}^{{\dagger}}{\hat{b}}_{l}\rangle \right]=\langle {\hat{\sigma }}_{k}^{x}{\hat{\sigma }}_{l}^{x}\rangle +\langle {\hat{\sigma }}_{k}^{y}{\hat{\sigma }}_{l}^{y}\rangle .\end{equation} \tag{ 27 }$$

Using laser beams addressing specific atoms and a global microwave pulse, it is possible to apply spin rotations around different axes to different Rydberg atoms. In this way anisotropic spin correlations can be measured which give access also to the imaginary part of first-order correlations between bosons at different lattice sites

$$\begin{equation}\mathrm{I}\mathrm{m}\left[\langle {\hat{b}}_{k}^{{\dagger}}{\hat{b}}_{l}\rangle \right]=\langle {\hat{\sigma }}_{k}^{y}{\hat{\sigma }}_{l}^{x}\rangle -\langle {\hat{\sigma }}_{k}^{x}{\hat{\sigma }}_{l}^{y}\rangle .\end{equation} \tag{ 28 }$$

As shown in the phase diagram figure 4 the non-trivial liquid phases II and III can be distinguished from the trivial phase I by first-order correlations of neighboring lattice sites between the sub-chains, $\langle {\hat{b}}_{j}^{{\dagger}}{\hat{b}}_{j+1}\rangle$, and along the sub-chains, $\langle {\hat{b}}_{j}^{{\dagger}}{\hat{b}}_{j+2}\rangle$. In phase I the real part of both is positive and there is no imaginary component, indicating the absence of currents, see top row of figure 15. In phase II the imaginary parts of both correlations feature an oscillating pattern, while the real parts remain positive with reduced amplitudes, see middle row of figure 15. The nearly homogeneous oscillation in the right column is a signature of the constant chiral currents along the upper and lower sub-chain. Finally in phase III the real part of the inter-chain correlation $\langle {\hat{b}}_{j}^{{\dagger}}{\hat{b}}_{j+1}\rangle$ flips its sign and the imaginary parts attain a more complex pattern reminiscent of vortex currents. To summarize, using the technique of global (and local) microwave pulses and subsequent space-resolved measurement of Rydberg excitations allows to experimentally confirm the non-trivial liquid regimes of the model.

Zoom In Zoom Out Reset image size

We studied the effects of a density-dependent, complex hopping of Rydberg excitations arising from second-order processes to the many-body ground state of a one-dimensional zig-zag-chain at half filling. The strength of the second-order hopping can be controlled by microscopic parameters of the Rydberg atoms and the system can be tuned from a regime where these processes are negligible to one where they dominate. Second order processes also lead to density–density interactions, which are modified by additional off-resonant van-der-Waals couplings between Rydberg atoms. The competition between direct and second-order hopping processes as well as density–density repulsion leads to a rich ground-state phase diagram. Most importantly the complex valuedness of the second-order hopping amplitudes breaks time-reversal symmetry and, as shown in the present paper, leads to two types of effective gauge fields. The first one is equivalent to a classical vector potential. It is worth noting that this field emerges directly from the microscopic physics of Rydberg atoms and is very different from other mechanisms that have been suggested for effective gauge fields for neutral atoms. This classical gauge field can be employed for example to realize topological band structures. Even more interesting is the second kind of gauge field which arises from the density-dependence of the second-order hopping. It is a quantum field linked to density modulations. Thus we here encounter a dynamical quantum gauge field. While based on the mechanism discussed here it is possible to implement true lattice gauge theories, the operator-valued field in the present context is a transversal vector potential. It is thus gauge invariant and does not lead to additional conserved local charges. Depending on the values of second order hopping and density–density interaction we identified four distinct phases. An insulating phase with long-range density order in a repulsion-dominated regime, a superfluid for weak interactions, and two nontrivial liquid phases where density repulsion is weak but second order hopping dominates over direct hopping. In the insulating phase long-range density correlations induce an oscillatory lattice of fluxes of the quantum gauge field, which in turn leads to a regular pattern of local vortex currents. Surprisingly, even in the complete absence of density–density interactions the quantum gauge field attains a non-vanishing value with liquid-like correlations. We must interpret this as a spontaneous creation of the gauge field by breaking of the lattice translational symmetry.

Our findings demonstrate that systems of Rydberg atoms can give rise to a broad spectrum of nontrivial interaction processes and offer new and experimentally accessible approaches to study dynamical gauge fields. Additionally, even topologically non-trivial states are realizable using this Rydberg system, as we demonstrate in a Haldane-type model on a honeycomb lattice [28].

SO, MK and MF gratefully acknowledge financial support from the DFG through SFB TR 185, Project Number 277625399. AB acknowledges financial support from the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 817482 (PASQuanS). AB and HPB acknowledge financial support from the French-German collaboration for joint projects in NLE Sciences funded by the Deutsche Forschungsgemeinschaft (DFG) and the Agence National de la Recherche (ANR PRCI, project RYBOTIN). The authors would like to thank David Petrosyan for very insightful discussions.

The data that support the findings of this study are available upon reasonable request from the authors.

In this section we provide a microscopic derivation of Hamiltonian (2). In this derivation, we assume that the dynamics of the atoms can be restricted to the three atomic states shown in figure 1(c), namely $\left\vert 1\right\rangle =\left\vert 60{P}_{3/2},{m}_{J}=-1/2\right\rangle$, $\left\vert +\right\rangle =\left\vert 60{P}_{3/2},{m}_{J}=+3/2\right\rangle$ as well as $\left\vert 0\right\rangle =\left\vert 60{S}_{1/2},{m}_{J}=+1/2\right\rangle$. The validity of this reduction to three states has been demonstrated in [18], we shortly summarize the arguments here: The states of the P_1/2 manifold can be neglected, as spontaneous decay into these states is negligible and they are not populated during preparation. If initially only the states $\left\vert 1\right\rangle$ and $\left\vert 0\right\rangle$ are populated, the only possible exchange process that couples out of the three-state manifold is the process $\left\vert 1\right\rangle \to \left\vert 60{S}_{1/2},{m}_{J}=-1/2\right\rangle$ and simultaneously $\left\vert 0\right\rangle \to \left\vert 60{P}_{3/2},{m}_{J}=+1/2\right\rangle$ in another atom. As mentioned in the appendix of [18], this transition is off-resonant and energetically suppressed.

Using this three-state approximation we derive the Hamiltonian (2) starting from the static dipole–dipole interaction between two atoms

$$\begin{equation}{\hat{V}}_{ij}=\frac{1}{4\pi {{\epsilon}}_{0}{R}_{ij}^{3}}\left[{\hat{d}}_{i}^{z}{\hat{d}}_{j}^{z}+\frac{1}{2}\left({\hat{d}}_{i}^{+}{\hat{d}}_{j}^{-}+{\hat{d}}_{i}^{-}{\hat{d}}_{j}^{+}\right)-\frac{3}{2}\left({\hat{d}}_{i}^{+}{\hat{d}}_{j}^{+}{e}^{-2i{\phi }_{ij}}+{\hat{d}}_{i}^{-}{\hat{d}}_{j}^{-}{e}^{2i{\phi }_{ij}}\right)\right],\end{equation} \tag{ A.1 }$$

where ${R}_{i,j}=\left\vert {\boldsymbol{r}}_{i,j}\right\vert$ and ϕ_ij refer to the interatomic distance and angle in polar coordinates, respectively. Additionally, the operators ${\hat{d}}_{i}^{\pm }=\mp \frac{1}{\sqrt{2}}({\hat{d}}_{i}^{x}\pm i{\hat{d}}_{i}^{y})$ represent ladder operators raising and lowering the angular momentum of the ith Rydberg atom. The total Hamiltonian of the three-site setup shown in figure 1 then reads

$$\begin{equation}\hat{H}=\sum\limits _{i=1}^{3}\left(\omega \left\vert 1\right\rangle {\left\langle 1\right\vert }_{i}+\left(\omega +{\Delta}\right)\left\vert +\right\rangle {\left\langle +\right\vert }_{i}\right)+\sum\limits _{i\ne j}{\hat{V}}_{ij}.\end{equation} \tag{ A.2 }$$

Here, ω refers to the energy difference between the states $\left\vert 1\right\rangle$ and $\left\vert 0\right\rangle$, while Δ denotes the detuning of the state $\left\vert +\right\rangle$ shown in figure 1. The subscript of the projectors indicates the respective atom. In deriving the spin excitation transport in the triangle of Rydberg atoms we only need to consider the case of one and two excitations in the system, where we define an atom to be excited if it occupies either the state $\left\vert 1\right\rangle$ or $\left\vert +\right\rangle$. This definition of occupation differs from the main text, but is convenient for the derivation.

For the case of a single excitation, we adiabatically eliminate the state $\left\vert +\right\rangle$, which is detuned by Δ and thus will only be virtually excited. This leads to a 3 × 3 matrix describing the effective dipole–dipole interaction of three two-level systems with a single excitation (in state $\left\vert 1\right\rangle$ or $\left\vert +\right\rangle$)

$$\begin{equation}{\hat{H}}_{\text{eff}}=\left(\begin{matrix}\hfill 0\hfill & \hfill \left\langle 100\right\vert \hat{H}\left\vert 010\right\rangle \hfill & \hfill \left\langle 100\right\vert \hat{H}\left\vert 001\right\rangle \hfill \\ \hfill \left\langle 010\right\vert \hat{H}\left\vert 100\right\rangle \hfill & \hfill 0\hfill & \hfill \left\langle 010\right\vert \hat{H}\left\vert 001\right\rangle \hfill \\ \hfill \left\langle 001\right\vert \hat{H}\left\vert 100\right\rangle \hfill & \hfill \left\langle 001\right\vert \hat{H}\left\vert 010\right\rangle \hfill & \hfill 0\hfill \end{matrix}\right)-\frac{1}{{\Delta}}\left(\begin{matrix}\hfill {U}_{1}\hfill & \hfill {h}_{1\to 3\to 2}\hfill & \hfill {h}_{1\to 2\to 3}\hfill \\ \hfill {h}_{2\to 3\to 1}\hfill & \hfill {U}_{2}\hfill & \hfill {h}_{2\to 1\to 3}\hfill \\ \hfill {h}_{3\to 2\to 1}\hfill & \hfill {h}_{3\to 1\to 2}\hfill & \hfill {U}_{3}\hfill \end{matrix}\right),\end{equation} \tag{ A.3 }$$

where |100⟩ denotes the state where the three atoms are in states |1⟩, $\left\vert 0\right\rangle$, $\left\vert 0\right\rangle$, respectively. In this calculation we neglected all terms of higher order in 1/Δ. In equation (A.3), the symbol h_i→j→k refers to the indirect hopping process going from the $\left\vert 1\right\rangle$ state at site i via the $\left\vert +\right\rangle$ state at site j to $\left\vert 1\right\rangle$ at site k:

$$\begin{equation}{h}_{1\to 2\to 3}=\left\langle 0,0,1\right\vert \hat{H}\left\vert 0,+,0\right\rangle \left\langle 0,+,0\right\vert \hat{H}\left\vert 1,0,0\right\rangle .\end{equation} \tag{ A.4 }$$

Therefore, h_i→j→k describes exactly the process shown in figure 2 in the main text. Similarly, the symbol U_i denotes the virtual processes where the initial site i and the final site k are identical, e.g. for i = 1:

$$\begin{equation}{U}_{1}={h}_{1\to 2\to 1}+{h}_{1\to 3\to 1}.\end{equation} \tag{ A.5 }$$

A key aspect is the fact that the latter 1/Δ term in (A.3) only occurs within the single-excitation manifold. E.g. if site 2 is occupied, the matrix element vanishes

$$\begin{equation}\left\langle 0,1,1\right\vert \hat{H}\left\vert 0,+,0\right\rangle \left\langle 0,+,0\right\vert \hat{H}\left\vert 1,1,0\right\rangle =0,\end{equation} \tag{ A.6 }$$

and only the direct dipole–dipole flopping process between site 1 and site 3 remains (first term in (A.3)). The identical consideration holds for the terms U_i . As a result, both the hopping amplitudes and the NN interaction terms are density-dependent. We can make this effect explicit by writing down the Hamiltonian using spin operators. We define the energy scale

$$\begin{equation}J=-\left\langle 0,1,0\right\vert \hat{H}\left\vert 0,0,1\right\rangle =-\left\langle 0,1,0\right\vert {\hat{V}}_{1,2}\left\vert 1,0,0\right\rangle ,\end{equation} \tag{ A.7 }$$

of the direct hopping process over the distance R as shown in figure 1. The explicit evaluation of these matrix elements can be done using Wigner-symbols. For this, we reduce the expectation values of the interaction Hamiltonian (A.7) to the dipole operators

$$\begin{equation}\left\langle 0,1,0\right\vert {\hat{V}}_{1,2}\left\vert 1,0,0\right\rangle \end{equation} \tag{ A.8 }$$

$$\begin{equation}=\frac{1}{8\pi {{\epsilon}}_{0}{R}_{1,2}^{3}}\left\langle 0,1,0\right\vert \left({\hat{d}}_{1}^{+}{\hat{d}}_{2}^{-}+{\hat{d}}_{1}^{-}{\hat{d}}_{2}^{+}\right)\left\vert 1,0,0\right\rangle \end{equation} \tag{ A.9 }$$

$$\begin{equation}=\frac{1}{8\pi {{\epsilon}}_{0}{R}^{3}}\left\langle 0,1,0\right\vert {\hat{d}}_{1}^{-}{\hat{d}}_{2}^{+}\left\vert 1,0,0\right\rangle =\frac{1}{8\pi {{\epsilon}}_{0}{R}^{3}}\left\langle 1\right\vert {\hat{d}}_{1}^{-}\left\vert 0\right\rangle \left\langle 0\right\vert {\hat{d}}_{2}^{+}\left\vert 1\right\rangle ,\end{equation} \tag{ A.10 }$$

In this way, all matrix elements can be reduced to obtaining a result for the dipole operator transition elements

$$\begin{equation}\left\langle 0\right\vert {\hat{d}}^{-}\left\vert +\right\rangle \quad \left\langle 1\right\vert {\hat{d}}^{-}\left\vert 0\right\rangle \quad \left\langle 0\right\vert {\hat{d}}^{+}\left\vert 1\right\rangle \quad \left\langle +\right\vert {\hat{d}}^{+}\left\vert 0\right\rangle .\end{equation} \tag{ A.11 }$$

We can evaluate these matrix elements by rewriting the dipole operator

$$\begin{equation}{\hat{d}}^{+}=-\frac{1}{\sqrt{2}}\left({\hat{d}}^{x}+i{\hat{d}}^{y}\right)=-\frac{q}{\sqrt{2}}\left(\hat{x}+i\hat{y}\right)=-{\boldsymbol{e}}_{+}\cdot \hat{\boldsymbol{d}}=-q{\boldsymbol{e}}_{+}\cdot \hat{\boldsymbol{r}}=-q{\boldsymbol{e}}_{+}\cdot {\hat{\boldsymbol{e}}}_{r}\cdot r,\end{equation} \tag{ A.12 }$$

where we have used ${\boldsymbol{e}}_{+}=\frac{1}{\sqrt{2}}\left({\boldsymbol{e}}_{x}+i{\boldsymbol{e}}_{y}\right)$. Similar expression can be derived for ${\hat{d}}^{-}$. Additionally, we can decompose the position unit vector using the spherical harmonics Y_l,m (see [29])

$$\begin{equation}{\boldsymbol{e}}_{r}=\sqrt{\frac{4\pi }{3}}\left[{Y}_{1,-1}\frac{{\boldsymbol{e}}_{x}+\mathrm{i}{\boldsymbol{e}}_{y}}{\sqrt{2}}+{Y}_{1,1}\frac{-{\boldsymbol{e}}_{x}+\mathrm{i}{\boldsymbol{e}}_{y}}{\sqrt{2}}+{Y}_{1,0}{\boldsymbol{e}}_{z}\right].\end{equation} \tag{ A.13 }$$

To calculate the transition matrix element between two states $\left\vert {\Psi}\right\rangle$ and $\left\vert {{\Psi}}^{\prime }\right\rangle$ we will consider the angular and radial parts separately:

$$\begin{equation}\left\langle {{\Psi}}^{\prime }\right\vert {\hat{d}}^{\pm }\left\vert {\Psi}\right\rangle =-q{\boldsymbol{e}}_{\pm }\left\langle {{\Psi}}^{\prime }\right\vert {\hat{\boldsymbol{e}}}_{r}\cdot r\left\vert {\Psi}\right\rangle ,\end{equation} \tag{ A.14 }$$

$$\begin{equation}=-q\left\langle {\psi }^{\prime }\right\vert {C}_{1,\pm 1}\left\vert \psi \right\rangle {\int }_{0}^{\infty }\mathrm{d}r\enspace {R}_{n,{l}^{\prime }}^{\prime }(r)r{R}_{n,l}(r).\end{equation} \tag{ A.15 }$$

In the above equation, $\left\vert \psi \right\rangle$ and $\left\vert {\psi }^{\prime }\right\rangle$ denote the angular and ${R}_{n,{l}^{\prime }}^{\prime },{R}_{n,l}$ the radial wavefunctions. The latter is equal to

$$\begin{equation}\xi \equiv {\int }_{0}^{\infty }\mathrm{d}r\enspace {R}_{n,{l}^{\prime }}^{\prime }(r)r{R}_{n,l}(r),\end{equation} \tag{ A.16 }$$

and is constant over all transitions relevant here, so we simply set it to ξ. For the angular part, we must evaluate the following expression

$$\begin{equation}\left\langle {\psi }^{\prime }\right\vert {\hat{d}}^{\pm }\left\vert \psi \right\rangle =q\left\langle {\psi }^{\prime }\right\vert {C}_{1,\pm 1}\left\vert \psi \right\rangle ,\end{equation} \tag{ A.17 }$$

$$\begin{equation}{C}_{l,m}=\sqrt{\frac{4\pi }{2l+1}}{Y}_{l,m},\end{equation} \tag{ A.18 }$$

where q = −e is the electron charge. Matrix elements of this kind can be evaluated using Wigner symbols (see [30])

$$\begin{equation}\left\langle {\psi }^{\prime }\right\vert {C}_{k,p}\left\vert \psi \right\rangle ={\delta }_{{S}^{\prime },S}\sqrt{(2{J}^{\prime }+1)(2J+1)(2{L}^{\prime }+1)(2L+1)}{(-1)}^{{M}^{\prime }-S}\end{equation} \tag{ A.19 }$$

$$\begin{align}\hfill & \quad \times \left(\begin{matrix}\hfill {J}^{\prime }\hfill & \hfill J\hfill & \hfill k\hfill \\ \hfill -{M}^{\prime }\hfill & \hfill M\hfill & \hfill p\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {L}^{\prime }\hfill & \hfill k\hfill & \hfill L\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right)\left\{\begin{matrix}\hfill {L}^{\prime }\hfill & \hfill L\hfill & \hfill k\hfill \\ \hfill J\hfill & \hfill {J}^{\prime }\hfill & \hfill S\hfill \end{matrix}\right\}.\hfill \end{align} \tag{ A.20 }$$

The brackets in the second line are referred to as Wigner-3j-symbols and Wigner-6j-symbols, respectively, and the variables contained in the above equation represent the quantum numbers of both states. Using the atomic states in question, we arrive at (the common factor of ξ is omitted)

$$\begin{equation}\left\langle 0\right\vert {\hat{d}}^{-}\left\vert +\right\rangle =\frac{e}{\sqrt{3}}\quad \left\langle 1\right\vert {\hat{d}}^{-}\left\vert 0\right\rangle \enspace \enspace =-\frac{e}{3},\end{equation} \tag{ A.21 }$$

$$\begin{equation}\left\langle 0\right\vert {\hat{d}}^{+}\left\vert 1\right\rangle =\frac{e}{3}\quad \enspace \enspace \enspace \left\langle +\right\vert {\hat{d}}^{+}\left\vert 0\right\rangle =-\frac{e}{\sqrt{3}}.\end{equation} \tag{ A.22 }$$

Now we can calculate explicitly the transition elements from above, starting from the unit of energy J:

$$\begin{equation}J=-\left\langle 0,1,0\right\vert {\hat{V}}_{1,2}\left\vert 1,0,0\right\rangle =-\frac{1}{8\pi {{\epsilon}}_{0}{R}^{3}}\left\langle 1\right\vert {\hat{d}}_{1}^{-}\left\vert 0\right\rangle \left\langle 0\right\vert {\hat{d}}_{2}^{+}\left\vert 1\right\rangle \end{equation} \tag{ A.23 }$$

$$\begin{equation}=\frac{{e}^{2}}{72\pi {{\epsilon}}_{0}{R}^{3}}.\end{equation} \tag{ A.24 }$$

Analogously, we calculate the indirect transport processes from above, h_1→2→3

$$\begin{equation}{h}_{1\to 2\to 3}=\left\langle 0,0,1\right\vert \hat{H}\left\vert 0,+,0\right\rangle \left\langle 0,+,0\right\vert \hat{H}\left\vert 1,0,0\right\rangle \end{equation} \tag{ A.25 }$$

$$\begin{equation}=\left\langle 0,0,1\right\vert {\hat{V}}_{2,3}\left\vert 0,+,0\right\rangle \left\langle 0,+,0\right\vert {\hat{V}}_{1,2}\left\vert 1,0,0\right\rangle \end{equation} \tag{ A.26 }$$

$$\begin{equation}={\left(\frac{1}{4\pi {{\epsilon}}_{0}{R}^{3}}\frac{3}{2}\right)}^{2}\left\langle 0,0,1\right\vert {\hat{d}}_{2}^{-}{\hat{d}}_{3}^{-}{e}^{2\mathrm{i}{\phi }_{2,3}}\left\vert 0,+,0\right\rangle \left\langle 0,+,0\right\vert {\hat{d}}_{1}^{+}{\hat{d}}_{2}^{+}{e}^{-2\mathrm{i}{\phi }_{1,2}}\left\vert 1,0,0\right\rangle \end{equation} \tag{ A.27 }$$

$$\begin{equation}={\left(\frac{1}{4\pi {{\epsilon}}_{0}{R}^{3}}\frac{3}{2}\right)}^{2}{e}^{-4\mathrm{i}\alpha }\frac{{e}^{4}}{27}\end{equation} \tag{ A.28 }$$

$$\begin{equation}=27{J}^{2}{e}^{-4\mathrm{i}\alpha }.\end{equation} \tag{ A.29 }$$

Similarly, one can evaluate the density–density terms U_i .

So far we have performed the derivation in the single-excitation manifold of the Hilbert space. In the case of two excitations on the triangle, all h_i→j→k vanish, whereas h_i→j→i remains non-zero if sites i and k are occupied.

Using this we arrive at the Hamiltonian in hard-core boson language as already shown in the main text

$$\begin{equation}\hat{H}=-J\sum\limits _{i\ne j=1}^{3}{\hat{b}}_{i}^{{\dagger}}{\hat{b}}_{j}+\text{h.c.}-2gJ\sum\limits _{i\ne j\ne k}^{3}{\hat{b}}_{i}^{{\dagger}}{\hat{b}}_{j}\left(1-{\hat{n}}_{k}\right){e}^{-4i{{\epsilon}}_{i,j,k}\alpha }+\text{h.c.}-2gJ\sum\limits _{i\ne k}{\hat{n}}_{i}\left(1-{\hat{n}}_{k}\right),\end{equation} \tag{ A.30 }$$

where we again introduce $g=\frac{27J}{2{\Delta}}$ as a dimensionless parameter.

We have seen in section 5.1 that in the density-ordered phase of large η an oscillatory regular flux pattern of the quantum gauge field is induced with long-range correlations. A measurable consequence of this is the formation of a vortex lattice of local currents as seen in figure 10(b). Although there is no net current in an insulating phase with density order and first-order NN correlations are suppressed, their amplitudes decays only algebraically with the density–density repulsion η. Thus local vortex currents can emerge even in the insulating phase. Their amplitude is expected to decrease as 1/η, which is verified by the very good agreement of the numerical data for the current I_j→j+2 with an 1/r fit in figure B1. We marked the transition for g = 2 (green curve) by the red dotted line.

Zoom In Zoom Out Reset image size