Superfluid fraction tensor of a two-dimensional supersolid

We investigate the superfluid fraction of crystalline stationary states within the framework of mean-field Gross–Pitaevskii theory. Our primary focus is on a two-dimensional Bose–Einstein condensate with a non-local soft-core interaction, where the superfluid fraction is described by a rank-2 tensor. We then calculate the superfluid fraction tensor for crystalline states exhibiting triangular, square, and stripe geometries across a broad range of interaction parameters. Factors leading to an anisotropic superfluid fraction tensor are also considered. We also refine the Leggett bounds for the superfluid fraction of the 2D system. We systematically compare these bounds to our full numerical results, and other results in the literature. This work is of direct relevance to other supersolid systems of current interest, such as supersolids produced using dipolar Bose–Einstein condensates.


I. INTRODUCTION
Superfluidity emerges as a distinctive state of matter within quantum many-body systems, manifesting under specific conditions and exhibits various traits, most notably the absence of viscosity.An enduring question has been whether superfluidity could coexist with crystalline order realizing a supersolid state of matter [1][2][3][4].In recent years several different types of experiment utilizing Bose-Einstein condensates (BECs) have successfully observed supersolidity [5][6][7][8][9].
Leggett's seminal work highlighted that the introduction of crystalline order, with the concomitant breakdown of translational invariance, leads to a reduction in the system's superfluid fraction even at zero temperature [3,10].Although the emergence of crystalline order has been directly witnessed in supersolids produced using BECs, the characterization of superfluid transport remains less straightforward.For dipolar BEC supersolids, excitation spectroscopy has been employed to infer superfluidity [11][12][13] (also see [14,15]).However, this does not directly quantify the superfluid fraction of the system.We also note that scissor mode excitations have furnished evidence of the diminished superfluid fraction [16].
Recently two groups have quantified the superfluid fraction in an artificial supersolid produced by applying an optical lattice to a BEC [17,18] (i.e. with imposed, rather than spontaneously formed density order).Their superfluid measurements were based on two methods: (i) a superfluid bound developed by Leggett, which requires the precise measurement of the density profile [18]; and (2) a ratio of the speeds of sound propagating parallel and perpendicular to the lattice used to obtain the effective mass (and thereby the superfluid fraction) along the lattice [17,18].It is currently not clear how to extend these ideas to two-dimensional (2D) or threedimensional (3D) supersolids, where the superfluid fraction is properly a rank-2 tensor.For example, the sound measurements do not generally apply to supersolids where spontaneous translational symmetry breaking leads to multiple gapless sound branches [19][20][21][22].The first experimental observation of 2D dipolar supersolids [23] underscores the need for a deeper understanding of superfluidity in higher dimensions.
In this work we consider theoretical methods to quantify the superfluid fraction for a 2D supersolid.We use meanfield model of a 2D Bose gas with soft-core interactions to illustrate our study.This model exhibits a transition to a triangular crystal when the dimensionless interaction parameter, Λ, surpasses a threshold value [24].Furthermore, our investigation encompasses metastable crystal configurations characterized by square or stripe geometries, thereby furnishing a diverse platform to implement and validate superfluid measures.The three main components of this work are: (1) We review the definition of the superfluid tensor via the nonclassical translational inertia.This describes the steady state response of the system to moving walls, which carries the normal fluid but leaves the superfluid at rest.For the crystalline ground state the superfluidity can be quantified by understanding the phase profile in a unit cell in response to the motion.The equations determining this phase profile were presented in Ref. [25], but are equivalent to a formulation by Saslow in 1976 [26].The superfluid tensor for the 2D soft-core system was previously computed using this approach in Ref. [27], however our results disagree with these, underscoring the complexity of the calculation.(2) We examine the effective mass of the BEC in a crystalline state, which also relates to the superfluid fraction [28].This approach has been used to determine the superfluid fraction tensor in a 2D case by solving the nonlinear Gross-Pitaevskii equation (GPE) for the condensate in motion [29].We discuss the connection to the phase profile formulation and show that the effective mass can be determined by a simpler linearized problem, which is equivalent to solving for the Hartree excitations.(3) The Leggett upper bound was used by Zhang et al. [30] to quantify the superfluid fraction of a 2D dipolar supersolid, although this involved minimising over angles, and no comparison was made to a direct calculation of the superfluid fraction to quantify the reliability of the bound.
Here we comprehensively study the use of the Leggett upper and lower bounds in application to a 2D supersolid.From this formulate a useful form of the bound involving a single unit cell without the need to consider various angles.Also, our general techniques allow us to explore conditions under which 2D crystalline states exhibit anisotropic superfluidity (cf.[29]), hence revealing situations where the tensorial nature of the superfluid fraction manifests.

