# Superfluidity of a Raman spin-orbit-coupled Bose gas at finite temperature

Xiao-Long Chen<EMAIL_ADDRESS>Department of Physics, Zhejiang Sci-
Tech University, Hangzhou 310018, China  Institute for Advanced Study,
Tsinghua University, Beijing 100084, China  Centre for Quantum Technology
Theory, Swinburne University of Technology, Melbourne 3122, Australia  Xia-Ji
Liu  Centre for Quantum Technology Theory, Swinburne University of Technology,
Melbourne 3122, Australia  Hui Hu  Centre for Quantum Technology Theory,
Swinburne University of Technology, Melbourne 3122, Australia

###### Abstract

We investigate the superfluidity of a three-dimensional weakly interacting
Bose gas with a one-dimensional Raman-type spin-orbit coupling at both zero
and finite temperatures. Using the imaginary-time Green’s function within the
Bogoliubov approximation, we explicitly derive analytic expressions of the
current-current response functions in the plane-wave and zero-momentum phases,
from which we extract the superfluid density in the limits of long wavelength
and zero frequency. At zero temperature, we check that the resultant
superfluid density agrees exactly with our previous analytic prediction
obtained from a phase-twist approach. Both results also satisfy a generalized
Josephson relation in the presence of spin-orbit coupling. At finite
temperature, we find a significant non-monotonic temperature dependence of
superfluid density near the transition from the plane-wave phase to the zero-
momentum phase. We show that this non-trivial behavior might be understood
from the sound velocity, which has a similar temperature dependence. The non-
monotonic temperature dependence is also shared by Landau critical velocity,
above which the spin-orbit-coupled Bose gas loses its superfluidity. Our
results would be useful for further theoretical and experimental studies of
superfluidity in exotic spin-orbit coupled quantum gases.

##  I Introduction

Spin-orbit coupling (SOC), linking a particle’s spin to its motion, can lead
to ubiquitous quantum effects in various areas of physics [1, 2]. It plays a
central role in the emergent exotic bosonic and fermionic states of matter [3,
4, 5, 6, 7, 8, 9], such as topological insulators featuring quantum spin Hall
effect [10, 11, 12], and topological Fermi superfluids and Majorana fermions
[13, 14, 15, 16, *liu2013topological, 18]. Owing to its versatility [19],
ultracold atomic gases have became a powerful platform to create synthetic
gauge fields and to simulate SOC [20, 21, 22, 23]. In the last decade, by
utilizing Raman laser beams to carefully modulate atom-light interactions,
experimentalists have successfully realized one-dimensional Raman-laser-
induced SOC [24, 25, 26] and two-dimensional SOC [27, 28, 29, 30] in both
ultracold Bose and Fermi gases. These crucial achievements have driven
tremendous investigations on the nontrivial role of SOC in quantum many-body
physics [31, 32, 33, 34, 35, 36, 37, 38, 39].

Superfluidity is a prominent phenomenon in quantum liquids, supporting
frictionless flow in a capillary without dissipating energy [40, 41]. This
property is closely related to Bose-Einstein condensation and has been
qualitatively understood by using the celebrated Landau criterion. For an
isotropic system that is invariant under Galilean transformation, a critical
velocity $v_{c}\equiv\mathrm{min}\left[\omega({\bf p})/|{\bf p}|\right]$ can
be formally defined with the elementary excitation spectrum $\omega({\bf p})$
of the superfluid. Below the critical velocity, an impurity’s motion can not
create excitations to induce dissipation in energy and hence to destroy the
superfluidity of the system. In conventional weakly-interacting Bose gases,
this critical velocity is dominated by the linear phonon mode $\omega({\bf
p})=c_{s}{\bf p}$ of the lowest-lying excitation spectrum in the long-
wavelength limit, so it coincides with the sound velocity $c_{s}$.

SOC is known to violates the Galilean invariance [42, 41]. Therefore, it
becomes an intriguing question to characterize the superfluidity of a weakly-
interacting Bose gas in the presence of, e.g., the Raman-induced SOC. Previous
works have shown that, this Raman-type SOC can not only lead to various
ground-state phases by tuning the Raman coupling (i.e., the Rabi frequency)
[31, 32, 33], i.e., the stripe (ST), plane-wave (PW), and zero-momentum (ZM)
phases, but also significantly changes the excitation spectrum with the
appearance of a roton-maxon structure in the PW phase and a quadratic
dispersion at the PW-ZM transition [32, 33, 34]. As a result, the critical
velocity is determined by the roton minimum instead of the phonon mode,
leading to anisotropic superfluidity. The critical velocity can even vanish at
the critical PW-ZM transition point [33, 43]. These novel properties in the
elementary excitation spectrum as well as the sound velocity have already been
demonstrated in cold-atom laboratories [36, 37, 38].

However, there are only a few works on the calculation of superfluid density,
a key quantity that directly characterizes the superfluidity of the system.
The superfluid density of a SOC Bose gas was first predicted for the PW and ZM
phases in Ref. [44] at _zero_ temperature by Zhang et al., using a sum-rule
analysis of current-current correlation functions. The study shows that the
gapped branch in the elementary excitation spectrum plays a crucial role in
inducing a nonzero normal density. Most recently, by employing a first-order
stripe Ansatz and a plane-wave Ansatz within a phase-twist method, the present
authors and collaborators derived analytic expressions of zero-temperature
superfluid density as a function of Rabi frequency in different phases [35],
consistent with the findings in Ref. [44]. In addition, the dependence of
superfluid density on the high-order harmonics has also been investigated
numerically in the exotic stripe phase. The high-order harmonics shift the
critical Rabi frequency of the ST-PW phase transition point and suppress the
superfluid density, compared to the results obtained by using a first-order
stripe Ansatz. Recently, by means of a perturbative approach at small Rabi
frequency, Martone and Stringari also addressed the superfluid density in the
stripe phase at zero temperature within the sum rule approach [45].

In this work, we aim to investigate the effect of thermal fluctuations on the
superfluidity of a three-dimensional weakly interacting Bose gas with Raman-
type SOC, by conducting a microscopic calculation of superfluid density using
current-current correlation functions at both zero temperature and _finite_
temperature. This direct calculation is crucial, since the previous sum-rule
approach [44] is no longer applicable at finite temperature. In more detail,
by means of the Gross-Pitaevskii theory within the Bogoliubov approximation,
we first calculate the condensate wave function from a plane-wave Ansatz and
the Bogoliubov quasiparticle wave functions with respect to the tunable Rabi
frequency. The current-current response function is then derived analytically
using the imaginary-time Green’s function formalism. The total and normal
density can be extracted from the longitudinal and transverse response
functions in the long-wavelength and low-frequency limits, respectively.

To check our direct calculation of the superfluid density, at zero temperature
we compare the obtained superfluid density fraction with the previous analytic
expression derived from a phase-twist method [35]. We also consider an
alternative examination by using the so-called Josephson relation, which
relates the superfluid density to condensate density. Josephson relation has
been studied in detail in conventional spinless weakly interacting Bose gases
[46, 47, 48, 49], and has also been considered recently in the multi-component
case [50]. Here, we generalize the Josephson relation to the case of a SOC
Bose gas, from which we determine the zero-temperature superfluid density
directly via the single-particle Green’s function and the condensate density.

At finite temperature, we find a non-monotonic temperature dependence of the
extracted superfluid density near the transition from the plane-wave phase to
the zero-momentum phase. We discuss the causes of this non-trivial behavior in
superfluid density, by considering the sound velocity and Landau critical
velocity at finite temperature.

The rest of the paper is organized as follows. The model Hamiltonian and the
theoretical framework are introduced in Sec. II. In Sec. III, we derive the
analytic expression of the current-current response functions, and show how to
calculate the superfluid density via the longitudinal and transverse
components of the response functions (see Fig. 1). We then apply the
generalized Josephson relation to spin-orbit coupled Bose gases to obtain the
superfluid density. In the plane-wave and zero-momentum phases, the
predictions of superfluid density at zero temperature from two approaches are
compared with our previous analytic result obtained from a phase-twist method
(see Fig. 2). We further calculate the current-current response functions at
finite temperature and study the temperature dependence of superfluid density.
To understand the non-trivial temperature dependence, we also show the sound
velocity and Landau critical velocity as a function of temperature (see Figs.
3 and 4). A summary and outlook are given in Sec. IV. Appendix A is devoted to
the technical details of the current-current response functions.

##  II Theoretical Frameworks

###  II.1 The model Hamiltonian

We consider a three-dimensional weakly interacting spin-1/2 Bose gas with a
one-dimensional Raman-induced spin-orbit coupling, the same as in our previous
works [34, 35]. The system can be described by the model Hamiltonian,
$\hat{H}=\hat{H}_{0}+\hat{H}_{\mathrm{int}}$, where the single-particle
Hamiltonian $\hat{H}_{0}$ and the interaction Hamiltonian
$\hat{H}_{\mathrm{int}}$ read, respectively ($\hbar=1$ and the superscript “*”
represents the Hermitian conjugate) [24, 5, 33]

$\displaystyle\hat{H}_{0}$ $\displaystyle=$ $\displaystyle\int d^{3}{\bf
r}\begin{pmatrix}\hat{\Psi}_{\uparrow}^{*}({\bf
r}),&\hat{\Psi}_{\downarrow}^{*}({\bf
r})\end{pmatrix}\mathcal{H}_{\mathrm{s}}(\hat{{\bf
p}})\begin{pmatrix}\hat{\Psi}_{\uparrow}({\bf r})\\\
\hat{\Psi}_{\downarrow}({\bf r})\end{pmatrix},$ (1a)
$\displaystyle\hat{H}_{\mathrm{int}}$ $\displaystyle=$ $\displaystyle\int
d^{3}{\bf
r}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}\frac{g_{\sigma\sigma^{\prime}}}{2}\hat{\Psi}_{\sigma}^{*}\hat{\Psi}_{\sigma^{\prime}}^{*}\hat{\Psi}_{\sigma^{\prime}}\hat{\Psi}_{\sigma}({\bf
r}).$ (1b)

Here, the single-particle part $\mathcal{H}_{\mathrm{s}}(\hat{{\bf p}})$ is
given by

$\mathcal{H}_{\mathrm{s}}(\hat{{\bf p}})=\frac{(\hat{{\bf
p}}-k_{\mathrm{r}}\hat{{\bf
e}}_{x}\sigma_{z})^{2}}{2m}+\frac{\Omega}{2}\sigma_{x}+\frac{\delta}{2}\sigma_{z},$
(2)

with the canonical momentum operator $\hat{{\bf p}}=-i\nabla$ and Pauli
matrices $\sigma_{x,z}$. $k_{\mathrm{r}}\hat{{\bf e}}_{x}$ is the recoil
momentum of the Raman lasers along the $x$ axis, with a recoil energy
$E_{\mathrm{r}}=k_{\mathrm{r}}^{2}/(2m)$. It is straightforward to see that
the momentum operator $\hat{{p}}_{x}$ is coupled to the spin via the one-
dimensional physical momentum term
$(\hat{{p}}_{x}-k_{\mathrm{r}}\sigma_{z})^{2}$. For simplicity, the detuning
of the Raman lasers is assumed to be zero $\delta=0$. The Rabi frequency
$\Omega$ can be flexibly tuned, in accord with the recent experiments [39,
38]. $g_{\sigma\sigma^{\prime}}=4\pi a_{\sigma\sigma^{\prime}}/m$ are
interaction strengths for intra- ($\sigma=\sigma^{\prime}$) and inter-species
($\sigma\neq\sigma^{\prime}$), and $a_{\sigma\sigma^{\prime}}$ are the
corresponding $s$-wave scattering lengths.

###  II.2 A self-consistent approach at finite temperature

In this work, we employ the quasiparticle formalism at the Bogoliubov level to
describe a weakly interacting dilute Bose gas with SOC at zero and finite
temperature [51, 52, 53, 54, 41]. Following the standard procedure as in our
previous works [55, 34, 35], the Bose field operator $\hat{\Psi}_{\sigma}({\bf
r},t)$ for spin component $\sigma=(\uparrow,\downarrow)$ can be rewritten as a
combination of the condensate wave function $\psi_{\sigma}$ and the
noncondensate fluctuation operator $\hat{\eta}_{\sigma}$ as

$\hat{\Psi}_{\sigma}({\bf r},t)=\psi_{\sigma}({\bf
r})+\hat{\eta}_{\sigma}({\bf r},t).$ (3)

Using a Bogoliubov transformation, the fluctuation operator
$\hat{\eta}_{\sigma}({\bf r},t)$ and its conjugate can be expanded as

$\displaystyle\hat{\eta}_{\sigma}$ $\displaystyle=$
$\displaystyle\underset{j}{\sum}\left[u_{j\sigma}({\bf
r})e^{-i\omega_{j}t}\hat{a}_{j}+v_{j\sigma}^{*}({\bf
r})e^{i\omega_{j}t}\hat{a}_{j}^{\dagger}\right],$ (4a)
$\displaystyle\hat{\eta}_{\sigma}^{*}$ $\displaystyle=$
$\displaystyle\underset{j}{\sum}\left[u_{j\sigma}^{*}({\bf
r})e^{i\omega_{j}t}\hat{a}_{j}^{\dagger}+v_{j\sigma}({\bf
r})e^{-i\omega_{j}t}\hat{a}_{j}\right],$ (4b)

