Supersolid phase of a spin-orbit-coupled Bose-Einstein condensate: a perturbation approach

The phase diagram of a Bose-Einstein condensate with Raman-induced spin-orbit coupling includes a stripe phase with supersolid features. In this work we develop a perturbation approach to study the ground state and the Bogoliubov modes of this phase, holding for small values of the Raman coupling. We obtain analytical predictions for the most relevant observables (including the periodicity of stripes, sound velocities, compressibility, and magnetic susceptibility) which are in excellent agreement with the exact (non perturbative) numerical results, obtained for significantly large values of the coupling. We further unveil the nature of the two gapless Bogoliubov modes in the long-wavelength limit. We find that the spin branch of the spectrum, corresponding in this limit to the dynamics of the relative phase between the two spin components, describes a translation of the fringes of the equilibrium density profile, thereby providing the crystal Goldstone mode typical of a supersolid configuration. Finally, using sum-rule arguments, we show that the superfluid density can be experimentally accessed by measuring the ratio of the sound velocities parallel and perpendicular to the direction of the spin-orbit coupling.


Introduction
Bose-Einstein condensates (BECs) with synthetic spin-orbit coupling have been the subject of intense investigations in the last decade (see the reviews [1][2][3][4][5][6][7][8] and references therein). They represent an ideal framework for observing new quantum phenomena, resulting from the interplay between the modified single-particle spectrum and the interaction. One of the most intriguing configurations appearing in the phase diagram of a spin-orbit-coupled BEC is the so-called stripe phase, where the atoms condense in a superposition of states with finite momenta. Its importance lies in being a crystalline structure with superfluid properties, meaning that it exhibits the phenomenon of supersolidity (see the review [9]). This feature prompted important efforts towards the experimental observation of the supersolid behavior, which was eventually first reported in a spin-orbit-coupled configuration [10] and in a BEC gas inside optical resonators [11,12]. A later study on the stripe phase of BECs with Ramaninduced spin-orbit coupling showed the presence of long-range coherence between the different spin components [13]. More recently, the supersolid phase has been identified in a series of experiments with dipolar Bose gases, in the form of linear [14][15][16] and two-dimensional [17,18] droplet arrays; in such systems the investigation of collective modes [19][20][21][22], rotational and superfluid properties [23][24][25], out-of-equilibrium phase coherence effects [26], and the role of temperature [27] represents a very active field of research. Although the above experiments were carried out in three dimensions, the appearance of a stripe phase was also predicted by Monte Carlo calculations in two-dimensional configurations, where the polarization direction of the dipoles is tilted with respect to the axis perpendicular to the plane in which they are confined; however, it is still a subject of debate whether such stripes are superfluid or not [28,29]. It is well theoretically understood that a common feature of the excitation spectra of supersolids is the presence of two or more Goldstone modes [30][31][32][33][34][35][36][37][38][39][40][41][42]. In general, one of them is associated with the superfluid motion of particles through the crystalline structure, which stays at rest. The other modes involve oscillations of the density fringes about the equilibrium position, which are typical of crystalline solid bodies. Actually, the superfluid and crystal characters completely separate only in the limit of infinite wavelength; far from this limit they are both present to some extent, leading to mode hybridization [19][20][21]. This physics is well understood in single-component systems, such as dipolar gases or BECs in optical resonators. Spin-orbit-coupled BECs in the stripe phase are of special interest because of their additional spin degree of freedom. Its role is well understood in Bose-Bose mixtures without spin-orbit coupling, whose Bogoliubov spectrum is made of a density and a spin branch [43,44]. In the infinite-wavelength limit the density and spin branches are both gapless and correspond to the dynamics of the total and relative phase between the two spin components, respectively. In the presence of spin-orbit coupling the nature of the spin excitations is instead far from being obvious. On the one hand, one expects that the spin modes are gapped as a consequence of the Raman coupling. On the other hand, a second gapless mode is predicted to occur in the stripe phase [36] as a result of the spontaneous breaking of translational invariance.
The main purpose of this paper is to gain insight into the Goldstone modes of a spinorbit-coupled BEC in the stripe phase, and in particular to point out the spin nature of the gapless crystal mode. This aspect has been recently numerically investigated in trapped configurations [45]. Here we consider homogeneous setups where analytic results can be obtained in an explicit way. For concreteness we will focus on BECs featuring a combination of the Rashba [46] and Dresselhaus [47] spin-orbit coupling with equal weights. Such systems can be realized by coupling a condensate to the light field generated by a pair of Raman lasers, detuned in frequency and with orthogonal polarizations. This setup was first employed by Ian Spielman's group at NIST [48]. When the Raman coupling term between the two spin components is zero the ground state of the system corresponds to an unpolarized BEC mixture. We will carry out our analysis in the limit where the Raman coupling is small, and study how this perturbation modifies the properties of the ground state and the Bogoliubov modes, particularly those in the phonon regime. Further information about the compressibility, magnetic susceptibility, and superfluid density will be obtained combining the perturbation approach with the method of sum rules.
This work is structured as follows. In Sec. 2 we briefly present the model describing a Raman-induced spin-orbit-coupled BEC and review its main properties, with a special focus on the stripe phase. Section 3 deals with the results of our perturbation approach for a BEC at equilibrium in the stripe phase, providing analytic formulas for several important observables. The perturbation method is subsequently extended to the study of the Bogoliubov modes in Sec. 4: we compute the velocities of sound waves and the corresponding amplitudes, the most relevant sum rules, and the superfluid density. We summarize in Sec. 5. The details of the calculations are reported in the Appendices.

The model
In this section we briefly introduce the model describing the system under consideration. Section 2.1 is devoted to the single-particle physics, while in Sec. 2.2 we summarize the properties of the many-body ground state obtained including interactions within the meanfield approach. More details about the stripe phase are discussed in Sec. 2.3.

Single-particle physics
A three-dimensional (3D) spin-1/2 BEC with Raman-induced spin-orbit coupling is described by the single-particle Hamiltonian Here m is the atom mass, p 2 ⊥ = p 2 y + p 2 z , and σ j (j = x, y, z) denote the usual 2 × 2 Pauli matrices. The Raman recoil momentum k R plays the role of the strength of the spin-orbit coupling; it also fixes the value of the Raman recoil energy, E R = 2 k 2 R /2m. The Raman coupling is quantified by Ω R . In writing Eq. (1) we have omitted a Zeeman shift term δ R σ z /2, where δ R is the detuning of the Raman light field from resonance. The latter can be set equal to zero by properly choosing the frequency detuning ∆ω L of the Raman lasers. Equation (1) actually represents the Hamiltonian of the system in a spin-rotated frame, that is connected to the laboratory frame by the space-and time-dependent unitary transformation [49] Hamiltonian (1) enjoys translational invariance as it commutes with the canonical momentum p. Hence, one can look for eigenstates in the form of plane waves. The energy spectrum includes two branches, a lower and an upper one. For Ω R < 4E R the lower branch exhibits two degenerate minima at p = ± k SP 1 e x , where e x denotes the unit vector along the x direction and Instead, when Ω R ≥ 4E R one has a single minimum located at p = 0.