II. 2D SOFT-CORE SUPERSOLID MODEL
Our model is a meanfield description of a 2D BEC, without any confining trap, and interacting via a soft-core interaction potential U (x) = U 0 θ(a sc − |x|), where a sc is the core radius, U 0 is the potential strength, θ is the Heaviside step function, and x = (x, y).The meanfield energy functional for this system is where the wavefunction Ψ is normalized to the total atom number N .Here we will be interested in the ground state for N atoms in a large plane of area A, in the limit that both of these are large, and with a mean areal density of n = N/A.We then choose to set The interaction term is efficiently evaluated by the convolution theorem, where the k-space form of the interaction potential is Ũ (k) = 2πU 0 a sc J 1 (k ρ a sc )/k.Adopting a sc as the unit of length and ℏω 0 = ℏ 2 /ma 2 sc as the unit of energy, it is convenient to define the dimensionless interaction parameter We look for spatially periodic (crystalline) stationary state solutions.We consider 2D periodic solutions on a unit cell with lattice vectors a 1 = ax, and a 2 = a(cos θx + sin θŷ), where we will consider θ = π/3 (triangular lattice) or θ = π/2 (square lattice).We also consider a 1D stripe phase with a single lattice vector a 1 = ax, with Ψ(x) = √ nψ(x).In both the 1D and 2D cases, the stationary states are characterized by a single cell size parameter a.
Numerically the stationary states can be calculated with spectral accuracy using a planewave presentation in reciprocal space associated with the unit cell we choose to consider (e.g.see [31]).We use the gradient flow method [32,33] to minimise the energy and obtain the ground state solution for the unit cell.We note that these solutions satisfy the timeindependent GPE µψ = L GP ψ, where is the GP operator and µ is the chemical potential.We denote the ground state energy per particle for a cell of size a as where uc denotes the unit cell, A uc = |a 1 × a 2 | is the unit cell area and For the uniform state the energy per particle is independent of the unit cell and is given by E u = 1 2 Λℏω 0 .For the modulated states there is an optimal value of cell size E(n) ≡ min a E(n; a), that minimises the energy per particle.
The phase diagram of this model is well-known (see [15,19,24,27,31,34]) and is summarised in Fig. 1.As Λ changes the ground state undergoes a first order phase transition from the uniform state to a triangular state at Λ m = 39.49[see Fig. 1(a)].The stripe and square states are metastable states, and undergo a second order transition to the uniform state at Λ rot = 46.298,where a roton excitation softens in the uniform state.The optimal unit cell size a for the modulated states is shown in Fig. 1(b).At the transition point the cell size of the square and stripe states matches the roton wave-length.In Fig. 1(c) we show the density contrast with |ψ| 2 max and |ψ| 2 min being the maximum and minimum values of |ψ| 2 in the unit cell, respectively.The density contrast is 0 in the uniform state, and thus serves as an order parameter for crystalline order.Having introduced our model we now turn to consider how to evaluate the superfluid fraction.We focus on the 2D case, relevant to our model, although our results generalize to 1D and 3D cases.

