Supersolidity of cnoidal waves in an ultracold Bose gas

A one-dimensional Bose-Einstein condensate may experience nonlinear periodic modulations known as ``cnoidal waves''. We argue that such structures represent promising candidates for the study of supersolidity-related phenomena in a non-equilibrium state. A mean-field treatment makes it possible to rederive Leggett's formula for the superfluid fraction of the system and to estimate it analytically. We determine the excitation spectrum, for which we obtain analytical results in the two opposite limiting cases of (i) a linearly modulated background and (ii) a train of dark solitons. The presence of two Goldstone (gapless) modes -- associated with the spontaneous breaking of $\mathrm{U}(1)$ symmetry and of continuous translational invariance -- at large wavelength is verified. We also calculate the static structure factor and the compressibility of cnoidal waves, which show a divergent behavior at the edges of each Brillouin zone.


I. INTRODUCTION
Supersolid phases of matter have attracted an increasing interest in the last few years. In these configurations two apparently conflicting properties, namely superfluidity and crystalline order, can coexist giving rise to novel features (see for instance the review [1]). The existence of such a phenomenon had initially been investigated, and apparently ruled out, by Penrose and Onsager in the 1950s [2]. It was reproposed shortly after by Gross. In Refs. [3,4] he considered a system of interacting bosons in the semiclassical limit, where the bosonic quantum field can be replaced by a classical one, which obeys a nonlinear field equation. The latter admits periodic solutions, describing a uniform background with a crystal lattice on top of it. In the subsequent decades the search for possible superfluid solid phases was extended and other scenarios in which supersolidity could occur were examined [5][6][7][8][9][10][11]. The main candidate has been for many years the solid phase of helium. However, the most recent experimental results and theory analyses seem to preclude superfluidity in bulk solid helium [12], and the attention turned to solid-helium two-dimensional films [13].
On the other hand, significant progress has been made with ultracold atomic gases starting from 2017, with the first observations of an incompressible supersolid state in bosonic systems coupled to two optical cavities [14] and of a superfluid smectic state in spin-orbit-coupled Bose-Einstein condensates [15]. Even more recently, coherent droplet arrays have been realized in dipolar quantum gases [16][17][18]. This has stimulated a large amount of further experimental work, shedding light on the spectrum and the collective modes [19][20][21][22], the superfluidity properties [23], and the out-of-equilibrium dynamics [24] of this exotic phase of matter.
A rather intriguing scenario for the occurrence of supersolidity is the one pointed out by Pitaevskii in 1984 [10]. He proved that a sample of superfluid 4 He flowing along a capillary with a velocity exceeding Landau's critical value develops a layered structure. This structure results from the condensation of excitations close to the roton momentum [25] and is at rest with respect to the walls of the capillary. Its excitation spectrum is deformed such that the system remains superfluid. These findings were later confirmed by numerical simulations based on a density functional approach [26]. The same physics can be observed in ultracold Bose gases as well, as found in recent times by Baym and Pethick [27]. In this reference the authors assumed a finite-range interaction between particles. This shifts the critical momentum at which the Landau instability can occur to a finite value. Similar to 4 He, a large number of excitations with momentum close to this value (called "levons") are created when crossing the Landau velocity, which represents the onset of the transition to the layered phase. The latter features a superfluid fraction smaller than one.
In general, supersolid-like configurations can have smaller energy than uniform ones only in special circumstances. Typically one needs to have either specific kinds of interparticle interaction (such as dipole-dipole, finite-range, or cavity-mediated), or a properly modified single-particle spectrum (as in the case of spin-orbitcoupled Bose-Einstein condensates). However, when none of these conditions is fulfilled one may still have a supersolid behavior in some excited state. This is the case of a standard quasi-one-dimensional dilute Bose gas with repulsive contact interaction, which is described by the Gross-Pitaevskii equation. This equation is known to have spatially periodic stationary solutions, which were studied by Tsuzuki in 1971 [28]. Korteweg and de Vries [29] had coined the term "cnoidal waves" for solutions of this type, because they can be expressed in terms of the Jacobi cosine amplitude function, denoted by cn [30]. In systems of bosons rotating in a ring trap, transitions between metastable uniform and cnoidal configurations have been predicted [31,32]. Very recently, cnoidal-wave-like solutions have been found for the extended Gross-Pitaevskii equation describing a self-trapped cigar-shaped Bose gas [33].
Cnoidal waves can be regarded as the equivalent of Pitaevskii's layered phase for an ultracold Bose gas. At variance with the case considered in Ref. [27], for repulsive contact interaction the Landau instability takes place at vanishing momentum, and there is no mechanism similar to levon condensation. Nevertheless, one can achieve a cnoidal structure by moving an obstacle at a suitable velocity through the condensate. The scope of this work is to highlight that these configurations exhibit typical features of supersolids in both their static and dynamic behavior. As such, they are new candidates for studying phenomena related to supersolidity within the most standard Bose-Einstein condensates. The latter do not suffer from the strong three-body losses typical of dilute ultracold systems in which the stabilization is due to beyond-mean-field equation of states, as for dipolar gases and quantum mixtures.
This article is organized as follows. In Sec. II we introduce the model to investigate and the equations governing it. In Section III we present the derivation of the cnoidalwave solution and illustrate some of its most important properties. The dynamic behavior of a cnoidal wave is discussed in Sec. IV, where we discuss the excitation spectrum, the static structure factor and the compressibility. We summarize in Sec. V. Finally, some technical details are presented in the Appendices: Appendix A presents the properties of the cnoidal wave in some limiting cases; the procedure for solving the Bogoliubov equations is explained in Appendix B; and Appendix C computes the lower branch of the spectrum of a train of dark solitons.