in terms of the quasiparticle wave functions $u(u^{*})$, $v(v^{*})$ and the
quasiparticle frequency $\omega_{j}$. In free space, the Bogoliubov wave
functions $u$, $v$ can be expanded as $u_{j\sigma}({\bf
r})=u^{(\lambda)}_{{\bf q}\sigma}e^{i{\bf qr}}/\sqrt{V}$ and $v_{j\sigma}({\bf
r})=v^{(\lambda)}_{{\bf q}\sigma}e^{i{\bf qr}}/\sqrt{V}$, where $V$ is the
system volume. Here, $j\equiv({\bf q},\lambda)$ is the index of the
quasiparticle energy level, with the momentum ${\bf q}$ and the branch index
$\lambda$. $\hat{a}^{\dagger}$ ($\hat{a}$) are respectively the creation
(annihilation) operators for quasiparticles, satisfying the bosonic
commutation relations $[\hat{a}_{i},\hat{a}^{\dagger}_{j}]=\delta_{ij}$,
$[\hat{a}_{i}^{\dagger},\hat{a}_{j}^{\dagger}]=[\hat{a}_{i},\hat{a}_{j}]=0$.
After substituting Eqs. (3) and (4) into the equations of motion

$i\partial_{t}\hat{\Psi}_{\sigma}({\bf
r},t)=\begin{bmatrix}\hat{\Psi}_{\sigma},\hat{H}\end{bmatrix},$ (5)

and applying the mean-field decoupling of the cubic terms in $\hat{\Psi}$ and
$\hat{\Psi}^{\dagger}$ [51], we obtain two coupled equations as in Refs. [34,
35].

The first equation is the modified Gross-Pitaevskii (GP) equation for the
condensate,

$\left[\mathcal{H}_{s}(\hat{{\bf
p}})+g_{{}_{\uparrow\downarrow}}n_{\mathrm{sf}}\sigma_{x}+\mathrm{diag}(\mathcal{L}_{\uparrow},\mathcal{L}_{\downarrow})\right]\psi=\mu\psi$
(6)

with the spinor
$\psi\equiv\left(\psi_{\uparrow},\psi_{\downarrow}\right)^{T}$, chemical
potential $\mu$, the spin-flip density term $n_{\mathrm{sf}}$ and the diagonal
element $\mathcal{L}_{\sigma}\equiv
g(n_{c\sigma}+2n_{t\sigma})+g_{{}_{\uparrow\downarrow}}n_{\bar{\sigma}}$ (here
$\bar{\sigma}\neq\sigma$). The second equation is the coupled Bogoliubov
equation for quasiparticles,

$\displaystyle\left[\mathcal{H}_{s}(\hat{\bf
p})-\mu+\mathcal{A}_{\uparrow}\right]U_{j}+\mathcal{B}V_{j}$ $\displaystyle=$
$\displaystyle\omega_{j}U_{j},$ (7a)
$\displaystyle-\mathcal{B}U^{*}_{j}-\left[\mathcal{H}_{s}(\hat{\bf
p})-\mu+\mathcal{A}^{*}_{\downarrow}\right]V^{*}_{j}$ $\displaystyle=$
$\displaystyle\omega_{j}V^{*}_{j},$ (7b)

where $U_{j}\equiv\left(u_{j\uparrow},u_{j\downarrow}\right)^{T}$,
$V_{j}\equiv\left(v_{j\uparrow},v_{j\downarrow}\right)^{T}$, and

$\mathcal{A}_{\sigma}\equiv\begin{pmatrix}2gn_{\uparrow}+g_{{}_{\uparrow\downarrow}}n_{\downarrow}&g_{{}_{\uparrow\downarrow}}\left(\psi_{\sigma}\psi_{\bar{\sigma}}^{*}+n_{\mathrm{sf}}\right)\\\
g_{{}_{\uparrow\downarrow}}\left(\psi_{\bar{\sigma}}\psi_{\sigma}^{*}+n_{\mathrm{sf}}\right)&2gn_{\downarrow}+g_{{}_{\uparrow\downarrow}}n_{\uparrow}\end{pmatrix},$
(8a)
$\mathcal{B}\equiv\begin{pmatrix}g\phi^{2}_{\uparrow}&g_{{}_{\uparrow\downarrow}}\psi_{\uparrow}\psi_{\downarrow}\\\
g_{{}_{\uparrow\downarrow}}\psi_{\uparrow}\psi_{\downarrow}&g\psi^{2}_{\downarrow}\end{pmatrix}.$
(8b)

In these equations, $n_{c\sigma}$ is the condensate density of spin components
$\sigma=\uparrow,\downarrow$ and we have applied the Popov approximation to
ensure a gapless theory [51, 55], i.e., omitting the anomalous densities
$\langle\hat{\eta}^{\dagger}\hat{\eta}^{\dagger}\rangle$ and
$\langle\hat{\eta}\hat{\eta}\rangle$. To take quantum and thermal fluctuations
into account, we have introduced a non-condensate density
$n_{t\sigma}\equiv\langle\hat{\eta}_{\sigma}^{\dagger}\hat{\eta}_{\sigma}\rangle=(1/V)\sum_{{\bf
q\lambda}}[(|u_{{\bf q}\sigma}^{(\lambda)}|^{2}+|v_{{\bf
q}\sigma}^{(\lambda)}|^{2})/(e^{\beta\omega_{{\bf q}\lambda}}-1)+|v_{{\bf
q}\sigma}^{(\lambda)}|^{2}]$ and a spin-flip density term
$n_{\mathrm{sf}}\equiv\langle\hat{\eta}_{\uparrow}^{\dagger}\hat{\eta}_{\downarrow}\rangle=(1/V)\sum_{{\bf
q\lambda}}[(u_{{\bf q}\uparrow}^{(\lambda)}u_{{\bf
q}\downarrow}^{(\lambda)}+v_{{\bf q}\uparrow}^{(\lambda)}v_{{\bf
q}\downarrow}^{(\lambda)})/(e^{\beta\omega_{{\bf q}\lambda}}-1)+v_{{\bf
q}\uparrow}^{(\lambda)}v_{{\bf q}\downarrow}^{(\lambda)}]$ at a temperature
$\beta\equiv 1/(k_{B}T)$. Thus, the total density of spin $\sigma$ is given as
$n_{\sigma}=n_{c\sigma}+n_{t\sigma}$.

At zero temperature, the thermal part of $n_{t\sigma}$ and the spin-flip term
$n_{\mathrm{sf}}$ vanish and we can safely neglect them. Nevertheless, the
condensate is still depleted by a small fraction of the total density, even at
zero temperature, due to quantum fluctuations. This is the so-called quantum
depletion, $n_{\mathrm{qd}}=\sum_{{\bf q}\lambda\sigma}|v_{{\bf
q}\sigma}^{(\lambda)}|^{2}/V$, involving typically about $1\%$ of the total
density in the weakly-interacting regime, i.e.,
$n_{c\uparrow}+n_{c\downarrow}\simeq n$. In our previous work, we have shown
that the SOC effect does not affect remarkably the quantum depletion fraction
in the condensate density (see Fig. 4 of Ref. [35]).

###  II.3 The Ansatz and ground-state phases

We aim to solve the modified GP equation and the coupled Bogoliubov equations
at finite temperature and then derive the current-current response functions.
Before that, for self-containedness let us briefly review the phase diagram at
zero temperature, where the model Hamiltonian of a Raman-type SOC Bose gas in
Sec. II.1 can be solved straightforwardly using a variational formalism.

To obtain the ground-state phases, a variational _first-order_ stripe Ansatz
is usually taken as [5, 56, 57, 35]

$\psi(\mathbf{r})=\sqrt{n}\left[C_{1}\begin{pmatrix}\sin\theta\\\
-\cos\theta\end{pmatrix}e^{-iP_{x}x}+C_{2}\begin{pmatrix}\cos\theta\\\
-\sin\theta\end{pmatrix}e^{iP_{x}x}\right],$ (9)

with the uniform average density $n=N/V$, and the variational angle $\theta$
in the range $[0,\pi/4]$ weighing the spin components of the condensate. By
substituting this trial wave function, Eq. (9), into the model Hamiltonian,
the total energy of the system can be written in terms of $\psi({\bf r})$ as

$\displaystyle E=$ $\displaystyle\int d^{3}{\bf
r}\left[\begin{pmatrix}\psi_{\uparrow}^{*}({\bf r}),\psi_{\downarrow}^{*}({\bf
r})\end{pmatrix}\mathcal{H}_{\mathrm{s}}(\hat{{\bf
p}})\begin{pmatrix}\psi_{\uparrow}({\bf r})\\\ \psi_{\downarrow}({\bf
r})\end{pmatrix}\right.$ (10)
$\displaystyle\left.+\frac{1}{2}g\left(|\psi_{\uparrow}({\bf
r})|^{4}+|\psi_{\downarrow}({\bf
r})|^{4}\right)+g_{\uparrow\downarrow}|\psi_{\uparrow}({\bf
r})|^{2}|\psi_{\downarrow}({\bf r})|^{2}\right].$

After minimizing the ground-state energy with respect to the variational
parameters $C_{1,2}$, $\theta$, and $P_{x}$, one can in general find three
exotic phases, i.e., the stripe (ST), plane-wave (PW), and zero-momentum (ZM)
phases, in the appropriate interaction regimes with
$g>g_{\uparrow\downarrow}$. This condition is necessary for the existence of
the exotic stripe phase [5]. In addition, by defining two interaction
parameters $G_{1}=(g+g_{\uparrow\downarrow})n/4$ and
$G_{2}=(g-g_{\uparrow\downarrow})n/4$, the critical Rabi frequency $\Omega$ of
three phases can be determined respectively by [5]

$\Omega_{c1}=2\left[(2E_{\mathrm{r}}+G_{1})(2E_{\mathrm{r}}-2G_{2})\frac{2G_{2}}{G_{1}+2G_{2}}\right]^{1/2}$
(11a) for the ST-PW phase transition, and $\Omega_{c2}=4E_{\mathrm{r}}-4G_{2}$
(11b)

for the PW-ZM phase transition. As illustrated in our previous work [34],
these critical positions can be notably altered by the quantum and thermal
fluctuations which favor the existence of PW phase. It is worth mentioning
that, the Ansatz of Eq. (9) is a superposition of two plane waves with momenta
$\pm P_{x}$ without including any higher-order terms. Higher-order harmonics
could be non-negligible when the interaction energies $G_{1,2}$ become
relatively large [35]. For example, these high-order contributions will shift
the ST-PW phase transition point $\Omega_{c1}$ to a larger Rabi frequency,
since the higher-order stripe Ansatz can host lower energy solutions than that
of the PW phase.

The consideration of the stripe Ansatz at finite temperature is much more
involved and is beyond the scope of this work. Here, we focus on the plane-
wave and zero-momentum phases, which are characterized by a single-plane-wave
Ansatz with a condensate density $n_{c}$ at a momentum $P_{x}$, and are
obtained by taking $C_{1}=0$ or $C_{2}=0$ in Eq. (9) [32, 33, 44, 34, 35]:

$\psi({\bf r})=\begin{pmatrix}\psi_{\uparrow}\\\
\psi_{\downarrow}\end{pmatrix}e^{iP_{x}x}=\sqrt{n_{c}}\left(\begin{array}[]{c}\cos\theta\\\
-\sin\theta\end{array}\right)e^{iP_{x}x}.$ (12)

With this trial Ansatz, the minimization of the mean-field energy gives rise
to two solutions with respective to Rabi frequency [5, 32, 33, 34]. At
$\Omega<\Omega_{c2}$, the magnetic plane-wave phase appears with a nonzero
momentum $P_{x}=\pm k_{\mathrm{r}}\sqrt{1-\Omega^{2}/\Omega_{c2}^{2}}$
occupied by the condensates, a nonzero magnetization
$\langle\sigma_{z}\rangle\neq 0$, and a typical roton-maxon structure in the
lowest-lying excitation spectrum. In contrast, the nonmagnetic zero-momentum
phase occurs for $\Omega\geq\Omega_{c2}$ where the condensate momentum and the
magnetization are both zero, i.e., $P_{x}=0$ and $\langle\sigma_{z}\rangle=0$,
and there exist only the conventional linear phonon modes in the long-
wavelength limit in the lowest-lying excitation spectrum.

###  II.4 The PW and ZM phases at finite temperature

At a certain finite temperature $T$ and a Rabi frequency $\Omega$, we can also
take this plane-wave Ansatz in terms of variational parameters
($\theta,P_{x}$), and solve self-consistently the GP and Bogoliubov equations
in Eqs. (6) and (7). Thus, the free energy of the system, i.e.,
$\mathcal{F}(\theta,P_{x})$ in Ref. [34], can be calculated straightforwardly.
The condensate wave function $\psi({\bf r})$, the Bogoliubov quasiparticle
wave functions $u_{j}(v_{j})$ and excitation spectrum $\omega_{j}$ can be then
determined by minimizing the free energy $\mathcal{F}$ with respect to two
variational parameters, namely, $(\partial\mathcal{F}/\partial\theta)_{N}=0$
and $(\partial\mathcal{F}/\partial P_{x})_{N}=0$.