Many-body physics in the mean-field approach
In the weakly interacting regime a spin-orbit-coupled BEC, with the single-particle Hamiltonian of the form (1), can be safely described by mean-field theory. Indeed, the addition of the spin-orbit coupling has limited effects on the quantum depletion of a three-dimensional condensate [50,51]; we checked that, for the set of parameters used in the present work, the quantum depletion remains on the order of a few percent. In the mean-field analysis the state of the system is represented by a two-component (spinor) wave function normalized to the total number of particles according to where V is the volume enclosing the atoms. We point out that, unlike N , the numbers of atoms in the up and down component, given by N ↑,↓ = V d 3 r |Ψ ↑,↓ | 2 , are not separately conserved because of the Raman coupling term in the spin-orbit Hamiltonian (1). The order parameter Ψ(r, t) relative to the ground-state many-body Bose-Einstein-condensed configuration evolves in time according to Ψ(r, t) = Ψ 0 (r)e −iµt/ , with µ the chemical potential [44]. The total energy, including both single-particle and interaction terms, reads Here n = Ψ † Ψ and s z = Ψ † σ z Ψ are the total and spin density, respectively; notice that these quantities are unaffected by the unitary transformation (2) and take the same value in both the laboratory and the spin-rotated frame. In writing Eq. (6) we have assumed that the atoms in the BEC interact via a two-body contact potential. The density-density and spinspin coupling constants are given by g dd = (g + g ↑↓ )/2 and g ss = (g − g ↑↓ )/2, respectively. In the following we will always assume equal intraspecies couplings, i.e., g ↑↑ = g ↓↓ ≡ g, together with the conditions g > 0 and g dd > 0 ensuring stability against collapse. In the more general case (not addressed in this work) of asymmetric intraspecies couplings one should replace g with (g ↑↑ + g ↓↓ )/2 and add a term g ds ns z , with g ds = (g ↑↑ − g ↓↓ )/4, inside the integral of Eq. (6). The couplings in each channel are related to the corresponding s-wave scattering lengths via g σσ = 4π 2 a σσ /m (σ, σ =↑, ↓). The zero-temperature phase diagram of the spin-orbit-coupled BEC under consideration was first investigated in Refs. [52,53]. In these papers an Ansatz for the ground-state wave function Ψ 0 was introduced. This Ansatz consists of a linear combination of two plane waves with opposite momenta along x; each plane wave is multiplied by a two-component spinor wave function. All the quantities characterizing the order parameter Ψ 0 (the wave vectors of the two plane waves, their relative weights, and the components of the corresponding spinor wave functions) are determined by a procedure of minimization of the energy (6). The results of this variational procedure depend in a nontrivial way on the value of the average density of the gasn = N/V . At relatively low densities and for g ss > 0 the ground state of the BEC is compatible with three distinct quantum phases: 1. At low values of the Raman coupling Ω R the system is in the stripe phase. In this phase the two momentum components in the wave function have equal weights, yielding a vanishing value of the spin polarization σ z = V d 3 r s z along the z axis. When calculating the condensate density, these two components interfere, giving rise to a periodically modulated density profile along the x direction. The appearance of density modulations signals the spontaneous breaking of translational symmetry, occurring simultaneously with the spontaneous breaking of U(1) invariance. These features are typical of supersolid configurations. For this reason we will often refer to the stripe phase as to the supersolid phase.
2. At intermediate Ω R the system enters the plane-wave phase, where the atoms condense in a state with finite momentum. Unlike the stripe phase, in the plane-wave phase the spin polarization σ z takes a nonzero value.
3. Finally, at large Ω R the gas is in the zero-momentum phase, where both the condensation momentum and the spin polarization along z vanish.
The transition between the stripe and the plane-wave phase has a first-order nature. In the limit of vanishingly smalln it occurs at a critical value of the Raman coupling given by [52] Ω ST−PW = 4E R 2g ss g dd + 2g ss .
At finite average density the value of Ω ST−PW can be computed imposing that the chemical potential and pressure be equal in the two phases. One finds that at the transition such phases are at equilibrium with different values of their densities. However, the density differences are typically extremely small [49]. Hence, in the present work, and particularly in Figs. 1, 2, 3, 4, and 7, we estimate Ω ST−PW without taking this effect into account and using the density as a fixed thermodynamic parameter. Another feature that is visible in the above figures is that, because of the first-order nature of the transition, the stripe phase can be found as a metastable configuration even above the critical Raman coupling Ω ST−PW , and up to the so-called spinodal point. The situation is different for the transition between the plane-wave and zero-momentum phases, which has a second-order nature and takes place at the critical coupling [53] The plane-wave and zero-momentum phases, as well as the corresponding phase transition at the critical point (8), also occur when g ss < 0. Instead, the stripe phase does not exist under this condition, as it favors immiscible configurations [53]. We conclude by mentioning that for g ss > 0 in the large-density regime the plane-wave phase is no longer energetically favored, and one has a direct first-order transition from the stripe to the zero-momentum phase [53,54].

The stripe phase. General remarks
In the rest of this paper we will focus almost exclusively on the stripe phase. As mentioned in the previous section, in the analysis of Refs. [52,53] the condensate wave function in this phase was taken as a superposition of two plane waves with opposite momenta. However, this is only an approximation to the exact wave function, that we write in the form [36] Ψ 0 (r) = √n withm = ±1, ±3, ±5, . . .. Here the spinor expansion coefficientsΨm = Ψm ,↑Ψm,↓ T obey the condition mΨ †mΨm = 1, which ensures the correct normalization of Ψ 0 . They also enjoy the symmetry propertyΨ −m = σ xΨ * m , 1 which implies the vanishing of the spin polarization σ z . Notice that the expansion (9) contains an infinite number of harmonic terms, each oscillating at an odd multiple wave vector of k 1 , this last being a parameter to be adjusted to find the ground state. The variational Ansatz employed in Refs. [52,53] can be recovered by truncating the expansion to the two terms withm = ±1. The presence of higher-order harmonics is due to the nonlinear terms in the Gross-Pitaevskii equation [see Eq. (12)].
The order parameter (9) can be put in the form of a Bloch wave, i.e., a plane wave with quasimomentum k 1 times a periodic function with the same periodicity of the stripes. One can easily see that the wave vector of the density modulations corresponds to the wave-vector difference 2k 1 between any two consecutive terms. As we will discuss below, it differs from the momentum difference of the two single-particle minima, given by 2k SP 1 [see Eq. (3)], because of interaction effects. These features are also reflected in the periodicity properties of the structure (9). It is antiperiodic in x with antiperiod π/k 1 , that is, Ψ 0 (x + π/k 1 , y, z) = −Ψ 0 (x, y, z), and thus it is also periodic with period 2π/k 1 .
One can prove that k 1 is related to theΨm's in a simple way. This can be done inserting the Ansatz (9) into the energy (6). After performing the integration over the spatial coordinates one finds that k 1 explicitly appears only in the kinetic energy term along x, Since Ψ 0 is the ground state of the system, the condition of stationarity ∂E ∂k 1 = ∂E kin,x ∂k 1 = 0 must hold, yielding This expression reproduces the result of Ref. [53] if one truncates Eq.
In the above expression we have introduced dimensionless coordinates defined by the following scaling transformations The main advantage of this choice is that Ψ 0 as a function of the scaled x coordinate has period 2π; this will prove useful when dealing with the Ω R -dependence of k 1 . For consistency we have to scale also the system volume, V → V /(k 1 k 2 R ) (with the consequent scaling of the densityn and the wave function Ψ 0 ), and the interaction strengths, g dd,ss → k 1 k 2 R g dd,ss . Notice that the products G dd =ng dd and G ss =ng ss remain unaffected by these scaling transformations.
Equation (12) represents the starting point of the analysis of the next section. In particular, we will develop a perturbation approach to solve it in the limit of small values of the Raman coupling.