II. THE MODEL
Let us consider a quasi-one-dimensional weakly interacting Bose-Einstein condensate at zero temperature. The condensate wave function ψ(x, t) obeys the Gross-Pitaevskii equation Here m is the mass of a particle, g > 0 the interaction strength, and µ the chemical potential. We use in Eq. (1) and throughout the paper the convention that subscripts denote derivatives with respect to x and t. We now perform the Madelung decomposition of the wave function, ψ = A e iΘ . Inserting it into the Gross-Pitaevskii equation (1) yields two coupled equations for the real quantities A ≥ 0 and Θ, equivalent to the collisionless hydrodynamic equations for a superfluid. The first one, expressing the particle number conservation, is the continuity equation, which reads with n = A 2 and we recall that v = Θ x /m is the velocity field. The second Madelung equation reads After taking the gradient on both sides, it becomes formally identical to the Euler equation for the potential flow of an inviscid fluid.

III. CNOIDAL-WAVE SOLUTION
The cnoidal-wave solution exhibited by the Gross-Pitaevskii equation (1) has been extensively studied in the literature. In this section we review its derivation in order to fix the notation, set the background for the subsequent calculations, and present some of its most relevant features.

A. Derivation of the cnoidal-wave solution
In order to find stationary solutions of the Gross-Pitaevskii equation (1) one has to set n t = 0 and Θ t = 0. This turns Eqs. (2) and (3) into ordinary differential equations in x. Following Ref. [28] we shall integrate these equations imposing the condition that the condensate density and velocity oscillate in space at a given wavelength Λ around fixed average valuesn andv. Integrating once Eq. (2) with respect to x one obtains where J denotes the constant value of the current density. We can use this result to eliminate Θ x from Eq. (3). This yields where Equation (5) has the same mathematical structure as the energy conservation of a classical particle having "position" A at "time" x [34,35]. The integration constant E plays the role of the energy and W that of the external potential.
In the following we assume that the current J verifies the inequality [36,37] which ensures that W (n > 0) has a local minimum W min and a local maximum W max , as illustrated in Fig. 1. The maximal value (7) of J is analogous to the Ginzburg-Landau critical current in a superconductor [38]. For given values of the parameters, the range of values that the solutions of Eq. (5) can take is fixed by the condition W (n) ≤ E. In particular, the roots of W (n) = E identify the extrema of n, where A x = 0. They are the analogous of the turning points of a classical system. In the case W min < E < W max considered in Fig. 1 there are 3 such roots which we denote as n 1 , n 2 and n 3 with 0 ≤ n 1 ≤ n 2 ≤ n 3 . Let us then rewrite Eq. (5) as Comparing the two rows of the above equation one immediately finds out how to express µ, E, and J in terms of n 1 , n 2 , and n 3 . The result reads µ = g(n 1 + n 2 + n 3 ) 2 , (9a) E = g(n 1 n 2 + n 2 n 3 + n 3 n 1 ) 2 , (9b) A bounded solution of Eq. (8) oscillates between n 1 and n 2 and thus is of the form Inserting this Ansatz into Eq. (8) yields (upon properly defining the spatial origin) where am is Jacobi's amplitude function [30] and The corresponding density and phase read n(x) = n 1 + (n 2 − n 1 ) sn 2 mg(n 3 − n 1 ) x m e , Θ(x) = ± n 2 n 3 n 1 (n 3 − n 1 ) Π (−n e ; ϕ(x)|m e ) .
Here sn(u|m e ) = sin(am(u|m e )) is the Jacobi sine amplitude function, and Π (−n e ; ϕ|m e ) denotes the incomplete elliptic integral of the third kind [30]. The quantity − n e = −(n 2 − n 1 )/n 1 is called the "characteristic". The condensate phase (14) was determined integrating Eq. (4) with respect to x, imposing Θ(0) = 0 for simplicity; the plus (minus) sign corresponds to a positive (negative) value of the current J . Equations (13) and (14) express the cnoidal-wave solution of the Gross-Pitaevskii equation (1). It was first investigated by Tsuzuki in Ref. [28], see also Ref. [39]. This solution depends on the three parameters n 1 ≤ n 2 ≤ n 3 . It represents a stationary layered structure, i.e., such that its density profile exhibits periodic spatial modulations; a fixed current J flows through the fringes. The oscillation wavelength and average density are computed in Sec. III B below and are given by Eqs. (15) and (16), respectively. The modulations correspond to a spontaneous breaking of continuous translational invariance. Because of the simultaneous presence of superfluid and crystal order, cnoidal waves are expected to exhibit a supersolid behavior in both their static and dynamic properties. These aspects will be elucidated in the following sections.

