Observation of anisotropic superfluid density in an artificial crystal

We experimentally and theoretically investigate the anisotropic speed of sound of an atomic superfluid (SF) Bose-Einstein condensate in a 1D optical lattice. Because the speed of sound derives from the SF density, this implies that the SF density is itself anisotropic. We find that the speed of sound is decreased by the optical lattice, and the SF density is concomitantly reduced. This reduction is accompanied by the appearance of a normal fluid in the purely Bose condensed phase. The reduction in SF density -- first predicted [A. J. Leggett, Phys. Rev. Lett. 1543--1546 (1970)] in the context of supersolidity -- results from the coexistence of superfluidity and density modulations, but is agnostic about the origin of the modulations. We additionally measure the moment of inertia of the system in a scissors mode experiment, demonstrating the existence of rotational flow. As such we shed light on some supersolid properties using imposed, rather than spontaneously formed, density-order.

Superfluidity and Bose-Einstein condensation (BEC) are deeply connected.In dilute atomic BECs, the superfluid (SF) and condensate densities are generally equal [1,2].By contrast, SF 4 He can be a nearly pure SF, with only about 14 % condensate fraction [3], and infinite 2D Berezinskii-Kosterlitz-Thouless (BKT) SFs have no condensate at all [4,5].In 1970 Tony Leggett showed that supersolids-systems spontaneously forming both SF and crystalline order (i.e., density modulations)exhibit the reverse behavior: SF density far below the condensate density [6].Here we observe this effect in a nearly pure atomic BEC with artificial crystal order imprinted by an optical lattice.
The complex-valued order parameter φ(r) = ρ sf exp[iϕ(r)], describing a SF with number density ρ sf and phase ϕ(r), gives rise to two hallmark SF properties: dissipationless supercurrents associated with spatial gradients in ϕ(r) and (Bogoliubov [2]) sound described by traveling waves in ϕ(r).Because supercurrents arise from phase gradients, they are locally irrotational; in liquid 4 He, the resulting non-classical rotational inertia appears below the SF transition temperature T c .Supersolids are more exotic systems spontaneously forming crystalline order while exhibiting SF transport properties.
Recent experiments with dipolar BECs of Dy and Er are suggestive of these properties [7,8].Leggett argued that the modulated density ρ(r) of a supersolid leads to an unavoidable reduction in ρ sf , and derived an upper bound for ρ sf [6].This reduction results from the 3D density distribution, and as such is masked in tight binding descriptions such as the Bose-Hubbard model, which makes the unrelated prediction of vanishing ρ sf at the superfluid to Mott insulator transition [9,10].
We created an artificial SF crystal by imprinting periodic density modulations into an atomic BEC using a 1D optical lattice as in Fig. 1(a).While these modulations do not form spontaneously, Leggett's result still applies, making this an ideal system for understanding crystalline SFs without the added complexity of spontaneously broken symmetries.We experimentally measured an anisotropic speed of sound via Bragg spectroscopy of the phonon mode.This implies the existence of an effective anisotropic superfluid density-which can be expressed as a second rank tensor ρ sf ij -and we find that it saturates Leggett's bound, in agreement with Gross-Pitaveskii equation (GPE) simulations.We also determined an associated anisotropic suppression of the moment of inertia in terms of the scissor-mode frequencies [11,12].
Even disregarding potential differences in ρ sf (r) and ρ(r), we argue that φ(r) is not simply equal to ψ(r) averaged over some scale large compared to a.The fundamental origin of this effect can be understood by considering a 1D system of size L with periodic boundary conditions in which both the condensate phase ϑ and SF phase ϕ wind by an integer multiple N of 2π [Fig.1(b,c)ii], yielding a metastable quantized supercurrent [13].
To satisfy the steady-state continuity equation, the microscopic current
v/v iv.

