Transport-induced melting of crystals of Rydberg dressed atoms in a one-dimensional lattice

Recently, many-body systems with non-local, power-law interactions have gained considerable interest, as their non-local coupling can give rise to quantum phases that do not exist for point-like interactions [1]. For example, in one spatial dimension repulsive dipole–dipole or van der Waals interactions of atoms can lead to a variety of crystalline ground state phases with less than unity filling in the presence of a commensurable periodic lattice [2–5]. For very strong lattice confinements, the filling fraction forms a so-called complete devil's staircase as a function of the chemical potential μ, i.e. every rational filling between 0 and 1 is stable for a finite interval of μ and these intervals form a dense staircase [6]. Power-law interactions arise, e.g., for dipolar atoms or molecules [1]. A very interesting alternative approach is Rydberg gases, for which there has been considerable experimental progress in recent years [7–13]. While most previous experiments implemented schemes where Rydberg excitations were created by continuous near-resonant laser driving, alternative proposals have been put forward to use far-detuned light fields to admix a small component of a Rydberg state to the ground state of atoms [3, 14, 15]. The potential advantage of this 'Rydberg-dressing' is that the interaction strength can be tuned to a certain extent and is reduced to an energy scale where the center-of-mass motion of atoms can become relevant. Furthermore, continuous driving by a near-resonant laser does not, in general, conserve the number of Rydberg excitations, while in the case of Rydberg-dressing the number of atoms is, to a good approximation, a conserved quantity. While the non-local repulsive interaction favors crystalline structures with non-unitary filling in lattice gases, hopping between lattice sites induced by the tunneling of the atoms can lead to a melting of the crystals [16, 17].

In this paper, we analyze the melting process induced by the hopping of Rydberg dressed atoms. In particular, we derive the phase boundaries of stable crystalline structures as a function of the hopping amplitude J both analytically and with exact numerical simulations. To this end we employ a second-order strong coupling expansion on the one hand and numerical density-matrix renormalization group (DMRG) simulations, adapted to long-range interactions, on the other. Both show excellent agreement up to second order in J. A corresponding strong-coupling analysis of the phase diagram has been discussed before for the case of dipolar 1/r³ interactions in [6]; however, no explicit expressions were given. We derive here analytic expressions that, in the case of fast decaying potentials such as the 1/r⁶ van der Waals potential, can to a very good approximation be written in a compact closed form. Furthermore, for the first time we provide a comparison to numerical data from DMRG simulations adapted to long-range interactions.

As a starting point, we consider the generalized Bose–Hubbard model in one dimension

$$\begin{equation} \skew3\hat{H} = \underbrace{\frac{1}{2}\sum_{i\ne j} V(|r_i-r_j|)\hat{b}^\dagger_i\hat{b}_i\hat{b}^\dagger_j\hat{b}_j}_{\skew3\hat{H}_0} - \underbrace{ J \sum_{i} \big(\hat{b}^\dagger_i\hat{b}_{i+1} + \hat{b}^\dagger_{i+1}\hat{b}_i\big)}_{\skew3\hat{H}_J}. \end{equation} \tag{ 1 }$$

For convenience we have set ℏ = a = 1, where a is the lattice constant. $\hat {b}_i,\hat {b}_i^\dagger$ are the annihilation and creation operators at site i, and J is the hopping amplitude between neighboring lattice sites describing the tunneling of atoms between adjacent lattice sites [18]. Most parts of our discussion are valid for any convex interaction potential, i.e.

$$\begin{equation} V(r+1)+V(r-1)\geqslant 2\,V(r). \end{equation} \tag{ 2 }$$

However, to be specific, we discuss the power-law potentials of the form

$$\begin{equation} V(r)= \frac{\tilde{C}_\beta}{|r|^\beta}, \end{equation} \tag{ 3 }$$

with $\beta \in \mathbb {N}\backslash \{1\}$ and $\tilde {C}_\beta >0$ being the interaction coefficients. β = 6 corresponds to a van der Waals-type interaction which results from virtual (i.e. off-resonant) long-wavelength electromagnetic transitions in the manifold of Rydberg states. If there is an accidental or engineered resonance, a so-called Förster resonance, the case β = 3 can also be of relevance. The model interaction potential (3) diverges for r → 0 and thus two particles cannot sit on the same lattice site. To a very good approximation the same is true if one considers the exact interaction potential between Rydberg atoms, which deviates from the above power law for small distances (see, e.g., [19, 20]). In fact, for the Rydberg-dressing scheme, one finds an effective two-body interaction potential between atoms at positions r_i and r_k [3, 14, 15],