Perturbation analysis of the ground state
The present section is devoted to a perturbative study of the ground state of a spin-orbitcoupled BEC in the stripe phase. We first provide the exact solution in the limit of zero Raman coupling (Sec. 3.1). Then, in Sec. 3.2 we introduce the perturbation approach and present the results for several observables up to second order in the Raman coupling. The explicit calculation of the ground-state wave function is illustrated in Appendix A.

Ground state at Ω R = 0
In the absence of Raman coupling (Ω R = 0) the order parameter minimizing the spin-orbit energy functional (6) takes the form withΨ This configuration corresponds to the ground state of an interacting quantum mixture with N ↑ = N ↓ , equal intraspecies interactions, and g ss > 0. The plane-wave dependence of the wave function, fixed by the value of k R , is the consequence of the unitary transformation (2) employed to write the spin-orbit Hamiltonian in the spin-rotated frame [see Eq. (1)]. We stress that the state described by the wave function (14) has not a supersolid character. Actually, this configuration has uniform density (n (0) =n), as a consequence of the orthogonality of the two spinorsΨ The physical meaning of the relative phase χ 0 between the two spin components can be better understood by reverting to the laboratory frame. In this frame the spin polarization vector ( σ x , σ y , σ z ) = N (cos(∆ω L t + χ 0 ), sin(∆ω L t + χ 0 ), 0) lies in the (x, y) plane, its direction at a given time t being fixed by χ 0 ; hence, this miscible configuration can be regarded as an easy-plane ferromagnet. The arbitrariness of the value of χ 0 reveals that spin-rotation invariance about the z axis is spontaneously broken in the miscible phase. Notice that these considerations no longer hold at finite Ω R , as the Raman coupling term in the Hamiltonian (1) explicitly breaks spin-rotation symmetry.

Perturbation results for the ground state
In the following we will consider the Ω R /4E R 1 case, corresponding to the deep doubleminimum regime of the single-particle spectrum. In this regime the Raman coupling term can be treated as a perturbation. The procedure used for determining Ψ 0 is very similar to that of Ref. [55], which considered a single-component BEC in a one-dimensional optical lattice. We perturbatively solve Eq. (12) expanding the wave function around the unperturbed configuration (14): Here the superscript "(n)" denotes the correction of order ( Ω R /4E R ) n to the wave function of the condensate. Inserting the expansion (15), as well the corresponding expansions for k 1 and for the chemical potential, into Eq. (12), and equating terms of the same order on both sides, one finds recurrence relations that can be solved at each order. The explicit form of these recurrence relations [which take explicitly into account the normalization of the wave function and the constraints imposed by the U(1) and translational symmetries of the model], as well as of their solutions up to second order in Ω R /4E R , is reported in Appendix A.
Here we provide the relevant results for the expansion of the physical quantities of major interest. In particular the contrast of stripes, defined by C = (n max −n min )/(n max +n min ) with n min (n max ) the minimum (maximum) value taken by the density modulations as a function of x, is found to obey the linear dependence on Ω R /4E R , higher order corrections being of order ( Ω R /4E R ) 3 . In general, C is antisymmetric under Ω R → −Ω R because this operation exchanges the locations of the maxima and minima of the density modulations. It is worth mentioning that the actual position of the fringes is the result of a spontaneous breaking mechanism of translational invariance. In the present formalism it is fixed by the value of the relative phase χ 0 between the two components of the Ω R = 0 wave function (14) and by the condition (68b) imposed on the perturbative corrections (see the related discussion in Appendix A.1). Other physical quantities, like the chemical potential and the energy per particle, as well as the wave vector k 1 fixing the periodicity of the density modulations have instead vanishing first-order contributions, the lowest corrections to the unperturbed values being quadratic in Ω R /4E R . The perturbation approach provides the following results: The above perturbation results can be compared with the exact values, obtained by numerically minimizing the energy (6) for a wave function of the form (9). Such a comparison is presented in Fig. 1, where the same interaction parameters as Ref. [36] have been taken. Notice that for all the considered quantities the agreement is excellent within a wide range of values of Ω R , extending up to values close to the transition to the plane-wave phase. This confirms the validity and usefulness of the perturbation approach developed in the present work. The above choice of the interaction parameters actually amplifies the range of values of Ω R characterizing the stripe-phase region, yielding a significant increase of the critical Raman coupling for the transition to the plane-wave phase, which occurs at Ω R = 2.70 E R (even though the stripe phase continues to exist as a metastable configuration up to the spinodal point at Ω R = 2.85 E R ). Unfortunately these values of the parameters do not correspond to the current values of the scattering lengths of alkali atoms. In particular, because of the smallness of the ratio g ss /g dd ≈ 10 −3 in 87 Rb the transition to the plane-wave phase takes place at a much lower critical coupling, whose value Ω R = 0.19 E R was measured in [48] and found to be in excellent agreement with the theoretical prediction (7). Consequently, the emerging supersolid effects, like the contrast of fringes, are very small in this case. However, a decrease of the effective interspecies coupling constant g ↑↓ (and hence an increase of the critical value of Ω R ) can be achieved by reducing the spatial overlap between the wave functions of the two spin [56] or pseudo-spin [57] components. This strategy, which already led to the first successful detection of the stripe phase employing the pseudo-spin states of a superlattice [10], could open realistic perspectives for the observability of the sizable effects predicted in the present work. Another promising strategy is provided by the use of different atomic species, like 39 K, where the availability of Feshbach resonances [58] is opening novel perspectives to increase the critical value of the Raman coupling giving the transition to the plane-wave phase [45].

