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

We discuss the many-body physics of an ensemble of Rydberg dressed atoms with van der Waals dipole-dipole interactions in a one-dimensional lattice. Using a strong coupling expansion and numerical density-matrix renormalisation group simulations, we calculate the many-body phase diagram. A devil's staircase structure emerges with Mott-insulating phases at any rational filling fraction. Closed analytic expressions are given for the phase boundaries in second order of the tunnelling amplitude and shown to agree very well with the numerical results. The transition point where the incompressible phases melt due to the kinetic energy term depends strongly on the denominator of the filling fraction and varies over many orders of magnitude between different phases.


Introduction
Recently many-body systems with non-local, power-law interactions gained considerable interest as the 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,3,4,5]. For very strong lattice confinements the filling fraction forms a so-called complete devils staircase as 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 are Rydberg gases for which there has been considerable experimental progress in the recent years [7,8,9,10,11,12,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 extend 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 in general not conserve the number of Rydberg excitations, while in the case of Rydberg-dressing the number of atoms are to 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 tunnelling of the atoms can lead to a melting of the crystals [16,17].
In the present paper we analyse the melting process induced by 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 renormalisation 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 3 interactions in [6], however no explicit expressions were given. We here derive analytic expressions which in the case of fast decaying potentials, such as the 1/r 6 van der Waals potential, can in very good approximation be written in a compact closed form. Furthermore we provide for the first time a comparison to numerical data from DMRG simulations adapted to long-range interactions.
As a starting point we consider the generalised Bose-Hubbard model in one dimensionĤ For convenience we have set = a = 1, where a is the lattice constant.b i ,b † i are the annihilation and creation operators at site i, and J is the hopping amplitude between neighbouring lattice sites describing the tunnelling of atoms between adjacent lattice sites [18]. Most parts of our discussion are valid for any convex interaction potential, i.e However to be specific we discuss in the following power-law potentials of the form with β ∈ N\{1} andC β > 0 being the interaction coefficients. β = 6 corresponds to a van der Waals type interaction which results from virtual (i.e. off-resonant) longwavelength electromagnetic transitions in the manifold of Rydberg states. If there is an accidental or engineered resonance, a so called Förster resonance, also the case β = 3 can be of relevance. The model interaction potential (3) diverges for r → 0 and thus two particles cannot sit on the same lattice site. To 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] Here C 6 is the Rydberg interaction coefficient and R c = (C 6 /(2∆)) 1/6 describes a characteristic length scale of the interaction potential, where ∆ |Ω| is an effective detuning of the off-resonant laser driving. η ∼ Ω 2 /∆ 2 1 accounts for the small admixture of the Rydberg state. Considering e.g. 87 Rb with a Rydberg state of principle quantum number around n = 60 and detunings of a few tens of MHz one finds a value of R c on the order of one 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 devided 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 reasonb i ,b † i are assumed to describe hard-core-bosons with commutation relations 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 β. 2. Ground state for J = 0 -classical limit a.marueg@web.de 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 minimised just 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 = m n ≤ 1, where m, n ∈ 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.
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, 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 ≤ 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 we will make use of their notation: X 0 i denotes the position of the i-th particle.
describes the distance between the p-th and the i-th particle where the particles are numbered from the left to the right. In the ground state ofĤ 0 Hubbard's solution [22] shows that all possible distances are given by where r p < p q < r p + 1. If p q ∈ N then X p i = r p = p q and all particles are equivalent. all nearest-neighbour separations are maximally close to the average separation, i.e. they are either 1 q or 1 q . The chemical potential for J = 0 is given by |q denotes the ground state ofĤ 0 for filling fractions q, whilst |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 = m n there will be n defects for |q ± , because it is energetically favourable to break up an extra particle (or a hole) into n defects each with fractional charge ± 1 n . Some of these defects are displayed in figure 2. In the thermodynamical limit they will be separated by arbitrarily large 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 = 1 n the chemical potential is given by [23] µ (0) We can go a step further and evaluate (12) for the power-law potential (3) analytically, using the digamma function Ψ(z) = Γ (z) Γ(z) and its l-th derivative Ψ (l) (z). To describe more general filling fractions q = m n with m = 1 as well, we only have to replace r p → r p + 1 for µ 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 devils staircase. The union of these intervals covers the full range of µ. It is clearly seen that filling fractions with small denominator n are the 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 = 1 2 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 µ (0) ≡ (µ As mentioned in the introduction the true interaction potential between Rydbergdressed 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 only the case if the filling q becomes large and in particular exceeds the ratio R c /a.