A. Supersolid in a slowly moving frame
To understand the superfluid response we consider our system with fictitious walls moving at a velocity of v with respect to the lab frame (see Fig. 2).The moving walls define the preferred frame for the normal component.We compute the ground state in the co-moving frame by solving the timeindependent GPE For small v we can make a perturbative expansion for the spatially periodic solution where ψ 0 is the stationary solution for v = 0 and ϕ 1 is the first order correction to the phase.We emphasize that changes in the density profile occur beyond first order and are neglected in our treatment.In the moving frame the current is where ρ 0 ≡ n|ψ 0 | 2 .The phase term is determined from the condition ∇ • j = 0, which ensures the density is stationary in the moving frame (i.e. the crystal is co-moving).In 1D this condition means that the current is constant.In higher spatial dimensions non-trivial flow patterns can emerge, determined by solving Saslow [26] originally obtained this result, albeit in reciprocal space, and used it to estimate the superfluid fraction of solid 4 He (also see [35]).From Eq. ( 8) the momentum in the lab frame is which can potentially be in a different direction to v (see Fig. 2).The energy of this state is where x nΦ(x)ρ 0 (x) denotes the interaction energy and E 0 is the energy 1 of ψ 0 .
If the integration region is taken over a whole number of unit cells, then using integration by parts and Eq. ( 10), we obtain Under these conditions Eq. ( 13) can be written as 1. Linear solution The linear system (10) can be recast by introducing the auxiliary vector function K = (K x , K y ) as The component functions satisfy [from Eq. ( 10)] 1 The kinetic energy of ψv is The first term is the energy of ψ 0 , the next term is the additional energy of ψv.
where we use the index notation i = {x, y}.In terms of the auxiliary vector function, the lab frame momentum (11) is For a given density ρ 0 , we can solve Eq. ( 16) for K by inverting this linear system2 .Because ρ 0 and K are periodic we can solve Eq. ( 16) using a single unit cell.Equation ( 16) was obtained by Josserand et al. [25,27,36,37] by applying the technique of homogenization to the GPE in the context of studying the superfluidity in soft-core models 3 .

Energy response to moving walls
The superfluidity can be identified from the lab frame energy.Writing this energy the form [3] where M = N m is the classical inertia of the system and N is the atom number.The term ∆E(v) describes the departure of the energy from the classical result ( 1 2 M v 2 ), associated with the superfluid component.From this we can define the superfluid fraction tensor f ij as Notably, when f ij → δ ij the system is a pure superfluid and E(v) is independent of v for sufficiently small velocities.Comparing Eqs. ( 14) and ( 18), we see that and using (17) we obtain

Momentum response to moving walls
We can also consider the momentum of the system in the lab frame, associating this momentum with the normal component that moves with the walls.In the low velocity limit this leads to the identification such that for f ij = 0 we have P = M v. Thus the superfluid fraction, using Eq. ( 17), is which is identical to result (21).

B. Superfluid phase twist and effective mass
We now discuss an alternative approach to identifying the superfluid tensor.To do this we consider the system in the frame where the crystal is stationary, and apply a phase twist to the wavefunction.For our 2D system we can specify the direction of phase twist as ∆θ = ∆θ x x + ∆θ y ŷ, over some macroscopic sample length of L. We assume this twist manifests as an average phase gradient that can be characterized by the quasimomentum q = ∆θ/L.For small q the energy, which we denote at E ′ (q), is (to second order in q) which defines the effective mass tensor m * ij .The relationship between the superfluid and effective mass tensors is as discussed in [18,28,29,38].Here, motivated by this definition, we explore an alternative method to compute the superfluid fraction and then discuss the equivalence of (25) to our earlier formulation.