B. Properties of the cnoidal-wave solution
We shall now examine some characteristic features of cnoidal waves. These include the average density, the contrast of the density modulations, the superfluid fraction, and the energy per particle.

Density profile and contrast of the fringes
The density profile (13) oscillates with a wavelength where K(m e ) is the complete elliptic integral of the first kind [30]. These oscillations occur around an average value given by [28] n = 1 Λ where E(m e ) is the complete elliptic integral of the second kind [30]. Using this average density we can define the healing length ξ = / √ mgn and the sound velocity c = gn/m. It is useful to rewrite n 1 , n 2 , and n 3 in terms ofn and of the two dimensionless parameters m e and η = n 3 − n 1 n .
From Eqs. (12) and (16) one gets where Γ(m e ) = E(m e )/K(m e ). One can easily check that the conditions 0 ≤ m e ≤ 1 and η ≥ 0 are sufficient to ensure that n 2 and n 3 are non-negative. Additional constraints come from the requirement n 1 ≥ 0. The latter is satisfied for any 0 ≤ m e ≤ 1 if 0 ≤ η ≤ 1; but, if η > 1, m e should not be larger than a threshold value m max e defined by Γ(m max e ) = (η − 1)/η. In the following we shall see that, when considered as functions of m e , the various observables have different behaviors, depending on whether η is smaller or larger than 1.
Making use of the average density Eq. (16) we can decompose the density (13) into a uniform and a modulated component as n(x) =n + ∆n(x), with In Fig. 2 we report a few density profiles of cnoidal waves for different values of m e and η. At small m e the oscillations have small amplitude and are practically sinusoidal, as discussed in Appendix A 2. Increasing m e at fixed η produces fringes with larger amplitude and wavelength, as well as significant deviations from the sinusoidal behavior. When m e is close to 1 the density profile takes the characteristic shape of a "soliton train", made by quasi-uniform regions separated by thin deep valleys.
A useful quantity to characterize the fringes is their contrast, At small m e the contrast behaves like C m e η/2, whereas beyond this regime two cases should be distinguished. When η ≤ 1 the parameter m e can vary between 0 and 1, the two extreme values corresponding to a uniform and a dark-soliton configuration, respectively (see Appendix A). Consequently the contrast increases from 0 to a value η/(2 − η) ≤ 1 at increasing m e [see Fig. 3

(a)].
In particular for m e = η = 1, which corresponds to a black soliton, one has C = 1. The situation is different for η > 1, where m e can only vary in a smaller range of values, as discussed earlier. As shown in Fig. 3(b), in this case the contrast always reaches its maximum value

Average velocity and superfluid fraction
The velocity field v = Θ x /m oscillates with the same wavelength Λ as the density. Its average value is where we define and Π (−n e |m e ) = Π (−n e ; π/2|m e ) is the complete elliptic integral of the third kind [30]. Notice that Eq. (21) can be rewritten in the natural form which indicates that f s is precisely the superfluid fraction of the system [40]. Our first equality in Eq. (22)  fraction for a supersolid introduced by Leggett [8,41]. Actually in these works the first row of Eq. (22) was shown to be an upper bound to the real superfluid fraction. It was derived using an Ansatz wave function that assumes all the particles in the superfluid to have the same phase. This assumption is weaker than the one we make in the present work using the Gross-Pitaevskii theory, in which all the atoms have the same wave function. This is why within this approximation Eq. (22) is found as an exact result.
Equations (21) small m e (an analogous relation was recently derived in Ref. [42] for a shallow sine-modulated supersolid). For η < 1 [ Fig. 3(c)] f s decreases at increasing m e down to a minimum, that is typically attained at some m e very close to 1; then, it undergoes a very steep ascent and goes back to 1 at m e = 1, where the cnoidal wave turns into a dark soliton (see Appendix A 2). Instead, when η = 1 the superfluid density continues to drop down to 0 as m e approaches 1. Also in the η > 1 regime [ Fig. 3(d)] f s monotonously decreases with m e from 1 to 0, the latter value being attained at m e = m max e , where the contrast of the fringes (20) is 1. Thus a cnoidal wave with strong modulations is very weakly superfluid, again in agreement with Leggett's arguments [8,41].
On the theory side, the situation encountered here is common also to the modulated configurations studied for dipolar Bose gases. Leggett's equation coincides with the superfluid density obtained from single-orbital density functional theory, a.k.a. extended Gross-Pitaevskii equation, and it becomes zero when the periodic structure has contrast C = 1 (see, e.g., Refs. [42][43][44]). Although a number of properties have been experimentally measured, the smallness of the sample and its short lifetime have precluded direct access to the superfluid density so far (see however Ref. [23] for a first try in this direction).
Let us also mention that in the stripe phase of spinorbit-coupled Bose gases the maximum achievable value of the contrast depends on the interaction strength in the various spin channels, and the deeply modulated regime with a small superfluid fraction is more challenging to reach [45].