$$\begin{equation} V(\mathbf{r}_i,\mathbf{r}_k) = \frac{\eta^2 C_6}{|\mathbf{r}_i-\mathbf{r}_k|^6 + R_{\mathrm{c}}^6} =\frac{{\tilde C}_6}{r^6 + R_{\mathrm{c}}^6},\quad r=|\mathbf{r}_i-\mathbf{r}_k|. \end{equation} \tag{ 4 }$$

Here C₆ is the Rydberg interaction coefficient and $R_{\mathrm {c}}=\left (C_6/(2\Delta )\right )^{1/6}$ describes a characteristic length scale of the interaction potential, where Δ ≫ |Ω| is an effective detuning of the off-resonant laser driving. η ∼ Ω²/Δ² ≪ 1 accounts for the small admixture of the Rydberg state. Considering, e.g., ⁸⁷Rb with a Rydberg state of principal quantum number around n = 60 and detunings of a few tens of MHz, one finds a value for R_c of the order of 1 μm to a few μm. In a recent experiment, atoms were loaded into optical lattices and excited to Rydberg states with n = 55–80 [21]. In this experiment, the lattice spacing was tunable between 0.5 μm and 13 μm. Thus the cut-off length R_c may become smaller than the average distance between particles given by the lattice constant divided by the average occupation number per lattice site. In this case, its effect can be disregarded and this is the parameter regime we are considering here. For that reason $\hat {b}_i,\hat {b}_i^\dagger$ are assumed to describe hard-core bosons with commutation relations

$$\begin{equation} [ \hat{b}_i, \hat{b}_j^\dagger]=[ \hat{b}_i, \hat{b}_j]=[ \hat{b}_i^\dagger, \hat{b}_j^\dagger] = 0, \quad i\ne j, \end{equation} \tag{ 5 }$$

$$\begin{equation} \{\hat{b}_i, \hat{b}_i\}=\{\hat{b}_i^\dagger, \hat{b}_i^\dagger\} = 0, \end{equation} \tag{ 6 }$$

$$\begin{equation} \{\hat{b}_i, \hat{b}_i^\dagger\}=1. \end{equation} \tag{ 7 }$$

We note that although we discuss here the particular case of the integer value β = 6, the qualitative structure of the phase diagram holds true for any positive value of β.

Let us first consider the ground state of Hamiltonian (1) in the limit of vanishing hopping, i.e. J → 0. In this limit, the Hamiltonian is diagonal in the number basis, where each lattice site contains exactly zero or exactly one particle, which corresponds to a classical situation. Consequently, the energy can be minimized by finding the configuration of classical particles that gives the smallest interaction energy. As the interaction potential is convex (see equation (2)), the minimum energy for commensurable filling fractions $q=\frac {m}{n}\leqslant 1$, where $m,n\in \mathbb {N}$ are relatively prime, is attained by a regular pattern with unit cells of size n [22]. For a given chemical potential the ground state is such a phase with rational q. For these values of q the corresponding phase is incompressible, i.e. particle as well as hole excitations require finite energy, i.e. μ⁺ > μ⁻, where

$$\begin{equation} \mu^\pm(q)=\left.\frac{\partial E}{\partial N}\right\vert_{\rho=q\pm0} = \pm\left(E^\pm(q)-E(q)\right). \end{equation} \tag{ 8 }$$

Here E(q) is the ground state energy for filling fraction q. E^±(q) is the corresponding energy where one particle has been added (E⁺(q)) or respectively removed (E⁻(q)). ρ is the average number of particles per lattice site. Because of particle–hole symmetry,

$$\begin{equation} \mu(q) = \Omega-\mu(1-q), \end{equation} \tag{ 9 }$$

where Ω is the energy per particle in the completely filled lattice (Ω is the Riemann zeta function ζ(β) in the case of (3)), it is sufficient to consider only filling fractions $q\leqslant \frac {1}{2}$.

The phase diagram in zeroth order of the hopping J has been discussed some time ago by Bak and Bruinsma [23]. Here, we will briefly review their results and make use of their notation: X⁰_i denotes the position of the ith particle. X^p_i = X⁰_p − X⁰_i describes the distance between the pth and the ith particle where the particles are numbered from left to right. In the ground state of $\skew3\hat {H}_0$, Hubbard's solution [22] shows that all possible distances are given by