Calculation of effective mass of a supersolid
The ground state with quasimomentum q is a Bloch wave solution of the GPE Similar to the analysis of Sec.III A 2, we find that to leading order in q the density of the ground state is unchanged, so in this limit we can set This solution obeys Bloch's theorem where ψ 0 (x)e iϕ1(x) is periodic on the unit cell.Thus we see that ϕ 1 (x) describes the response to the quasimomentum phase gradient arising from the crystalline density.In 1D this can lead to the overall phase profile having a staircase pattern (see Ref. [17]).Also, the unchanged density profile means that we can set Thus ψ q can be determined by solving the linearized GP equation which coincides with the single particle Hartree excitations for the lowest band.Comparing the full energy functional (1) to the GPE operator (3) we see the energy these Bloch states can be evaluated using the result with E int,0 being the interaction energy calculated with ψ 0 (noting that E int,0 is independent of q).Since methods for solving the Schrödinger equation are well developed, it is relatively straightforward to numerically evaluate µ q and hence determine the effective mass by computing finite difference second derivatives of N µ q with small changes in q.

Relation to moving frame treatment
It is worth considering the relationship to the analysis presented in Sec.III A. In defining the effective mass we are working in the (co-moving) frame where the crystal is stationary, and the superfluid moves with an average superfluid velocity of v s = ℏq/m.Equivalently, the crystal is moving at a velocity of v = −ℏq/m relative to the superfluid (lab) frame.First, this observation reveals that ψ q , as given in (27), in just result (8) boosted by v to the crystal frame.This also confirms that the function ϕ 1 appearing in both treatments is identical.
The previous result we obtained for the lab frame energy E(v) [Eq.(14)] and E ′ (q), are thus related by the Galilean transformation Here primed coordinates are quantities in the crystal (comoving) frame and unprimed coordinates are in the lab frame.
Using results (18) and (20) we obtain where ∆E defined in (20), and we have changed variable v → q.From this it is easily verified that the effective mass definition (25) of the superfluid fraction tensor is identical to Eq. ( 21).

IV. SUPERFLUID FRACTION RESULTS
In this section we evaluate the superfluid fraction for the various modulated states of the 2D soft-core gas.We perform these calculations by solving for the auxiliary vector function and then using (23) to determine f ij .We find identical results using the effective mass approach.The numerical calculations are preformed in a single unit cell using the spectral approach discussed in Sec.II.

A. Stripe state
The stripe state realizes a 1D supersolid, in that crystalline order is along one spatial dimension, here taken to be x.For this case the solution of Eq. ( 16) is analytic.Notably, K y is zero, K x reduces to a function of x, and Eq. ( 16) simplifies to This has the solution , where the integration constant is c = [1/a a 2 − a 2 dx ′ /ρ 0 (x ′ )] −1 to ensure that K(x) is periodic.Using Eq. ( 23) we obtain the tensor components to be with (33) being the Leggett upper bound for the superfluid fraction [3,10], which is known to be exact for 1D cases [37].
Examples of the stripe state density profile and the K x functions are shown in Fig. 3 for a low contrast state close to the phase transition [subplots (a) and (b)] and a high contrast state deep in the crystalline phase [subplots (c) and (d)].We also show the momentum for the unit cell in terms of the moving frame velocity component v x .In the high contrast state we observe that K x ≈ x in the center of the unit cell where the wavefunction has high density.This leads to the state having a high momentum (i.e.majority of the system moving with the walls) and hence a low superfluid fraction along x.  4 .In all triangular and square lattice cases we find the superfluid fraction tensor to be diagonal and of the isotropic (scalar) form This observation of isotropy was also found for the triangular state of the 2D soft-core model in Ref. [27].For the case with Λ = 100 the auxiliary function is noticeably sharper [similar to the stripe result in Fig. 3(d)] and we observe that K x ≈ x and K y ≈ y in the central part of the unit cell where the density is high.The superfluid fraction tensor for this case is f s = 0.02.