Energy per particle
The energy per particle is given by where the prefactor accounts for the number of particles in each layer, N Λ =nΛ. The evaluation of the integral in the above expression can be simplified using Eqs. (4), (5), (9), and (18). The final result is ε gn We have checked that the minimization of ε with respect to m e and η at fixed average densityn and velocity |v| > c always gives a uniform state. Hence, unlike in the case of superfluid helium [10] and of Bose gases with finite-range interaction [27], here there is no spontaneous transition from a uniform to a layered structure when the fluid velocity crosses the critical one (equal to c in our case). For this reason cnoidal waves should be regarded as (nonlinear) excited states of the system.

IV. DYNAMIC PROPERTIES
This section is devoted to the study of the quantities characterizing the dynamic behavior of a cnoidal wave. We first derive the Bogoliubov equations (Sec. IV A). Then, in Sec. IV B we compute and discuss the excitation spectrum, whereas in Sec. IV C we study the dynamic structure factor, its moments and the sum rules they obey. We note here that the same problem has been studied by the mathematical physics community (see, e.g., Refs. [46][47][48][49] and references therein), with a focus different from ours, and mainly dedicated to the dynamical stability of the cnoidal wave.

A. Bogoliubov equations
We shall now use the Bogoliubov approach [50][51][52] to study small oscillations about the equilibrium configuration derived in Sec. III A. In the present context it is convenient to describe the collective modes in terms of the fluctuations of the density and the phase. To this aim we decompose the total density and phase as n(x) + δn(x, t) and Θ(x) + δΘ(x, t), respectively. At first order in δn and δΘ Eqs.
We look for solutions oscillating in time of the form This turns Eqs. (26) into an eigenvalue problem, which enables one to determine the frequency ω and the complex amplitudes δñ and δΘ. The latter obey the normalization condition [52] i For each solution δñ and δΘ with frequency ω there exists another one, δñ * and δΘ * , having frequency −ω [50]. The integral of Eq. (28) evaluates to −1 (instead of 1) for the latter solution. Both solutions correspond to the same physical oscillation, as clear from the structure of Eqs. (27). In order to avoid this redundancy we shall only consider solutions having positive norm. This choice is customary because, in a second-quantization framework, it is naturally associated to the usual boson commutation relation.