Bogoliubov modes in the stripe phase
The perturbation approach developed in the previous section can be naturally extended to the study of the Bogoliubov modes of our striped condensate. As a preliminary step, in Sec. 4.1 we briefly review the fundamentals of Bogoliubov theory and of its formulation in the stripe phase. In Sec. 4.2 we summarize the properties of the Bogoliubov modes in the absence of Raman coupling while in Sec. 4.3 we discuss the results for the most relevant features of the Bogoliubov solutions holding up to first and second order in the perturbation parameter Ω R /4E R . Special focus is given on the Goldstone modes exhibited by the stripe phase. Subsequently we address the computation of the sum rules in the long-wavelength limit (Sec. 4.4) and of the superfluid density (Sec. 4.5). In all these sections we systematically compare the perturbation approach with the exact solution of the Bogoliubov equations. The detailed perturbation calculation of the Bogoliubov modes is illustrated in Appendix B.

Basics of Bogoliubov theory. Results in the stripe phase
The starting point of our analysis is the time-dependent Gross-Pitaevskii equation whose linearized solutions, that we write in the form are known to be equivalent to the formulation of Bogoliubov theory [43,44,59]. In the above equation Ψ 0 (r) is the equilibrium configuration studied in the previous section, while the Figure  fluctuation term δΨ obeys the linearized Bogoliubov equation Here we have introduced the two matrices where I 2 denotes the 2 × 2 identity matrix. We look for solutions of Eq. (24) in the form Here ω is the frequency at which the wave function oscillates around Ψ 0 , with U and V the corresponding two-component (spinor) Bogoliubov amplitudes, and λ is a complex multiplicative factor such that |λ| 1. Inserting Eq. (26) into Eq. (24) one finds that the Bogoliubov frequencies and amplitudes are solutions of the eigenvalue problem Here, for notational compactness, we have defined the two 4 × 4 matrices and as well as the four-component Bogoliubov column vectors The column vector W is normalized to For further details concerning the Bogoliubov formalism applied to the spinor configuration see Appendix B.
If the condensate is in the stripe phase described by the mean-field wave function (9), the solutions of Eq. (27) can be expressed as Bloch waves of the form [36,60] with ω b,k the corresponding frequency. HereŨ b,k,m andṼ b,k,m are expansion coefficients, and k is a quasimomentum. We take the x component of k in the first Brillouin zone, i.e., 0 ≤ k x < 2, while k y and k z are unbounded. Notice that, for consistency with Eq. (13), the scaling has been performed, such that k is dimensionless. For a fixed value of k one finds an infinite number of solutions of Eq. (27), which implies that the Bogoliubov spectrum has a band structure; we use the additional index b = 1, 2, . . . to distinguish between the various bands (see also [36,60]). Here we mention that for each value of k the spectrum exhibits two gapless branches, whose frequency vanishes at the edge of the Brillouin zone (see Fig. 6 in Appendix B). The occurrence of these Goldstone modes is the consequence of the spontaneous breaking of U(1) and translational symmetry. The Bogoliubov formalism is particularly useful for studying the dynamics of the system when a weak perturbation, of the form V pert = − Oe −iωt + H.c. with | | 1 and O a given operator, is added to the single-particle Hamiltonian (1). Typical choices for O are the qcomponent of the density operator, ρ q = e iq·r , and of the spin-density operator, s z,q = σ z e iq·r . A fundamental object to compute is the dynamic structure factor [44] Here 0|O|b, k denotes the matrix element of O between the ground state and a given excited mode. Its square modulus is called the strength of O in the mode (b, k); a simple calculation yields Notice that for the above-mentioned operators ρ q and s z,q , as well as for the transverse current operator introduced in Eq. (58) below, the strength vanishes unless the difference q − k equals an integer multiple of 2k 1êx . We define the p-th moment of the dynamic structure factor as These moments are known to obey important sum rules [44]. For instance, the p = 0 moment gives the static structure factor S(O) = m 0 (O)/N . In this work we will make large use of the identity relating the p = −1 moment to the static response function χ(O).
In the following we implement a perturbation approach to study the Bogoliubov modes and analytically calculate several relevant quantities in the limit of small Raman coupling.

Bogoliubov spectrum at Ω R = 0
As discussed in Sec. 3.1, in the Ω R → 0 limit the condensate wave function approaches that of a uniform unpolarized mixture. Furthermore, because of the vanishing of the Raman coupling term, the spin-orbit Hamiltonian (1) commutes with the physical momentum q = p − σ zêx . Thus, we can label the excited states using the eigenvalues of the q operator. The Bogoliubov spectrum as a function of q is made of two branches, one featuring density oscillations, the other spin oscillations [43,44], with frequencies respectively. In the above equations Ω q = E R q 2 and we employed, for the components of q, the same scaling transformations used for the quasimomentum [see Eq. (33)]. The corresponding Bogoliubov amplitudes are given by where Ψ (0) 0 is the ground-state wave function [see Eq. (14)], and with = d, s the branch index. Notice that the dispersion laws (38) are isotropic, i.e., they do not depend on the direction of q. In the low-q limit they are linear in q, ω ,q ∼ c (0) q, and the sound velocities for density and spin waves, using the unscaled momentum q, are given by As soon as Ω R = 0 the physical momentum q is no longer a good quantum number, and has to be replaced by the quasimomentum k (see Sec. 4.1). Actually one can define k also at vanishing Ω R . In this case the two original branches (38) of the Bogoliubov spectrum split into an infinite number of bands, each identified by the band index b and defined in the first Brillouin zone, as discussed in Appendix B. As shown in Fig. 5, these bands can cross each other, and the crossing points (marked in black in the figure) correspond to Bogoliubov modes that are degenerate at Ω R = 0. This degeneracy is lifted as soon as the Raman coupling is turned on, leading to the opening of gaps between the excitation bands. For simplicity we will not deal with this mechanism, whose theoretical treatment requires the use of an involved degenerate perturbation method. On the other hand, we stress that the simpler nondegenerate perturbation approach developed in this work is well suited for studying the long-wavelength limit of the lowest-lying bands. This will be the core of the discussion of the next section.

