# Quantum crystals in a trapped Rydberg-dressed Bose-Einstein condensate

## Abstract

Spontaneously crystalline ground states, called quantum crystals, of a trapped Rydberg-dressed Bose-Einstein condensate are numerically investigated. As a result described by a mean-field order parameter, such states simultaneously possess crystalline and superfluid properties. A hexagonal droplet lattice is observed in a quasi-two-dimensional system when dressing interaction is sufficiently strong. Onset of these states is characterized by a drastic drop of the non-classical rotational inertia proposed by Leggett [Phys. Rev. Lett. 25, 1543 (1970)]. In addition, an AB stacking bilayer lattice can also be attained. Due to an anisotropic interaction possibly induced by an external electric field, transition from a hexagonal to a nearly square droplet lattice is also observed.

Superfluidity which implies a long-range phase coherence is a crucial property at low temperatures of many quantum liquids or gases such as liquid Helium or Bose-Einstein condensates (BECs), whereas crystallization implies a long-range configurational order. Superfluidity and crystallization are generally two conflicting properties. Penrose and Onsager (1) were the first to consider a BEC in a solid and concluded that such a supersolid state simultaneously possessing crystalline and superfluid properties was impossible. Since then, this question has been revisited by a number of authors (2); (3); (4) and has been a matter of large speculation for the last forty years. Recently, the observation of a supersolid phase in He systems (5) revitalized this fundamental interest.

An alternative and excellent candidate to study supersolidity is in atomic
BECs which provide a clean and experimentally controllable system. The
crystal structure in solid helium can be replaced by the modulated density
in BEC. Density modulated BECs are already formed by the imposition of an
external potential, creating the so-called optical lattices. In these
systems, by varying the properties of the optical lattice, the condensate
was shown to exhibit a Mott insulator/superfluid phase transition (6); (7); (8); (9). More recently, it has been shown that supersolidity might
be present for Rydberg atoms in the dipole/van der Waals (vdW) blockade
regime (10); (11); (12); (13). Cinti et al.(11) considered a dipole-dipole interaction softening at short
distance, allowing for a ground-state computation that happens to display
the properties of supersolidity. It proved that a quantum system of
interacting particles can exhibit both crystalline structure and
superfluidity property. Similar results were obtained by Saccani *et al*.(12) by using a Heaviside-function interaction. Based on a
mean-field treatment, Henkel et al.(13) proposed that a BEC of
particles interacting through an isotropically repulsive vdW interaction
with a softened core might support a density modulation. They found that the
Fourier transform of such interaction has a partial attraction in momentum
space, which gives rise to a transition from a homogeneous BEC to a
supersolid phase because of the roton instability (see also Refs. (14); (15)).

Based on the Gross-Pitaevskii (GP) treatment, this Letter aims to study the ground states of a trapped Rydberg-dressed BEC. Comparing to other GP work that did not consider the effect of trapping (13), we exactly solve the nonlocal GPE with a trap. In particular, we focus on the quasi-2D geometry and the strong interacting regime that lead to various fascinating ground-state structures of Rydberg-dressed BEC. It will be shown that in the supersolid phase, periodic structures of Rydberg-dressed BEC can undergo a transition from concentric rings to a lattice (or crystalline) if the interaction is above some critical value. The lattices are formed in terms of crystalline superfluid droplets, called quantum crystals, whose onset is characterized by a drastic drop of the non-classical rotational inertia fraction (NCRIF) (4). The crystalline structure appears to be a hexagonal lattice in a quasi-two-dimensional (quasi-2D) geometry which can turn into a nearly square lattice if interaction acquires an anisotropic component in the presence of an external electric field (Stark effect). Moreover, multilayer crystal structure such as an AB stacking bilayer is also obtained when the frozen axis is relaxed or particle number is increased. Crystalline in the case of a quasi-one-dimensional (quasi-1D) geometry will also be presented.

We start from a nonlocal Gross-Pitaevskii equation (GPE)

(1) | |||||

where is an isotropically repulsive vdW interaction between Rydberg-dressed ground-state atoms (16) with and the effective coupling constant and blockade radius respectively (we will return to the interaction later). Here , normalized as ( is the atom number), denotes the condensate wave function of dressed atoms. denotes the atomic mass and describes the strength of local interaction due to -wave scattering with scattering length . In the cylindrical coordinates, (), the harmonic trapping potential with radius frequency and aspect ratio is included for simulating real experiments.

By introducing useful length () and time () scales, Eq. (1) can be rewritten as

(2) |

where we have redefined the normalized wave function , the strength of the radius potential , and the interaction constants and . To obtain ground-state wave functions, we computed the governing Eq. (2) with imaginary time propagation till the convergence of the normalized wave function with error less than . Moreover, we have used the method of lines with spatial discretization by the Fourier pseudospectral method. The time integration in Eq. (2) is done by the adaptive Runge-Kutta method of order 2 and 3 (RK23), which is more time efficient due to an adjustable time step.