B. Excitation spectrum
The procedure for solving this eigenvalue problem is similar to that employed in the previous works [53,54], and is detailed in Appendix B. Since the coefficients of the linear coupled equations (26) are periodic in x, we can look for solutions δñ and δΘ in the form of Bloch waves [55]. They are given by a plane wave, with wave vector q, times a periodic function with period Λ [see Eqs. (B3)]. To any fixed value of q there correspond infinitely many solutions, with different amplitudes and frequencies. This is at the origin of the band structure exhibited by the Bogoliubov spectrum. This structure is clearly visible in Fig. 4, where we plot the lowest three branches of the spectrum of elementary excitations of two given cnoidal-wave solutions. To distinguish between the various Bogoliubov modes we make use of two subscripts, the quasimomentum q and the band index = 1, 2, . . .. The spectrum is periodic in q, with period Q = 2π/Λ equal to the wave vector of the density modulations. Each range of values of q enclosed between consecutive integer multiples of Q defines a Brillouin zone. Notice that the frequencies ω ,q are not invariant under inversion of q into −q; this reflects the fact that cnoidal-wave solutions do not enjoy parity and time-reversal symmetry separately when the current J they carry is not zero. For the sake of comparison, in each panel of Fig. 4 we also plot (dashed curve) the spectrum of a uniform Bose gas having the same average densityn and velocityv as the cnoidal wave considered in the panel.
The main feature of the spectra of Fig. 4 is that the two lowest bands ( = 1, 2) are gapless and have linear dispersion close to the edges of each Brillouin zone. The higher mode ( = 2) at small positive q and the lower one ( = 1) at negative q are already present in a uniform system (dashed curve), whereas the other two branches are specific of cnoidal waves. 1 The presence of two gapless Goldstone modes is a feature expected for the spectrum of a one-dimensional supersolid intended as a system which breaks both U(1) and continuous translational invariance (see, e.g., Ref. [56] for a detailed discussion). Such an increase of gapless modes has been indeed theoretically discussed for the supersolid phase of solid Helium [6], of soft-core Bose gases [57][58][59][60], of dipolar Bose gases, as well as for the stripe phase in spin-orbit-coupled Bose gases [53,61]. It is under very active experimental investigation for dipolar gases [19][20][21][22].
Different from other supersolids considered in literature, here the lowest gapless mode has negative frequency. This means that cnoidal waves in one-dimensional Bose gases with repulsive contact interaction are energetically unstable. This agrees with the fact that they are excited states of the system, as discussed in Sec. III B 3. existence reveals that these waves tend to decay to a lower-energy state in the presence of external perturbations. However, if this decay takes place over sufficiently long timescales, measurements of the dispersion relation based on Bragg spectroscopy techniques would still be feasible. It is worth mentioning that the situation is not very different from that of dipolar gases. Indeed, due to three-body losses (energetic instability) the lifetime of the supersolid phase in those systems is very short (few to tens milliseconds) but many measurements, from phase coherence to collective excitations to Bragg spectroscopy, have been performed [16][17][18][19][20][21][22][23][24]. This phenomenon is often referred to as "transient supersolidity".
The nature of the two gapless modes can be understood looking at the limit where their frequencies vanish. This happens when q lies on an edge of a Brillouin zone. Setting δn t = 0 and δΘ t = 0 in Eqs. (26) one finds two kinds have been found in soliton trains trapped in ring geometries [32]. of solutions. The first one is δn = 0 and δΘ equal to a constant; it corresponds to an infinitesimal U(1) transformation of the phase of the condensate wave function. The second solution is δn = n x δx 0 and δΘ = Θ x δx 0 , which performs a translation of the wave function by an infinitesimal displacement δx 0 . This finding further reinforces the idea that the appearance of these modes is a result of the spontaneous breaking of U(1) and continuous translational symmetry.
For modes with nonzero frequency one can still distinguish between a crystal and phase character. The former involves small oscillations of the density peaks about their equilibrium positions; the latter features a superfluid current of particles tunneling from one peak to another [19,21,58,60]. However, hybridization can occur and both characters can be present in a single mode when ω ,q = 0. In cnoidal waves with small m e , the = 1 branch at q 0 and the = 2 one at q 0 have a dominant phase character, which is explained by their closeness to the corresponding modes of a uniform gas [dashed curve of Fig. 4(a)]; conversely, the branches that do not appear in uniform systems are mainly crystal modes. This change of behavior when crossing q = 0 becomes less and less pronounced at increasing m e because of stronger hybridization. When m e is large [close to 1 for η ≤ 1 and to m max e for η > 1, as in Fig. 4(b)] we find that the = 1 branch becomes a dominantly crystal mode for both positive and negative q. In particular, in the η ≤ 1 and m e 1 case the frequency of this band is almost 0 at all q, and in the m e → 1 limit it reduces to the zero-frequency mode of the excitation spectrum of a dark soliton. For η > 1 and m e close to m max e the phase character of both gapless modes is further suppressed because of the strong reduction of the superfluid fraction pointed out in Sec. III B 2.
As we mentioned in Sec. III B 1, at η ≤ 1, m e can reach values close to unity, and in this regime cnoidal waves can be regarded as chains of dark solitons. It is proven in Appendix C that in this case the dispersion relation of the lowest mode has the following analytic expression: We have checked that this expression reproduces very accurately the lower branch of the spectrum in the regime where m e → 1 and η − 1 2 ln(1 − m e ), see Appendix C. Finally, we examined the regions of the spectrum where two different bands approach each other and tried to determine whether they cross or not. Our numerical results suggest that there is no gap separating any couple of adjacent bands, and thus we are in the presence of a phenomenon of level crossing. Hence, the usual argument of gap opening because of Bragg scattering at the boundary of the Brillouin zone [55] does not seem to apply here, presumably because cnoidal waves do not scatter linear excitations.