C. Superfluid fraction results
In Fig. 5(a) we show the superfluid faction f x x for a range of Λ values and for the three different lattice geometries we consider.The triangular state undergoes a first order transition from the uniform state at Λ = Λ m , and here the superfluid fraction suddenly drops from unity.In contrast, the stripe and square states undergo a continuous transition at Λ = Λ rot and the superfluidity changes continuously.The stripe state has a f xx superfluid fraction than the square state, consistent with the stripe state having a lower density contrast relative to the square state [see Fig. 1(c)], and hence more tunnelling between unit cells.
For comparison we show the superfluid fraction results (also for f xx ) reported in Fig. 3 of Sepúlveda et al. [27].These results are also for the triangular state and were computed using the auxiliary vector function approach we employ, but are seen to compare poorly to our results.In that work the ground states were calculated on a large square region containing many unit cells, but not commensurate with the unit cells.This could strain the resulting crystal and effect the superfluidity (see discussion in Sec.IV D).Also, as our results show, the auxiliary vector functions exhibit quite sharp features for large Λ [e.g.see Figs. 4(b2) and (b3)], which make them difficult to calculate accurately.In Ref. [27] a finite element methods was used to calculate the auxiliary vector functions, whereas all our calculations are performed using a spectral method (i.e.planewave expansion in reciprocal lattice vectors) that is specific to the unit cell.We check that our results are converged with respect to the number of reciprocal lattice vectors (typically we use > 10 3 planewaves in the expansion).Also, we have validated details of our calculations (e.g. using predictions of the melting point, excitation spectra, optimal lattice parameters) against more recent calculations on the 2D soft-core model [31,34].
In Fig. 5(a) we also show the lower bound for the superfluid fraction by Aftalion et al. [39].Their work was developed for soft-core models, with the 2D result we show here being Our results in Fig. 5(a) confirm that this is a lower bound, but is more than two orders of magnitude smaller than our calculated superfluid fraction for the interaction parameter regime we consider.

D. Anisotropic 2D Superfluidity
Our earlier results showed that for the 2D crystalline square and triangular states states have an isotropic superfluid response (36), whereas the 1D stripe phase is anisotropic.We investigate here the conditions where 2D crystals become anisotropic.To do this we can extend our consideration to unit cells of the form a 1 = a 1 x, and a 2 = a 2 (cos θx + sin θŷ), i.e. allowing the direct lattice vectors to be of different length and for θ to vary continuously.To analyse these cases we consider the superfluid fraction tensor, and here we denote the matrix representing this as f .This is a real symmetric matrix and can be taken to a diagonal form by the similarity transfor- is the rotation matrix with angle φ.From this we identify the eigenvalues ( f1 , f2 ) of f as the superfluid fractions along the principal axes, which are orthogonal axes rotated at an angle φ with respect to x and y.In Fig. 6(a) we consider a crystalline ground state with the lengths of the unit cell fixed (a 1 = a 2 ), but allow the angle θ to vary.The results show that the isotropic condition ( f1 = f2 ) occurs for triangular (θ = π/3, 2π/3) and square (θ = π/2) lattices, but is anisotropic at other angles.For these results the rotation angle defining the principal axes is φ = θ/2.
In Fig. 6(b) we fix the angle to the rectangular case (θ = π/2) and vary the length of a 2 relative to a 1 .In general a difference in length causes the superfluid tensor to become anisotropic, although in this case the principal axes remain aligned to x and y (i.e.φ = 0).In addition to the square case (a 2 = a 1 ), we also observe that an isotropic response occurs at a 2 − a 1 ≈ −0.144a sc .

V. LEGGETT BOUNDS FOR AN ISOTROPIC 2D CRYSTAL
Leggett has developed bounds for the superfluid fraction based on integrals of the ground state density [3,10] (also see [38]).Here we specialize these bounds to the 2D case and restrict our attention to cases of isotropic superfluidity (i.e.triangular and square crystals).We first investigate them in application to a macroscopic region (with many unit cells), in the spirit of Leggett's original formulation.We use this to investigate the sensitivity to size and orientation of the region.These results motivate us to specialize the bounds to a form that can be evaluated, without any adjustable parameters, using the system density on a single unit cell.