For the Rydberg-dressed BEC trapped in a harmonic potential, one can define useful characteristic lengths, and , corresponding to radial and axial potential, respectively. The spectrum and the onset of instability are tunable by varying the particle number or the confining potential. By varying the two ratios and , one can effectively have quasi-1D ( and ) or quasi-2D ( and ) limits.

A quasi-2D system with and is first studied. Fig. 1 shows the condensate density profile varying with the strength of dressing-induced interaction . When is small, the system displays superfluidity, and the ground-state density profile exhibits a central peak [see Fig. 1(a)]. As increases, the central density is too high to be stable and thus starts to modulate. Fig. 1(b) shows a cratered condensate due to the central instability. As increases further, owing to the roton instability occurring in the modulated density (see later), condensate wrinkles violently and forms a ring structure [see Fig. 1(c)]. When is increased above a critical value , condensate eventually forms a droplet lattice. Fig. 1(d) shows the quasi-2D Rydberg-dressed BEC forming a hexagonal droplet lattice.

It is important to check whether the parameter regime discussed above is actually experimentally accessible. The vdW interaction can be generated from the off-resonant dressing of ground-state atoms with high-lying Rydberg states. Considering, for example, ground-state Rb atoms coupled to excited Rydberg state Rb atoms with via a Rabi frequency and a red laser detunning , it will admix a small fraction of Rydberg character into the ground-state atoms for weak dressing (). For two far-distant atoms, it leads to an effective interaction (), arising from the strong vdW interaction ( a.u.) (17). At shorter distances, the two atoms enter the vdW blockade regime (18) to which the effective interaction saturates. Altogether, as given earlier, the effective potential between Rydberg-dressed ground-state atoms is with the blockade radius (13). As shown, the above off-resonant scheme produces only repulsive nonlinearities for alkaline atoms (19). Using typical value of Rabi frequency KHz and a red detunning MHz, blockade radius m. Moreover, the effective lifetime of dressed atoms, , is as large as several seconds with Rydberg state decaying rate ms (16). To achieve the largest coupling constant that we are considering, , based on the above parameters one needs the total atom number which corresponds to for the number of excited Rydberg atoms. By counting the number of droplets in Fig. 1(d), an average of one excited Rydberg atom together with ground-state atoms is within a single droplet. This justifies the validity of the GP treatment with the two-body dressing interaction (20).

A qualitative understanding of forming the ring-like structures [Fig. 1(c)] in a quasi-2D system is given in the uniform 2D limit. Considering that the system is strongly confined in -direction by a harmonic potential, but free moving in the plane, superfluid BEC density can be approximated by . The 2D excitation spectrum calculated from the corresponding Bogoliubov-de Gennes equations is thus

(3) |

where (in units of ) is momentum component in the plane and with and () being Fourier transforms of the product of and the scaled interaction in Eq. (2) (21). When deriving Eq. (3) and throughout this Letter, contact interaction is set to zero. Including will not affect the behaviors of roton instability if is sufficiently large, nor will affect the phonon behavior if is small to which leading term of the expansion in is a positive constant, i.e., nonzero only modifies the sound velocity. has asymptotically a phonon and a free-particle character at small and large , respectively. However, with having a negative minimum at some finite momentum, drops near that particular momentum, and eventually becomes imaginary when increasing the strength [see Fig. 1(e)]. It is estimated that when , the assumed uniform superfluid state is unstable towards formation of nonuniform (periodic ring-like) structures.

To characterize the transition from concentric rings to a crystalline hexagonal lattice, we study the non-classical rotational inertia fraction (NCRIF), defined by . Here is moment of inertia of the superfluid system under study and is its corresponding classical value (5). As proposed by Leggett (4), NCRIF of the superfluid system can be calculated under a small rotation. In the rotating frame, free energy of a rotating BEC with rotation velocity about axis is

(4) | |||||

where is the -component angular momentum operator and is the free energy of the system without rotation. When , can be expanded as with , taken to be real, being the ground state of . Since classical moment of inertia is given by , we obtain for ,

(5) |

where is the ground state of . In arriving (5), we have assumed for . Therefore NCRIF can be obtained by computing Eq. (5) with the solved and . Fig. 1(f) plots the calculated NCRIF as a function of the strength . It is evident that onset of crystallization of the BEC droplets is characterized by the drastic drop of NCRIF, occurring at . No similar drop appears for the onset of the ring-like supersolid phase occurring at smaller as a result of full rotational symmetry.