C. Dynamic structure factor and sum rules
The dynamic structure factor provides important information on the dynamic behavior of the system. It is given by where the sum is extended over all the bands, and δρ q is the q-component of the density fluctuation operator. Its matrix operator between the ground state |0 and the -th excited band | can be easily computed: For a given integer p one defines the p-th moment of the dynamic structure factor as [52] We first consider the p = 0 moment m 0 (q) = N Λ S(q), where we have introduced the static structure factor We plot S(q) in Fig. 5 for the same values of the parameters as Fig. 4. We also plot the contributions of the two gapless Bogoliubov modes (sometimes referred to as the "strengths" of δρ q ). These contributions are not symmetric under inversion of q into −q for the same reason the spectra of Fig. 4 are not; however, as was shown in Ref. [62], the full static structure factor is symmetric as a consequence of its definition, regardless of the properties of the ground state. The contributions of the gapless modes to S(q) are dominant at small |q|. For a shallow cnoidal wave (small m e ) the lower gapless mode ( = 1) exhausts S(q) at q < 0 and the upper one ( = 2) at q > 0, as visible in Fig. 5(a). This behavior has the same explanation as that of the excitation spectrum (see Sec. IV B), namely, it stems from the closeness of these shallow waves to uniform gases. It is also shared by all the moments m p (q) with p = 0. As m e increases, the strength of the lower mode grows significantly and eventually, at high m e , it dominates S(q) in a wide range of q (both positive and negative), as shown in Fig. 5(b).
Another remarkable feature is that the strengths of both gapless modes [and consequently S(q) itself] diverge when q equals an integer multiple of Q, i.e., at the edges of each Brillouin zone (except at q = 0). An analogous behavior occurs for the supersolid phase of dipolar gases [63], as well as for the stripe phase of spin-orbit-coupled Bose-Einstein condensates [53], where, using sum-rule arguments, it was shown that the existence of a nonzero crystalline order parameter causes a |q − Q| −1 divergence of S(q) at the boundary of the first Brillouin zone.
The p = −1 moment is related to the compressibility χ(q) by The behavior of the compressibility, as well as that of the contributions of the two gapless modes, is displayed in Fig. 6. Notice that χ(q) is dominated by the lowest negative-frequency mode for a wide range of values about q = 0, and is thus itself negative in this range. Like the static structure factor, it diverges at the edges of the Brillouin zones except q = 0, again in agreement with the findings of Ref. [53].
Finally, we have checked that the p = 1 moment satisfies the f -sum rule m 1 (q) + m 1 (−q) = N Λ 2 q 2 /m [52]. Different from the sum rules discussed previously, for large m e and small |q| the f -sum rule is dominated mainly by the upper gapless mode. This is because the lower mode, despite having bigger strength, has much smaller frequency (in absolute value) than the upper one in this regime.

V. CONCLUSION
We have studied several relevant features of an ultracold Bose gas in a cnoidal-wave state. The equilibrium wave function is characterized by periodic spatial density modulations described in terms of Jacobi's elliptic functions. Cnoidal waves spontaneously break both U(1) and continuous translational symmetry, thus exhibiting typical supersolid features. Besides, as argued by Leggett [8,41], their superfluid fraction is depleted even at zero temperature, and gets smaller and smaller as the contrast of the fringes increases. A further signature of supersolidity is represented by the behavior of the excitation spectrum, featuring a band structure with two gapless modes. The latter exhibit a mixed phase and crystal character. The presence of a crystalline structure causes the divergence of the static structure factor and the compressibility at the edges of the Brillouin zones.
Our results open new perspectives for the study of supersolidity in ultracold atomic gases. Cnoidal waves are excited states that could be realized, for instance, moving an obstacle into the gas at an appropriate speed [64][65][66][67][68][69]. The density modulations can be probed either in situ or after time of flight. The excitation spectrum and the dynamic structure factor can then be accessed using twophoton Bragg spectroscopy.
The spectrum of a cnoidal wave demonstrates dynamical stability (all the eigenvalues are real) and, more im-portantly, energetical instability (some of the eigenvalues are negative). This aspect may be relevant in the context of analog gravity: it has been shown [35,69] that in some circumstances, an obstacle moving at supersonic speed in a Bose-Einstein condensate may give rise to an upstream cnoidal wave and a downstream supersonic flat density pattern, both stationary in the reference frame of the obstacle. It would then be of great interest to study the analogous Hawking radiation in this realistic, and experimentally relevant setting, where negative-norm modes exist on both sides of the acoustic pseudo-horizon [70].
We finally note that cnoidal waves have already been experimentally realized in the framework of nonlinear optics, both for repulsive and attractive interaction, in nonlinear whispering gallery modes resonators [71][72][73][74], but also in two-dimensional photo-refractive media [75,76] and in optic fibers [77]. One could as well experimentally study their excitation spectrum and other supersolidity phenomena, either in the above quoted contexts, or also by considering the propagation of a flow of microcavity polaritons. As discussed in Ref. [28], the cnoidal-wave solution admits several important limiting cases. In this appendix we shall focus on the uniform and linear-wave limit at small m e , as well as on the dark-soliton limit corresponding to m e → 1.