The recent experiments have successfully utilized ultracold 87Rb atoms to
realize this Raman-type SOC [24, 39, 38]. For instance in Ref. [38], the
typical interaction energy is $gn=0.38E_{\mathrm{r}}$ with a peak density
$n=0.46k_{\mathrm{r}}^{3}$ in harmonic traps, and the ratio
$g_{{}_{\uparrow\downarrow}}/g=100.99/101.20$ between the inter- and intra-
species interactions. Accordingly, the two critical Rabi frequencies
characterizing the ST-PW and PW-ZM phase transitions at zero temperature, can
be determined respectively from Eqs. (11a) and (11b), i.e.,
$\Omega_{c1}=0.2E_{\mathrm{r}}$ and $\Omega_{c2}=4.0E_{\mathrm{r}}$. Previous
works have shown that the regime of the stripe phase (i.e., $\Omega_{c1}$) can
be actually tuned by the difference between intra- and inter-species
interactions [5, 56]. As we are concentrating on the plane-wave and zero-
momentum phases in this work, we will set
$g_{{}_{\uparrow\uparrow}}=g_{{}_{\downarrow\downarrow}}=g\geq
g_{{}_{\uparrow\downarrow}}$ in the following calculations. Also, we will
consider a relatively small difference between inter- and intra-species
interaction strengths, i.e., zero or relatively small $G_{2}$, in order to
enlarge the window for the plane-wave and zero-momentum phases in the phase
diagram.

##  III Results and Discussions

We are now ready to perform our derivations of the current-current response
functions and extract the superfluid density at finite temperature. To address
the effect of the one-dimensional Raman-type SOC, we will focus only on the
SOC direction (i.e., in the $x$-axis) in the following. In the perpendicular
plane, the SOC does not affect the elementary excitations in the long-
wavelength limit and the superfluid density fraction remains unity at zero
temperature (see Refs. [44, 35]).

###  III.1 Current-current response functions

In this subsection, we study the current-current response functions that are
closely associated with the superfluidity of the SOC Bose gas. Technically,
the superfluid density can be calculated by the difference between the
longitudinal and transverse current-current response functions in the long-
wavelength and low-frequency limits [58, 41, 59].

We start by finding the expression for the current operator in the presence of
the Raman-type SOC. From the equation of motion in Eq. (5) for the Bose field
operators, we obtain the continuity equation [60, 61]

$\displaystyle\partial_{t}n+\nabla\cdot(J^{0}-\frac{2k_{\mathrm{r}}}{m}s_{z}\hat{\bf
e}_{x})$ $\displaystyle=$ $\displaystyle 0.$ (13)

Here $J^{0}=J^{0}_{\uparrow\uparrow}+J^{0}_{\downarrow\downarrow}$ is the
conventional total-density current with the current operators
$\hat{J}^{0}_{\mu\nu}\equiv\frac{1}{2mi}(\hat{\Psi}_{\mu}^{*}\nabla\hat{\Psi}_{\nu}-\nabla\hat{\Psi}_{\mu}^{*}\hat{\Psi}_{\nu})$
in the absence of SOC, and
$\hat{s}_{z}\equiv\frac{1}{2}\sum_{\mu,\nu}\hat{\Psi}_{\mu}^{*}\hat{\sigma}_{z}^{\mu\nu}\hat{\Psi}_{\nu}$
is the spin density. It is clear that due to SOC the total-density current is
coupled to the spin density. Therefore, it is necessary to consider the
modified current operators for each spin component
$\sigma=(\uparrow,\downarrow)$,

$\displaystyle\hat{j}_{\uparrow\uparrow}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{2mi}(\hat{\Psi}_{\uparrow}^{*}\nabla_{x}\hat{\Psi}_{\uparrow}-\nabla_{x}\hat{\Psi}_{\uparrow}^{*}\hat{\Psi}_{\uparrow})-\frac{k_{\mathrm{r}}}{m}\hat{\Psi}_{\uparrow}^{*}\hat{\Psi}_{\uparrow},$
(14a) $\displaystyle\hat{j}_{\downarrow\downarrow}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{2mi}(\hat{\Psi}_{\downarrow}^{*}\nabla_{x}\hat{\Psi}_{\downarrow}-\nabla_{x}\hat{\Psi}_{\downarrow}^{*}\hat{\Psi}_{\downarrow})+\frac{k_{\mathrm{r}}}{m}\hat{\Psi}_{\downarrow}^{*}\hat{\Psi}_{\downarrow},$
(14b)

such that the total-density current $\hat{j}_{d}({\bf
r})\equiv\hat{j}_{\uparrow\uparrow}+\hat{j}_{\downarrow\downarrow}$ reads

$\displaystyle\hat{j}_{d}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{2mi}\sum_{\mu}(\hat{\Psi}_{\mu}^{*}\nabla_{x}\hat{\Psi}_{\mu}-\nabla_{x}\hat{\Psi}_{\mu}^{*}\hat{\Psi}_{\mu})-\frac{2k_{\mathrm{r}}}{m}\hat{s}_{z}.$
(15)

We expand the spin-component current operators using the well-defined Bose
field operator in Eq. (3), and the operators are divided into three parts

$\displaystyle\hat{j}_{\sigma\sigma}({\bf r},t)$ $\displaystyle\equiv$
$\displaystyle j_{0\sigma}({\bf r})+\delta\hat{j}_{\sigma}({\bf
r},t)+\tilde{j}_{\sigma}({\bf r},t),$ (16)

where $j_{0\sigma}$, $\delta\hat{j}_{\sigma}$, and $\tilde{j}_{\sigma}$
represent the zeroth-, first-, and second-order terms of the fluctuation
operator $\hat{\eta}$, respectively, being