FIG. 1. Concept. (a)
A BEC is confined in a harmonic trap superimposed with a 1D optical lattice (along ex, green), spatially modulating the condensate density (red).The dashed and dotted lines call out a region of nominally constant mean density and the left and right columns indicate the (b) state of the condensate and (c) SF in the presence of a current.These were computed for a 5Er deep lattice and plot: i. density (red), ii.current (green), iii.phase (orange), and iv.local velocity (blue).The red dashed line plots the mean density ρ.
mode and the SF order parameter and integrating over a unit cell yields Leggett's equation [6] This implies that ρ sf ≤ ρ, where ρ is the spatial average of the condensate density over a UC, and as we discuss below the remaining density ρ n = ρ − ρ sf behaves as a pseudo-normal fluid.In the more general context where the GPE is inapplicable, the Leggett expression for ρ sf is an upper bound for the SF density in systems with crystalline order [6].In a 3D system, the current ] derives from a SF density tensor.For systems with rectangular symmetry [14] ρ sf ij is diagonal, and the analogs to Eq. ( 1) for each of the three elements use a 1D density integrated along the transverse directions.In our experiments this implies that the superfluid density is only reduced along the direction of the optical lattice, so ρ sf yy = ρ sf zz = ρ.Experiment-We used 87 Rb BECs with N ≈ 2 × 10 5 atoms in the |F = 1, m F = 1 hyperfine ground state.A 1064 nm trapping laser with an elliptical cross-section, traveling along e x provided strong vertical confinement with frequency ω z /(2π) = 220 Hz; the in-plane frequencies, from ω x,y /(2π) = (34, 51) Hz to (56, 36) Hz, were optimized for our different experiments.We created a 1D optical lattice using a retro-reflected λ = 532 nm laser traveling along e x , giving an a = 266 nm lattice pe-riod, comparable to the ξ = 170 (20) nm minimum healing length.The optical lattice was linearly ramped on in 100 ms to a final depth ≤ 10 E r , with single photon recoil energy and momentum E r = 2 k 2 r /(2m), and k r = 2π /λ respectively [15].For Bragg experiments the final state was measured using resonant absorption imaging after a 15 ms time of flight (TOF); scissors mode measurements were performed in-situ using partial transfer absorption imaging [16].
Anisotropic speed of sound-The speed of sound for diagonal ρ sf ij is predicted to result from c 2 i = f sf ii /(κm) in terms of the superfluid fractions f sf ii = ρ sf ij /ρ, the compressibility κ = ρ−1 (∂ ρ/∂µ), and the chemical potential µ.This reduces to the well-known value c 2 = µ/m for an isotropic homogeneous system (See [17] for the full dispersion beyond the linear approximation).The sound speed ratio provides direct access to the different components of the superfluid density [see [17] for a Josephson sum rule [18] argument].Because the density is y-independent, Eq. ( 1) implies ρ sf yy = ρ.We Bragg-scattered the BEC off a weak sinusoidal potential with reciprocal lattice vector δk slowly moving with velocity v by patterning a laser beam with a digital micro-mirror device (DMD [19]) and measured the scattered fraction p.This results from what are effectively two interfering laser beams driving two-photon transitions with difference-wavevector δk and angular frequency δω = δk v.We applied this potential for ≈ 5 ms.Bragg transitions ensued when the difference energy and momentum were resonant with the BEC's Bogoliubov dispersion, and Fig. 2(a) shows data in the linear regime.The width of this spectral feature is limited by our BEC's inhomogeneous density profile; the resonance (vertical dashed line) obtained from a Lorentzian fit (solid curve) therefore reflects an average speed of sound [20].
A series of such fits lead to phonon dispersion relations with Bragg-lattice period from 2.25 µm to 8.5 µm.Representative dispersions taken along e x and e y are shown in Fig. 2(b), and we obtain the phonon speed of sound using linear fits.Figure 2(c) summarizes these data showing the speed of sound decreasing along the lattice direction e x , but slightly increasing along e y (resulting from the increased atomic density in the individual lattice sites).Finally Fig. 2(d) shows our main result: the normalized superfluid density obtained from these data using Eq. ( 2) decreases as a function of U 0 .
We compared these data to GPE simulations in two ways, we: (1) used the Bogoliubov-de Gennes (BdG) equations to obtain c x and c y and (2) directly evaluated Eq. ( 1) from the GPE ground state density.The solid curves in Fig. 2(c (no free parameters [17]).(d) SF density obtained from speed of sound measurements (blue markers, error bars mark singlesigma statistical uncertainties).We compare with two models: the red dashed curve plots a homogeneous gas BdG calculation, and the solid black curve plots the result of Eq. ( 1).The simulations used our calibrated experimental parameters.
solving the 1D BdG [21], and the red dashed curve in (d) is the ratio of these speeds.To compare with Leggett's prediction, we found the ground state of the 2D GPE for our experimental parameters and evaluated Eq. ( 1) throughout our inhomogeneous system.The black curve in Fig. 2 plots the resulting weighted average.Remarkably the BdG results are in near-perfect agreement with Leggett's expression.
Scissors mode-The single-valued nature of the SF order parameter greatly affects rotational properties such as the moment of inertia I.For highly anisotropic traps, the scissors mode [11] describes a fixed density distribution pivoting by a small angle θ about the trap center with frequency ω sc .Scissors mode experiments are reminiscent of torsional balance experiments in 4 He [22] which give access to the non-classical rotational inertia [6].
It is suggestive to quantify these dynamics in terms of the Lagrangian L = I θ2 /2 − V (θ), for moment of inertia I and potential energy V (θ).For small θ the potential can be expanded as V (θ) ≈ Iω 2 sc θ 2 /2 with in terms of the classical moment of inertia I c and in agreement with Ref. [11] for isotropic superfluids.Although this interpretation is highly intuitive, it does not survive careful consideration.The anisotropic superfluid density couples radial and azimuthal flow and as a result a single parameter Lagrangian is insufficient to describe rotational dynamics.Instead the superfluid hydrodynamic equations predict a moment of inertia scaled by a factor of (f sf xx ω 2 x − f sf yy ω 2 y )/(ω 2 x − ω 2 y ) (see [17]) compared to Eq. ( 3).Therefore we expect ω sc , in conjunction with the superfluid density will give I/I c as a function of lattice depth.
The inset to Fig. 3(a) plots the dipole mode frequencies ω x,d and ω y,d for a trap with frequencies (54, 36) Hz.The frequency reduction is also related to ρ sf via f sf = (ω x,d /ω x ) 2 along the lattice direction [17].This ratio can also be expressed in terms of an increased effective mass that converges to the predictions of single-particle band structure [23] when the lattice period falls below the healing length; in our case the value computed perturbatively from the GPE differs by about 20 % from the band structure prediction.The result of this modeling is shown by the solid curves.
We excited the scissors mode using our DMD to tilt the harmonic potential by 50 to 140 mrad for ≈ 1 ms (shorter than the trap periods) and let the BEC evolve in the original trap for a variable time.We measured the resulting dynamics in-situ and extracted the angle by fitting the resulting density profile to a rotated Gaussian.Figure 3(a) shows scissor mode frequency normalized to the expected frequency [24] of ω 2 sc,0 = f sf xx ω 2 x + f sf yy ω 2 y for a trap elongated either along e x [with frequencies (56, 36) Hz, blue] or along e y [with frequencies (36, 50) Hz, green].In both cases ω sc is about 5 % in excess of the simple prediction, perhaps from finite temperature or anharmonicities in the ODT.
We combine these observations in Fig. 3(b) to obtain I/I c ; the data (symbols) and our 2D GPE simulations (curves) are in agreement.For traps elongated along e x (green) I/I c unexpectedly changes sign when ω x,d = ω y,d .To understand the physical origin of this effect we now turn our attention to rotating systems.
Rotation-Thus far we focused exclusively on the superfluid density, while avoiding questions about any associated normal fluid.We can deduce the existence of a normal fluid component by considering two thought experiments in a 1D ring geometry (with radius R) and quantify both in terms of the resulting angular momentum [25].In case (i), we consider an Aharonov-Bohm geometry and slowly thread the ring with a single quanta of magnetic flux (see Ref. [26] for an artificial gauge field proposal).The process is equivalent to imprinting a 2π phase winding (of the type discussed on page 1), giving angular velocity Ω = /(mR 2 ) and angular momentum L z / = 2πRρ sf .In case (ii), we consider a complimentary experiment in which the lattice is very slowly accelerated to a final angular velocity Ω; this is best understood by transforming into the frame co-rotating with the lattice.This leads to a lab frame angular momentum L z / = 2πR(ρ−ρ sf ) which we interpret as resulting from the normal fluid co-moving with the lattice.
With this insight we extended our 2D numerical simulations to analogous cases for rotating harmonically trapped systems where : (i) the lattice is static in the lab frame (as in scissors mode experiments) or (ii) it corotates with the confining potential.In both cases we use the coarse graining defined in Eq. ( 1) to obtain the superfluid density and phase.In this way we compute the total moment of inertia I from ψ(r, t), the superfluid component I sf from φ(r, t), and we define the normal component as the difference I n = I − I sf .
Case (i): as in our 1D thought experiment only the SF component responds.Then although ∇ϕ is manifestly irrotational, because ρ sf xx = ρ sf yy the superfluid current can be rotational.In this case, the relative magnitude of the co-and counter-rotating contributions vary with the lattice depth, leading to regions of negative angular momentum density L(r) along the BEC's semi-minor axis [Fig.4(a)].The superfluid moment of inertia computed from these simulations [Fig.4(c)] is in full agreement with the scissor mode simulation, and as expected for a static lattice I sf = I (no normal flow).
When the lattice is along the semi-minor axis, as pictured in (a) and the green curve in (c), the counterrotating contribution increases with U 0 , until the dipole mode frequencies along e x and e y invert, after which point, I/I c becomes negative.The reverse is the case when the lattice is along the semi-major axis and I/I c increases monotonically.This novel observation confirms the negative kinetic energy resulting from θ.
Case (ii): In contrast, the angular momentum density is strictly positive [Fig.4(b)] for both lattice orientations and I/I c increases with lattice depth [Fig.4(d)].In this case the normal fluid to co-rotate with the trap giving the current J n = (−ρ n xx y, ρ n yy x) θ.The total I/I c is then the sum of the superfluid [17] and normal contribution This result, along with our 2D GPE simulations, are plotted in Fig. 4(d).The dashed curve plots the superfluid contribution to I sf /I c in agreement with the coarsegrained GPE (crosses).The solid curve and the triangles plot the corresponding total moment of inertia, in excess of the SF contribution.This implies the appearance of normal fluid flow.This agreement confirms that the superfluid contribution derives from gradients of the coarse-grained phase ϕ, while the normal contribution stems from variations of ϑ within each lattice site.
Discussion and outlook-Our inability to obtain I/I c from scissors mode measurements without detailed modeling reinforces similar conclusions in dipolar gases [27].In both cases the simple argument fails because θ couples to more internal degrees of freedom than L z alone.In this context Ref. [27] concluded that the scissors mode does yield the moment of inertia when 1D density modulations comove with the oscillatory motion: this is consistent with our findings comparing motion in static and rotating lattices.Our GPE simulations indicate that the analytical relations generalize to lattices with period in excess of the healing length.
Although we conclude that a normal fluid exists, it is inseparable from the optical lattice and lacks any internal dynamics of its own, i.e., it is not described by a dynamical equation of motion.In contrast, both the superstripe phase in spin-orbit coupled BECs [28][29][30][31][32] and supersolid phases of dipolar gases [33][34][35], support dynamical density modulations.Leggett's expression applies to both of these systems implying a reduced superfluid density, which in this case could exhibit dynamics, as expected for a system described by a two-fluid model [11,12].
This leaves open questions regarding nature the normal fluid of spin-orbit coupled systems where an interplay between single-particle physics and interactions govern supersolid-like properties.In addition, ρ sf is expected to be reduced outside of the superstripe phase [31,32] where the density is uniform (making Leggett's expression inapplicable), but the BEC's spin vector is spatially periodic.
Note: During the early stages of manuscript preparation we become aware of a related work, using a long period 1D lattice applied to a homogeneously confined 2D BEC.

FIG. 2 .
FIG. 2. Bragg spectroscopy.Black and red symbols mark excitations created along ex and ey respectively.(a) Transferred population fraction p as a function of frequency difference δω with wavevetor δk/2π = 0.26 µm −1 and lattice depth U0 = 5.7Er.The solid curve is a Lorentzian fit giving the resonance frequency marked by the vertical dashed line.(b) Phonon dispersion obtained from Bragg spectra.The bold symbols resulted from (a) and the linear fit (with zero intercept) gives the speed of sound.(c) Anisotropic speed of sound.The bold symbols are derived from (b) and the solid curves are from BdG simulations (no free parameters[17]).(d) SF density obtained from speed of sound measurements (blue markers, error bars mark singlesigma statistical uncertainties).We compare with two models: the red dashed curve plots a homogeneous gas BdG calculation, and the solid black curve plots the result of Eq. (1).The simulations used our calibrated experimental parameters.

FIG. 3 .
FIG. 3. Moment of inertia from scissors mode.(a-inset) Measured dipole mode frequencies (markers) along with fits (curves) where the frequency at U0 is the only free parameter for each curve.(a) Scissors mode frequency.Blue and green correspond to U = 0 trap frequencies (34, 51) Hz and (54, 36) Hz respectively.(b) Moment of inertia in units of Ic.Symbols are the data computed as described in the text, and the solid curves are GPE predictions using I = ∂Ω Lz , with angular frequency Ω.

1 FIG. 4 .
FIG. 4.Moment of inertia in rotating systems computed using 2D GPE simulations.The left column indicates simulations in which the lattice is static while in the right column the lattice co-rotates with the confining potential.(a,b) Angular momentum density for trap frequencies 2π×(56,36) and U0 = 10Er.(c, d) Total momentum of inertia in traps with frequencies 2π×(56,36) (top, green) and 2π × (36, 56) Hz (bottom, blue).Dashed curves plot I sf /Ic and the solid curve plots I/Ic.