Uniform limit
Let us first look at the m e = 0 case. In this situation the amplitude of the density oscillations vanishes. Besides, the velocity (21) simplifies tov = ±c √ 1 + η. This yields η = η 0 with η 0 =v 2 /c 2 − 1, and implies that the flow is supersonic, i.e., |v| ≥ c (recall that η > 0 by definition). The chemical potential (9a) takes the standard form µ = gn + mv 2 /2 (the same happens in the linear-wave and dark-soliton limits discussed below).

Linear-wave limit
At first order in m e one can approximate Γ(m e ) − 1 −m e /2 and replace the sn function with an ordinary sine in Eq. (19). As a consequence, in this limit (nonlinear) cnoidal waves reduce to (linear) sinusoidal waves, where Q 0 = 2 √ η 0 /ξ = 2m √v 2 − c 2 / . It is interesting to study these waves in a frame where they travel with a given phase velocity V . In such a frame, at every point in space the density and the velocity field oscillate in time with frequency Ω 0 = Q 0 V . Let us set v =v − V , wherev is the average velocity in the new frame. One can invert the definition of Q 0 given just above to express V , and hence Ω 0 , as a function ofv and Q 0 . One finds the two values These are the two branches of the Bogoliubov spectrum of the uniform condensate discussed in the previous section (which flows with velocityv in the new frame). Besides, the ratio between the amplitudes of the phase and density modulations is ∓ (1 + η 0 )/(4η 0 ) = ∓m|Ω 0 − Q 0v |/ Q 2 0 , in agreement with the prediction of the Bogoliubov theory [52]. We thus conclude that in the small-m e limit traveling cnoidal waves describe standard Bogoliubov modes.

Dark-soliton limit
In the m e → 1 limit one has Γ(m e ) → 0 and the sn function approaches the hyperbolic tangent. The velocity (21) becomesv = ±c √ 1 − η. This gives η = 1−v 2 /c 2 , together with the condition of subsonic flow |v| ≤ c. Thus, the density and phase turn into those of a dark soliton [51,52], written in the frame where the latter is at rest and the background has velocityv: Θ(x) = mvx + arctan tanh(cos θ x/ξ) tan θ . (A5) Here we have set sin θ =v/c and cos θ = √ η, with θ ∈ [−π/2, π/2]. In particular, whenv = 0 one obtains a black soliton, whose density vanishes at the center. In this appendix we present the method we have used for solving the Bogoliubov equations (26). Using the Ansatz (27), and equating the coefficients of e −iωt and e iωt on both sides, they can be cast in matrix form as follows: Here we have defined the four operators Because of the periodicity of the coefficients entering these operators, the solutions of Eq. (B1) can be expressed as Bloch waves, δΘ ,q (x) = e iqx ν∈Z δΘ ,q,ν e iνQx .
Here q denotes the quasimomentum, δñ ,q,ν and δΘ ,q,ν are Fourier expansion coefficient, and the sums run over all integer ν. As explained in Sec. IV B, the band index is needed to distinguish between different solutions at fixed q.
The normalization condition for δñ ,q and δΘ ,q follows from Eq. (28) and reads iΛ δñ † ,q δΘ ,q − δΘ † ,q δñ ,q = 1. The above procedure leads to the eigenvalue equation where the B's are infinite-dimensional matrices with en-tries B (nn) ν1ν2 (q) = Jñ Numerically solving Eq. (B6) one recovers all the results of Sec. IV. Of course, in order to reduce the problem to a finite-dimensional one it is necessary to fix a cutoff ν max , and truncate all the above Fourier expansions retaining only the terms with −ν max ≤ ν ≤ ν max . The choice of the best value of ν max depends on m e and η. At fixed η and small m e , where cnoidal waves do not significantly deviate from linear waves, taking ν max equal to 5 or 6 can be sufficient to achieve good accuracy in the results. Conversely, increasing m e one needs larger and larger values of ν max . These can even exceed 100 when the cnoidal wave is close to the soliton limit (for η ≤ 1) or its contrast is close to 1 (for η > 1).