$\displaystyle j_{0\sigma}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{2mi}(\psi_{\sigma}^{*}\nabla_{x}\psi_{\sigma}-\nabla_{x}\psi_{\sigma}^{*}\psi_{\sigma})\mp\frac{k_{\mathrm{r}}}{m}\psi_{\sigma}^{*}\psi_{\sigma},$
(17a) $\displaystyle\delta\hat{j}_{\sigma}({\bf r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{2mi}(\psi_{\sigma}^{*}\nabla_{x}\hat{\eta}_{\sigma}+\nabla_{x}\psi_{\sigma}\hat{\eta}_{\sigma}^{*}-\nabla_{x}\psi_{\sigma}^{*}\hat{\eta}_{\sigma}-\psi_{\sigma}\nabla_{x}\hat{\eta}_{\sigma}^{*})$
(17b)
$\displaystyle\mp\frac{k_{\mathrm{r}}}{m}(\psi_{\sigma}^{*}\hat{\eta}_{\sigma}+\hat{\eta}_{\sigma}^{*}\psi_{\sigma}),$
$\displaystyle\tilde{j}_{\sigma}({\bf r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{2mi}(\hat{\eta}_{\sigma}^{*}\nabla_{x}\hat{\eta}_{\sigma}-\nabla_{x}\hat{\eta}_{\sigma}^{*}\hat{\eta}_{\sigma})\mp\frac{k_{\mathrm{r}}}{m}\hat{\eta}_{\sigma}^{*}\hat{\eta}_{\sigma}.$
(17c)

Here the sign ’$\mp$’ is $-$ ($+$) for spin component $\sigma=\uparrow$
($\downarrow$).

Within the Bogoliubov approximation, we are interested only in the current
fluctuation (i.e., $\delta\hat{j}$) or the linear terms of the fluctuation
operator $\hat{\eta}$ [58]. Therefore, within the imaginary-time (i.e.,
$\tau=it$) Green’s function method [62, 63, 64], the time-ordered correlation
function $\langle T_{\tau}\delta\hat{j}_{\sigma}({\bf
r},\tau)\delta\hat{j}_{\sigma^{\prime}}({\bf r^{\prime}})\rangle$ for the
current fluctuation density operator can be derived explicitly in terms of the
condensate wave function $\psi$ and the Bogoliubov quasiparticle wave
functions $u$($v$). The details can be seen in Appendix A. We summarize the
main results as follows.

We take a plane-wave Ansatz in Eq. (12) for the classic condensate wave
function and the time-ordered correlation function of the current fluctuation
density can be consequently rewritten as

$\langle T_{\tau}\delta\hat{j}_{\sigma}({\bf
r},\tau)\delta\hat{j}_{\sigma^{\prime}}({\bf
r^{\prime}})\rangle\equiv\mathcal{A}_{1}+\mathcal{A}_{2}+\mathcal{A}_{3}+\mathcal{A}_{4},$
(18)

in terms of four correlation functions $\mathcal{A}_{m}$ $(m=1,2,3,4)$ for the
fluctuation operator $\hat{\eta}$, which are given explicitly in Appendix A.
Thus, the response functions of the elements in $\mathcal{A}_{m}$ can be
calculated by

$\chi^{\mu\nu}_{m}({\bf r},{\bf r^{\prime}};i\omega_{n})=\int_{0}^{\beta}d\tau
e^{i\omega_{n}\tau}\frac{\theta(\tau)}{V}\mathcal{A}_{m}^{\mu\nu}({\bf r},{\bf
r^{\prime}};\tau)$ (19)

with a step function $\theta(\tau)$, the system volume $V$ and the Matsubara
frequencies $i\omega_{n}\equiv 2n\pi i/\beta$ for bosons. The current response
function is then calculated collectively by using Eq. (17b) as

$\displaystyle\chi_{JJ}({\bf r},{\bf
r^{\prime}};i\omega_{n})\equiv\chi^{\uparrow\uparrow}+\chi^{\uparrow\downarrow}+\chi^{\downarrow\uparrow}+\chi^{\downarrow\downarrow}$
(20) $\displaystyle=$ $\displaystyle\int_{0}^{\beta}d\tau
e^{i\omega_{n}\tau}\theta(\tau)\langle T_{\tau}\delta\hat{j}({\bf
r},\tau)\delta\hat{j}({\bf r^{\prime}})\rangle.$

By employing the analytic continuation $(i\omega_{n}\rightarrow\omega+i\eta)$
with $\eta=0^{+}$, setting $\omega=0$ and then taking the Fourier transform to
momentum space, the ${\bf q}$-component static response function is obtained
in terms of the elements $\chi^{\sigma\sigma^{\prime}}({\bf q};\omega=0)$,
i.e.,

$\displaystyle\chi^{\uparrow\uparrow}$ $\displaystyle=$
$\displaystyle\frac{\psi_{\uparrow}^{2}}{m^{2}V}\sum_{j}\frac{C_{1}u^{2}_{\uparrow}+C_{2}v^{2}_{\uparrow}+C_{3}u_{\uparrow}v_{\uparrow}}{\omega_{j}},$
(21a) $\displaystyle\chi^{\downarrow\downarrow}$ $\displaystyle=$
$\displaystyle\frac{\psi_{\downarrow}^{2}}{m^{2}V}\sum_{j}\frac{C_{4}u^{2}_{\downarrow}+C_{5}v^{2}_{\downarrow}+C_{6}u_{\downarrow}v_{\downarrow}}{\omega_{j}},$
(21b) $\displaystyle\chi^{\uparrow\downarrow}$ $\displaystyle=$
$\displaystyle\chi^{\downarrow\uparrow}=\frac{\psi_{\uparrow}\psi_{\downarrow}}{m^{2}V}\cdot\sum_{j}$
(21c)
$\displaystyle\frac{C_{7}u_{\uparrow}u_{\downarrow}+C_{8}v_{\uparrow}v_{\downarrow}+C_{9}u_{\uparrow}v_{\downarrow}+C_{10}v_{\uparrow}u_{\downarrow}}{\omega_{j}},$

where the index $j$ runs over all the possible single-particle states. Without
confusion, we have made the index $j$ implicit in the quasiparticle wave-
functions $u$ and $v$, i.e., $u_{\sigma}\equiv u_{j,\sigma}$ and
$v_{\sigma}\equiv v_{j,\sigma}$ are calculated at each momentum ${\bf q}$ for
the energy level $\omega_{j}$ (see the notations in Sec. II.2). The
coefficients $C_{n}$ are given by

$\displaystyle C_{1/2}$ $\displaystyle=$ $\displaystyle\frac{1}{2}(2P_{x}\pm
q_{x}-2k_{\mathrm{r}})^{2},$ (22a) $\displaystyle C_{3/6}$ $\displaystyle=$
$\displaystyle(2P_{x}-q_{x}\mp 2k_{\mathrm{r}})(2P_{x}+q_{x}\mp
2k_{\mathrm{r}}),$ (22b) $\displaystyle C_{4/5}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(2P_{x}\pm q_{x}+2k_{\mathrm{r}})^{2},$ (22c)
$\displaystyle C_{7/8}$ $\displaystyle=$ $\displaystyle\frac{1}{2}(2P_{x}\pm
q_{x})^{2}-2k_{\mathrm{r}}^{2},$ (22d) $\displaystyle C_{9/10}$
$\displaystyle=$ $\displaystyle 2P_{x}^{2}-\frac{1}{2}(q_{x}\mp
2k_{\mathrm{r}})^{2}.$ (22e)

The explicit derivation can be seen in Appendix A. Particularly in the zero-
momentum phase with the condensates located at $P_{x}=0$, the coefficients
reduce to

$\displaystyle C^{(\mathrm{ZM})}_{1/2}$ $\displaystyle=$ $\displaystyle
C^{(\mathrm{ZM})}_{5/4}=-C^{(\mathrm{ZM})}_{9/10}=\frac{1}{2}(q_{x}\mp
2k_{\mathrm{r}})^{2},$ $\displaystyle C^{(\mathrm{ZM})}_{3}$ $\displaystyle=$
$\displaystyle
C^{(\mathrm{ZM})}_{6}=-2C^{(\mathrm{ZM})}_{7}=-2C^{(\mathrm{ZM})}_{8}=4k_{\mathrm{r}}^{2}-q_{x}^{2}.$

Finally, the ${\bf q}$-component static current-current response function is
then obtained by the summation

$\chi_{JJ}({\bf
q};\omega=0)=\sum_{\sigma,\sigma^{\prime}}\chi^{\sigma\sigma^{\prime}}({\bf
q};\omega=0).$ (23)

Figure 1: The longitudinal and transverse (i.e., dotted red and dashed blue
lines, respectively) static current-current response functions
$\frac{m}{N}\chi^{L,T}_{JJ}({\bf q};0)$ in the ${\bf q}\to 0$ limit, as a
function of Rabi frequency $\Omega$. Their difference (i.e., solid black line)
reveals the superfluid density fraction. Here we take the parameters as in the
experiment [38], with a total density $n=0.46k_{\mathrm{r}}^{3}$,
$gn=0.38E_{\mathrm{r}}$ and $g_{{}_{\uparrow\downarrow}}/g=100.99/101.20$,
i.e., $G_{1}=0.19E_{\mathrm{r}}$ and $G_{2}\approx 0$.

We consider only the SOC direction (i.e., $x$ axis) and the corresponding
total and normal densities can be then calculated via the longitudinal and
transverse components of the total-current response function $\chi_{JJ}({\bf
q};0)$ in the long wavelength limit [48, 41]. Therefore, the superfluid mass
density along the SOC direction is given by

$\rho_{s}=\rho-\rho_{n}=\rho\frac{m}{N}(\lim_{q_{x}\to 0}\lim_{q_{\perp}\to
0}\chi_{JJ}-\lim_{q_{\perp}\to 0}\lim_{q_{x}\to 0}\chi_{JJ}).$ (24)

An example of the current-current response function of a zero-temperature SOC
87Rb Bose gas in the static and low-wavelength limits is illustrated in Fig.
1, with a total density $n=0.46k_{\mathrm{r}}^{3}$, the interaction energy
$gn=0.38E_{\mathrm{r}}$ and the ratio of inter-, intra-species interaction
strength $g_{{}_{\uparrow\downarrow}}/g=100.99/101.20$ taken from a recent
experiment [38]. The longitudinal static current-current response function or
the total density (i.e., dotted red line) remains unity while the transverse
component or the normal density fraction (i.e., dashed blue line) changes non-
monotonically. Their difference (i.e., solid black line) reveals the
superfluid density fraction, which is significantly affected by the SOC
effect. The zero-temperature superfluid fraction in Fig. 1, explicitly
calculated using the current-current response functions, agrees with the
earlier result obtained via a sum-rule approach [44].

###  III.2 Generalized Josephson relation

Superfluid density is intimately related to Bose-Einstein condensate density.
Nonetheless, they are not exactly equivalent even at zero temperature, for
instance in liquid helium with strong interactions [65], or in weakly-
interacting Bose gases with SOC [44, 35]. In 1965, B.D. Josephson and other
researchers have derived the so-called Josephson relation or Josephson sum
rule for a conventional weakly-interacting Bose gas, which associates the
superfluid mass density $\rho_{s}$ with the condensate density $n_{0}$ at zero
or finite temperatures as [46, 47, 48, 49]

$\rho_{s}=-\lim_{{\bf q}\to 0}\frac{n_{0}m^{2}}{{\bf q}^{2}G_{11}({\bf
q},0)},$ (25)

with the atomic mass $m$ and the single-particle Green’s function at momentum
${\bf q}$ and zero frequency. It is straightforward to show that, for the
conventional weakly-interacting Bose gas at zero temperature, within the
Bogoliubov approximation the superfluid mass density is connected to the
condensate density via the Josephson relation [58]. Here, we are interested in
examining the Josephson relation within the same Bogoliubov approximation at
exactly _zero_ temperature. The situation at nonzero temperature is more
involved, as accurate calculations of the superfluid density and condensate
fraction may require different levels of approximation. Therefore, we defer to
a future study for the nonzero-temperature case.

For the examination at zero temperature, let us first discuss the Josephson
relation in the presence of SOC. In 2018, Yi-Cai Zhang generalizes the
Josephson relation to multi-component Bose gases using the Green’s function
within the linear response theory [50]. In this work, we will follow his
procedure and consider a two-component SOC Bose gas. We focus only on the SOC
direction, i.e., the $x$-axis. Thus, in terms of the order parameter
$\langle\hat{\Psi}\rangle$ and Green’s function matrix in the long-wavelength
limit, the expression for the superfluid density in the SOC direction is given
by [50]

$\rho_{\mathrm{s}}^{\mathrm{(SOC)}}=-\lim_{{\bf q}\rightarrow
0}\frac{1}{q^{2}}\begin{pmatrix}\langle\hat{\Psi}\rangle,\langle\hat{\Psi}^{*}\rangle\end{pmatrix}\sigma^{\mathbb{Z}}\mathcal{G}^{-1}\sigma^{\mathbb{Z}}\begin{pmatrix}\langle\hat{\Psi}\rangle\\\
\langle\hat{\Psi}^{*}\rangle\end{pmatrix}.$ (26)

Here we have introduced a $4\times 4$ matrix
$\sigma^{\mathbb{Z}}\equiv\left(\begin{array}[]{cc}\mathbb{I}&0\\\
0&-\mathbb{I}\end{array}\right)$ and a $2\times 2$ identity matrix
$\mathbb{I}$. The Green’s matrix is defined as

$\displaystyle\mathcal{G}$ $\displaystyle\equiv$
$\displaystyle\begin{pmatrix}G({\bf q},0)&F^{T}({\bf q},0)\\\ (F^{T}({\bf
q},0))^{*}&G^{*}({\bf q},0)\end{pmatrix},$ (27)

and

$\displaystyle G({\bf q},0)$ $\displaystyle\equiv$
$\displaystyle\begin{pmatrix}G_{\uparrow\uparrow}({\bf
q},0)&G_{\uparrow\downarrow}({\bf q},0)\\\ G_{\downarrow\uparrow}({\bf
q},0)&G_{\downarrow\downarrow}({\bf q},0)\end{pmatrix},$ (28a) $\displaystyle
F({\bf q},0)$ $\displaystyle\equiv$
$\displaystyle\begin{pmatrix}F_{\uparrow\uparrow}({\bf
q},0)&F_{\uparrow\downarrow}({\bf q},0)\\\ F_{\downarrow\uparrow}({\bf
q},0)&F_{\downarrow\downarrow}({\bf q},0)\end{pmatrix}.$ (28b)

Here $G_{\sigma\sigma^{\prime}}$ and $F_{\sigma\sigma^{\prime}}$ are normal
and anomalous Green’s function matrix elements. It is worth noting that the
above derivations are assumed for homogeneous systems, which possess the
translational symmetry.

Let us now derive the expression for Green functions. In the presence of the
Raman-type SOC, one can start with a plane-wave Ansatz for the condensate wave
function in the plane-wave and zero-momentum phases [32, 33, 34, 35], i.e.,
$\psi({\bf r})=\begin{pmatrix}\psi_{\uparrow}\\\
\psi_{\downarrow}\end{pmatrix}e^{iP_{x}x}$ as in Eq. (12), satisfying the
Gross-Pitaevskii equation with the condensation momentum $P_{x}\geq 0$ and the
average density $n=N/V$. To the next order, i.e., within the Bogoliubov
approximation, the Hamiltonian for the fluctuation operators is written as
[33]

$\displaystyle\hat{H}_{\mathrm{bogo}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{{\bf q}}\hat{\Phi}^{\dagger}_{{\bf
q}}H_{\mathrm{B}}({\bf q})\hat{\Phi}_{{\bf q}}+\epsilon_{\mathrm{shift}},$
(29)

with $\epsilon_{\mathrm{shift}}=-\frac{1}{2}\left[\xi_{{}_{P_{x}-{\bf
q},\uparrow}}+\xi_{{}_{P_{x}-{\bf
q},\uparrow}}-2\mu+2gn+g_{{}_{\uparrow\downarrow}}n\right]$ being an energy
shift arising from the bosonic commutation relations and $\xi_{\bf
k}\equiv{\bf k}^{2}/(2m)$. Here $H_{\mathrm{B}}$ is the Bogoliubov Hamiltonian
for the $4$-component Nambu spinor $\hat{\Phi}_{{\bf
q}}\equiv[\hat{\phi}_{P_{x}+{\bf q},\uparrow},\hat{\phi}_{P_{x}+{\bf
q},\downarrow},\hat{\phi}^{\dagger}_{P_{x}-{\bf
q},\uparrow},\hat{\phi}^{\dagger}_{P_{x}-{\bf q},\downarrow}]^{T}$ as

$H_{\mathrm{B}}({\bf q})\equiv\begin{bmatrix}K_{0}(P_{x}+{\bf
q})+\Sigma_{\mathrm{N}}&\Sigma_{\mathrm{A}}\\\
\Sigma_{\mathrm{A}}&K_{0}(P_{x}-{\bf q})+\Sigma_{\mathrm{N}}\end{bmatrix},$
(30)

where we have defined three matrices

$\displaystyle K_{0}({\bf q})$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\frac{(q_{x}-k_{\mathrm{r}})^{2}+q_{\perp}^{2}}{2m}-\mu&\frac{\Omega}{2}\\\
\frac{\Omega}{2}&\frac{(q_{x}+k_{\mathrm{r}})^{2}+q_{\perp}^{2}}{2m}-\mu\end{bmatrix},$
(31a) $\displaystyle\Sigma_{\mathrm{N}}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}2gn_{\uparrow}+g_{{}_{\uparrow\downarrow}}n_{\downarrow}&g_{{}_{\uparrow\downarrow}}\psi_{\uparrow}\psi_{\downarrow}\\\
g_{{}_{\uparrow\downarrow}}\psi_{\uparrow}\psi_{\downarrow}&2gn_{\downarrow}+g_{{}_{\uparrow\downarrow}}n_{\uparrow}\end{bmatrix},$
(31b) $\displaystyle\Sigma_{\mathrm{A}}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}gn_{\uparrow}&g_{{}_{\uparrow\downarrow}}\psi_{\uparrow}\psi_{\downarrow}\\\
g_{{}_{\uparrow\downarrow}}\psi_{\uparrow}\psi_{\downarrow}&gn_{\downarrow}\end{bmatrix},$
(31c)

with the spin density $n_{\sigma}\equiv|\psi_{\sigma}|^{2}$ and
$q_{\perp}^{2}\equiv q_{y}^{2}+q_{z}^{2}$.

The gapless condition
$\mathrm{Det}[K_{0}(P_{x})+\Sigma_{\mathrm{N}}-\Sigma_{\mathrm{A}}]=0$ at
${\bf q}=0$ is ensured by the GP equation, with which the chemical potential
can be solved as
$\mu=(P_{x}^{2}+k_{\mathrm{r}}^{2})/2m+(g+g_{{}_{\uparrow\downarrow}})n/2-\sqrt{\left(-P_{x}k_{\mathrm{r}}/m+(g-g_{{}_{\uparrow\downarrow}})(n_{\uparrow}-n_{\downarrow})/2\right)^{2}+\Omega^{2}/4}$.
Hence, the diagonal elements of $K_{0}(P_{x}+{\bf q})$ can be rewritten in
terms of the explicit chemical potential and we can write straightforwardly
the inverse single-particle Green’s function matrix with SOC as

$\displaystyle G^{-1}_{\mathrm{B}}({\bf q},i\omega_{n})=$
$\displaystyle\begin{bmatrix}i\omega_{n}\mathbb{I}-K_{0}(P_{x}+{\bf
q})-\Sigma_{\mathrm{N}}&-\Sigma_{\mathrm{A}}\\\
-\Sigma_{\mathrm{A}}&-i\omega_{n}\mathbb{I}-K_{0}(P_{x}-{\bf
q})-\Sigma_{\mathrm{N}}\end{bmatrix}_{4\times 4}.$ (32)

In particular, in the zero-momentum phase with $P_{x}=0$ and
$\psi_{\uparrow}=-\psi_{\downarrow}=\sqrt{n/2}$, the chemical potential
reduces to
$\mu^{\mathrm{ZM}}=k_{\mathrm{r}}^{2}/2m+(g+g_{{}_{\uparrow\downarrow}})n/2-\Omega/2$,
and the relevant matrices become $K^{\mathrm{ZM}}_{0}({\pm}{\bf
q})=\begin{bmatrix}\frac{{\pm}2k_{\mathrm{r}}q_{x}+{\bf
q}^{2}}{2m}-2G_{1}+\frac{\Omega}{2}&\frac{\Omega}{2}\\\
\frac{\Omega}{2}&\frac{{\mp}2k_{\mathrm{r}}q_{x}+{\bf
q}^{2}}{2m}-2G_{1}+\frac{\Omega}{2}\end{bmatrix}$ and

$\displaystyle\Sigma^{\mathrm{ZM}}_{\mathrm{N}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\begin{bmatrix}2gn+g_{{}_{\uparrow\downarrow}}n&-g_{{}_{\uparrow\downarrow}}n\\\
-g_{{}_{\uparrow\downarrow}}n&2gn+g_{{}_{\uparrow\downarrow}}n\end{bmatrix},$
(33a) $\displaystyle\Sigma^{\mathrm{ZM}}_{\mathrm{A}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\begin{bmatrix}gn&-g_{{}_{\uparrow\downarrow}}n\\\
-g_{{}_{\uparrow\downarrow}}n&gn\end{bmatrix}.$ (33b)

The full elements of inverse Green’s function at zero frequency (after
analytic continuation) then becomes

$G^{-1}_{\mathrm{ZM}}({\bf q},0)=-\begin{bmatrix}K^{\mathrm{ZM}}_{0}({\bf
q})+\Sigma^{\mathrm{ZM}}_{\mathrm{N}}&\Sigma^{\mathrm{ZM}}_{\mathrm{A}}\\\
\Sigma^{\mathrm{ZM}}_{\mathrm{A}}&K^{\mathrm{ZM}}_{0}(-{\bf
q})+\Sigma^{\mathrm{ZM}}_{\mathrm{N}}\end{bmatrix}.$ (34)

Therefore, we can straightforwardly find out the single-particle Green’s
function matrix $G_{\mathrm{B}}({\bf q},0)$ as well as the elements $G$ and
$F$ defined in Eq. (28). Thus, the new matrix $\mathcal{G}$ in the plane-wave
and zero-momentum phases can be then obtained and substituted into Eq. (26) to
calculate the superfluid density. We will verify the generalized Josephson
relation at zero temperature by comparing the resulting superfluid density
fraction with these obtained from other approaches in the following subsection
(see Fig. 2).

###  III.3 Zero-temperature superfluid density from different approaches

Figure 2: Superfluid density fraction
$\rho^{(\mathrm{SOC})}_{\mathrm{s}}/\rho$ in the SOC direction as a function
of Rabi frequency at (a) equal and (b) unequal intra-, inter-spin
interactions. Here the predictions given by the current-current response
function and by the generalized Josephson relation are denoted by red circles
and blue crosses, respectively. The black lines indicate the analytic
prediction using a phase-twist approach in Ref. [35]. The vertical dotted and
dashed curves indicate the critical $\Omega_{c1}$ and $\Omega_{c2}$,
respectively, i.e., Eqs. (11a) and (11b). Here, we take the interaction
energies [$G_{1}$, $G_{2}$]$/E_{\mathrm{r}}$= (a) [0.8, 0] and (b) [0.7, 0.1].

Here, we discuss the superfluid density in the SOC direction at zero
temperature obtained from the current-current response function and from the
generalized Josephson relation. In Fig. 2, the superfluid density faction
$\rho^{(\mathrm{SOC})}_{\mathrm{s}}/\rho$ in the uniform plane-wave and zero-
momentum phase regimes is presented as a function of Rabi frequency at equal
and unequal intra-, inter-spin interaction strengths.

It is clearly seen that the superfluid density
$\rho^{(\mathrm{SOC})}_{\mathrm{s}}/\rho$ along the SOC direction calculated
via Eq. (24) from the current-current response function is identical to the
one obtained via Eq. (26) from the generalized Josephson relation (see the red
circles and blue crosses, respectively). Meanwhile, they agree with our
previous analytic prediction in Eqs. (30-31) of Ref. [35] found using a phase-
twist approach, i.e.,
$\rho^{(x,\mathrm{PW})}_{\mathrm{s}}/\rho=1-\frac{E_{\mathrm{r}}}{(E_{\mathrm{r}}-G_{2})\Omega_{c2}^{2}/\Omega^{2}+G_{2}}$
and
$\rho^{(x,\mathrm{ZM})}_{\mathrm{s}}/\rho=1-\frac{4E_{\mathrm{r}}}{\Omega+4G_{2}}$,
as shown in solid black lines.

The superfluid density exhibits a non-monotonic behavior with respect to Rabi
frequency, decreasing smoothly in the PW phase towards zero at $\Omega_{c2}$
(vertical dashed line), and then rising back in the ZM phase. This non-trivial
behavior at the PW-ZM phase transition point $\Omega_{c2}$ could be understood
from its diverging effective mass, which suppresses the super-flow in the
system [33], and from the vanishing sound velocity or critical velocity
favoring the creation of excitations to destroy its superfluidity [32, 34,
35]. It is worth noting that in Fig. 2(b), the superfluid density in the
stripe phase at $\Omega<\Omega_{c1}$ is merely calculated from the phase-twist
approach in Ref. [35] and there’s a discontinuity at $\Omega_{c1}$ due to the
first-order nature of the ST-PW transition. In contrast, the critical Rabi
frequency $\Omega_{c1}$ in Fig. 2(a) shrinks to zero owing to the vanishing
$G_{2}$.

Finally, we note that the calculation of the current-current response function
and the generalization of the Josephson relation in the stripe phase would be
much more involved. We will consider these two interesting issues in future
works.

###  III.4 Superfluid density, sound velocity and Landau critical velocity at
finite temperature

Figure 3: Superfluid density fraction
$\rho^{(\mathrm{SOC})}_{\mathrm{s}}/\rho$ in the SOC direction as a function
of Rabi frequency at various values of temperature $T/T_{0}=10^{-4}$, $0.1$,
$0.3$, and $0.5$ (i.e., lines) for a SOC 87Rb gas. Here $T_{0}$ is the
critical BEC temperature of an ideal spinless Bose gas with density $n$, i.e.,
$T_{0}=2\pi\hbar^{2}[n/\zeta(3/2)]^{2/3}/(mk_{B})$, and the circles denote the
zero-temperature analytic result from Ref. [35]. Other parameters are the same
as in Fig. 1.

We now turn to discuss the superfluid density at finite temperature using
merely the current-current response function described in Sec. III.1, and
address its relation to sound velocity and Landau critical velocity.

In Fig. 3, we present the behaviour of superfluid density
$\rho_{\mathrm{s}}^{\mathrm{(SOC)}}$ in the SOC direction as a function of
Rabi frequency in a SOC 87Rb Bose gas, at four typical temperatures
$T/T_{0}=10^{-4}$, $0.1$, $0.3$, and $0.5$ (i.e., lines) with $T_{0}$ being
the critical BEC temperature of an ideal spinless Bose gas. We find that the
superfluid density at a tiny temperature $T=10^{-4}T_{0}$ in solid black line
is slightly different from the zero-temperature one denoted by the circles
near the PW-ZM transition point, i.e., $\Omega_{c2}=4E_{\mathrm{r}}$. This is
due to the quantum fluctuations taken into account in the current calculation,
which can slightly shift the critical point and the minimum of sound velocity
as discussed in Ref. [34].

At zero temperature SOC plays a significant role on the superfluid density,
leading to a non-monotonic behaviour with respect to Rabi frequency [44, 35]
as we already discussed in the last subsection. At nonzero temperature,
thermal fluctuations start to play a non-negligible role, by suppressing the
superfluidity and exhausting the condensate in the system. As a result,
generally the superfluid density fraction reduces correspondingly with
increasing temperature, as shown by the dashed, dotted, dashed-dotted colorful
lines in Fig. 3. In particular, as shown by the dotted and dashed lines in
Fig. 4(a), the superfluid density for a wide range of Rabi frequency exhibits
a monotonic decreasing dependence on temperature.

However, near the PW-ZM transition point at $\Omega_{c2}=4E_{\mathrm{r}}$,
this consensus is no longer true. There is a striking non-monotonic behavior
of the superfluid density with respect to temperature, as shown by the solid
blue line in Fig. 4(a). As the temperature increases, the superfluid fraction
increases first from a tiny value to the maximum about 0.1, and then
decreases.

We attribute this interesting non-monotonic temperature dependence to the non-
trivial interplay between the spin-orbit coupling and thermal fluctuations. At
zero temperature, at the critical Rabi frequency the SOC effect strongly
suppresses the superfluidity of the system with a vanishing sound velocity
[36, 44, 35]. While at finite temperature, thermal fluctuations shift the
critical point to larger Rabi frequency favoring the PW phase. As a result, we
anticipate that the sound velocity may become nonzero, which instead restores
the superfluidity of the system near the critical point [38, 34]. The
anticipated temperature dependence of the sound velocity is examined in Fig.
4(b). Indeed, near the transition point we find the very similar non-monotonic
dependence as in the superfluid fraction. In contrast, away from the
transition point in the deep PW and ZM phases (see, i.e., the dashed red and
dotted blue lines for $\Omega/E_{\mathrm{r}}=3$ and $5$, respectively), the
monotonic decreasing dependence on temperature is recovered in the sound
velocity. This is consistent with the monotonic temperature dependence of the
superfluid fraction at the corresponding Rabi frequency.

Figure 4: (a) Superfluid density fraction
$\rho^{(\mathrm{SOC})}_{\mathrm{s}}/\rho$ in the SOC direction, (b) sound
velocities $c_{s}^{\pm}$ and (c) critical velocities $v_{c}^{\pm}$ in the $\pm
x$ directions, as a function of temperature at three sets of Rabi frequency
$\Omega/E_{\mathrm{r}}=3$, $4$, and $5$ for a SOC 87Rb gas. Other parameters
are the same as in Fig. 1.

It is worth mentioning that, the unexpected, counterintuitive effect of
thermal fluctuations near the PW-ZM transition point can persist up to a
relatively large temperature and is responsible to the maximum superfluid
fraction at $T\sim 0.6T_{0}$. As the temperature increases further, thermal
fluctuations eventually depletes the superfluidity of the system. The
superfluid fraction then decreases until reaching the BEC transition
temperature.

We note also that, in the PW phase we have two sound velocities along and
opposite to the SOC direction, owing to the emergence of the roton-maxon
structure in the elementary excitation spectrum [32]. The non-trivial roton
excitation spectrum also gives rise to two distinct Landau critical velocities
$v_{c}^{\pm}$ in a PW phase [43], as shown in Fig. 4(c), which provides
additional information on the superfluidity. It is readily seen that, in the
deep PW regime ($\Omega/E_{\mathrm{r}}=3$), one of the critical velocities in
the dashed red line shows a significant difference with respect to the
corresponding sound velocity due to the pronounced roton structure and a
softened energy gap [34]. Near the transition point
($\Omega/E_{\mathrm{r}}=4$), the critical velocity in the solid blue line
starts to deviate from the sound velocity at certain value of temperature
$T/T_{0}\sim 0.2$. This deviation becomes even larger as the temperature
rises, revealing that thermal fluctuations favor the emergence of the roton
structure near the transition point. The overall temperature dependence of
Landau critical velocities are consistent with that of the superfluid
fraction.

##  IV Conclusions and Outlooks

In conclusions, we have studied the current-current response functions and the
Josephson relation in a three-dimensional weakly interacting Bose gas with
one-dimensional Raman-type spin-orbit coupling. The analytic expression of the
current-current response functions and single-particle Green’s functions are
derived explicitly within the Bogoliubov approximation. The superfluid density
fraction along the spin-orbit coupling direction is calculated using both
approaches at zero temperature. The resulting superfluid density fraction
agrees with previous analytic prediction obtained from a phase-twist approach.
The zero-temperature superfluid density exhibits an intriguing behaviour as
Rabi frequency rises, decreasing smoothly in the plane-wave phase, becoming
zero at the critical point, and then rising back again in the zero-momentum
phase.

At finite temperature, we have calculated the superfluid density by using the
current-current response functions, and have discussed its relation to sound
velocity and Landau critical velocity. A significant non-monotonic temperature
dependence is seen in these quantities near the transition point from the
plane-wave to zero-momentum phases. Future works remain to be done, in order
to calculate the current-current response functions and to establish the
Josephson relation in the more interesting stripe phase.

###### Acknowledgements.

X.L.C. acknowledges the insightful discussions with Zeng-Qiang Yu and the
support from Swinburne University of Technology through an HDR Publication
Award. Our research was supported by the Science Foundation of Zhejiang Sci-
Tech University (ZSTU), Grant No. 21062339-Y (X.L.C.), and China Postdoctoral
Science Foundation, Grant No. 2020M680495 (X.L.C.), the Australian Research
Council’s (ARC) Discovery Project, Grant No. DP180102018 (X.J.L.).

##  Appendix A The current-current response functions at finite temperature

In this Appendix, we give an explicit derivation of the current-current
response functions and other results in Sec. III.1, in terms of the condensate
wave function and the Bogoliubov wave functions at finite temperature.

By decoupling the Bose field operator
$\hat{\Psi}_{\sigma}=\psi_{\sigma}+\hat{\eta}_{\sigma}$ in Eq. (3), we first
expand the spin-component current operators in Eq. (14) for spin component
$\sigma$ in the presence of the Raman-type spin-orbit coupling as
$\hat{j}_{\sigma\sigma}({\bf r},t)=j_{0\sigma}({\bf
r})+\delta\hat{j}_{\sigma}({\bf r},t)+\tilde{j}_{\sigma}({\bf r},t)$, where

$\displaystyle j_{0\sigma}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{2mi}(\psi_{\sigma}^{*}\nabla_{x}\psi_{\sigma}-\nabla_{x}\psi_{\sigma}^{*}\psi_{\sigma})\mp\frac{k_{\mathrm{r}}}{m}\psi_{\sigma}^{*}\psi_{\sigma}$
(35a) $\displaystyle\delta\hat{j}_{\sigma}({\bf r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{2mi}(\psi_{\sigma}^{*}\nabla_{x}\hat{\eta}_{\sigma}+\nabla_{x}\psi_{\sigma}\hat{\eta}_{\sigma}^{*}-\nabla_{x}\psi_{\sigma}^{*}\hat{\eta}_{\sigma}-\psi_{\sigma}\nabla_{x}\hat{\eta}_{\sigma}^{*})\mp\frac{k_{\mathrm{r}}}{m}(\psi_{\sigma}^{*}\hat{\eta}_{\sigma}+\hat{\eta}_{\sigma}^{*}\psi_{\sigma})$
(35b) $\displaystyle\tilde{j}_{\sigma}({\bf r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{2mi}(\hat{\eta}_{\sigma}^{*}\nabla_{x}\hat{\eta}_{\sigma}-\nabla_{x}\hat{\eta}_{\sigma}^{*}\hat{\eta}_{\sigma})\mp\frac{k_{\mathrm{r}}}{m}\hat{\eta}_{\sigma}^{*}\hat{\eta}_{\sigma},$
(35c)

with the sign ’$\mp$’ is $-$ ($+$) for spin component $\sigma=\uparrow$
($\downarrow$).

Within the Bogoliubov approximation, we are interested in the current
fluctuation (i.e., $\delta\hat{j}$) only and the time-ordered correlation
function for the current fluctuation density operator is given by

$\displaystyle\langle T_{t}\delta\hat{j}_{\sigma}({\bf
r},t)\delta\hat{j}_{\sigma^{\prime}}({\bf r^{\prime}})\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{m^{2}}\begin{pmatrix}\psi^{*}_{\sigma}({\bf
r})&\psi_{\sigma}({\bf r})&\nabla_{x}\psi^{*}_{\sigma}({\bf
r})&\nabla_{x}\psi_{\sigma}({\bf
r})\end{pmatrix}\begin{pmatrix}\mathcal{A}&\mathcal{B}\\\
\mathcal{C}&\mathcal{D}\end{pmatrix}\begin{pmatrix}\psi^{*}_{\sigma^{\prime}}({\bf
r^{\prime}})\\\ \psi_{\sigma^{\prime}}({\bf r^{\prime}})\\\
\nabla_{x}\psi^{*}_{\sigma^{\prime}}({\bf r^{\prime}})\\\
\nabla_{x}\psi_{\sigma^{\prime}}({\bf r^{\prime}})\end{pmatrix},$ (36)

where we have introduced $\mathrm{sgn}(\sigma)=-1$ (1) for $\sigma=\uparrow$
($\downarrow$), and the matrices read

$\displaystyle\mathcal{A}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\begin{pmatrix}-\langle\nabla_{x}\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\nabla_{x}\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}+\frac{\mathrm{sgn}(\sigma^{\prime})k_{\mathrm{r}}}{2i}\begin{pmatrix}\langle\nabla_{x}\hat{\eta}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\nabla_{x}\hat{\eta}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
-\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}$
(37a)
$\displaystyle+\frac{\mathrm{sgn}(\sigma)k_{\mathrm{r}}}{2i}\begin{pmatrix}\langle\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
\langle\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}+\mathrm{sgn}(\sigma)\mathrm{sgn}(\sigma^{\prime})k_{\mathrm{r}}^{2}\begin{pmatrix}\langle\hat{\eta}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\hat{\eta}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix},$
$\displaystyle\mathcal{B}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\begin{pmatrix}\langle\nabla_{x}\hat{\eta}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\nabla_{x}\hat{\eta}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
-\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}+\frac{\mathrm{sgn}(\sigma)k_{\mathrm{r}}}{2i}\begin{pmatrix}-\langle\hat{\eta}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\hat{\eta}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
-\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix},$
(37b) $\displaystyle\mathcal{C}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\begin{pmatrix}\langle\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
-\langle\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}+\frac{\mathrm{sgn}(\sigma^{\prime})k_{\mathrm{r}}}{2i}\begin{pmatrix}-\langle\hat{\eta}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\hat{\eta}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix},$
(37c) $\displaystyle\mathcal{D}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\begin{pmatrix}-\langle\hat{\eta}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\hat{\eta}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}.$
(37d)