Perturbation results for the Bogoliubov modes
We shall now perform a perturbation study of the Bogoliubov modes in the Ω R /4E R 1 limit, analogous to that of the ground state made in Sec. 3. We start by expanding the Bogoliubov frequencies and amplitudes, as well as the matrix B of Eq. (29), where the scaling transformations (13) are implicitly assumed. The formalism of the expansion procedure is discussed in details in Appendix B, where we derive explicit expressions for the frequencies and the Bogoliubov amplitudes of the various bands up to quadratic terms in the perturbation parameter Ω R /4E R .
Here we provide explicit results for the most interesting lowest-energy branches in the longwavelength limit (corresponding to = s, d,K = 0, and k → 0 in the notation of Appendix B), where they reveal their Goldstone nature. These modes, which are characterized by both small quasimomentum and small frequency, have a phonon-like dispersion law, ω ,0,k ∼ c (θ k )k. Here c is the sound velocity in the branch ; due to the anisotropy of our system, which is caused by the spin-orbit coupling, it generally depends on the angle θ k between k and the positive x direction [36,49] (see Appendix B).
The second-order expansion of the sound velocities relative to the two bands, calculated along the x direction (θ k = 0), yields the results The expansion of the sound velocities calculated along the directions perpendicular to x (θ k = π/2) instead gives The velocities of phonons propagating parallel and perpendicular to the x axis are plotted in the upper and lower panel of Fig. 2, respectively. Remarkably, the perturbative estimates (dashed lines) turn out to be very close to the exact numerical solutions of the Bogoliubov equations (solid lines) for a wide range of values of the Raman coupling. Equations (44) and (45) show that the dependence of the sound velocities on the interaction parameters is considerably more involved than in binary mixtures without spin-orbit coupling. Notice that, although in the figures of this paper we take G dd > G ss , our analytical formulas also hold in the opposite case. An important feature of the sound velocities plotted in Fig. 2 is that they remain positive even beyond the transition to the plane-wave phase. More generally, the stripe phase is dynamically and energetically stable (i.e., the excitation spectrum is real and positive) up to the spinodal point as a consequence of its metastability, see Sec. 2.2. As shown in Fig. 2(a), the velocity of spin waves along x vanishes at the spinodal point, where the stripe phase becomes dynamically unstable; this effect, which is associated with the divergence of the magnetic susceptibility (see Sec. 4.4), in trapped systems manifests itself in the softening of the spin-dipole frequency as the spinodal point is approached [45]. We point out that in the (Ω R , g ss ) plane the transition from the stripe to the plane-wave phase continuously connects to the miscible-immiscible transition, occurring at Ω R = 0 when g ss is tuned from positive to negative values. However, no metastability window exists across the miscible-immiscible phase transition, as the miscible phase becomes dynamically unstable precisely at the critical point g ss = 0, following the vanishing of the spin sound velocity c We next look at the low-k behavior of the Bogoliubov amplitudes of the two gapless bands, which can be related to the ground-state wave function Ψ 0 as follows: (in this section and in the next two ones we use dimensional coordinates and momenta), where the dependence on the Raman coupling Ω R is implicitly contained in Ψ 0 . One can immediately see that Eqs. (46) reduce to results (39) at zero Raman coupling, as a consequence of the identity Equations (46) have been evaluated including linear terms in Ω R /4E R . More general results accounting for higher-order corrections are presented in Appendix B.4. The full condensate wave function (23)-(26) for these two modes, including both the equilibrium and the fluctuation terms, is readily evaluated, The previous equations give direct insight on the physical nature of the two phonon modes in the macroscopic limit k → 0. The former result represent an infinitesimal U(1) transformation of the wave function, that changes its phase by δϕ = Im λ 4 4mG dd / 2 k 2 N 2 , where λ accounts for the size of the fluctuation term δΨ [see Eq. (26)]; this is the standard Goldstone mode associated with the superfluid character of a Bose-Einstein condensate, featuring, to leading order, a simple change of the phase. Equation (48b) instead corresponds to an infinitesimal translation of the wave function by a displacement δx = Im λ 4 4mG ss / 2 k 2 N 2 k −1 R along x; in particular the total density is characterized by shifted fringes compared to the equilibrium profile n(r). This is the crystal Goldstone mode that typically appears in solid configurations. The above discussion emphasizes in an explicit way the deeply different nature of the spin Goldstone mode in the presence of Raman coupling as compared to the spin mode in mixtures with vanishing Raman coupling. In the latter case the identity (47) holds, and employing it in Eq. (48b) one finds meaning that the spin mode at Ω R = 0 represents a mere shift by δχ = k R δx of the relative phase between the two spin components of the wave function. It rotates the direction of the spin polarization vector (see Sec. 3.1) by δχ but has no effect on the (uniform) density of the system. Result (48b) is particularly remarkable because it shows that in the long-wavelength limit the spin collective mode corresponds to the translational motion of stripes, reflecting the crucial interplay between the spin and density degrees of freedom caused by the spin-orbit coupling. This effect has been recently pointed out numerically in a harmonically trapped spin-orbit-coupled Bose gas by observing that the translational motion of the stripes can be actually induced by the sudden release of a uniform spin perturbation [45].

Sum rules at long wavelengths
Further information about our system can be obtained combining our perturbation approach with the sum-rule method (see Sec. 4.1). We focus on external perturbations proportional to the q-component of either the density operator, ρ q = e iq·r , or the spin-density operator, s z,q = σ z e iq·r . The corresponding moments can be deduced from Eq. (36) with O = ρ q , s z,q . Besides the Bogoliubov spectrum ω b,k , this requires the knowledge of the strengths (35); the perturbative expressions of the latter can be straightforwardly obtained from those of the ground-state wave function Ψ 0 and the Bogoliubov amplitudes W b,k (see Appendices A and B). It is possible to evaluate any moment at arbitrary q, but here for simplicity we restrict ourselves to the small-q limit. As shown in Eq. (37), the static response is directly related to the inverse-energy-weighted moment. We recall that the q → 0 value of the static density response coincides with the compressibility, whose second-order expansion takes a simple expression, Notice that the relation holds, reflecting the fact that the phonon mode propagating perpendicular to x exhausts both the inverse-energy-weighted sum rule and the f -sum rule [44] Instead, the phonon mode along x exhausts the inverse-energy-weighted sum rule but not the f -sum rule. In particular we have checked that the ratio between the part of the f -sum rule accounted for by the phonon mode as q → 0 and the compressibility sum rule (51) gives direct access to the square of the velocity (44a) of sound waves propagating along the x direction. The static spin response χ(s z,q ) approaches the magnetic susceptibility as q → 0. At second order in Ω R /4E R one has The fact that there is no counterpart to Eq. (52) in the spin channel, i.e., χ −1 M = mc 2 s,⊥ , has a simple explanation in terms of sum rules. Indeed, the energy-weighted spin sum rule with σ x the spin polarization along x, is gapped at q → 0 because of the Raman coupling. Thus, it cannot be exhausted by the phonon mode even in the directions perpendicular to x. We plot the inverse compressibility and magnetic susceptibility as functions of Ω R in the upper and lower panel of Fig. 3, respectively. As for the results of the previous sections, also here the perturbative formulas match the exact numerical estimates up to relatively large values of the Raman coupling. Notice, in particular, that χ −1 M vanishes at the spinodal point, pointing out the divergent behavior of the magnetic susceptibility and the occurrence of the dynamic instability discussed in Sec. 4.3. In the bottom panel we additionally compare the obtained magnetic susceptibility with the estimate This relation was derived in [61] using a two-harmonic Ansatz for the ground-state wave function, which is expected to be accurate in the limit of small densities. Notice that Eq. (56) diverges at the critical Raman coupling (7), reflecting the fact that the transition and the spinodal point tend to coincide at low densities. In addition, when the two conditions G dd,ss E R and Ω R 4E R are simultaneously satisfied, Eqs. (54) and (56)