A. Application of a macroscopic region
Here we consider evaluating the Leggett bounds using the system density profile in a square macroscopic region M of size L × L, with L ≫ a (i.e enclosing many unit cells).Such a region is indicated in Fig. 7(a) 5 .The expressions are the upper and lower bounds, respectively.We vary the orientation of the region M by rotating the square region by an angle of θ with respect to the standard axes [see Fig. 7(a)].
As the angle changes relative to the crystal axes (defined by a 1 and a 2 ) the value of f ± M changes.We notice that there are minima (maxima) in f + M (f − M ) when the axis of the first integration (i.e.y ′ -axis) is parallel to the direct lattice vectors [e.g. it is parallel or anti-parallel to a 1 when θ = ± π 2 , and similarly for a 2 when θ = ± π 6 , etc., see Fig. 7(b)] 6 .Thus, we can find the strictest bounds by selecting θ to align to these directions.In Figs.7(c) and (d) we consider how f ± M , evaluated at θ = π/2, changes with increasing size of the integration region.The results oscillate, but converge as L increases.These oscillations occur because the integration region is not commensurate with the unit cells, and suggest that we would be able to extract the value from a careful analysis of a single unit cell.

B. Application to a unit cell
For the 2D cases we considered in this paper, the primitive unit cell defined by a 1 and a 2 is a rhombus with side length a, angle θ and vertical height a y = a sin θ.We specialize the Leggett bound definitions to the unit cell as follows Here the integral aj refers to integration parallel to the direct lattice vector a j , i.e. along the horizontal dotted blue lines or diagonal magenta lines in Fig. 8 for a 1 and a 2 , respectively.This turns to be be convenient numerically, as the Fourier based sampling of the wavefunction is on a grid of the type represented in that figure.Hence these integrals can be numerically calculated with spectral accuracy using modest resources.As verified in Figs.7(c) and (d), the macroscopic application of the Leggett bounds at the optimal angles converge to f ± as L → ∞.
To make the integrals appearing in Eqs. ( 41) and ( 42) unambiguous we note that the unit cell can be mapped to a rectangular region of dimensions a × a y with the affine transformation (a shear in the x-direction) of Applying the affine transform we define the coordinates u = A t x, with u = (u x , u y ).The primitive unit cell becomes the rectangular region − 1 2 a ≤ u x ≤ 1 2 a and − 1 2 a y ≤ u y ≤ 1 2 a y .Noting that the Jacobian determinant of this transformation is unity, the integrals can be written as, for example

C. Comparison to exact calculations
In Fig. 5(b) we compare f ± to our earlier results for f xx superfluidity for the triangular case.For larger Λ values f + can be an order of magnitude larger than f − , so the bounds are not that tight, but still provide a useful estimate.In contrast for the square lattice states [inset to Fig. 5(b)] the bounds are much tighter.Here our results reveal that the bounds estimate the superfluid fraction with an absolute error |∆f | = |f ± − f xx | ≲ 10 −3 , similar to results for a 1D dipolar supersolid in a tube found in Ref. [40].The bounds are more effective in this case because the wavefunction is close to being separable.For the stripe phase the well-known 1D version of the bound is exact (also see Sec.IV A).