By relaxing the originally frozen direction potential and/or increasing the particle number, density starts to modulate in direction and eventually forms a multilayer structure. Figs. 2(a) and 2(b) compare formations of a monolayer and a bilayer lattice. The inset in Fig. 2(b) indicates clearly that such a bilayer structure is an AB stack.

It is also interesting to note that when Rydberg dressing interaction becomes anisotropic, hexagonal lattice can shift to a nearly square lattice due to distortion of the interaction. This can occur when an external static electric field is applied to the system (Stark effect) to which two-photon mechanism will acquire an anisotropic component for the interaction, as compared to the purely isotropic case with (22). Fig. 3 shows a nearly (though not perfectly) square lattice by using an interaction of the Heaviside-function form instead of as in Eq. (2). Here corresponds to the anisotropic ratio. The static electric field is considered applied along the axis. Without losing the generality, while Heaviside-function interaction make the simulation of anisotropy more conveniently, it does capture a roton minimum in the excitation spectrum.

The condensate ground state in a quasi-1D system is also investigated with and . Fig. 4 shows the modulation of condensate density of a quasi-1D system by varying the strength of dressing-induced interaction . When is small, the system displays superfluidity, and the density has a central peak [see Fig. 4(a)]. As increases, the density starts to modulate and has multiple peaks. Fig. 4(b) shows that there are five peaks in the condensate, along the axial direction. With a sufficiently large , the condensate spontaneously crystallizes in the axial direction [see Fig. 4(c)]. As increases further, the condensate starts to modulate in the frozen direction, forming a crystalline structure with a central hole [see Figs. 4(d) and 4(e)]. If is extremely large, the condensate starts to cluster in the originally frozen direction, and these clusters form a gyroidal chain.

In summary, based on the Gross-Pitaevskii treatment, spontaneously crystalline ground states, called quantum crystals, are numerically studied for a trapped Rydberg-dressed Bose-Einstein condensate. In a quasi-2D system, a hexagonal droplet lattice characterized by a drastic drop of the non-classical rotational inertia is shown when dressing interaction is sufficiently large. By relaxing the originally frozen axis, an AB stacking bilayer lattice is observed. We also show that by applying a static electric field to make the interaction anisotropic, a nearly square droplet lattice can be obtained.

We are grateful to Yu-Ching Tsai for many helpful discussions. This work was supported by National Science Council of Taiwan (Grant No. 98-2112-M-018-001-MY2). We also acknowledge the support from the National Center for Theoretical Sciences, Taiwan.

### References

- O. Penrose and L. Onsager, Phys. Rev. 104, 576 (1956).
- A. F. Andreev and I. M. Lifshitz, Sov. Phys. JETP 29, 1107 (1969).
- G. V. Chester, Phys. Rev A 2, 256 (1970).
- A. J. Leggett, Phys. Rev. Lett. 25, 1543 (1970).
- E. Kim and M. H. W. Chan, Nature (London) 427, 225 (2004); Science 305, 1921 (2004).
- M. Greiner, O. Mandel, T. Esskinger, T. W. Hänsch, and I. Bloch, Nature (London) 415, 39 (2002).
- O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006).
- K. Goral, L. Santos, and M. Lewenstein, Phys. Rev. Lett. 88, 170406 (2002).
- T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Rep. Prog. Phys. 72, 126401 (2009).
- G. Pupillo, A. Micheli, M. Boninsegni, I. Lesanovsky, and P. Zoller, Phys. Rev. Lett. 104, 223002 (2010).
- F. Cinti, P. Jain, M. Boninsegni, A. Micheli, P. Zoller, and G. Pupillo, Phys. Rev. Lett. 105, 135301 (2010).
- S. Saccani, S. Moroni, and M. Boninsegni, Phys. Rev. B 83, 092506 (2011).
- N. Henkel, R. Nath, and T. Pohl, Phys. Rev. Lett. 104, 195302 (2010).
- Y. Pomeau and S. Rica, Phys. Rev. Lett. 72, 2426 (1994).
- X. Li, W. V. Liu, and C. Lin, Phys. Rev A 83, 021602(R) (2011).
- J. E. Johnson and S. L. Rolston, Phys. Rev. A 82, 033412 (2010).
- K. Singer et al., J. Phys. B 38, S295 (2005).
- D. Jaksch et al., Phys. Rev. Lett. 85, 2208 (2000).
- F. Maucher et al., Phys. Rev. Lett. 106, 170401 (2011).
- J. Honer, H. Weimer, T. Pfau, and H. P. Büchler, Phys. Rev. Lett. 105, 160404 (2010).
- is Fourier transform of . As the softened range of is much longer than the width of , i.e., , the integration is well approximated by . Consequently with the Meijer G function.
- A. V. Gorshkov et al., Phys. Rev. Lett. 101, 073201 (2008).