Superfluid density
The results derived in the previous sections concerning the sound velocities and the sum rules can be usefully employed to calculate the superfluid density of a spin-orbit-coupled Bose gas in the stripe (supersolid) phase. We will follow the procedure developed in [62], based on Baym's approach [63] to the normal (non superfluid) density in terms of the macroscopic limit of the static response to the transverse current field: whereρ = mn is the average mass density and is the transverse current operator along the x direction (here q ⊥ is taken perpendicular to x). A crucial feature of spin-orbit-coupled gases is that the current along the x direction contains a novel spin-dependent term which is responsible for the violation of Galilean invariance. For this reason, even for uniform-density configurations, like the plane-wave and the zeromomentum phases, the superfluid density ρ s =ρ − ρ n , calculated at zero temperature, does not coincide with the total density [62] and is expected to be an anisotropic tensor since the y and z components of the current are not affected by the spin term.
Since the transverse current operator does not excite the gapless phonon mode, which is of longitudinal nature, the only contribution to Eq. (57) comes from the gapped branch. On the other hand, the q → 0 contribution associated with the gapped modes can be safely calculated by replacing q ⊥ with q xêx , a procedure that would not be allowed in the calculation of the q → 0 contributions of the gapless excitations. Making explicit use of the equation of continuity to relate the static longitudinal current response to the energy-weighted fsum rule (53), one then finds that the normal density is directly related to the contribution provided by the gapped excitations to the f -sum rule. By taking the difference ρ s =ρ − ρ n and noticing that the phonon modes exhaust the inverse-energy-weighted moment, fixed by the compressibility κ of the system (see Sec. 4.4), one finds the following useful relationship between the superfluid density, the velocity of sound, and the compressibility holding at zero temperature: Analogously one finds the result ρ y s /ρ = ρ z s /ρ = mc 2 d,⊥ κ for the superfluid density calculated for the motion parallel to the stripes (y and z directions). It reveals, as expected, that the superfluid density is not isotropic. As shown in Sec. 4.3, the sound velocities along the x and y (or z) directions are actually different. Since, as already pointed out in the previous section, the phonon mode propagating along the directions perpendicular to x exhausts the corresponding f -sum rule, the equation κ −1 = mc 2 d,⊥ holds and one then concludes that the value of the superfluid density can be determined from the ratio of the sound velocities propagating parallel and perpendicular to the x direction, according to Result (60) holds independently of the value of the Raman coupling Ω R , in the stripe as well as in the other quantum phases exhibited by spin-orbit-coupled BECs, provided that the symmetric intraspecies coupling case (g ↑↑ = g ↓↓ ) is considered. In Fig. 4 we report the value of the ratio (60) calculated either employing the ( Ω R /4E R ) 2 expansion using the results (44a)-(45a) for the sound velocities, and the full results employing the exact values for c 2 d,x and c 2 d,⊥ reported in Fig. 2. Result (60) is also well suited to determine experimentally the superfluid density of spin-orbit-coupled Bose gases through the direct measurements of the sound velocities. In the figure we also show the results for the superfluid density calculated in the plane-wave and zero-momentum phases, using the findings of [49] for the sound velocities and clearly revealing the first-order nature of the phase transition between the stripe and the plane-wave phase. Notice that in the plane-wave phase the sound velocity c 2 d,x should be actually replaced by the product c x,+ c x,− of the sound velocities along the +x and −x directions. The general behavior of the superfluid density as a function of the Raman coupling Ω R is qualitatively similar to the results derived in [51], where the phase twist method was employed using a variational Ansatz for the order parameter in the stripe phase; first studies based on Monte Carlo methods have shown that quantum fluctuations do not significantly affect the value of the superfluid density, which in three-dimensional dilute systems remains close to the mean-field prediction [54]. As pointed out in [62] the effects of the spin-orbit coupling have a dramatic consequence near the second-order transition between the plane-wave and zeromomentum phases, characterized by the divergence of the magnetic susceptibility χ M . In the same paper the theoretical predictions for the superfluid density were successfully compared with experiments in both the plane-wave and zero-momentum phases. The present approach generalizes the study of the superfluid density to the supersolid stripe phase, suggesting a direct procedure for its experimental measurement.
We finally compare the above results with those obtained from the Leggett criterion [64,65]. According to this criterion, the superfluid fraction of a system exhibiting density modulations along x, with periodicity π/k 1 , is bounded from above by the quantity with ρ = mn the local mass density (here assumed to be a function of the sole variable x). Leggett deduced this upper bound using a class of Ansatz many-body wave functions where all the particles in the superfluid have the same phase. The mean-field Ansatz belongs to this class, as it makes the stronger assumption that all the particles have identical wave functions. Consequently, within the mean-field approximation the upper bound ρ x s,L is always saturated, i.e., it coincides with the real superfluid density ρ x s . However, our system represents an exception to this rule because the single-particle Hamiltonian (1) lacks Galilean invariance, which is another requirement for the derivation of Leggett's criterion. Hence, Eqs. (60) and (61) can give different results in the case of spin-orbit-coupled Bose gases. This is obvious in the uniform phases, where Leggett's criterion predicts that the superfluid density coincides with the total density, ρ x s,L /ρ = 1, in stark contrast with the correct value reported in [62] and shown in Fig. 4. A similar effect takes place in the stripe phase, where ρ x s,L /ρ is not one but still significantly larger than the correct superfluid density ρ x s /ρ, see the blue dash-dot line in Fig. 4. This line actually corresponds to the exact numerical evaluation of the integral (61), but one can also compute it within the perturbative framework of the present work. The final results up to second order in ( Ω R /4E R ) 2 can be cast in the form where the contrast C is given by Eq. (18). This simple relation between Leggett's upper bound and the contrast of fringes is a common feature of shallow one-dimensional supersolids; it was deduced by general arguments in Ref. [66], and directly proven in the case of cnoidal waves [41]. Our perturbative approach enabled us to check its validity in spin-orbit-coupled Bose gases. Besides the spinodal behavior of the stripe phase we previously pointed out, here we see that the plane-wave phase can be found as a metastable state down to Ω R = 2.52 E R , where the roton gap in its excitation spectrum vanishes (see Ref. [49]). The interaction parameters are the same as Fig. 1.

