Supercurrent-carrying supersolid in spin-orbit-coupled Bose-Einstein condensates

In the recent years, synthetic spin-orbit coupling in ultracold atoms has been recognized as a powerful tool to simulate many-body physics and explore exotic superfluids [1, 2, 3]. The spin-orbit coupling can be artificially introduced in neutral atoms by coupling two hyperfine states via a pair of Raman lasers, where the two states are treated as pseudospin [4, 5, 6]. The celebrated spin-orbit coupled Bose-Einstein condensates (BECs) possess unconventional and rich ground-state phase diagrams  [7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. Outstanding exotic ground states of such a system include the plane-wave phase and the stripe phase. The plane-wave phase breaks the time-reversal symmetry so that it is spin-polarized [17, 8]. The fundamental feature is that the plane-wave phase supports a phonon-maxon-roton structure in its elementary excitation spectrum [18, 19]. This phase provides a unique route to generate rotons [20, 21, 22]. The stripe phase simultaneously breaks the continuous translational symmetry and the gauge symmetry and accordingly possesses supersolidity [23, 24, 25, 26, 27]. In the stripe phase, the density crystalline structure offers an intriguing approach for experimental observation, while the elementary excitation spectrum has Bloch band-gap structures, which includes two gapless Nambu-Goldstone modes in the long wave-length regime. The theoretical discovery [24, 23, 25] and experimental observation [27] of the supersolid stripe phase constitute a main achievement in the field of spin-orbit-coupled BECs. Rich physics relevant to supersolid stripes have been revealed, including dynamics and excitations in these exotic states [28, 29, 30, 31, 32, 33]. On the other hand, atomic BECs with long-range interactions are active platforms to support roton excitations and to generate supersolid phases [34, 35, 36]. Superfluid properties of long-range-interaction-induced supersolids have attracted a great deal of research attention [37, 38, 39, 40]. The spin-orbit-coupling-induced supersolid stripes provide an alternative means to detect the exotic superfluidity.

The fundamental property of superfluids is that they can carry current flow without dissipation. The supercurrent-carrying ability is limited by two different instabilities [41]. Energetic instability happens when elementary excitations of a supercurrent-carrying superfluid have negative energies. The instability physically relates to the Laudau’s criterion of superfluidity. When the velocity of the superfluid is larger than a critical velocity, the superfluid becomes energetically unfavourable and then loses the supercurrent-carrying ability. Dynamical instability emerges when imaginary energies appear in elementary excitations [42]. The instability causes the exponential growth of excitation perturbations and therefore breaks superfluidity. Immediately after the first experimental realization of spin-orbit-coupled BECs [4], the existence of supercurrent-carrying plane-wave states are proposed [7, 20, 43]. Dynamical instability of these states is caused by their negative effective mass, and the energetic instability is shown to relate to the negative-energy phonon or negative-energy roton [20]. So far, all studied supersolid stripe phases in spin-orbit-coupled BECs in literature do not carry supercurrent. Whether there exist supercurrent-carrying supersolid stripes becomes an intriguing question. The main intuitive concern is that the density crystalline structure of stripes is not in favour of current flow. If they exist, the instabilities of their superfluidity naturally arise as a further question.

In this paper, we address the above two questions by constructing supercurrent-carrying supersolid stripes and analyzing their stabilities. The construction is inspired by Bloch waves in optical lattices [44, 45, 46, 47]. In a one-dimensional optical lattice, the Bloch wave function is $\exp(ikx)\phi(x)$, where $k$ is the quasimomentum, and $\phi(x)$ is a periodic function. It is well known that the Bloch waves carry current flow if $k\neq 0$. Therefore, the plane-wave prefactor $\exp(ikx)$ becomes the source of the current carried by the Bloch waves. We assume the wave function of supersolids to have a similar form as Bloch waves. The obvious difference is that the periods of the supersolid wave functions should be determined by minimizing the associated energy functional, while the period of the Bloch waves in optical lattice has a fixed value related to the wave length of the lattice beams. The supersolid family we found has the same density period and carries current flow that is proportional to the quasimomentum. The features of the supersolid family, including contrasts and spin polarizations, are identified. Energetic and dynamical instabilities are revealed, which can lead to the breakdown of their supercurrent-carrying ability. We calculate the relevant elementary excitations of the supersolid family, from which the energetic and dynamical instabilities are identified. The energetic instability relates to the well-known Landau’s criterion of superfluidity. The dynamical instability is qualitatively understood from the pure plane-wave phase. In spin-orbit-coupled BEC experiments, the detuning is a tunable and important parameter. We also study the supercurrent-carrying supersolid stripes in the presence of the detuning.

The outline of the paper is as follows. In Sec. II, we present our theoretical frame for investigating supercurrent-carrying supersolids and their instabilities. In Sec. III, we reveal the existence of the supercurrent-carrying supersolid family in the absence of detuning. Important features of the supersolid family are identified. The reasons of their energetic and dynamical instabilities are explained. The supercurrent-carrying supersolid stripes in the presence of detuning are considered in Sec. IV. In Sec. V, we discuss experimental accessibility of the novel supercurrent-carrying supersolids. Finally, the conclusions follow in Sec. VI.

We consider a spin-orbit-coupled spin-1/2 BEC. The spin-orbit coupling can be induced between two hyperfine states of the atoms via a two-photon transition process induced by a pair of Raman beams [4]. The system is described by the Gross-Pitaevskii (GP) equations,

$\displaystyle i\frac{\partial\psi}{\partial t}=\left(H_{\mathrm{SOC}}+H_{\mathrm{int}}\right)\psi,$

(1)

where $\psi=(\psi_{1},\psi_{2})^{T}$ is the two-component wave function. The GP equations are dimensionless. We choose the units of momentum, length, and energy as $\hbar k_{\mathrm{Ram}}$, $1/k_{\mathrm{Ram}}$, and $\hbar^{2}k^{2}_{\mathrm{Ram}}/m$, respectively. Here, $k_{\mathrm{Ram}}=2\pi/\lambda_{\mathrm{Ram}}$ is the wave vector of the Raman lasers with $\lambda_{\mathrm{Ram}}$ being the corresponding wave length, and $m$ is the atom mass. The spin-orbit-coupled single-particle Hamiltonian is

$$H_{\mathrm{SOC}}=-\frac{1}{2}\frac{\partial^{2}}{\partial x^{2}}-i\frac{\partial}{\partial x}\sigma_{z}+\frac{\Omega}{2}\sigma_{x}+\frac{\delta}{2}\sigma_{z},$$

(2)

with $\sigma_{x,z}$ being Pauli matrices. $\Omega$ is the Rabi frequency which depends on Raman laser intensities, and $\delta$ is the detuning of the two laser beams with respect to the energy difference between the two hyperfine levels. $H_{\mathrm{int}}$ denotes the contact interactions,

$\displaystyle H_{\mathrm{int}}[\psi]=\left(\begin{matrix}g|\psi_{1}|^{2}+g_{12}|\psi_{2}|^{2}&0\\ 0&g|\psi_{2}|^{2}+g_{12}|\psi_{1}|^{2}\end{matrix}\right).$

$g$ and $g_{12}$ are nonlinear coefficients for intra and intercomponent interactions respectively, which are proportional to the corresponding $s$-wave scatting lengths.

At the first stage, we search supercurrent-carrying supersolid solutions to the GP equations. The wave function of supersolids can be constructed as

	$\displaystyle\psi(x,t)$	$\displaystyle=\sqrt{n_{0}}e^{ikx-i\mu t}\left[\begin{array}[]{c}\phi_{1}(x)\\ \phi_{2}(x)\end{array}\right],$		(5)
	$\displaystyle\phi_{1,2}(x)$	$\displaystyle=\sum^{L}_{j=-L}e^{ij\xi x}\phi^{(j)}_{1,2},$		(6)

where $n_{0}$ denotes the mean atom density and $\mu$ is the chemical potential. The plane-wave prefactor $\exp(ikx)$ is introduced for the generation of current and $k$ is the quasimomentum. $\phi_{1,2}(x)$ are periodic functions and are expressed as a superposition of plane-wave modes $\exp(ij\xi x)$ with coefficients $\phi^{(j)}_{1,2}$. $L$ is the cutoff of the mode number and $\xi$ determines the period of supersolids. The coefficients $\phi^{(j)}_{1,2}$ satisfy the normalization condition,

$\displaystyle\sum^{L}_{j=-L}\left[\left|\phi^{(j)}_{1}\right|^{2}+\left|\phi^{(j)}_{2}\right|^{2}\right]=1.$

For a given quasimomentum $k$, the constructed wave functions are fully determined by minimizing the energy functional associated to the GP equations,

	$\displaystyle\mathcal{E}[\psi]$	$\displaystyle=\int dx\psi^{\dagger}H_{\mathrm{SOC}}\psi+\frac{g}{2}\int dx\left(|\psi_{1}|^{4}+|\psi_{2}|^{4}\right)$
		$\displaystyle\phantom{={}}+g_{12}\int dx|\psi_{1}|^{2}|\psi_{2}|^{2}.$		(7)

After the minimizing procedure, the unknown quantities, i.e., $\xi$ and $\phi^{(j)}_{1,2}$, can be obtained. The dependence of the resultant minimized energy $\mathcal{E}$ on the quasimomentum $k$, i.e., $\mathcal{E}(k)$, constitutes the dispersion relationship of supersolid stripe family. Then substituting the resultant wave function into the GP equations, we can calculate the corresponding chemical potential,

	$\displaystyle\mu$	$\displaystyle=\int dx\psi^{\dagger}H_{\mathrm{SOC}}\psi+g\int dx\left(|\psi_{1}|^{4}+|\psi_{2}|^{4}\right)$
		$\displaystyle\phantom{={}}+g_{12}\int dx|\psi_{1}|^{2}|\psi_{2}|^{2}.$		(8)

Whether the states carry supercurrent can be checked by calculating the current density $J$,

	$\displaystyle J$	$\displaystyle=\langle\psi|\hat{v}|\psi\rangle$
		$\displaystyle=n_{0}k+n_{0}\sum^{L}_{j=-L}\left[j\xi\left(|\phi^{(j)}_{1}|^{2}+|\phi^{(j)}_{2}|^{2}\right)\right.$
		$\displaystyle\phantom{={}}\left.+\left(|\phi^{(j)}_{1}|^{2}-|\phi^{(j)}_{2}|^{2}\right)\right],$		(9)

where $\hat{v}=-i\partial_{x}+\sigma_{z}$ is the velocity operator resulting from the the spin-orbit coupling [48, 49, 50, 51]. We can see that the current density is the product of the velocity and the density. The nonzero quasimomentum $k$ plays an important role in the generation of current flow.

The elementary excitation spectrum can be calculated once we know the wave function $\psi$ and associated chemical potential $\mu$. When the system is perturbed, the total wave functions can be written as the sum of the ground-state supersolid stripe $\psi$ and perturbations $\delta\psi$,

$$\Psi(x,t)=\sqrt{n_{0}}e^{ikx-i\mu t}\left\{\left[\begin{array}[]{c}\phi_{1}(x)\\ \phi_{2}(x)\end{array}\right]+\delta\psi(x,t)\right\},$$

(10)

with

$$\delta\psi(x,t)=\begin{bmatrix}u_{1}(x)e^{-i\omega t}+v^{\ast}_{1}(x)e^{i\omega^{*}t}\\ u_{2}(x)e^{-i\omega t}+v^{\ast}_{2}(x)e^{i\omega^{*}t}\end{bmatrix}.$$

(11)

Here, $\omega$ is the excitation energy, and $u_{1,2}(x)$ and $v_{1,2}(x)$ are perturbation amplitudes, satisfying the normalization condition,

$\displaystyle\sum_{j=1,2}\int dx\left[|u_{j}(x)|^{2}-|v_{j}(x)|^{2}\right]=1.$

By substituting the total wave functions into the GP equations [Eq. (1)] and keeping the linear terms with respect to perturbation amplitudes $u_{1,2}$ and $v_{1,2}$, we obtain the Bogoliubov–de Gennes (BdG) equation [52],

$\displaystyle\mathcal{M}\Phi=\omega\Phi,$

(12)

with

$\displaystyle\Phi=(u_{1},u_{2},v_{1},v_{2})^{T},~{}\mathcal{M}=\left(\begin{array}[]{cc}\mathcal{A}&\mathcal{B}\\ -\mathcal{B}^{\ast}&-\mathcal{A}^{*}\end{array}\right).$

(15)

The matrices $\mathcal{A}$ and $\mathcal{B}$ are

	$\displaystyle\mathcal{A}$	$\displaystyle=\frac{\Omega}{2}\sigma_{x}+\left(\begin{array}[]{cc}\mathcal{H}_{1}-i(\partial_{x}+ik)&n_{0}g_{12}\phi_{1}\phi^{\ast}_{2}\\ n_{0}g_{12}\phi^{\ast}_{1}\phi_{2}&\mathcal{H}_{2}+i(\partial_{x}+ik)\end{array}\right),$		(18)
	$\displaystyle\mathcal{B}$	$\displaystyle=n_{0}\left(\begin{array}[]{cc}g\phi^{2}_{1}&g_{12}\phi_{1}\phi_{2}\\ g_{12}\phi_{1}\phi_{2}&g\phi^{2}_{2}\end{array}\right),$		(21)

with

	$\displaystyle\mathcal{H}_{1}$	$\displaystyle=-\frac{1}{2}(\frac{\partial}{\partial x}+ik)^{2}+\frac{\delta}{2}+2n_{0}g|\phi_{1}|^{2}+n_{0}g_{12}|\phi_{2}|^{2}-\mu,$
	$\displaystyle\mathcal{H}_{2}$	$\displaystyle=-\frac{1}{2}(\frac{\partial}{\partial x}+ik)^{2}-\frac{\delta}{2}+2n_{0}g|\phi_{2}|^{2}+n_{0}g_{12}|\phi_{1}|^{2}-\mu.$

By diagonalizing the BdG equation, we can get the excitation energy $\omega$. Since the wave functions $\phi_{1,2}(x)$ have the period of $2\pi/\xi$, the matrices $\mathcal{M}$ depending on quadratic terms of $\phi_{1,2}$ are also periodic with the period $\pi/\xi$. Therefore, the excitation spectrum will have Bloch band-gap structures and perturbation amplitudes are in the form of Bloch waves. To proceed the diagonalization of the BdG equation, we set the perturbation amplitudes as Bloch waves and expand them by the plane-wave basis,

	$\displaystyle u_{1,2}(x)$	$\displaystyle=e^{iqx}\sum^{L}_{j=-L}U^{(j)}_{1,2}e^{i(2j+1)\xi x},$
	$\displaystyle v_{1,2}(x)$	$\displaystyle=e^{iqx}\sum^{L}_{j=-L}V^{(j)}_{1,2}e^{i(2j+1)\xi x},$		(22)

where $q$ is the perturbation quasimomentum, and $U^{(j)}_{1,2}$ and $V^{(j)}_{1,2}$ are expansion coefficients. The excitation spectrum $\omega(q)$ can be straightforwardly obtained by substituting Eqs. (22) into Eq. (12) and diagonalizing the resultant BdG equation. It was shown that the excitation spectrum has two gapless Nambu-Goldstone modes in the long wave-length regime ($q\rightarrow 0$) for a ground-state supersolid stripe phase and their appearance characterizes the superfluidity of the stripe phase [25]. Based on the elementary excitations, the properties such as quantum depletion [26] and structure factors [25, 53] can be analysed.

The matrix $\mathcal{M}$ is non-Hermitian, allowing for the existence of imaginary excitation energies. In the case of imaginary energies, the perturbations in Eq. (11) grow exponentially, so that the associated supersolid stripe is dynamically unstable. Such dynamical instability will destroy the associated supersolid stripe. The ground-state supersolid stripes are always dynamically stable. Recently, it is proposed that dynamical instability can be induced in a metastable supersolid stripe in spin-orbit-coupled BECs [54].

The energetic instability is associated with the excitation spectrum corresponding to $\tau\mathcal{M}$ [42], with

$\displaystyle\tau=\left(\begin{matrix}I&\\ &-I\end{matrix}\right),$

where $I$ is a $2\times 2$ identical matrix. The eigenvalue equation can then be written as

$\displaystyle\tau\mathcal{M}\Phi=\omega^{\prime}\Phi.$

(23)

The matrix $\tau\mathcal{M}$ is Hermitian so that it always has a real spectrum. The energetic instability is featured with negative values of the energy spectrum. The appearance of negative energy modes in $\omega^{\prime}$ means that the associated BEC state is energetically unfavourable. It was shown that the negative energies may happen around phonon and roton excitations for a supercurrent-carrying plane-wave state in spin-orbit-coupled BECs [20].

In spin-orbit-coupled BEC experiments, the detuning $\delta$ is an important tunable parameter. With or without the detuning, the symmetries of $H_{\mathrm{SOC}}$ changes dramatically, which will lead to a profound impact on relevant physics. We first study the existence and stability of supercurrent-carrying supersolids in the absence of the detuning. It is already known that the ground-state supersolid stripe phase which does not carry supercurrent only exists in the miscible interactions $g>g_{12}$ [9, 55, 27]. This motivates us to study the supercurrent-carrying analogs in this interaction regime with an intuitive expectation that they may share the same stability with the ground states if they exist. We choose $n_{0}g=0.5$, $n_{0}g_{12}=0.2$, and $\Omega=0.4$ in all the numerical calculations, and thus the supersolid stripe at $k=0$ is the ground state of the system.

By applying the trial wave functions in Eq. (6) and the minimization procedure illustrated in the previous section, we get the supercurrent-carrying supersolid states, whose dispersion relation $\mathcal{E}(k)$ is shown in Fig. 1(a). The dispersion of the stripe family has a parabolic-like shape. The supersolid stripe with $k=0$ has the lowest energy, corresponding to the ground state of the system, which is widely studied in literature. Its density crystalline structures are demonstrated in Fig. 1(b). The two components have the same periodic structure, i.e., $|\psi_{1}(x)|^{2}=|\psi_{2}(x)|^{2}$. Therefore, its spin polarization is $\langle\sigma_{z}\rangle=\int dx(|\psi_{1}|^{2}-|\psi_{2}|^{2})=0$. The period of crystalline densities is $\pi/\xi$ with $\xi=0.98$. Such value of $\xi$ is very close to the spin-orbit-coupling strength which is $1$ in our dimensionless units. The supersolid family demonstrated by the solid line in Fig. 1(a) has the same period with the $k=0$ supersolid. Two typical density profiles at $k=0.4$ and $0.7$ are depicted in Figs. 1(c) and 1(d), respectively. Compared with the $k=0$ supersolid, the family has three distinguished features.

(1) The supersolid family does carry current flow. The calculated current density from Eq. (II) for the associated supersolids with different $k$ is shown in Fig. 1(e). Only the $k=0$ supersolid has a zero current density, and all others are nonzero. The sign of the current density is same as the quasimomentum $k$. For $k>0$ states, the current flows in the positive direction. In contrast, for $k<0$ states, the current density is negative, which means that it flows towards the negative direction. The nonzero current witnesses that it is essentially important to introduce the plane-wave prefactor $\exp(ikx)$ in the construction of supersolid wave functions.

(2) The contrasts of crystalline densities are spin-dependent and vary obviously depending on the quasimomentum $k$. The contrast is defined as

$\displaystyle C=\frac{n_{\mathrm{max}}-n_{\mathrm{min}}}{n_{\mathrm{max}}+n_{\mathrm{min}}},$

(24)

where $n_{\mathrm{max}}$ ($n_{\mathrm{min}}$) is the maximum (minimum) of the density. From the density profiles in Figs. 1(c) and 1(d), we can see that two components still overlap with each other, while they have apparently different contrasts. The contrasts of the two components as a function of the quasimomentum $k$ are shown in Fig. 1(f). For $k>0$, the contrast of the first component decreases (the red-solid line), and that of the second component increases (the blue-dashed line) by raising $k$. However, the dependence of $C$ on $k$ changes oppositely for $k<0$. When $|k|$ reaches a critical value, the contrast of the one of two components completely disappears, and this component turns into a plane wave. A bit beyond the critical $|k|$, the supersolid family can not exist and the solutions become the plane-wave phase which is represented by the purple-dotted line in Fig. 1(a). A typical plane-wave density distribution is plotted in Fig. 1(h). The obvious feature of the plane-wave phase is its spin polarization, which is a result of the spontaneous breakdown of time-reversal symmetry.

(3) Most of the supercurrent-carrying supersolid states are spin-balanced. At first glance on Figs. 1(c) and 1(d), two components may be spin-imbalanced. However, the density oscillating centers of the two components are not coincident, resulting in the spin balance between the two components. The calculated spin polarization $\langle\sigma_{z}\rangle$ and $\langle\sigma_{x}\rangle$ are demonstrated in Fig. 1(g). It is shown that the two components are spin-balanced ($\langle\sigma_{z}\rangle\approx 0$) in the dominated region of the quasimomenta ($|k|<0.7$). When $|k|>0.7$, the supersolids change dramatically to adjust distributions including contrast and spin polarization for the preparation of transition to the plane-wave phase.

The three features are intrinsically correlated, which can be understood in the following way. We start from the GP equation and neglect the nonlinear interactions, $H_{\mathrm{SOC}}|\psi\rangle=\mathcal{E}|\psi\rangle$, with $|\psi\rangle$ being the supersolid stripe state and the associated energy $\mathcal{E}$. The stripe wave function is Bloch-wave-like, i.e., $|\psi\rangle=e^{ikx}|\phi\rangle$ with $|\phi\rangle$ being a periodic function. Then the GP equation becomes $H^{\prime}|\phi\rangle=\mathcal{E}|\phi\rangle$. Here, $H^{\prime}=e^{-ikx}H_{\mathrm{SOC}}e^{ikx}=-(\partial_{x}+ik)^{2}/2-i(\partial_{x}+ik)\sigma_{z}+\Omega\sigma_{x}/2$. By applying the Hellman-Feynman theorem, we get $\partial\mathcal{E}/\partial k=\langle\phi|\partial H^{\prime}/\partial k|\phi\rangle=k+\langle\phi|(-i\partial_{x})|\phi\rangle+\langle\phi|\sigma_{z}|\phi\rangle=J$ [56]. Meanwhile, the spin polarization of the supersolids is zero, $\langle\phi|\sigma_{z}|\phi\rangle=0$, as shown in Fig. 1(g), and $\langle\phi|(-i\partial_{x})|\phi\rangle\approx 0$. We then have $J\approx k$ and $\partial\mathcal{E}/\partial k\approx k$. Therefore, the current density linearly depends on the quasimomentum $k$ and its sign is relevant to the sign of the latter. This provides an insight for the understanding of the results in Fig. 1(e). $\partial\mathcal{E}/\partial k\approx k$ leads to $\mathcal{E}\propto k^{2}$ which explains the parabolic-like shape of dispersion relation shown in Fig. 1(a). On the other hand, when $k=0$, the Hamiltonian $H^{\prime}$ has the time-reversal symmetry, i.e., $\mathcal{T}H^{\prime}\mathcal{T}^{-1}=H^{\prime}$, with $\mathcal{T}=\mathcal{K}\sigma_{x}$ and $\mathcal{K}$ being the complex conjugate operator. The $k=0$ supersolid inherits the time-reversal symmetry, giving rise to $|\psi_{1}(x)|^{2}=|\psi_{2}(x)|^{2}$. Therefore the two components have the same contrast. However, the presence of nonzero quasimomentum breaks the time-reversal symmetry. Consequently, the two components in the supercurrent-carrying supersolids have different contrasts.

Given that the supercurrent-carrying supersolid stripes can exist, we study their elementary excitations by solving the BdG equations [Eq. (12)] with the purpose for examining their dynamical instabilities, and we also calculate Eq. (23) to examine energetic instabilities. The elementary excitations $\omega(q)$ for the supersolid states with $k=0,0.2,0.5,0.7$ are shown in Figs. 2(a1)-2(a4). As the supersolids are periodic in densities with the period of $\pi/\xi$, the elementary excitation spectra have Bloch band-gap structures with the Brillouin zone size $2\xi$. For the supersolid with $k=0$, the excitation spectrum is symmetric with respect to the Brillouin zone centered at $q=0$ [Fig. 2(a)]. The typical feature is that the lowest two bands have two gapless Nambu-Goldstone modes in the long wave-length region ($q\rightarrow 0$). The time-reversal symmetry is broken by the nonzero quasimomentum $k$. Consequently, the excitation spectrum losses the symmetry with respect to the Brillouin center. The excitation spectrum of the $k=0.2$ supersolid in Fig. 2(a2) shows such an asymmetry. The energy of two gapless modes softens in the $q<0$ region (i.e., around $q/\xi=2$ in the figure). There is no imaginary-energy excitation in both $k=0$ and $k=0.2$ supersolids, which means that these two supersolids are dynamically stable. For the supersolid with $k=0.5$, the softening of the gapless modes around $q/\xi=2$ leads to the lowest mode becoming energy-negative [see Fig. 2(a3)]. It was proved that the energetic instability can also be indicated from the existence of negative-energy modes in elementary excitation based on the BdG equation instead of calculating Eq. (23) [57]. The appearance of negative-energy modes shown in Fig. 2(a3) is a signal that the corresponding $k=0.5$ supersolid is energetically unstable. However, this supersolid is dynamically stable since there is no imaginary energy excitation. The excitation spectrum demonstrated in Fig. 2(a4) for the $k=0.7$ supersolid is very different from others in Figs. 2(a1)-2(a3). There are imaginary-energy excitations (represented by dashed lines in the figure), which indicates the $k=0.7$ supersolid dynamically unstable.

The imaginary-energy and negative-energy excitations emerge in the lowest two bands of the elementary excitation spectra, and they are confined in a finite quasimomentum $q$ region. We systematically identify the region of the quasimomentum $q$ where imaginary-energy and negative-energy excitations happen by calculating Eqs. (12) and (23), respectively. The results are shown in the ($k,q$) plane in Fig. 2(b). In this figure, the white region represents that there is no imaginary-energy and negative-energy excitation, there are negative-energy excitations in the light-green region, and in the dark-green region there are imaginary-energy excitations. The structures in Fig. 2(b) have an exact symmetry $(k,q)\rightarrow(-k,-q)$. This is because the BdG equations in Eq. (23) possess the symmetry $\mathcal{M}(k,q)=\mathcal{M}(-k,-q)$. Appearance of imaginary-energy (negative-energy) excitations means that the associated supersolids are dynamically (energetically) unstable. From Fig. 2(b), we can demarcate instability regions: the supersolids in $|k|<0.38$ are stable, in $0.38\leq|k|<0.58$ are energetically unstable, and in $0.58\leq|k|<0.78$ are dynamically and energetically unstable. We also incorporate the instability regions by different colored lines in Fig. 1(a).

The above identified energetic instability can be understood straightforwardly. The velocity of carried supercurrent in supersolids is proportional to the current density which is approximately $J\approx n_{0}k$. As the quasimomentum $|k|$ increasing from zero, the supercurrent velocity also increases. When $|k|$ reaches a critical value, the supercurrent-carrying ability is broken according to Landau’s criterion of superfluidity, which leads to energetic instability. Therefore, the supercurrent-carrying supersolids are energetically stable only within a certain quasimomentum $k$ region around $k=0$.

The dynamical instability of the supercurrent-carrying supersolids identified above is of interest. We provide an insightful picture to qualitatively understand its emergence. It is noticed that the supersolids are composed of plane waves $\exp[i(k+j\xi)x]$ with $j$ being an integer [see Eq. (6)]. This motivates us to study the pure plane-wave solutions firstly. The wave function of plane-wave phase is [20]

$$\psi=\sqrt{n_{0}}e^{ik^{\prime}x-i\mu^{\prime}t}\begin{pmatrix}\varphi_{1}\\ \varphi_{2}\end{pmatrix}.$$

(25)

Here the spinor $(\varphi_{1},\varphi_{2})^{T}$ is spatially independent and satisfies $|\varphi_{1}|^{2}+|\varphi_{2}|^{2}=1$. The spinor together with the chemical potential $\mu^{\prime}$ can be determined by minimizing the associated energy functional. The resultant energy $\mathcal{E}^{\prime}(k^{\prime})$ of the plane-wave phase is shown in Fig. 3(a), which has a similar structure to the single-particle spin-orbit-coupled energy. The analysis of dynamical instability of the plane-wave phase can be performed by solving the associated BdG equations similar to Eq. (12). The result is that the plane waves in the region of $|k^{\prime}|<0.51$ are dynamically unstable [represented by the red-dashed line in Fig. 3(a)], and others are dynamically stable. Such dynamical instability of the plane-wave phase was revealed to relate with the negative effective mass [20]. We interpret the dynamical instability of the supercurrent-carrying supersolids from the instability of the plane-wave phase. In the wave function of the supersolids $\sum_{j}\exp[i(k+j\xi)x](\phi_{1}^{(j)},\phi_{2}^{(j)})^{T}$, the $j=\pm 1$ plane waves have a dominant occupation, which can be seen from the numerically calculated wave function. When $k\geq 0.47$, the $j=-1$ plane wave $\exp[i(k-\xi)x]$ enters into the dynamically unstable region of the pure plane-wave phase (since $\xi=0.98$) [see the green dots in Fig. 3(a)], the consequence of which is that the whole wave function of the supersolids may become dynamically unstable. In the same way, when $k\leq-0.47$, the $j=1$ plane wave $\exp[i(k+\xi)x]$ accesses to the dynamically unstable region of the pure plane-wave phase. Therefore, the supersolids are expected to be dynamically unstable when $|k|\geq 0.47$. This predicted region is different from the previous numerical result ($|k|\geq 0.58$). This indicates that the interpretation of the dynamical instability of the supersolids by that of the pure plane-wave phase is not exact. After all, the supersolids are an assembly of many plane waves. Nevertheless, it provides a profound insight to understand the dynamical instability of the supercurrent-carrying supersolids.

Finally, we emphasize that it is the plane-wave factor $\exp(ikx)$ that induces supercurrent and instability in supersolid stripes. The stripe at $k=0$ is the ground state of the system, which does not carry current and is always dynamically and energetically stable. The energetic instability indicates that the associated supersolid stripe is energetically unfavourable, and the system actually prefers the plane-wave phase [3]. The dynamical instability breaks the associated supersolid stripe to give rise to complicate and random density patterns. The fate of such an unstable stripe is similar to that of a dynamically unstable Bloch wave in optical lattices [58, 59], where the loss of atoms and random density patterns have been experimentally observed [60, 61]. Meanwhile, the dynamical instability of the Bloch wave is used as a signal for the existence of solitonic solutions [62, 63].

In the presence of the detuning term $(\delta/2)\sigma_{z}$, $H_{\mathrm{SOC}}$ does not have the time-reversal symmetry, and the system prefers spin-polarized solutions with nonzero $\langle\sigma_{z}\rangle$. The experimentally tunable detuning brings novel physics in spin-orbit-coupled BECs [10, 11]. We will show that supercurrent-carrying supersolid stripes can still exist in the presence of the detuning. Using the same constructed wave functions in Eq. (6) and following the procedure of minimizing the associated energy functional, we obtain the supersolid stripes in the parameter regime of $\delta=0.1$, $\Omega=0.4$, $n_{0}g=0.5$, and $n_{0}g_{12}=0.2$.

The dispersion relation $\mathcal{E}(k)$ is demonstrated in Fig. 4(a). The whole family has the same density period $\pi/\xi$ with $\xi\approx 0.98$, which is nearly the same as the zero-detuning case in the previous section. The density distributions for the $k=0$ and $k=0.5$ supersolids are described in Figs. 4(b) and 4(c), respectively. The outstanding feature is that the two components are spin-imbalanced. The second component has a larger occupation when $\delta>0$. The family carries supercurrent and the current density is shown in Fig. 4(d). It is obvious that the current density is still linearly dependent on the quasimomentum, $J\propto k$. The contrast of each component as a function of the quasimomentum is demonstrated in Fig. 4(e). In the positive quasimomentum region, the contrast of the first component decreases with $k$ increasing, and it completely disappears after $k$ reaching a critical value, indicating that the first component turns to be a plane wave. Beyond the critical quasimomentum, the supersolids can not exist. Similarly, in the negative quasimomentum region, the contrast of the second component gradually disappears. Therefore, the supersolids can only exist in a finite quasimomentum region, which is $-0.72\leqslant k\leqslant 0.85$. Such an asymmetric existence region with respect to $k=0$ becomes a main difference compared with the zero-detuning case shown in Fig. 1(a). In the main part of the existence region, the spin polarizations $\langle\sigma_{z}\rangle$ and $\langle\sigma_{x}\rangle$ are almost constant and negative [see Fig. 4(f)]. Around the existence region boundaries, the spin polarizations change dramatically.

The breakdown of supercurrent-carrying ability may be induced by energetic and dynamical instabilities. In order to analyze the instability, we calculate elementary excitations of supercurrent-carrying supersolid stripes. The elementary excitation spectrum $\omega(q)$ for the supersolids with $k=0,0.2,0.5,0.7$ are shown in Figs. 5(a1)-(a4), respectively. All the spectra show the asymmetry with respect to $q=0$ due to the lack of the time-reversal symmetry, which is a result of the nonzero detuning. The structures of the spectra are similar to these in the zero-detuning cases shown in Figs. 2(a1)-2(a4). The negative-energy excitations [Fig. 5(a3)] and the imaginary-energy excitations [Fig. 5(a4)] happen in the lowest two bands. From the elementary excitation spectrum, we capture instability regions and show them in the $(k,q)$ plane in Fig. 5(b). The instability regions (shadowed areas in the figure) lose the symmetry of $(k,q)\rightarrow(-k,-q)$ since in the BdG equations $\mathcal{M}(k,q)\neq\mathcal{M}(-k,-q)$ in presence of the detuning. From boundaries of light-green regions (for energetic instability) and dark-green regions (for dynamical instability), we conclude that in the region $-0.33<k<0.35$ [shown by the green line in Fig. 1(a)] the supersolids are stable, in $-0.52\leq k\leq-0.33$ and $0.35\leq k\leq 0.58$ [the blue line in Fig. 1(a)] they are energetically unstable, and in $k>0.58$ and $k<-0.52$ [the red line in Fig. 1(a)] they are both dynamically and energetically unstable.

The onsets of the dynamical instability are at $k=0.58$ and $k=-0.52$, which are not symmetric with respect to $k=0$. The occurring of the dynamical instability as well as its asymmetric parameter region can be explained from the pure plane-wave phase. The energy of the plane-wave phase shown in Fig. 3(b) still has a double-well structure. However, in comparison with the zero-detuning case shown in Fig. 3(a), the nonzero detuning biases the double wells, so that the dispersion loses the symmetry with respect to $k^{\prime}=0$. Furthermore, the region of the dynamical instability in the plane-wave phase also loses the symmetry with respect to $k^{\prime}=0$, and it is $-0.56<k^{\prime}<0.46$ [represented by the red-dashed line in Fig. 3(b)]. Even in the presence of the detuning, the $j=\pm 1$ plane waves in the wave function of supercurrent-carrying supersolids $\sum_{j}\exp[i(k+j\xi)x](\phi_{1}^{(j)},\phi_{2}^{(j)})^{T}$ still have a dominant occupation. Considering $\xi\approx 0.98$ in the supersolid family, when $k>0.42$, the $j=-1$ plane wave $\exp[i(k-\xi)x]$ enters into the dynamically unstable region of the pure plane-wave phase [see the green dots in Fig. 3(b)]. Similarly, when $k<-0.52$, the $j=1$ plane wave $\exp[i(k+\xi)x]$ enters into the dynamically unstable region of the pure plane-wave phase. Therefore, intuitively following the dynamical instability of the pure plane-wave phase, we predict the supersolid family becomes dynamically unstable in the region of $k>0.42$ and $k<-0.52$. The predicted onsets are at $k=0.42$ and $k=-0.52$, which are not symmetric with respect to $k=0$. The reason of such asymmetry is that the dispersion relation of the pure plane-wave phase does not have the symmetry with respect to $k^{\prime}=0$ as shown in Fig. 3(b) due to the detuning. It is noticed that the predicted onsets are not exact to the numerically calculated values.

We have uncovered the existence of supercurrent-carrying supersolid stripes. Now we discuss their experimental accessibility. Experiments may start from a two-component BEC with miscible interactions in an elongated trap. Dressing the two-component BEC along the longitudinal direction by adiabatically ramping up two Raman lasers can prepare a supersolid stripe [27]. Such an implementation of supersolid stripes corresponds to the ground state, i.e., the $k=0$ supersolid, which does not carry supercurrent. Then the $k=0$ supersolid may be accelerated by applying an external force $F$ along the longitudinal direction, such as the gravity-induced force [64, 65]. The quasimomentum of the supersolid will increase linearly with the time $t$, i.e., $k=Ft/\hbar$ [66]. By controlling the duration of the applied force, we can prepare the supersolid with a desired quasimomentum. Once the supercurrent-carrying supersolids are prepared, their dynamical instability can be observed by holding the supersolids for a certain time and measuring the loss of condensed atoms [12, 60]. In the spin-orbit-coupled supersolid stripe experiment in Ref. [27], the excitations in the transverse directions are completely irrelevant to the observed results. We expect that the transverse degrees of freedom may not qualitatively change the predicted instability results of the supercurrent-carrying supersolids, which are reminiscent of Bloch states in a longitudinal optical lattice. It has been shown that the transverse degrees of freedom do not qualitatively modify the instabilities of the Bloch states existing along the longitudinal direction [67].

By constructing Bloch-wave-like supersolid wave functions, we revealed the existence of supercurrent-carrying supersolid stripes in spin-orbit-coupled BECs. The family of supersolid stripes has the same density period and possesses a parabolic-like dispersion. The current carried by these supersolids is proportional to the quasimomentum. Energetic and dynamical instabilities emerge with a large current, which limit the supercurrent-carrying ability of the supersolids and even destroy them. The energetic instability relates to the well-known Landau’s criterion of superfluidity. The dynamical instability is interpreted as that the supersolid stripes involve a dynamically unstable plane wave. These results provide a possible route to exploring exotic supersolids carrying supercurrents and their instabilities.

We are grateful to Prof. Tomoki Ozawa for providing us the insight into the explanation of the dynamical instability of the supercurrent-carrying supersolids from that of the pure plane-wave phase. We also acknowledge valuable discussions with Prof. Xi-Wang Luo. This work is supported by National Natural Science Foundation of China with Grants No. 12374247 and No. 11974235 and Shanghai Municipal Science and Technology Major Project (Grant No. 2019SHZDZX01-ZX04). L.H. also acknowledges support from Okinawa Institute of Science and Technology Graduate University. Q.Z. is supported by NKRDPC (Grant No. 2022YFA1405304), NSFC (Grant No. 12004118), and Guangdong Basic and Applied Basic Research Foundation (Grants No. 2020A1515110228 and No. 2021A1515010212).