Here $\langle\cdots\rangle$ is the thermal average in the statistical
equilibrium and
$\hat{\eta}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\equiv\hat{\eta}_{\sigma}({\bf
r},t)\hat{\eta}_{\sigma^{\prime}}({\bf r^{\prime}})$. After taking a plane-
wave Ansatz $\psi_{\sigma}({\bf r})=\psi_{\sigma}e^{iP_{x}x}$ for the
condensate wave function, we can further rewrite the fluctuation correlation
function in the imaginary time ($\tau=it$) as

$\displaystyle\langle T_{\tau}\delta\hat{j}_{\sigma}({\bf
r},\tau)\delta\hat{j}_{\sigma^{\prime}}({\bf r^{\prime}})\rangle$
$\displaystyle=$
$\displaystyle\mathcal{A}_{1}+\mathcal{A}_{2}+\mathcal{A}_{3}+\mathcal{A}_{4}$
(38)

with

$\displaystyle\mathcal{A}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{4m^{2}}\begin{pmatrix}\psi^{*}_{\sigma}({\bf
r})&\psi_{\sigma}({\bf
r})\end{pmatrix}\begin{pmatrix}-\langle\nabla_{x}\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\nabla_{x}\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}\begin{pmatrix}\psi^{*}_{\sigma^{\prime}}({\bf
r^{\prime}})\\\ \psi_{\sigma^{\prime}}({\bf r^{\prime}})\end{pmatrix},$ (39a)
$\displaystyle\mathcal{A}_{2}$ $\displaystyle=$
$\displaystyle\frac{P_{x}+2bk_{\mathrm{r}}}{4m^{2}i}\begin{pmatrix}\psi^{*}_{\sigma}({\bf
r})&\psi_{\sigma}({\bf
r})\end{pmatrix}\begin{pmatrix}\langle\nabla_{x}\hat{\eta}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\nabla_{x}\hat{\eta}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
-\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\nabla_{x}\hat{\eta}^{*}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}\begin{pmatrix}\psi^{*}_{\sigma^{\prime}}({\bf
r^{\prime}})\\\ \psi_{\sigma^{\prime}}({\bf r^{\prime}})\end{pmatrix},$ (39b)
$\displaystyle\mathcal{A}_{3}$ $\displaystyle=$
$\displaystyle\frac{P_{x}+2ak_{\mathrm{r}}}{4m^{2}i}\begin{pmatrix}\psi^{*}_{\sigma}({\bf
r})&\psi_{\sigma}({\bf
r})\end{pmatrix}\begin{pmatrix}\langle\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\hat{\eta}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
\langle\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}_{\sigma^{\prime}}\rangle&-\langle\hat{\eta}^{*}_{\sigma}\nabla_{x}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}\begin{pmatrix}\psi^{*}_{\sigma^{\prime}}({\bf
r^{\prime}})\\\ \psi_{\sigma^{\prime}}({\bf r^{\prime}})\end{pmatrix},$ (39c)
$\displaystyle\mathcal{A}_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{m^{2}}\left[abk_{\mathrm{r}}^{2}+\frac{k_{\mathrm{r}}P_{x}}{2}(a+b)+\frac{P_{x}^{2}}{4}\right]\begin{pmatrix}\psi^{*}_{\sigma}({\bf
r})&\psi_{\sigma}({\bf
r})\end{pmatrix}\begin{pmatrix}\langle\hat{\eta}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\hat{\eta}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\\\
\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}_{\sigma^{\prime}}\rangle&\langle\hat{\eta}^{*}_{\sigma}\hat{\eta}^{*}_{\sigma^{\prime}}\rangle\end{pmatrix}\begin{pmatrix}\psi^{*}_{\sigma^{\prime}}({\bf
r^{\prime}})\\\ \psi_{\sigma^{\prime}}({\bf r^{\prime}})\end{pmatrix}$ (39d)