Perturbation theory in J
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 3 [6], but no compact analytic expression was given.
Let us start with first order processes in J. Then, the chemical potential reads µ (1) ± (q) = ± q ± |Ĥ J |q ± − q|Ĥ J |q . 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 localised defects are degenerated with respect toĤ 0 , we have to apply a) b) Figure 4. Hopping of defects for q = 1 2 , q = 1 3 and q = 2 5 . The black arrow indicates the hopping of one particle and the red one corresponds to the resulting hopping of the defect. The blue ones are the alternative hopping-possibility, which would result in a hopping of the defect in the opposite direction. a) b) Figure 5. All possible second order processes for |q + (a) and |q − (b) at filling q = 1 n , here q = 1 3 . (i) generates an effective chemical potential for single particles that depends on the distance to the defect, i.e., the defect polarises the background. (ii) refers to a virtual deformation of the defect. (iii) corresponds to an effective hopping of the defects over two unit cells.
degenerate perturbation theory. Therefore, we look for a basis of these states in whicĥ H J is diagonal. This is given by states where all the defects are delocalised over many lattice sites (yet well separated from each other) with some quasi-momentum k. The energy is minimised for a state with quasi-momentum zero. The resulting first order correction to the phase boundary thus takes the simple form µ (1) For small but nonzero J gaps open between the incompressible phases of rational filling, because the kinetic energy gained by delocalisation of defects favours 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 favourable 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 devils staircase.a.marueg@web.de In second order perturbation the energy corrections are given by where |Ψ denotes the ground state ofĤ 0 and {|Φ } is a complete set of orthonormal states. In the following we consider only filling fractions of the form q = 1 n . In principle it is posa.marueg@web.desible to extend the calculations to other fractions. However, finding a general formula which describes all energy differences in the denominator for all hopping processes is a far more tricky task. In figure 5 all possible processes for q ± = 1 n ± are shown. As has been seen in first order, the states with localised defects are degenerate with respect toĤ 0 . In degenerate second order perturbation theory the energy correction is given by q ± |Ĥ JPĤJ |q ± /(E (0) whereP projects out the subspace of ground states ofĤ 0 . The state |q ± has to be an eigenstate ofĤ J restricted to the degenerate subspace. Therefore, intermediate states 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 where N is the number of particles in the lattice.
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 In order to obtain the (finite) chemical potential we have to calculate also the energy corrections E (2) (q ± ) for the states with one extra particle or hole |q ± in second order perturbation theory. To this end let us have a closer look 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 where ∆E (0)± +j,− denotes the energy difference between the ground state |q ± and the same state, where the j-th particle right of the defect moves right (second subscript +) or left (second subscript −). For the meaning of position "+0" see figure 5. We have assumed here that the defect is in the centre of the system and thus the summation limit is N ± 0 = 1/2[(N ± 1)q − 2]. As |q ± has to be an eigenstate ofĤ J in the degenerate subspace the defects are completely delocalised. As in first order it turns out that the ground state has quasi-momentum k = 0. Then, the second order correction to the chemical potential can be written as 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 less particles in |q ± taking part in process (i) then 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 (0) , the chemical potential for filling fractions q = 1 n can finally be written in second order in J as where the terms S ± j (n) are given by If the power β of the interaction potential is sufficiently large and if the effect of the cutoff length R c can be disregarded, i.e. if R c < a/q, the entire sum j S ± 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 µ (2)

Exact numerical calculation
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 renormalisation group algorithm [24]. It is a variational technique that minimises the energy of the full HamiltonianĤ 0 +Ĥ J using a matrix product state (MPS) ansatz [25]. Although limited in the amount of entanglement along the lattice that they can capture, MPS are known [26] to approach the true ground state of one-dimensional 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 priory, which comes in handy to calculate 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 = m n can be estimated to be 2n 2 , 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 a 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/∆ 6 . 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.  Figure 6. Phase diagram for a van der Waals interaction potential. The dashed lines correspond to perturbation theory up to second order in J/∆ 6 , with ∆ 6 =C 6 /a 6 , continuous lines are calculated with the DMRG method using MPSs of dimension 32. Coloured lines are infinite size extrapolations for L 24 (n = 2, 3 only), L 60, and L 120 lattice sites. No infinite size extrapolation is applied for n = 6 because we decided not to make the effort to calculate large enough system sizes. For all cases interactions over distances larger than r = 7 lattice sites. a) Double logarithmic plot. The finite size results are shown in grey for q = 1 3 and q = 2 5 for illustration. b) Linear plot. Note the vast difference in size for different n (13), which is much more drastic than in the β = 3 case [6].
The continuous lines in the figure end at arbitrary values, while the phase boundary ends where these lines make contact.
As can be seen from fig.6 the melting points of the different incompressible phases differ by many orders of magnitude in the normalized detuning J/∆ 6 . This raises the question if the corresponding timescales are compatible with the finite lifetime of the dressed Rydberg atoms. It should be noted however that the relevant time scale is determined by ∆ 6 = C 6 /a 6 and thus depends on the lattice constant a. As mentioned in the introduction the r −6 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 c q. (For the same reason the smearing out of the interaction potential due to the convolution with the Wannier functions has little effect.) This means ∆ 6 and thus the relevant time scale can be modified by adjusting the lattice spacing according to a ∼ R c q resulting in a scaling ∆ 6 ∼ q −6 . For example taking a = R c q one finds that the melting points of the q = 1/2 and q = 1/6 phases are J 1/2 ≈ 6 × C 6 /R 6 c and J 1/6 ≈ 0.5 × C 6 /R 6 c . Since C c /R 6 c ∼ Ω 2 /∆, where Ω and ∆ are the effective Rabi-frequencies and detunings of the Rydberg dressing, the common energy scale of the meltig 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.

Summary
We have calculated the µ − J phase diagram for Rydberg dressed atoms in a deep one-dimensional lattice potential. For vanishing hopping strength, J = 0, the chemical potential µ forms a complete devils staircase as function of the filling fraction corresponding to stable crystalline phases at any rational filling. Following a similar analysis of reference [23] we derived a closed analytic expression for any power-law interaction potential. We then analysed the melting of the Rydberg crystals with increasing hopping. To this end we performed a second order strong-coupling expansion and found excellent agreement to exact numerical simulations based on DMRG adapted to long-range interactions. For the case of stronger localised 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.