Appendix C: The lower excitation branch of a train of dark solitons
In the regime where the period Λ of the cnoidal wave is large compared to the width ξ/ √ η of a dark soliton, the cnoidal wave can be considered as a train of regularly spaced identical solitons. From Eq. (15) this occurs when K(m e ) 1, i.e., when m e is close to unity. Since solitons are essentially classical objects, it is natural to expect that, in this regime, the lowest branch of the spectrum should be described as an excitation of an array of classical particles connected by springs. Denoting by Ω the resonant angular frequency associated to these springs, the corresponding spectrum is of the form [55] The value of Ω depends on the interaction between two solitons and on their inertial mass. It can be determined by means of Manton's method [80,81] as explained now. For studying the interaction between two solitons, one considers a configuration where the solitons are stationary, in a background with subsonic velocityv and otherwise uniform densityn. It is convenient to single out the velocity of the background and to write ψ(x, t) = φ(x, t) exp(ikx) where k = mv/ . Then, one has in (1) µ = 2 k 2 /2m + gn and φ is solution of: An Ansatz describing two identical stationary solitons separated by a distance ∆ is of the form where with θ ∈ [−π/2, π/2]. √n Φ(x) describes a stationary isolated soliton, solution of Eq. (C2). The soliton is stationary because its velocity V = −c sin θ is exactly opposed to the velocityv = c sin θ of the background. Notice that, up to a global phase factor, one has √n Φ(x) exp(ikx) = √ n e iΘ , with n(x) and Θ(x) given by Eqs. (A4) and (A5), respectively. As regards the twosoliton case, of course the Ansatz (C3) is not an exact solution of the Gross-Pitaevskii equation (C2), but it is expected to be a reasonable approximation if 3 ∆ cos θ ξ. The Lagrangian density associated to the Gross-Pitaevskii equation (C2) is Note the unfamiliar multiplicative term (1 −n|φ| −2 ) in the first term of the above expression. It corresponds to adding to the usual Lagrangian density a total derivative which does not affect the form of the Gross-Pitaevskii equation (C2), but yields the correct physical momentum of a soliton [52,[82][83][84], namely for a soliton of type (C4). Considering two points a and b located around the soliton centered at ∆ (a < ∆ < b), one has where is the momentum density and the stress tensor, both associated to the Lagrangian density (C5). If ∆ − a and b − ∆ are large compared to ξ/ cos θ, then the left-hand side of Eq. (C7) can be identified with the time derivative dP/dt = 4 nθ cos 2 θ of the momentum (C6) of the soliton centered around ∆. Manton's method amounts to identifying, in the righthand side of Eq. (C7), the contribution due to the soliton centered at the origin, from which one can infer the force exerted by one soliton onto the other. Retaining only the leading order of this contribution and discarding all the other contributions leads to θ −8 gn cos 4 θ exp(−2∆ cos θ/ξ) .
In this formulaθ can be expressed in term of the time derivative of the velocity V = −c sin θ of the soliton with respect to the background. The fact that V changes means that the soliton under scrutiny does not remain stationary and moves with an acceleration∆ =V = −cθ cos θ = f (∆), where f -which is easily evaluated from Eq. (C9)is the ratio of the force experienced by the soliton centered around ∆ to its inertial mass (m I = −4m ξn cos θ, see Ref. [85]). Both the force and the mass are negative and this results in a repulsive interaction between the solitons. Once the interaction between two dark solitons has been determined, it is easy to turn to the case of a chain of solitons, considered as a one-dimensional lattice of classical particles. For determining the elementary excitations of such a system, one writes the spacing between two successive solitons as ∆(t) = Λ + X(t) and the angular frequency Ω of the equivalent spring is just which, together with Eq. (C1), yields the result (29) in the regime where the cnoidal wave becomes a chain of well-separated solitons. In this regime, the spacing between nearest solitons being large, the intensity of their interaction is weak and the lowest branch has a decreasing amplitude: in the dark-soliton limit of Sec. A 3 it becomes a zero mode corresponding to the translational degree of freedom of an isolated soliton. The accuracy of the approximation (29) is illustrated in Fig. 7 in the case of a cnoidal wave with m e = 0.99, for several values of η, ranging from 0.1 to 0.9. As one can see, the agreement is excellent for low values of η and becomes less accurate when η increases. This could be considered as strange because we repeatedly stated that the validity of our approximation should only rely on the fact that the separation Λ between two successive solitons is large compared to the soliton's width ξ/ cos θ, and the ratio of these two quantities only depends on m e , not on η [see Eq. (15)]. This conundrum is solved by inspecting the two-soliton Ansatz (C3). This Ansatz is valid for evaluating the interaction between two nearest solitons inasmuch as the ground state of the train of solitons itself can be described by an approximate wave function of the type From the density profiles plotted in Fig. 8, it is clear that the validity of expression (C11) decreases for increasing values of η, since the density n 2 of the flat region between two solitons significantly exceedsn as η increases, contrarily to the situation depicted by Eq. (C11). From expression (18b) one sees that, in the regime 0 < 1 − m e 1, n 2 remains close ton when the additional condition