and the coefficients $a\equiv\mathrm{sgn}(\sigma)$,
$b\equiv\mathrm{sgn}(\sigma^{\prime})$.

At finite temperature $\beta\equiv 1/k_{B}T$, the current response function is
then calculated collectively by using Eq. (17b) as

$\displaystyle\chi_{JJ}({\bf r},{\bf
r^{\prime}};i\omega_{n})\equiv\chi^{\uparrow\uparrow}+\chi^{\uparrow\downarrow}+\chi^{\downarrow\uparrow}+\chi^{\downarrow\downarrow}=\int_{0}^{\beta}d\tau
e^{i\omega_{n}\tau}\frac{\theta(\tau)}{V}\langle T_{\tau}\delta\hat{j}({\bf
r},\tau)\delta\hat{j}({\bf r^{\prime}})\rangle.$ (40)

We can take an example for illustration of calculating $\langle
T_{\tau}\delta\hat{j}_{\sigma}({\bf
r},\tau)\delta\hat{j}_{\sigma^{\prime}}({\bf r^{\prime}})\rangle$ in Eq. (38).
For instance,

$\displaystyle\mathcal{A}_{1}^{11}({\bf r},{\bf
r^{\prime}};\tau)=\frac{1}{4m^{2}}\psi^{*}_{\sigma}({\bf
r})(-)\langle\nabla_{x}\hat{\eta}_{\sigma}({\bf
r},\tau)\nabla_{x}\hat{\eta}_{\sigma^{\prime}}({\bf
r^{\prime}})\rangle\psi^{*}_{\sigma^{\prime}}({\bf r^{\prime}})$ (41)
$\displaystyle=$ $\displaystyle\frac{1}{4m^{2}}\psi^{*}_{\sigma}({\bf
r})(-)\psi^{*}_{\sigma^{\prime}}({\bf r^{\prime}})\sum_{j\equiv({\bf
q}_{1},\lambda_{1}),k\equiv({\bf
q}_{2},\lambda_{2})}\left\langle\nabla_{x}\left[e^{iP_{x}x}\left(u^{(\lambda_{1})}_{{\bf
q}_{1}\sigma}e^{i{\bf
q_{1}r}}\hat{\alpha}_{j}e^{-\omega_{j}\tau}+(v^{(\lambda_{1})}_{{\bf
q}_{1}\sigma})^{*}e^{-i{\bf
q_{1}r}}\hat{\alpha}^{\dagger}_{j}e^{\omega_{j}\tau}\right)\right]\right.$
$\displaystyle\left.\nabla_{x}\left[e^{iP_{x}x}\left(u^{(\lambda_{2})}_{{\bf
q}_{2}\sigma^{\prime}}e^{i{\bf q}_{2}{\bf
r^{\prime}}}\hat{\alpha}_{k}+(v^{(\lambda_{2})}_{{\bf
q}_{2}\sigma^{\prime}})^{*}e^{-i{\bf q}_{2}{\bf
r^{\prime}}}\hat{\alpha}^{\dagger}_{k}\right)\right]\right\rangle$
$\displaystyle=$
$\displaystyle\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}}e^{i{\bf
q(r-r^{\prime})}}\sum_{j}(P_{x}^{2}-q_{x}^{2})\left[u_{\sigma}v_{\sigma^{\prime}}\left(1+f_{B}(\omega_{j})\right)e^{-\omega_{j}\tau}+v_{\sigma}u_{\sigma^{\prime}}f_{\mathrm{B}}(\omega_{j})e^{\omega_{j}\tau}\right],$

with the expanded fluctuation operator $\hat{\eta}_{\sigma}({\bf
r},t)=e^{iP_{x}x}\sum_{j\equiv({\bf q},\lambda)}\left(u^{(\lambda)}_{{\bf
q}\sigma}e^{i{\bf qr}}\hat{\alpha}_{j}e^{-i\omega_{j}t}+(v^{(\lambda)}_{{\bf
q}\sigma})^{*}e^{-i{\bf qr}}\hat{\alpha}^{\dagger}_{j}e^{i\omega_{j}t}\right)$
in a quasiparticle basis and $u_{\sigma}\equiv u^{(\lambda)}_{{\bf q}\sigma}$
and $v_{\sigma}\equiv v^{(\lambda)}_{{\bf q}\sigma}$ defined here and
hereafter.
$f_{\mathrm{B}}(\omega_{j})\equiv\langle\hat{\alpha}^{\dagger}_{j}\hat{\alpha}_{j}\rangle$
is the occupation number for $j$-th quasiparticle with energy $\omega_{j}$
satisfying the bosonic statistical distribution function as
$f_{\mathrm{B}}(\omega_{j})=1/(\mathrm{exp}(\beta\omega_{j})-1)$.

After integrating
$\int_{0}^{\beta}e^{(i\omega_{n}-\omega)\tau}\theta(\tau)d\tau=\frac{e^{\beta(i\omega_{n}-\omega)}-1}{i\omega_{n}-\omega}$
with the Matsubara frequencies $i\omega_{n}\equiv 2n\pi i/\beta$ for bosons
and $e^{\beta i\omega_{n}}=1$, the response function becomes

$\displaystyle\chi^{11}_{1}({\bf r},{\bf r^{\prime}};i\omega_{n})$
$\displaystyle=$ $\displaystyle\int_{0}^{\beta}d\tau
e^{i\omega_{n}\tau}\frac{\theta(\tau)}{V}\mathcal{A}_{1}^{11}({\bf r},{\bf
r^{\prime}};\tau)$ (42) $\displaystyle=$
$\displaystyle\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}e^{i{\bf
q(r-r^{\prime})}}\sum_{j}(P_{x}^{2}-q_{x}^{2})\int_{0}^{\beta}d\tau\theta(\tau)\left[u_{\sigma}v_{\sigma^{\prime}}\left(1+f_{B}(\omega_{j})\right)e^{(i\omega_{n}-\omega_{j})\tau}+v_{\sigma}u_{\sigma^{\prime}}f_{\mathrm{B}}(\omega_{j})e^{(i\omega_{n}+\omega_{j})\tau}\right]$
$\displaystyle=$
$\displaystyle\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}e^{i{\bf
q(r-r^{\prime})}}\sum_{j}(P_{x}^{2}-q_{x}^{2})\left[\frac{v_{\sigma}u_{\sigma^{\prime}}}{i\omega_{n}+\omega_{j}}-\frac{u_{\sigma}v_{\sigma^{\prime}}}{i\omega_{n}-\omega_{j}}\right],$