$$\begin{equation} X_i^p = r_p \quad {\rm or}\quad X_i^p = r_p +1, \end{equation} \tag{ 10 }$$

where $r_p < \frac {p}{q} < r_p + 1$. If $\frac {p}{q} \in \mathbb {N}$, then $X_i^p = r_p = \frac {p}{q}$ and all particles are equivalent. Examples of ground states are shown in figure 1. Due to the convexity of the interaction, all nearest-neighbor separations are maximally close to the average separation, i.e. they are either $\lfloor \frac 1q\rfloor$ or $\lceil \frac 1q\rceil$.

Zoom In Zoom Out Reset image size

The chemical potential for J = 0 is given by

$$\begin{equation} \mu^{(0)}_\pm(q)= \pm ( \langle q^\pm| \skew3\hat{H}_0 | q^\pm\rangle - \langle q| \skew3\hat{H}_0 | q\rangle), \end{equation} \tag{ 11 }$$

where |q〉 denotes the ground state of $\skew3\hat {H}_0$ for filling fractions q, while |q^±〉 denotes the corresponding ground states where one particle has been added, |q⁺〉, respectively removed, |q⁻〉. Although both expectation values in (11) are infinite in the thermodynamic limit, their difference is finite and can be calculated by summing the change in interaction energy particle by particle. The same is true for the first- and second-order corrections in the hopping amplitude J discussed below. For filling fractions $q=\frac {m}{n}$ there will be n defects for |q^±〉, because it is energetically favorable to break up an extra particle (or a hole) into n defects each with fractional charge $\pm \frac 1n$. Some of these defects are displayed in figure 2. In the thermodynamical limit they will be separated by arbitrarily long distances so that we can assume without approximation that they do not interact with each other at all. Note that |q^±〉 is not uniquely defined, but the position of defects is arbitrary, as long as they are well separated. For filling fractions of the form $q=\frac {1}{n}$ the chemical potential is given by [23]

$$\begin{equation} \mu_\pm^{(0)} = \pm\sum_{p=1}^\infty ( r_p\ V(r_p\mp1)- (r_p\mp1) V(r_p)). \end{equation} \tag{ 12 }$$

We can go a step further and evaluate (12) for the power-law potential (3) analytically,

$$\begin{equation} \fl \mu_\pm^{(0)}/\tilde{C}_\beta = \frac{(-1)^\beta \Psi^{(\beta-1)}(\frac{n\mp1}{n})}{n^\beta(\beta-1)! } + \frac{\zeta(\beta)}{n^\beta} \pm \left(\frac{(-1)^{(\beta-1)} \Psi^{(\beta-2)}(\frac{n\mp1}{n})}{n^{(\beta-1)}(\beta-2)! } - \frac{\zeta(\beta-1)}{n^{(\beta-1)}}\right), \end{equation}\noindent \tag{ 13 }$$