Conclusion
We carried out a detailed analysis of a spin-orbit-coupled Bose-Einstein condensate in the stripe phase. By developing a suitable perturbation approach we obtained analytical expressions for a number of observables of interest. In particular, for systems at equilibrium we computed the periodicity and contrast of the stripes, as well as the chemical potential and energy per particle. Although our formulas are derived in the regime of small Raman coupling, they turn out to be quantitatively accurate also for significantly large values of the coupling. The extension of the perturbation method to the Bogoliubov modes enabled us to compute with similar accuracy the dispersion relation of the elementary excitations. We focused on the long-wavelength limit, where we proved that the two lowest-lying bands are gapless despite the presence of the Raman coupling term. We evaluated the velocities of sound waves propagating perpendicular and parallel to the stripes, as well as the corresponding Bogoliubov amplitudes, that reveal the Goldstone nature of these modes. In particular, while in the long-wavelength limit the density mode is dominated by a mere rotation of the global phase of the order parameter, the spin mode corresponds to a translation of the density fringes with respect to the equilibrium position. The simultaneous presence of a superfluid and a crystal mode is a typical signature of supersolidity. It also marks a striking difference with respect to binary mixtures without spin-orbit coupling, where the spin mode represents a simple rotation of the spin polarization.
We subsequently calculated the most relevant sum rules in the long-wavelength regime, yielding the compressibility, the magnetic susceptibility, and the superfluid density. The latter was directly related to the ratio of the sound velocities along and perpendicular to the direction of the spin-orbit coupling. This provides a viable way to measure the superfluid density of spin-orbit-coupled Bose-Einstein condensates using currently available experimental setups. Comparing the results of this method to those of the Leggett criterion we found a significant discrepancy, that highlights the lack of Galilean invariance in our system.

Acknowledgements
Many useful discussions and collaborations with Y. Li are acknowledged.

A Perturbation calculation of the ground-state wave function
In this appendix we delineate the general procedure for determining the ground-state wave function at arbitrary order in the Raman coupling (Sec. A.1). Subsequently, we explicitly compute the first-(Sec. A.2) and second-order (Sec. A.3) corrections.
A.1 General procedure and recurrence relations for the ground state As mentioned in Sec. 3.2, after inserting the expansions (15), (16), and (17) into the Gross-Pitaevskii equation (12) and equating terms of the same order on both sides, one deduces for n ≥ 1 the recurrence relation On the left-hand side of this equation one has the two matrices where σ ± = (σ x ± iσ y )/2, σ ↑,↓ = (I 2 ± σ z )/2, and I 2 is the 2 × 2 identity matrix. For simplicity, in the calculations of this Appendix and all the next ones we take χ 0 = 0 in the unperturbed wave function (14). On the right-hand side of Eq. (63) one has a source term J (n) that only depends on quantities determined up to order n − 1. It reads In the next two sections we shall consider the cases n = 1 and n = 2, where this involved formula reduces to much simpler expressions. Notice that, before solving Eq. (63) for a given n, one has first to determine µ (n) . This can be done by multiplying both sides of Eq. (63) by Ψ (0) † 0 and integrating with respect to the spatial coordinates. The resulting expression will contain terms depending on Ψ (n) 0 , that can be eliminated using the relation that follows from the normalization condition (5). Since Eq. (63) is linear, any of its solutions can be expressed as the sum of a particular solution and of the general integral of the associated homogeneous equation. This general integral reads Here the A's and B's constitute a set of 8 real arbitrary constants. In order to fix their values we first require that the wave function of the stripe phase be 2π-periodic in x at any order in Ω R /4E R . This implies that the coefficients of the nonperiodic terms in Eq. (67) have to vanish, i.e., A s,1 = 0. Regarding the periodic part of the general integral (67), corresponding to its last row, one can easily verify that the first term merely changes the global phase of the wave function, while the second term produces a shift of the offset of the density fringes. We recall that, because of the U(1) and translation symmetry of the model, if Ψ 0 (x, y, z) is a solution of Eq. (12), so is e iϕ 0 Ψ 0 (x − x 0 , y, z) for arbitrary ϕ 0 and x 0 . We choose the values of these two parameters by requiring the additional constraints to be fulfilled at all perturbation orders. This corresponds to taking B

A.2 First-order corrections to the ground state
Using the prescriptions of the previous section one can easily show that the first-order correction to the chemical potential is vanishing, µ (1) = 0. Besides, the source term (65) with n = 1 is simply The solution of Eq. (63) with n = 1 that satisfies the conditions (68) has the form Here the coefficients can be determined by inserting Eq. (70) into (63) and equating the terms on both sides with the same oscillating behavior. This yields

A.3 Second-order corrections to the ground state
With the results of Sec. A.2 at hand, the second-order correction to the chemical potential is easily evaluated and found equal to the second term on the right-hand side of Eq. (19) of the main text. The source term (65) with n = 2 takes the form where The solution of Eq. (63) with n = 2 has the structure The procedure for determining the coefficients is the same as the previous section, and yields These results enable us to determine the second-order correction to k 1 [from Eq. (11)], and to the energy per particle [from Eq. (6)], which are provided in Eq. (21) and (20) of the main text, respectively. We point out that the perturbative corrections given in Eqs. (70) and (74) feature harmonic components that do not correspond to minima of the single-particle dispersion, their appearance being uniquely caused by the interaction. In general, the n-th-order correction to the condensate wave function contains harmonic terms with wave vectors ranging from ±1 to ±(2n + 1) (in units of k 1 ). It is natural to expect that the perturbation series for Ψ 0 converges to the Bloch-wave structure (9). Although evidence for the presence of such higher-order harmonics was first found numerically in Ref. [36], the perturbation approach employed in the present work has the advantage of providing a better insight into their origin.