where we have used the relations
$(e^{-\beta\omega_{j}}-1)\left(1+f_{B}(\omega_{j})\right)=-1$ and
$(e^{\beta\omega_{j}}-1)f_{\mathrm{B}}(\omega_{j})=1$. The real-time response
function can be then obtained by employing the analytic continuation
$(i\omega_{n}\rightarrow\omega+i\eta)$ with $\eta=0^{+}$ as

$\chi^{11}_{1}({\bf r},{\bf
r^{\prime}};\omega+i\eta)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}e^{i{\bf
q(r-r^{\prime})}}\sum_{j}(P_{x}^{2}-q_{x}^{2})\left[\frac{v_{\sigma}u_{\sigma^{\prime}}}{\omega+i\eta+\omega_{j}}-\frac{u_{\sigma}v_{\sigma^{\prime}}}{\omega+i\eta-\omega_{j}}\right].$
(43)

Similarly, another three response functions can be obtained as

$\displaystyle\chi^{12}_{1}({\bf r},{\bf r^{\prime}};\omega+i\eta)$
$\displaystyle=$
$\displaystyle\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}e^{i{\bf
q(r-r^{\prime})}}\sum_{j}\left[(P_{x}-q_{x})^{2}\frac{v_{\sigma}v_{\sigma^{\prime}}}{\omega+i\eta+\omega_{j}}-(P_{x}+q_{x})^{2}\frac{u_{\sigma}u_{\sigma^{\prime}}}{\omega+i\eta-\omega_{j}}\right],$
(44a) $\displaystyle\chi^{21}_{1}({\bf r},{\bf r^{\prime}};\omega+i\eta)$
$\displaystyle=$
$\displaystyle\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}e^{i{\bf
q(r-r^{\prime})}}\sum_{j}\left[(P_{x}+q_{x})^{2}\frac{u_{\sigma}u_{\sigma^{\prime}}}{\omega+i\eta+\omega_{j}}-(P_{x}-q_{x})^{2}\frac{v_{\sigma}v_{\sigma^{\prime}}}{\omega+i\eta-\omega_{j}}\right],$
(44b) $\displaystyle\chi^{22}_{1}({\bf r},{\bf r^{\prime}};\omega+i\eta)$
$\displaystyle=$
$\displaystyle\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}e^{i{\bf
q(r-r^{\prime})}}\sum_{j}(P_{x}^{2}-q_{x}^{2})\left[\frac{u_{\sigma}v_{\sigma^{\prime}}}{\omega+i\eta+\omega_{j}}-\frac{v_{\sigma}u_{\sigma^{\prime}}}{\omega+i\eta-\omega_{j}}\right].$
(44c)

The static response for ${\bf q}$-component can be then obtained by setting
$\omega=0$ and taking a Fourier transform to momentum space, and the equations
now read

$\displaystyle\chi^{11}_{1}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{22}_{1}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}\sum_{j}(P_{x}^{2}-q_{x}^{2})\left[\frac{u_{\sigma}v_{\sigma^{\prime}}+v_{\sigma}u_{\sigma^{\prime}}}{\omega_{j}}\right],$
(45a) $\displaystyle\chi^{12}_{1}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{21}_{1}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}\sum_{j}\left[\frac{(P_{x}+q_{x})^{2}u_{\sigma}u_{\sigma^{\prime}}+(P_{x}-q_{x})^{2}v_{\sigma}v_{\sigma^{\prime}}}{\omega_{j}}\right].$
(45b)

Similarly, the remaining static response functions for $\mathcal{A}_{2,3,4}$
can be obtained as

$\displaystyle\chi^{11}_{2}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{22}_{2}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}\sum_{j}(P_{x}+2bk_{\mathrm{r}})\left[\frac{(P_{x}+q_{x})u_{\sigma}v_{\sigma^{\prime}}+(P_{x}-q_{x})v_{\sigma}u_{\sigma^{\prime}}}{\omega_{j}}\right],$
(46a) $\displaystyle\chi^{12}_{2}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{21}_{2}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}\sum_{j}(P_{x}+2bk_{\mathrm{r}})\left[\frac{(P_{x}+q_{x})u_{\sigma}u_{\sigma^{\prime}}+(P_{x}-q_{x})v_{\sigma}v_{\sigma^{\prime}}}{\omega_{j}}\right],$
(46b) $\displaystyle\chi^{11}_{3}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{22}_{3}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}\sum_{j}(P_{x}+2ak_{\mathrm{r}})\left[\frac{(P_{x}-q_{x})u_{\sigma}v_{\sigma^{\prime}}+(P_{x}+q_{x})v_{\sigma}u_{\sigma^{\prime}}}{\omega_{j}}\right],$
(46c) $\displaystyle\chi^{12}_{3}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{21}_{3}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}\sum_{j}(P_{x}+2ak_{\mathrm{r}})\left[\frac{(P_{x}+q_{x})u_{\sigma}u_{\sigma^{\prime}}+(P_{x}-q_{x})v_{\sigma}v_{\sigma^{\prime}}}{\omega_{j}}\right],$
(46d) $\displaystyle\chi^{11}_{4}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{22}_{4}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{m^{2}V}\sum_{j}\left[abk_{\mathrm{r}}^{2}+\frac{k_{\mathrm{r}}P_{x}}{2}(a+b)+\frac{P_{x}^{2}}{4}\right]\left[\frac{u_{\sigma}v_{\sigma^{\prime}}+v_{\sigma}u_{\sigma^{\prime}}}{\omega_{j}}\right],$
(46e) $\displaystyle\chi^{12}_{4}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{21}_{4}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{m^{2}V}\sum_{j}\left[abk_{\mathrm{r}}^{2}+\frac{k_{\mathrm{r}}P_{x}}{2}(a+b)+\frac{P_{x}^{2}}{4}\right]\left[\frac{u_{\sigma}u_{\sigma^{\prime}}+v_{\sigma}v_{\sigma^{\prime}}}{\omega_{j}}\right].$
(46f)

Therefore, the static response functions are derived by summarizing
$\chi^{\sigma\sigma^{\prime}}({\bf
q};0)=\sum_{m,\mu\nu}\left[\chi^{\mu\nu}_{m}({\bf
q};0)\right]_{\sigma\sigma^{\prime}}$ for spin components
$\sigma\sigma^{\prime}$ as

$\displaystyle\chi^{11}_{\sigma\sigma^{\prime}}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{22}_{\sigma\sigma^{\prime}}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}\sum_{j}\frac{B_{1}u_{\sigma}v_{\sigma^{\prime}}+B_{2}v_{\sigma}u_{\sigma^{\prime}}}{\omega_{j}},$
(47) $\displaystyle\chi^{12}_{\sigma\sigma^{\prime}}({\bf q};0)$
$\displaystyle=$ $\displaystyle\chi^{21}_{\sigma\sigma^{\prime}}({\bf
q};0)=\frac{\psi_{\sigma}\psi_{\sigma^{\prime}}}{4m^{2}V}\sum_{j}\frac{B_{3}u_{\sigma}u_{\sigma^{\prime}}+B_{4}v_{\sigma}v_{\sigma^{\prime}}}{\omega_{j}},$
(48)

with $B_{1}=(2ak_{\mathrm{r}}+2P_{x}+q_{x})(2bk_{\mathrm{r}}+2P_{x}-q_{x})$,
$B_{2}=(2ak_{\mathrm{r}}+2P_{x}-q_{x})(2bk_{\mathrm{r}}+2P_{x}+q_{x})$,
$B_{3}=(2ak_{\mathrm{r}}+2P_{x}+q_{x})(2bk_{\mathrm{r}}+2P_{x}+q_{x})$ and
$B_{4}=(2ak_{\mathrm{r}}+2P_{x}-q_{x})(2bk_{\mathrm{r}}+2P_{x}-q_{x})$. Thus,
for spin index
$\sigma\sigma^{\prime}=(\uparrow\uparrow,\downarrow\downarrow,\uparrow\downarrow,\downarrow\uparrow)$,
we obtain the explicit expressions for Eq. (21) in the main text as

$\displaystyle\chi^{\uparrow\uparrow}({\bf q};0)$ $\displaystyle=$
$\displaystyle\sum_{ij}\chi^{ij}_{\uparrow\uparrow}=\frac{\psi_{\uparrow}^{2}}{m^{2}V}\sum_{j}\frac{C_{1}u_{\uparrow}u_{\uparrow}+C_{2}v_{\uparrow}v_{\uparrow}+C_{3}u_{\uparrow}v_{\uparrow}}{\omega_{j}},$
(49a) $\displaystyle\chi^{\downarrow\downarrow}({\bf q};0)$ $\displaystyle=$
$\displaystyle\sum_{ij}\chi^{ij}_{\downarrow\downarrow}=\frac{\psi_{\downarrow}^{2}}{m^{2}V}\sum_{j}\frac{C_{4}u_{\downarrow}u_{\downarrow}+C_{5}v_{\downarrow}v_{\downarrow}+C_{6}u_{\downarrow}v_{\downarrow}}{\omega_{j}},$
(49b) $\displaystyle\chi^{\uparrow\downarrow}({\bf q};0)$ $\displaystyle=$
$\displaystyle\chi^{\downarrow\uparrow}({\bf
q};0)=\sum_{ij}\chi^{ij}_{\uparrow\downarrow}=\frac{\psi_{\uparrow}\psi_{\downarrow}}{m^{2}V}\sum_{j}\frac{C_{7}u_{\uparrow}u_{\downarrow}+C_{8}v_{\uparrow}v_{\downarrow}+C_{9}u_{\uparrow}v_{\downarrow}+C_{10}v_{\uparrow}u_{\downarrow}}{\omega_{j}},$
(49c)

with

$\displaystyle C_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(2P_{x}+q_{x}-2k_{\mathrm{r}})^{2},\quad
C_{2}=\frac{1}{2}(2P_{x}-q_{x}-2k_{\mathrm{r}})^{2},$ (50a) $\displaystyle
C_{3}$ $\displaystyle=$
$\displaystyle(2P_{x}-q_{x}-2k_{\mathrm{r}})(2P_{x}+q_{x}-2k_{\mathrm{r}}),$
(50b) $\displaystyle C_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(2P_{x}+q_{x}+2k_{\mathrm{r}})^{2},\quad
C_{5}=\frac{1}{2}(2P_{x}-q_{x}+2k_{\mathrm{r}})^{2},$ (50c) $\displaystyle
C_{6}$ $\displaystyle=$
$\displaystyle(2P_{x}-q_{x}+2k_{\mathrm{r}})(2P_{x}+q_{x}+2k_{\mathrm{r}}),$
(50d) $\displaystyle C_{7}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(2P_{x}+q_{x})^{2}-2k_{\mathrm{r}}^{2},\quad
C_{8}=\frac{1}{2}(2P_{x}-q_{x})^{2}-2k_{\mathrm{r}}^{2},$ (50e) $\displaystyle
C_{9}$ $\displaystyle=$ $\displaystyle
2P_{x}^{2}-\frac{1}{2}(q_{x}-2k_{\mathrm{r}})^{2},\quad
C_{10}=2P_{x}^{2}-\frac{1}{2}(q_{x}+2k_{\mathrm{r}})^{2}.$ (50f)

In the zero-momentum phase with $P_{x}=0$, the coefficients $C_{n}$ reduces to
$C_{1}=C_{5}=-C_{9}=\frac{1}{2}(q_{x}-2k_{\mathrm{r}})^{2}$,
$C_{2}=C_{4}=-C_{10}=\frac{1}{2}(q_{x}+2k_{\mathrm{r}})^{2}$ and
$C_{3}=C_{6}=-2C_{7}=-2C_{8}=4k_{\mathrm{r}}^{2}-q_{x}^{2}$.

Eventually, the static current-current response function, i.e., Eq. (20) in
the main text, is given by

$\chi_{JJ}({\bf q};0)=\chi^{\uparrow\uparrow}({\bf
q};0)+\chi^{\downarrow\downarrow}({\bf q};0)+\chi^{\uparrow\downarrow}({\bf
q};0)+\chi^{\downarrow\uparrow}({\bf q};0).$ (51)