VI. CONCLUSIONS
We have discussed several computational approaches for determining the superfluid fraction tensor of the meanfield system.One approach involves determining the phase response to moving walls by solving two Poisson-like equations for the vector auxiliary equations (16).From these solutions the superfluid density can be evaluated using Eq. ( 21).The effective mass approach instead involves solving the linear Schrödinger problem (28), but requires a finite difference in q to calculate the superfluid tensor.We have discussed why these approaches are formally equivalent.Our results show that the ground state triangular supersolid and the metastable square supersolid both exhibit isotropic superfluidity.We have investigated how changes in the lattice vector orientation and length can cause the superfluity to become anisotropic.We also considered how the Leggett bounds can be applied to 2D supersolids with isotropic superfluidity, finding that the tightest bounds from the thermodynamic-limit of their application can be deduced from a calculation on a single unit cell.
Despite the theoretical definitions of superfluidity being well-known for many decades, there are few explicit results for higher dimensional supersolids.Indeed our results disagree with previously published results for the 2D soft-core model, which forms one of the most basic models for describing higher dimensional supersolidity.This lack of accurate numerical results has likely hindered a better understanding of the superfluid response in this system, and how to apply the Leggett bounds in 2D and 3D crystals.In this absence of understanding many proxy-measures of superfluidity have been used in the recent literature without justification (e.g.various forms of density contrast etc.) While our work here has focused on the thermodynamic limit, our results will have implications for finite trapped systems of the type realized in experiments.The auxiliary vector function calculation only requires a density profile within a unit cell, which could allow this approach to estimate the local superfluid fraction (associated with that cell) in a finite system.Similarly the Leggett bounds could also be applied to this case.Since recent experiments have been able to apply the Leggett upper bound to BECs in 1D optical lattices [17,18], we expect this avenue of investigation may be practical for experiments considering superfluidity of higher dimensional supersolids (or to BECs inhigher dimensional optical lattices).It would be interesting to compare the distribution of local superfluid behavior to global measurements of the trapped case, e.g.via the calculation of the rotational inertia (e.g.see [16,[41][42][43][44][45][46][47]).A future direction would be to consider the extension of these ideas to consider superfluidity in the stripe phase of spin-orbit coupled BECs (c.f.[48][49][50][51]), in which Galilean invariance is also broken.It might also be necessary to explore scenarios in which the supersolid state may not necessarily be the ground state, but rather exhibits specific phase excitations [52] (also see [53]) or lattice distortions.In the latter case, our results already suggest this will lead to anisotropic effects in the superfluid response.

FIG. 1 .
FIG. 1. Crystallisation phase transition of a 2D soft-core BEC.(a) The energy per particle relative to the uniform ground state for the triangular (black), square (blue), stripe (red) and uniform (magenta) states.(b) The energy minimising unit cell size of the modulated states.(c) The density contrast.Vertical dotted lines indicate the melting point for the triangular state and the critical point for the square and stripe states, which coincides with the roton softening of the uniform state.In (b) the roton wavelength is indicated for reference.

FIG. 2 .
FIG.2.Schematic of moving frame used to define superfluidity.The walls are taken to move at a velocity of v relative to the lab frame, and in steady-state the crystalline structure will move with the walls.The superfluid remains stationary in the lab frame.The momentum P of the moving system can be in a different direction to v due to the tensorial character of the superfluidity.

5 FIG. 3 .
FIG. 3. Ground state profiles (a),(c) and Kx auxiliary functions (b), (d) for the stripe state.Dashed line in (d) shows that Kx ≈ x near the center of the cell where the crystal density is high.

3 FIG. 5 .
FIG. 5. (a) Superfluid fraction fxx computed using (21) for the various geometry stationary states (solid lines) as a function of Λ.The data from Sepùlveda et al. [27] for the triangular state presented for comparison (black line squares) and the relevant analytic bound of Aftalion et al. [39].(b) The Leggett bounds f ± shown relative to the triangular state results from (a).The Leggett bounds are tighter for the square state and are shown in the inset.

FIG. 7 .
FIG. 7. Leggett bounds.(a) Our macroscopic sample M is specified by square grid of dimensions L × L indicated by the mesh.We have introduced the coordinates (x ′ , y ′ ), which are rotated by θ relative to the default axes.The crystalline density over this region is integrated along x ′ and y ′ to evaluate the Leggett bounds f ± M .(b) The Leggett upper (black line) and lower (blue line) bounds for L = 9.6asc as a function of θ.The local minima of f + M and local maxima of f − M occur when the rotated axes are parallel to the direct lattice vectors.The (c) f + M minima (least upper bound) and (d) f − M maxima (largest lower bound) at θ = π/2 as a function of L, relative to the unit cell values [from Sec.V B] f + = 0.4283 and f − = 0.2794, respectively.The triangular ground state for Λ = 48 is used for the results presented in this figure.