using the digamma function $\Psi (z)= \frac {\Gamma '(z)}{\Gamma (z)}$ and its lth derivative Ψ^(l)(z). To describe more general filling fractions $q=\frac {m}{n}$ with m ≠ 1 as well, we only have to replace r_p → r_p + 1 for μ⁽⁰⁾₊ in (12) if $r_p\not \in \mathbb {N}$. Plotting the filling fraction as a function of the chemical potential gives rise to a complete devil's staircase (see figure 3). Every rational filling fraction between 0 and 1 has a finite stable range with respect to the chemical potential and the functional dependence q(μ) consists of a dense set of steps between these stable plateaus, a so-called devil's staircase. The union of these intervals covers the full range of μ. It is clearly seen that filling fractions with a small denominator n are most stable, as can be seen directly from equation (13). These intervals of stability depend crucially on the exponent β: for large β, the interval of the stable phase $q=\frac 12$ overgrows by far the size of all the other phases. Note that for n big enough (n ± 1)/n is almost constant and hence the mean chemical potential μ⁽⁰⁾ ≡ (μ⁽⁰⁾₊ − μ⁽⁰⁾₋)/2 scales as n^−β. Thus, below half-filling the fraction q scales as $(\mu ^{(0)}/\tilde C_\beta )^{-\beta }$ as indicated in the inset of figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As mentioned in the introduction, the true interaction potential between Rydberg dressed atoms becomes flat below a certain distance. As a consequence, the atoms can be treated as hard-core bosons only for sufficiently small energies, i.e. sufficiently small chemical potentials. The finite cut-off R_c will lead to a qualitative change of the phase diagram in parameter regions where the separation between atoms becomes smaller than R_c. This is, however, the case only if the filling q becomes large and, in particular, exceeds the ratio R_c/a.

We are now interested in the melting of a crystal phase with increasing hopping rate J. Therefore we perform perturbation theory up to second order in J. Similar calculations have been done by Burnell et al for dipolar interactions V ∼1/r³ [6], but no compact analytic expression was given.

Let us start with first-order processes in J. Then, the chemical potential reads as $\mu ^{(1)}_\pm (q)= \pm ( \langle q^\pm | \skew3\hat {H}_J | q^\pm \rangle - \langle q| \skew3\hat {H}_J | q\rangle)$. Hopping of any single particle in state |q〉 will not contribute to any energy correction, because the resulting state has no overlap with the J = 0 ground state. The same is true for hopping of any but the 2n particles of the state |q^±〉 which sit at the left and right borders of any of the n defects. Hopping of these particles will lead to a hopping of the respective defect by n-sites; see figure 4. As the states with localized defects are degenerated with respect to $\skew3\hat {H}_0$, we have to apply degenerate perturbation theory. Therefore, we look for a basis of these states in which $\skew3\hat {H}_J$ is diagonal. This is given by states where all the defects are delocalized over many lattice sites (yet well separated from each other) with some quasi-momentum k. The energy is minimized for a state with quasi-momentum zero. The resulting first-order correction to the phase boundary thus takes the simple form

$$\begin{equation} \mu^{(1)}_\pm(q)= \mp 2 n J. \end{equation} \tag{ 14 }$$

Zoom In Zoom Out Reset image size

For small but nonzero J, gaps open between the incompressible phases of rational filling, because the kinetic energy gained by delocalization of defects favors a finite density of defects. At the tip of the phase J becomes large enough such that the creation of particle and hole defects becomes favorable even at commensurate filling fractions and the crystalline phase melts. This is completely analogous to the ordinary Hubbard model. However, for any arbitrarily small but finite value of J almost all phases with rational filling become unstable, leaving only a finite number (those with the smallest denominator n) and destroying the devil's staircase.

In second-order perturbation the energy corrections are given by $E^{(2)} = \sum _{\Phi \ne \Psi }\frac {\langle \Psi | \skew3\hat {H}_J| \Phi \rangle \langle \Phi | \skew3\hat {H}_J | \Psi \rangle }{E^{(0)}_\Psi -E^{(0)}_\Phi }$, where |Ψ〉 denotes the ground state of $\skew3\hat {H}_0$ and {|Φ〉} is a complete set of orthonormal states. In the following, we consider only filling fractions of the form $q=\frac {1}{n}$. In principle, it is possible to extend the calculations to other fractions. However, finding a general formula that describes all energy differences in the denominator for all hopping processes is a far more difficult task. In figure 5, all possible processes for $q^\pm =\frac {1}{n}^\pm$ are shown. As has been seen in first order, the states with localized defects are degenerate with respect to $\skew3\hat {H}_0$. In degenerate second-order perturbation theory the energy correction is given by $\langle q^\pm |\skew3\hat {H}_J \skew3\hat P \skew3\hat {H}_J|q^\pm \rangle /\big (E^{(0)}_{q^\pm }-E^{(0)}_\Phi \big )$, where $\skew3\hat P$ projects out the subspace of ground states of $\skew3\hat {H}_0$. The state |q^±〉 has to be an eigenstate of $\skew3\hat {H}_J$ restricted to the degenerate subspace. Therefore, intermediate states |Φ〉 with E⁽⁰⁾_Φ = E⁽⁰⁾_q^± do not contribute to the energy corrections. For states |q〉, only process (i) in figure 4 contribute. For the commensurate state |q〉 one finds for the second-order correction to the energy

$$\begin{equation} E^{(2)}(q) = 2J^2 N \frac{1}{\Delta E^{(0)}(q)}, \end{equation} \tag{ 15 }$$

where N is the number of particles in the lattice. ΔE⁽⁰⁾(q) = E⁽⁰⁾(q) − E⁽⁰⁾₊(q) is the energy difference between the ground state |q〉 and the same state but with one particle hopped by one site. With the notation introduced above, it can be evaluated to

$$\begin{eqnarray} \Delta E^{(0)}(q) &= \sum_{p=1}^{\infty}(2\,V(r_p)-V(r_p-1) - V(r_p +1))\\ &=\frac{2\zeta(\beta)}{n^\beta}-\frac{(-1)^\beta}{n^\beta(\beta-1)!}\! \left[\Psi^{(\beta-1)}\! \left(\frac{n-1}{n}\right)\! +\Psi^{(\beta-1)}\! \left(\frac{n+1}{n}\right)\right]\! .\nonumber \end{eqnarray} \tag{ 16 }$$

In order to obtain the (finite) chemical potential we also have to calculate the energy corrections E⁽²⁾(q^±) for the states with one extra particle or hole |q^±〉 in second-order perturbation theory. To this end, let us look more closely at all possible hopping processes displayed in figure 5. As has been done in first order, we will concentrate on a single defect multiplying the resulting contribution with its number n. Process (ii) describes a virtual deformation of the defect while process (iii) corresponds to an effective hopping of the defects over two unit cells. Both processes have only to be taken into account for one particle of the defect and one particle on its right side and therefore contribute only to a finite energy value. In the thermodynamical limit, process (i) takes place for an infinite number of particles. For this reason, its contributions to the energy corrections will diverge. However, for a decaying potential the contribution to the energy correction of a particle far away from the defect will be the same as the contribution of a particle of the commensurate state |q〉. Subtracting now the energy corrections for states |q〉 and |q^±〉, these contributions of process (i) will mostly cancel. We find that

$$\begin{equation} \fl E^{(2)}(q^\pm) = 2 n J^2 \left\{\sum_{i=1}^{N_0^\pm} \left(\frac{1}{\Delta E^{(0)\pm}_{+i,-}} +\frac{1}{\Delta E^{(0)\pm}_{+i,+}}\right) + \frac{1-\frac12(\delta_{n,2}\pm\delta_{n,2})}{\Delta E^{(0)\pm}_{+0,\mp}} +\cos(2 k) \frac{1}{\Delta E^{(0)\pm}_{+1,\pm}}\right\} , \end{equation} \tag{ 17 }$$

where ΔE^(0)±_+j,− denotes the energy difference between the ground state |q^±〉 and the same state, where the jth particle right of the defect moves right (the second subscript +) or left (the second subscript −). For the meaning of position '+0' see figure 5. We have assumed here that the defect is in the center of the system and thus the summation limit is N^±₀ = 1/2[(N ± 1)q − 2]. As |q^±〉 has to be an eigenstate of $\skew3\hat {H}_J$ in the degenerate subspace the defects are completely delocalized. As in first order, the ground state has quasi-momentum k = 0. Then, the second-order correction to the chemical potential can be written as

$$\begin{equation} \mu^{(2)}_\pm = \pm2 n J^2 \Bigg\{ \sum_{j=1}^{\infty}\Big[\frac{1}{\Delta E^{(0)\pm}_{+j,-}}+\frac{1}{\Delta E^{(0)\pm}_{+j,+}}-\frac{2}{\Delta E^{(0)}}\Big]-\left(1\mp\frac{1}{2n}\right)\frac{2}{\Delta E^{(0)}} \end{equation} \tag{ 18 }$$

$$\begin{equation} \hphantom{\mu^{(2)}_\pm = } +\frac{1-\frac12(\delta_{n,2}\pm\delta_{n,2})}{\Delta E^{(0)\pm}_{+0,\mp}} \end{equation} \tag{ 19 }$$

$$\begin{equation} \hphantom{\mu^{(2)}_\pm = } +\frac{1}{\Delta E^{(0)\pm}_{+1,\pm} }\Bigg\}. \end{equation} \tag{ 20 }$$

Part (18) takes into account all processes of type (i). Because of symmetry it is only necessary to sum over particles on the right side of the defect. The last term in this line is due to the fact that there are fewer particles in |q^±〉 taking part in process (i) than there are particles in state |q〉. Part (19) describes process (ii) and part (20) process (iii). Calculating all the energy differences analogous to the case of ΔE⁽⁰⁾, the chemical potential for filling fractions $q=\frac {1}{n}$ can finally be written in second order in J as

$$\begin{eqnarray} &&\fl \mu^{(2)}_{\pm} = \pm\, 2 n \frac{J^2}{\tilde{C}_\beta}\, \Bigg\{ \sum_{j=1}^\infty S_j^\pm(n) -\frac{2\mp\frac{1}{n}}{\frac{2\zeta(\beta)}{n^\beta}-\frac{(-1)^\beta}{ n^\beta(\beta-1)!}[\Psi^{(\beta-1)}(\frac{n-1}{n})+\Psi^{(\beta-1)}(\frac{n+1}{n})]}\nonumber\\ &&+\frac{1-\frac12(\delta_{n,2}\pm\delta_{n,2})}{\frac{\zeta(\beta)}{n^\beta}+\frac{(-1)^\beta}{ n^\beta(\beta-1)!}[\Psi^{(\beta-1)}(\frac{n\mp1}{n})-\Psi^{(\beta-1)}(\frac{n\pm1}{n})-\Psi^{(\beta-1)}(\frac{n\mp2}{n})]}\nonumber\\ &&+\frac{1}{\frac{2}{n^\beta}-\frac{1}{(n+1)^\beta}-\frac{1}{(n-1)^\beta}} \Bigg\}, \end{eqnarray} \noindent \tag{ 21 }$$

where the terms S^±_j(n) are given by

$$\begin{eqnarray} &&\fl S_j^\pm(n) =\Bigg\{\sum_{p=1}^j\left[\frac{2}{(pn)^\beta}-\frac{1}{(pn+1)^\beta}-\frac{1}{(pn-1)^\beta}\right]\Bigg\}^{-1}\nonumber\\ &&+\Bigg\{ \sum_{p=1}^j\!\Big[\frac{1}{(pn)^\beta}\!-\!\frac{1}{(pn+1)^\beta}\!-\!\frac{1}{(pn-1)^\beta}\Big]\!+\!\frac{\zeta(\beta)}{n^\beta}\!+\! \frac{(-1)^\beta}{ n^\beta(\beta-1)!}\nonumber\\ &&\times\Big[\Psi^{(\beta-1)}\left(\frac{n\mp1}{n}\!+\!j\right)-\Psi^{(\beta-1)}\left(\frac{n\pm1}{n}\!+\!j\right)-\Psi^{(\beta-1)}\left(\frac{n\mp2}{n}\!+\!j\right)\Big]\Bigg\}^{-1}\nonumber\\ &&-\frac{2}{\frac{2\zeta(\beta)}{n^\beta}-\frac{(-1)^\beta}{ n^\beta(\beta-1)!}[\Psi^{(\beta-1)}(\frac{n-1}{n})+\Psi^{(\beta-1)}(\frac{n+1}{n})]}. \end{eqnarray} \noindent \tag{ 22 }$$

If the power β of the interaction potential is sufficiently large and if the effect of the cut-off length R_c can be disregarded, i.e. if R_c < a/q, the entire sum $\sum _j\!S^\pm _j(n)$ contributes only very little. For the case of a van der Waals potential, i.e β = 6, we have numerically verified that the sum contributes less than a few per cent to the final result. Thus we can to a good approximation ignore the sum, which gives, for the chemical potential of the q = 1/n phases in a system with van der Waals interactions,

$$\begin{eqnarray} &&\fl \mu^{(2)}_{\pm} = \pm\, 2 n \frac{J^2}{\tilde{C}_\beta}\left\{ -\frac{7560 n^6 (2\mp\frac{1}{n})}{16\pi^3-63[\Psi^{(5)}(\frac{n-1}{n})+\Psi^{(5)}(\frac{n+1}{n})]}\right.\nonumber\\ &&\left.+\frac{7560 n^6(1-\frac12(\delta_{n,2}\pm\delta_{n,2}))}{8\pi^3+63[\Psi^{(5)}(\frac{n\mp1}{n})-\Psi^{(5)}(\frac{n\pm1}{n})-\Psi^{(5)}(\frac{n\mp2}{n})]}+\frac{1}{\frac{2}{n^6}-\frac{1}{(n+1)^6}-\frac{1}{(n-1)^6}} \right\}. \end{eqnarray}\noindent \tag{ 23 }$$

Zoom In Zoom Out Reset image size

Complementary to the perturbation analysis, we perform numerical calculations for the case of a van der Waals potential. For this we employ the density-matrix renormalization group algorithm [24]. It is a variational technique that minimizes the energy of the full Hamiltonian $\skew3\hat {H}_0+\skew3\hat {H}_J$ using a matrix product state (MPS) ansatz [25]. Although limited in the amount of entanglement along the lattice that they can capture, MPSs are known [26] to approach the true ground state of one-dimensional (1D) systems quickly with growing matrix dimension if interactions are of finite range. Although the interactions decay only polynomially in our model, we find that we can safely cut off the interaction at a finite distance of r lattice sites, and the ground state energy will be for all practical means independent of r as long as we restrict the filling fractions to denominators n that are small compared to r. In order to make our DMRG implementation more efficient, we group the interaction terms as originally introduced for computations in momentum space [27]. We also take advantage of the particle number conservation, which allows us to use MPS with the proper symmetry and fix the total particle number a priori, which is helpful in calculating the phase diagram.

To calculate the phase boundaries for a given filling fraction, we make direct use of (8) by computing E(q) and E^±(q) for a finite system². The minimum system size L to capture the physics in the thermodynamic limit for a given filling fraction $q=\frac {m}{n}$ can be estimated to be 2n², because it must fit n defects of size about n separated from each other and the boundaries by at least one unit cell of size n. For not too large n we can perform an infinite size extrapolation by employing different system sizes and assuming a 1/L scaling of the finite size error. In figure 6 we plot the resulting phase diagram for a van der Waals potential. The dashed lines display results from perturbation theory including up to the second order and the continuous lines are DMRG results. As expected the agreement is very good for small J/Δ₆. Perturbation theory, however, fails close to the critical points corresponding to the tips of the lobes of incompressible phases. Also using DMRG the exact position of the critical point is hard to compute due to strong finite-size effects. Accordingly, we have not attempted to extract these values. The continuous lines in the figure end at arbitrary values, while the phase boundary ends where these lines make contact.

Zoom In Zoom Out Reset image size

As can be seen from figure 6 the melting points of the different incompressible phases differ by many orders of magnitude in the normalized detuning J/Δ₆. This raises the question of whether the corresponding time scales are compatible with the finite lifetime of the dressed Rydberg atoms. It should be noted, however, that the relevant time scale is determined by $\Delta _6 = {\widetilde C}_6/a^6$ and thus depends on the lattice constant a. As mentioned in the introduction, the r⁻⁶ interaction potential is only valid beyond a cut-off length R_c which sets a constraint on a. On the other hand, the properties of the incompressible phases with filling q are only determined by the tails of the interaction potential at distances a/q and thus for our model to hold it is sufficient that a ⩾ R_cq. (For the same reason the smearing out of the interaction potential due to the convolution with the Wannier functions has little effect.) This means Δ₆ and thus the relevant time scale can be modified by adjusting the lattice spacing according to a ∼ R_cq, resulting in a scaling Δ₆ ∼ q⁻⁶. For example, taking a = R_cq, one finds that the melting points of the q = 1/2 and q = 1/6 phases are $J_{1/2} \approx 6\times {\widetilde C}_6/R_{\mathrm {c}}^6$ and $J_{1/6} \approx 0.5 \times {\widetilde C}_6/R_{\mathrm {c}}^6$. Since ${\widetilde C}_c/R_{\mathrm {c}}^6\sim \Omega ^2/\Delta$, where Ω and Δ are the effective Rabi frequencies and detunings of the Rydberg dressing, the common energy scale of the melting points can be tuned and can be in the kHz to MHz range. Thus melting can be observed also for small fillings well within the lifetime of Rydberg dressed states.

We have calculated the μ–J phase diagram for Rydberg dressed atoms in a deep 1D lattice potential. For vanishing hopping strength, J = 0, the chemical potential μ forms a complete devil's staircase as a function of the filling fraction corresponding to stable crystalline phases at any rational filling. Following a similar analysis of Bak and Bruinsma [23], we derived a closed analytic expression for any power-law interaction potential. We then analyzed the melting of the Rydberg crystals with increasing hopping. To this end, we performed a second-order strong-coupling expansion and found excellent agreement with exact numerical simulations based on DMRG adapted to long-range interactions. For the case of stronger localized interactions, such as a van der Waals coupling, a compact analytic expression for the phase boundaries of phases with filling q = 1/n was found.

We thank J Otterbach for inspiring discussions. This work was supported by the SFB TR49 of the Deutsche Forschungsgemeinschaft and the Graduate School of Excellence MAINZ.