## References

  * Galitski and Spielman [2013] V. Galitski and I. B. Spielman, Spin-orbit coupling in quantum gases, Nature (London) 494, 49 (2013).
  * Zhai [2015] H. Zhai, Degenerate quantum gases with spin-orbit coupling: a review, Reports on Progress in Physics 78, 026001 (2015).
  * Wang  _et al._ [2010] C. Wang, C. Gao, C.-M. Jian, and H. Zhai, Spin-orbit coupled spinor bose-einstein condensates, Phys. Rev. Lett. 105, 160403 (2010).
  * Ho and Zhang [2011] T.-L. Ho and S. Zhang, Bose-einstein condensates with spin-orbit interaction, Phys. Rev. Lett. 107, 150403 (2011).
  * Li  _et al._ [2012a] Y. Li, L. P. Pitaevskii, and S. Stringari, Quantum tricriticality and phase transitions in spin-orbit coupled bose-einstein condensates, Phys. Rev. Lett. 108, 225301 (2012a).
  * Hu  _et al._ [2012] H. Hu, B. Ramachandhran, H. Pu, and X.-J. Liu, Spin-orbit coupled weakly interacting bose-einstein condensates in harmonic traps, Phys. Rev. Lett. 108, 010402 (2012).
  * Gong  _et al._ [2011] M. Gong, S. Tewari, and C. Zhang, Bcs-bec crossover and topological phase transition in 3d spin-orbit coupled degenerate fermi gases, Phys. Rev. Lett. 107, 195303 (2011).
  * Hu  _et al._ [2011] H. Hu, L. Jiang, X.-J. Liu, and H. Pu, Probing anisotropic superfluidity in atomic fermi gases with rashba spin-orbit coupling, Phys. Rev. Lett. 107, 195304 (2011).
  * Yu and Zhai [2011] Z.-Q. Yu and H. Zhai, Spin-orbit coupled fermi gases across a feshbach resonance, Phys. Rev. Lett. 107, 195305 (2011).
  * Zhu  _et al._ [2006] S.-L. Zhu, H. Fu, C.-J. Wu, S.-C. Zhang, and L.-M. Duan, Spin hall effects for cold atoms in a light-induced gauge potential, Phys. Rev. Lett. 97, 240401 (2006).
  * Liu  _et al._ [2007] X.-J. Liu, X. Liu, L. C. Kwek, and C. H. Oh, Optically induced spin-hall effect in atoms, Phys. Rev. Lett. 98, 026602 (2007).
  * Beeler  _et al._ [2013] M. C. Beeler, R. A. Williams, L. J. Jimenez-Garcia, K.and LeBlanc, A. R. Perry, and I. B. Spielman, The spin hall effect in a quantum gas, Nature (London) 498, 201 (2013).
  * Sato  _et al._ [2009] M. Sato, Y. Takahashi, and S. Fujimoto, Non-abelian topological order in $s$-wave superfluids of ultracold fermionic atoms, Phys. Rev. Lett. 103, 020401 (2009).
  * Jiang  _et al._ [2011] L. Jiang, T. Kitagawa, J. Alicea, A. R. Akhmerov, D. Pekker, G. Refael, J. I. Cirac, E. Demler, M. D. Lukin, and P. Zoller, Majorana fermions in equilibrium and in driven cold-atom quantum wires, Phys. Rev. Lett. 106, 220402 (2011).
  * Liu  _et al._ [2012] X.-J. Liu, L. Jiang, H. Pu, and H. Hu, Probing majorana fermions in spin-orbit-coupled atomic fermi gases, Phys. Rev. A 85, 021603 (2012).
  * Liu and Hu [2012] X.-J. Liu and H. Hu, Topological superfluid in one-dimensional spin-orbit-coupled atomic fermi gases, Phys. Rev. A 85, 033622 (2012).
  * Liu and Hu [2013] X.-J. Liu and H. Hu, Topological fulde-ferrell superfluid in spin-orbit-coupled atomic fermi gases, Phys. Rev. A 88, 023622 (2013).
  * Liu  _et al._ [2014] X.-J. Liu, K. T. Law, and T. K. Ng, Realization of 2d spin-orbit interaction and exotic topological orders in cold atoms, Phys. Rev. Lett. 112, 086401 (2014).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Liu  _et al._ [2009] X.-J. Liu, M. F. Borunda, X. Liu, and J. Sinova, Effect of induced spin-orbit coupling for atoms via laser fields, Phys. Rev. Lett. 102, 046402 (2009).
  * Spielman [2009] I. B. Spielman, Raman processes and effective gauge potentials, Phys. Rev. A 79, 063613 (2009).
  * Dalibard  _et al._ [2011] J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Colloquium: Artificial gauge potentials for neutral atoms, Rev. Mod. Phys. 83, 1523 (2011).
  * Goldman  _et al._ [2014] N. Goldman, G. Juzeliūnas, P. Öhberg, and I. B. Spielman, Light-induced gauge fields for ultracold atoms, Reports on Progress in Physics 77, 126401 (2014).
  * Lin  _et al._ [2011] Y.-J. Lin, K. Jimenez-Garcia, and I. B. Spielman, Spin-orbit-coupled bose-einstein condensates, Nature (London) 471, 83 (2011).
  * Wang  _et al._ [2012] P. Wang, Z.-Q. Yu, Z. Fu, J. Miao, L. Huang, S. Chai, H. Zhai, and J. Zhang, Spin-orbit coupled degenerate fermi gases, Phys. Rev. Lett. 109, 095301 (2012).
  * Cheuk  _et al._ [2012] L. W. Cheuk, A. T. Sommer, Z. Hadzibabic, T. Yefsah, W. S. Bakr, and M. W. Zwierlein, Spin-injection spectroscopy of a spin-orbit coupled fermi gas, Phys. Rev. Lett. 109, 095302 (2012).
  * Wu  _et al._ [2016] Z. Wu, L. Zhang, W. Sun, X.-T. Xu, B.-Z. Wang, S.-C. Ji, Y. Deng, S. Chen, X.-J. Liu, and J.-W. Pan, Realization of two-dimensional spin-orbit coupling for bose-einstein condensates, Science 354, 83 (2016).
  * Huang  _et al._ [2016] L. Huang, Z. Meng, P. Wang, P. Peng, S.-L. Zhang, L. Chen, D. Li, Q. Zhou, and J. Zhang, Experimental realization of two-dimensional synthetic spin-orbit coupling in ultracold fermi gases, Nat Phys 12, 540 (2016).
  * Chen  _et al._ [2018a] H.-R. Chen, K.-Y. Lin, P.-K. Chen, N.-C. Chiu, J.-B. Wang, C.-A. Chen, P.-P. Huang, S.-K. Yip, Y. Kawaguchi, and Y.-J. Lin, Spin–orbital-angular-momentum coupled bose-einstein condensates, Phys. Rev. Lett. 121, 113204 (2018a).
  * Zhang  _et al._ [2019] D. Zhang, T. Gao, P. Zou, L. Kong, R. Li, X. Shen, X.-L. Chen, S.-G. Peng, M. Zhan, H. Pu, and K. Jiang, Ground-state phase diagram of a spin-orbital-angular-momentum coupled bose-einstein condensate, Phys. Rev. Lett. 122, 110402 (2019).
  * Li  _et al._ [2012b] Y. Li, G. I. Martone, and S. Stringari, Sum rules, dipole oscillation and spin polarizability of a spin-orbit coupled quantum gas, EPL (Europhysics Letters) 99, 56008 (2012b).
  * Martone  _et al._ [2012] G. I. Martone, Y. Li, L. P. Pitaevskii, and S. Stringari, Anisotropic dynamics of a spin-orbit-coupled bose-einstein condensate, Phys. Rev. A 86, 063621 (2012).
  * Zheng  _et al._ [2013] W. Zheng, Z.-Q. Yu, X. Cui, and H. Zhai, Properties of bose gases with the raman-induced spin–orbit coupling, Journal of Physics B: Atomic, Molecular and Optical Physics 46, 134007 (2013).
  * Chen  _et al._ [2017] X.-L. Chen, X.-J. Liu, and H. Hu, Quantum and thermal fluctuations in a raman spin-orbit-coupled bose gas, Phys. Rev. A 96, 013625 (2017).
  * Chen  _et al._ [2018b] X.-L. Chen, J. Wang, Y. Li, X.-J. Liu, and H. Hu, Quantum depletion and superfluid density of a supersolid in raman spin-orbit-coupled bose gases, Phys. Rev. A 98, 013614 (2018b).
  * Zhang  _et al._ [2012] J.-Y. Zhang, S.-C. Ji, Z. Chen, L. Zhang, Z.-D. Du, B. Yan, G.-S. Pan, B. Zhao, Y.-J. Deng, H. Zhai, S. Chen, and J.-W. Pan, Collective dipole oscillations of a spin-orbit coupled bose-einstein condensate, Phys. Rev. Lett. 109, 115301 (2012).
  * Khamehchi  _et al._ [2014] M. A. Khamehchi, Y. Zhang, C. Hamner, T. Busch, and P. Engels, Measurement of collective excitations in a spin-orbit-coupled bose-einstein condensate, Phys. Rev. A 90, 063624 (2014).
  * Ji  _et al._ [2015] S.-C. Ji, L. Zhang, X.-T. Xu, Z. Wu, Y. Deng, S. Chen, and J.-W. Pan, Softening of roton and phonon modes in a bose-einstein condensate with spin-orbit coupling, Phys. Rev. Lett. 114, 105301 (2015).
  * Ji  _et al._ [2014] S.-C. Ji, J.-Y. Zhang, L. Zhang, Z.-D. Du, W. Zheng, Y.-J. Deng, H. Zhai, S. Chen, and J.-W. Pan, Experimental determination of the finite-temperature phase diagram of a spin-orbit coupled bose gas, Nat Phys 10, 314 (2014).
  * Lifshitz and Pitaevskii [2013] E. Lifshitz and L. Pitaevskii,  _Statistical physics: theory of the condensed state_, Course of Theoretical Physics No. pt. 2 (Elsevier Science, 2013) Chap. III.
  * Pitaevskii and Stringari [2016] L. Pitaevskii and S. Stringari,  _Bose-Einstein Condensation and Superfluidity_, International Series of Monographs on Physics, Vol. 164 (Oxford University Press, Oxford, 2016).
  * Zhu  _et al._ [2012] Q. Zhu, C. Zhang, and B. Wu, Exotic superfluidity in spin-orbit coupled bose-einstein condensates, EPL (Europhysics Letters) 100, 50003 (2012).
  * Yu [2017] Z.-Q. Yu, Landau criterion for an anisotropic bose-einstein condensate, Phys. Rev. A 95, 033618 (2017).
  * Zhang  _et al._ [2016] Y.-C. Zhang, Z.-Q. Yu, T. K. Ng, S. Zhang, L. Pitaevskii, and S. Stringari, Superfluid density of a spin-orbit-coupled bose gas, Phys. Rev. A 94, 033635 (2016).
  * Martone and Stringari [2021] G. I. Martone and S. Stringari, Supersolid phase of a spin-orbit-coupled bose-einstein condensate: A perturbation approach, SciPost Phys. 11, 92 (2021).
  * Hohenberg and Martin [1965] P. Hohenberg and P. Martin, Microscopic theory of superfluid helium, Annals of Physics 34, 291 (1965).
  * Josephson [1966] B. Josephson, Relation between the superfluid density and order parameter for superfluid he near tc, Physics Letters 21, 608 (1966).
  * Baym [1968] G. Baym, The microscopic description of superfluidity, in  _Mathematical Methods in Solid State and Superfluid Theory_ (Springer, 1968) pp. 121–156.
  * Holzmann and Baym [2007] M. Holzmann and G. Baym, Condensate superfluidity and infrared structure of the single-particle green’s function: The josephson relation, Phys. Rev. B 76, 092502 (2007).
  * Zhang [2018] Y.-C. Zhang, Generalized josephson relation for conserved charges in multicomponent bosons, Phys. Rev. A 98, 033611 (2018).
  * Griffin [1996] A. Griffin, Conserving and gapless approximations for an inhomogeneous bose gas at finite temperatures, Phys. Rev. B 53, 9341 (1996).
  * Dodd  _et al._ [1998] R. J. Dodd, M. Edwards, C. W. Clark, and K. Burnett, Collective excitations of bose-einstein-condensed gases at finite temperatures, Phys. Rev. A 57, R32 (1998).
  * Buljan  _et al._ [2005] H. Buljan, M. Segev, and A. Vardi, Incoherent matter-wave solitons and pairing instability in an attractively interacting bose-einstein condensate, Phys. Rev. Lett. 95, 180401 (2005).
  * Dalfovo  _et al._ [1999] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of bose-einstein condensation in trapped gases, Rev. Mod. Phys. 71, 463 (1999).
  * Chen  _et al._ [2015] X.-L. Chen, Y. Li, and H. Hu, Collective modes of a harmonically trapped one-dimensional bose gas: The effects of finite particle number and nonzero temperature, Phys. Rev. A 91, 063631 (2015).
  * Martone  _et al._ [2014] G. I. Martone, Y. Li, and S. Stringari, Approach for making visible and stable stripes in a spin-orbit-coupled bose-einstein superfluid, Phys. Rev. A 90, 041604 (2014).
  * Yu [2014] Z.-Q. Yu, Equation of state and phase transition in spin-orbit-coupled bose gases at finite temperature: A perturbation approach, Phys. Rev. A 90, 053608 (2014).
  * Griffin [1993] A. Griffin,  _Excitations in a Bose-condensed Liquid_, Cambridge Studies in Low Temperature Physics (Cambridge University Press, 1993).
  * Nozieres and Pines [2018] P. Nozieres and D. Pines,  _Theory Of Quantum Liquids: Superfluid Bose Liquids_ (CRC Press, 2018) Chap. 4.
  * Stringari [2017] S. Stringari, Diffused vorticity and moment of inertia of a spin-orbit coupled bose-einstein condensate, Phys. Rev. Lett. 118, 145302 (2017).
  * Tang and Zhang [2018] W. H. Tang and S. Zhang, Quantum spin dynamics in a normal bose gas with spin-orbit coupling, Phys. Rev. Lett. 121, 120403 (2018).
  * Fetter and Walecka [2003] A. L. Fetter and J. D. Walecka,  _Quantum theory of many-particle systems_ (Courier Corporation, 2003).
  * Mahan [2013] G. D. Mahan,  _Many-particle physics_ (Springer Science & Business Media, 2013).
  * Liu  _et al._ [2004] X.-J. Liu, H. Hu, A. Minguzzi, and M. P. Tosi, Collective oscillations of a confined bose gas at finite temperature in the random-phase approximation, Phys. Rev. A 69, 043605 (2004).
  * Penrose and Onsager [1956] O. Penrose and L. Onsager, Bose-einstein condensation and liquid helium, Phys. Rev. 104, 576 (1956).