B Perturbation calculation of the Bogoliubov modes
The purpose of this appendix is to present the perturbation scheme for studying the Bogoliubov modes in the stripe phase. We start by introducing the band structure of the spectrum at vanishing Raman coupling (Sec. B.1) and deriving the recurrence relations (Sec. B.2). Then, we compute the perturbative corrections to the Bogoliubov frequencies and amplitudes at first (Sec. B.3) and second (Sec. B.4) order in the Raman coupling.
B.1 Band structure at Ω R = 0 As we mentioned at the end of Sec. 4.2, the Bogoliubov modes at Ω R = 0 can be labeled either by the physical momentum q or the quasimomentum k. Although the latter possibility may look unnatural, it is an obvious choice in view of the perturbation study of the Bogoliubov solutions at finite Raman coupling. One can define k at vanishing Ω R by rewriting the components of q as q x = k x +K , q y = k y , q z = k z .
Note that the x component of the quasimomentum, k x = 2{q x /2}, coincides with q x up to a reciprocal lattice vector,K = 2[q x /2]. Here we have used the standard notation {· · · } and [· · · ] for the fractional and integer part of a real number, respectively. On the one hand, these definitions ensure that k x is confined in the first Brillouin zone, i.e., 0 ≤ k x < 2 for any q x ; on the other hand, one has thatK is always an integer multiple of the size of the Brillouin zone, equal to 2 when expressed in units of k 1 . Let us insert Eq. (76) into Eqs. (38), (39), and (40), and label all quantities byK and k instead of q. In this new representation the Bogoliubov spectrum (38) exhibits an infinite number of branches, each labeled by two indices = d, s andK, and restricted to the first Brillouin zone. In Fig. 5 we plot some of the lowest-lying branches as functions of k x (we take k y = k z = 0). Notice that a band structure with two gapless modes is already visible in this figure, although there are no gaps separating the various bands because Ω R is zero.
As b,k . In performing the above procedure one has to take into account that the ω (0) ,K,k 's can cross one another (see Fig. 5), meaning that how to arrange them depends on k. For instance, for the two lowest bands shown in Fig. 5 one has where k c2 and k c3 are the crossing points marked in the figure. 2 We conclude this section by mentioning that the Bogoliubov amplitudes W as well as the completeness relation b,k We recall that, because of the symmetry properties of the operator B [see Eq. (29)], if W (0) b,k is solution of Eq. (27) at Ω R = 0 with frequency ω [59]. From Eq. (78a) we see that the two solutions have opposite norm. However, they actually correspond to the same physical oscillation of the system. Hence, without loss of generality we will restrict our perturbation analysis to the sole positive-norm modes. Nevertheless, the contribution of the negative-norm modes is crucial when writing the completeness relation (79). These properties will play a fundamental role in the formulation of the perturbation approach for the Bogoliubov modes.

B.2 General procedure and recurrence relations for the Bogoliubov modes
Let us now insert the expansions (42) and (43) into Eq. (27). To order n one obtains This recurrence relation is the starting point for our perturbative calculation of the Bogoliubov spectrum. As we will show in the next sections, an important building block of our formalism is represented by the integrals and analogous expressions for B b,k in terms of W ,K,k andW ,K,k according to the above prescriptions. One finds that only a limited subset of these integrals are nonzero; some of those that are needed for the calculations of the present work are reported in Appendix C. Notice that they vanish when k = k [for B The knowledge of the integrals (81) is crucial to compute the perturbative corrections to the Bogoliubov frequencies and amplitudes. For the latter we will use the completeness relation (79) to express W

B.3 First-order corrections to the Bogoliubov modes
Let us consider Eq. (80) with n = 1, The entries of B (1) are easily evaluated noting that h (1) SO = Ω R σ x /2, µ (1) = 0 (see Sec. A.2), and The first-order correction to the Bogoliubov frequency is vanishing. To prove this, we first multiply both sides of Eq. (83) by W (0) † b,k , and use the equality W ,k η to eliminate the terms containing W (1) b,k . Then, we integrate with respect to the spatial coordinates, and use the normalization condition (78a), yielding In order to calculate the first-order correction to the Bogoliubov amplitudes we first need to evaluate These results were obtained multiplying both sides of Eq. (83) by W (0) † b ,k , integrating over space, and using the orthonormalization conditions (78). Equation (86a) cannot be used when (b , k ) = (b, k); to address this case, we first notice that expanding the normalization condition (31) one finds at first order Besides, the imaginary part of this integral can be set equal to zero with an appropriate choice of the phase of W b,k . Hence, one eventually has V d 3 r W These relations, combined with formulas (100) and (101) in Appendix C, enable one to compute the first-order correction to the amplitude of any Bogoliubov mode.

B.4 Second-order corrections to the Bogoliubov modes
Equation (80) with n = 2 reads Here the operator B (2) can be evaluated recalling that 1 is given by the second term of Eq. (21), µ (2) by that of Eq. (19), and The second-order correction to the frequency can be calculated starting from Eq. (89) and performing the same kind of manipulations as those leading to Eq. (85). One eventually finds where we made use of Eqs. (85) and (88). This expression can be explicitly evaluated using the formulas in Appendix C. In general, for a given Bogoliubov mode the correction to the frequency (92) contains up to nine nonzero terms. In Fig. 6 we compare the numerically computed spectrum with the results of our perturbation approach up to second order. We consider excitations propagating along x and show the lowest four bands for a sufficiently small value of the Raman coupling. The agreement is excellent at almost any k x , except in the proximity of the modes where two bands cross at Ω R = 0 (black dots of Fig. 5), where Eq. (92) predicts an unphysical divergent behavior (see also the discussion at the end of Sec. 4.2). This effect becomes stronger at increasing Ω R .
For the computation of the two sound velocities we take the two lowest-lying bands (b = 1, 2) in the k → 0 limit; recall that at Ω R = 0 these bands correspond to = s, d and  44) and (45) of the main text, respectively. Notice that the scaling (33) of k requires an analogous transformation of the sound velocity, c → c /(k 2 1 cos 2 θ k +k 2 R sin 2 θ k ) 1/2 , so to leave ω ,0,k unaffected. In writing Eqs. (44) and (45) we inverted this transformation so to obtain an undistorted value of the c 's.
The second-order correction to the Bogoliubov amplitudes is given by Eq. (82) with n = 2 and the expansion coefficients The derivation of these formulas is analogous to that of Eqs. (86) for the first-order corrections. As in the latter case, when (b , k ) = (b, k) one has to replace Eq. (94a) by the relation which follows from the expansion of the normalization condition (31) up to second order [with the phase of W b,k again chosen such that the coefficient (95) be real].
In Sec. 4.3 we studied the low-k behavior of the amplitudes of the two gapless Bogoliubov bands up to first order in Ω R /4E R , see Eqs. (46). The general form of these amplitudes at arbitrary Raman coupling is (96b) Here we use dimensional coordinates and momenta. The coefficients α and β are real and depend on θ k , Ω R , and the interaction parameters. The structures (96) show that, because of the spin-orbit coupling, the density and spin modes do not exhibit a pure superfluid or crystal behavior, rather both characters are present even at low k. The only exception is for excitations propagating perpendicular to the x axis (θ k = π/2), for which β = 0 and the hybridization mechanism does not occur. At second order in Ω R /4E R the coefficients entering Eqs. (96) have the following expressions: α (θ k ) = α ,x cos 2 θ k + α ,⊥ sin 2 θ k , β (θ k ) = β ,x cos θ k , with α d,⊥ = 1, and The values of the α's and the β's as functions of Ω R are plotted in Fig. 7, where only excitations propagating parallel or perpendicular to x are considered. The above perturbative estimates are in excellent agreement with the exact numerical values up to significantly large values of the Raman coupling. In particular, in Fig. 7(c) it is shown that the property α d (θ k = π/2) = 1 holds at arbitrary Ω R , even beyond the perturbative regime.    bb,kk entering the second-order correction to the frequency (92): ,k ω s,K,k ΩK ,k . (102b)